Adding gradient of a vector to the Fourier collections class

[qiUip]

Added a method to the FourierCollections class to calculate the gradient of a vector field, returning a matrix. This is a general gradient operator on a vector field (see issue #122).
The unit tests added follow the same structure and functions used for the gradient and divergence unit tests.

In the current code base, this method does not have any use. It will form a component of subsequent solvers such as the ElasticSolver (see issue #120) and potentially more in the future (such as a Maxwell solver).

[github-actions]

Cyclomatic Complexity Check

L 97 (C 22-41): The McCabe cyclomatic complexity of 'initialiseSensorData' is 19.
L 153 (C 24-39): The McCabe cyclomatic complexity of 'recordSensorData' is 6.
L 159 (C 24-37): The McCabe cyclomatic complexity of 'AcousticSolver' is 5.
L 189 (C 29-40): The McCabe cyclomatic complexity of 'get.pressure' is 1.
L 192 (C 33-48): The McCabe cyclomatic complexity of 'get.densitySplit' is 1.
L 195 (C 29-40): The McCabe cyclomatic complexity of 'get.velocity' is 1.
L 38 (C 21-39): The McCabe cyclomatic complexity of 'autoComputeTimeStep' is 3.
L 34 (C 10-24): The McCabe cyclomatic complexity of 'executeTimeStep' is 25.
L 60 (C 7): The McCabe cyclomatic complexity of the anonymous function on line 60 is 1.
L 61 (C 9): The McCabe cyclomatic complexity of the anonymous function on line 61 is 1.
L 62 (C 12): The McCabe cyclomatic complexity of the anonymous function on line 62 is 1.
L 63 (C 14): The McCabe cyclomatic complexity of the anonymous function on line 63 is 1.
L 67 (C 54-77): The McCabe cyclomatic complexity of 'buildAbsorptionOperators' is 2.
L 33 (C 23-42): The McCabe cyclomatic complexity of 'kappaSplitCorrection' is 2.
L 46 (C 10-34): The McCabe cyclomatic complexity of 'setAbsorptionCoefficients' is 5.

L 59 (C 14-17): The McCabe cyclomatic complexity of 'curl' is 2.
L 64 (C 14-23): The McCabe cyclomatic complexity of 'divergence' is 1.
L 66 (C 14-28): The McCabe cyclomatic complexity of 'divergenceSplit' is 16.
L 58 (C 16-28): The McCabe cyclomatic complexity of 'fracLaplacian' is 4.
L 64 (C 15-22): The McCabe cyclomatic complexity of 'gradient' is 15.
L 69 (C 15-28): The McCabe cyclomatic complexity of 'gradientVector' is 16.
L 60 (C 17-23): The McCabe cyclomatic complexity of 'stagger' is 17.
L 52 (C 10-21): The McCabe cyclomatic complexity of 'validateSize' is 7.
L 68 (C 24-43): The McCabe cyclomatic complexity of 'validateNonNegScalar' is 2.
L 84 (C 24-41): The McCabe cyclomatic complexity of 'validateIndexField' is 6.
L 99 (C 24-32): The McCabe cyclomatic complexity of 'Materials' is 1.
L 124 (C 32-42): The McCabe cyclomatic complexity of 'addMaterial' is 20.
L 222 (C 22-34): The McCabe cyclomatic complexity of 'listMaterials' is 4.
L 254 (C 22-40): The McCabe cyclomatic complexity of 'listMaterialIndices' is 1.
L 263 (C 26-41): The McCabe cyclomatic complexity of 'getMaterialNames' is 1.
L 270 (C 45-76): The McCabe cyclomatic complexity of 'getMaterialNamesAndIndicesSorted' is 2.
L 128 (C 24-29): The McCabe cyclomatic complexity of 'Medium' is 1.
L 141 (C 24-31): The McCabe cyclomatic complexity of 'subsasgn' is 5.
L 171 (C 24-30): The McCabe cyclomatic complexity of 'refresh' is 1.
L 179 (C 18-40): The McCabe cyclomatic complexity of 'set.soundSpeedReference' is 2.
L 185 (C 18-39): The McCabe cyclomatic complexity of 'set.diffusionReference' is 2.
L 196 (C 24-46): The McCabe cyclomatic complexity of 'updateDerivedPaddedMaps' is 10.
L 246 (C 25-35): The McCabe cyclomatic complexity of 'mapProperty' is 6.
L 284 (C 22-40): The McCabe cyclomatic complexity of 'get.absorptionPower' is 4.
L 57 (C 36-56): The McCabe cyclomatic complexity of 'get.totalSensorPoints' is 1.
L 67 (C 33-49): The McCabe cyclomatic complexity of 'ProcessSensorData' is 4.
L 87 (C 24-29): The McCabe cyclomatic complexity of 'Solver' is 5.
L 81 (C 24-39): The McCabe cyclomatic complexity of 'TimeDomainSolver' is 1.
L 110 (C 18-37): The McCabe cyclomatic complexity of 'setInitialConditions' is 1.
L 57 (C 10-12): The McCabe cyclomatic complexity of 'run' is 5.
L 1 (C 1): The McCabe cyclomatic complexity is 0.

[codecov]

Codecov Report

❌ Patch coverage is 78.78788% with 14 lines in your changes missing coverage. Please review.
✅ Project coverage is 98.71%. Comparing base (c474ac2) to head (10715ab).
⚠️ Report is 28 commits behind head on develop.

Files with missing lines	Patch %	Lines
...wave/+toolbox/@FourierCollocation/gradientVector.m	78.78%	14 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop     #123      +/-   ##
===========================================
- Coverage   100.00%   98.71%   -1.29%     
===========================================
  Files           59       48      -11     
  Lines          911     1093     +182     
===========================================
+ Hits           911     1079     +168     
- Misses           0       14      +14     

Flag	Coverage Δ
unit	98.71% <78.78%> (-1.29%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[qiUip]

I realise the case with a N-dimensional kappa is not being tested for this method. However, I cannot see where it is tested for the gradient or divergence methods either. If it is tested, I could follow a similar approach for using the N-dimensional k-space correction so please point me to the right place in the code for this.

[btreeby]

@qiUip Some initial comments after discussions with @bencox

[btreeby]

If the operator is defined as the outer product between [dx; dy; dz] * [fx, fy, fz] the output would be the transpose of this. Is there a code reason for having it this way around, for example, when you combine with other tensors?

[qiUip]

Not really. I normally use this notation and when looking through the literature it also seemed like a correct approach see this for instance. The numerator layout (as above) is commonly used for Jacobians, which I guess is what I had in the back of my mind.

Happy to use whatever convention you prefer and I will make sure the convention is consistent in the solver as well.

[btreeby]

Depending on what we do with the staggering flags (see later comment), we should make some comments about the staggering scheme that is used. Currently, it's not generic, but connected to the scheme we are using in the elastic code (similarly for the divergenceSplitTensor function).

[qiUip]

Agree, lets see what we decide regarding the staggering in general and how we want to generalise the API, then we will address this as appropriate.

[btreeby]

Aside from staggering, this could be implemented by calling the gradient function on each of the vector components. However, this would mean we would need to have a more general way of defining the staggering in the gradient function.

If we could think of a nice API for this, it would potentially allow us to make the functions more general, including, the staggering used in divergenceSplitTensor (which currently has the elastic code staggering built in). Noting that in 3D, if you allow forward, backward, and none, there are 27 combinations, so just a look up table would be cumbersome.

The nice thing about re-using gradient is that any low level optimisations we do on this function would then be automatically used by gradientVector.

[qiUip]

Maybe something like 'staggerXPos' and 'staggerXNeg' as options (yes, there would be a lot of options this way, but it would be rather transparent?) could work, where XPos means YNeg and ZNeg and XNeg would mean YPos and ZPos? I lookup table is more elegant programmatically, but less readable / transparent to the user so I'm not sure what the 'best' option is.

I understand the issue with the staggering here though being elastic-solver specific, but given that calling a gradient of a vector would return a tensor, which I "assume" would be used for a stress component also in other solves, would it not be the case that this staggering strategy, which mathematically not being general, would always be used in practice by other solvers using this operator?

Might there also be some advantages to having the gradient and gradientVector operators separated, such as having the ability to apply different low-level optimisation strategies directly on the 5D array rather than the components of it? From the API perspective, it could also help end users knowing that from a gradient operator you always get a vector, and from a gradientVector operator you always get a tensor (I personally do not like that and would prefer to have a single operator with the output dependent on the input, but it might be worth thinking what would be better for end users as well).

I will let @btreeby and @bencox have a think and whatever the final decision is will implement it in the next development phase.

[qiUip]

As there seems to be some discussion which will have to take place, should I mark this as draft in the meanwhile?