Finding edge dofs based on values in integration points

[ettaka]

Assuming we know edge element dofs we can calculate values at integration points: $a\mathrm{ip}=R\mathrm{w} a_\mathrm{dofs}$

or

$a\mathrm{ip}^T=a\mathrm{dofs} R_\mathrm{w}^T$

Explicitly in code: a_ip(1:3) = MATMUL(a_dofs(nd-np), wbasis(nd-np, 1:3))

Assuming we know values at integration points, we need to do the reverse, solve the linear system: $a\mathrm{dofs} = R\mathrm{w}^{-1} a_\mathrm{ip}$

However, $R_\mathrm{w}$ is not a square matrix, it is

$3\times (n{d}-n{p})$, ($3\times6$ in first Whitney form).

This system is under determined. Though, we can assemble all the integration points to the same matrix, then we have a matrix size

$3n\mathrm{ip}\times (n\mathrm{d}-n_\mathrm{p})$, ($12\times6$ in the first Whitney form)

So to solve this, one can do

$a\mathrm{dofs} = (R^T R\mathrm{w})^{-1} R^T a_\mathrm{ip}$

We thus need a subroutine SUBROUTINE InvWbasis(Wbasis, a_ip) that can generate $(R^T R_\mathrm{w})^{-1} R^T$ based on curl basis and values of variable at integration points.

[ettaka]

@mmalinen what do you think about this?

[mmalinen]

I believe creating directly a square (N x N)-matrix, with N the count of edge element DOFs, would be the most natural approach. If one proceeds in a point-wise and coordinate-wise manner to obtain new rows, one could just stop adding new equations when N constraints have been obtained. The pattern of points where the matching has been done may then look unsymmetric, but I believe it shouldn't be a major problem ?

[ettaka]

Maybe I made a mistake earlier when creating NxN, but I had trouble finding the inverse. I made this in a new branch: Edge dofs linear system built here at the moment: https://github.com/ElmerCSC/elmerfem/blob/sr_decomposition/fem/src/modules/MagnetoDynamics/WhitneyAVSolver.F90#L1927

Test case: https://github.com/ElmerCSC/elmerfem/blob/sr_decomposition/fem/tests/mgdyn_sr_decomp/sif/box.sif

I get garbage out at the moment, If you see an obvious mistake, let me know.

Thanks!

[mmalinen]

To obtain a local interpolant having a guaranteed error estimate, I'd use a different approach and start by writing an additional integration loop in order to integrate the local mass matrix and evaluate suitable linear forms for the local interpolation. That is, I'd solve

$$ (Wi,\sum\limits{j=1}^N W_j) d_j = (W_i,A) $$

where $W_i$ are the edge basis functions and $A$ is the function to approximate. After this computation of the degrees of freedom $d_k$, one could evaluate the value of the approximation within another loop over integration points. The evaluation could of course be done for any point in the element.

[raback]

This may be a little similar as InterpolateMeshToMesh.F90 -> Ip2DgFieldInElement() The purpose of which is to allow on-the-fly computation of dg fields for postprocessing.

[ettaka]

Thanks for the info @mmalinen @raback !

Does this look about right to you https://github.com/ElmerCSC/elmerfem/blob/sr_decomposition/fem/src/modules/MagnetoDynamics/WhitneyAVSolver.F90#L1928

I am still getting garbage btw.

[ettaka]

I am trying to debug the field by directly outputting the field.

I get the variable https://github.com/ElmerCSC/elmerfem/blob/e401ba524d76df5afb559954d24818bf02ab412c/fem/src/modules/MagnetoDynamics/WhitneyAVSolver.F90#L539

But when I try to set a vector [1,2,3] everywhere (as a test) https://github.com/ElmerCSC/elmerfem/blob/e401ba524d76df5afb559954d24818bf02ab412c/fem/src/modules/MagnetoDynamics/WhitneyAVSolver.F90#L1981

I don't get exactly correct result in Paraview (see the small blue dots, it should be completely red)

Please advice. Thanks!

[ettaka]

Coming back to this. The previous seems to be due to face elements showing at the same time as body elements. But if we filter only the body elements, the field seems fine for debugging purposes.

I am now showing here the vector potential (z component, as all others are 0) from which we calculate the magnetic flux density: So this should be equivalent to a homogeneic magnetic flux density in y-direction (and other components 0)

This is the resulting magnetic flux density: x direction y direction z direction

This is clearly not what we expect.

So it seems there is some issue with the algorithm to compute the magnetic flux density.

Does the block at https://github.com/ElmerCSC/elmerfem/blob/4c822dc8aac012d1e71a7e05469b37560a86cf84/fem/src/modules/MagnetoDynamics/WhitneyAVSolver.F90#L1935 look correct to you @raback?

[mmalinen]

It seems that the DOFs (the content of the array Asloc after the first BLOCK) are computed correctly. I verified this with some simple source fields which should be contained in the FE space so that an accurate FE representation should be possible. These tests included the case of constant fields and the field $A = (z,0,-x)$, which produces constant curl along the y-axis.

However, the loop

DO k = 1,AsVar % DOFs AsVar % Values( AsVar % DOFs*(ind(t)-1)+k) = SRDAatIp(k) END DO

makes valgrind to say "Invalid read of size 4". I commented this out, since I checked only the first BLOCK which doesn't need AsVar % Values after writing it.

[mmalinen]

It seems that your code at the moment (https://github.com/ElmerCSC/elmerfem/commit/4c822dc8aac012d1e71a7e05469b37560a86cf84) supposes the variables AsVar and BsVar to have their value arrays suitable for representing the variable at integration points. However, they are defined to be elementwise variables and therefore picking indices as Ind(t) is not right when you write the arrays AsVar % Values and BsVar % Values.

[raback]

I was going through old issues and looked into this. My recommendation would be to use the Elemental routinies here since the initial field may be more accurate directly on IP points. Does not really make sense to go through the nodal points. Also the resulting fields may be saved directly as IP fields and let the code take care of the interpolation.

Unfortunately I do not be able to get the projection to work any better than you. For my eye it should work. The initial vector potential looks good as an IP field. But it fitted version and the resulting rotor does not.