[Getdp] OT: Theoretical Questions/Gauge Conditions

[Patrick Dular]

Dear Bernhard,

Here are answers to some of your questions...

Kubicek Bernhard wrote:

>Hello,
>I have already been reading (little)/collecting papers on gauge conditions in numerical electrodynamics, and honestly, there are some things which I would like to ask some professionals:
>
>The edge-Whitney elements are often defined by the node-Whitney elements n_i (which are also sometimes called hat functions): e_ij=w_i grad w_j-w_j grad w_i When one plots these (see http://kariert.org/files/edge.pdf), they are not really intuitive. However, by adding/subtracting them pairwise some structure appears. The added ones show homogenous fields, while the subtraction of two looks like carousel rot-fields with a origin that is centered in the remaining edge. a)Is the reason for defining the edge Whitney elements in this way, that they are divergence free and/or that they are rot-constant within the element ? If so, is there another edge-Whitney element that is rot-free and div-constant?
>  
>
Other edge elements, or more generally curl-conform elements exist, for 
example for higher order interpolations.

An edge element that is rot-free (curl-free) is obtained as the curl of 
a nodal element (adding conditions for considering possible 
non-single-valued scalar potentials).

>If one ones wants to calculate the mag. vec. potential (A) from a given current distribution, the "normal" way seems to define A at the edges. Then however, there would be be many possible solutions that would lead to zero residual. In [1] there is stated that it is still possible to solve this "ungauged" A using iterative methods. When I tried this with getDP, the preconditioning ILUTP complained about zero rows iirc (which seems logical since the Matrix is singulary). 
>
You are right. An iterative method (CG method) for solving the system 
should work with ungauged problems. The problem is here with the 
preconditioning ILUTP.

>One way of solving this problem seems to be to take the graph of nodes/edges and leave away as many edges, so that all nodes are still connected. Only on those edges the equations for A are solved while for the "loop-creating" ones A is set to zero.
>  
>
This is exactly what is done with the gauge condition fixing a zero 
circulation for A along the edges of a tree.

>b) Why is it relevant to tell getDP from which surfaces this tree "starts" ? Or does this states that all edges on the declared surfaces and all of a random cotree have zero A ?
>  
>
The way, the tree is complete on these surfaces (first built on them 
before being extended in the volumes). This is to avoid 
incompatibilities between the gauge condition (zero circulation fixed on 
the edges of the whole tree) and the boundary or interface conditions 
(essential conditions, i.e directly fixing degrees of freedom). As an 
example, think of a tree built in a wrong way:  1. all the nodes of a 
face on a boundary are only connected with edges coming from the 
adjacent volume;
2. This means that the path of tree edges connecting two of these nodes 
is contained in the volume; 3. Because the circulation is fixed to zero 
along all the edges of this path (the gauge condition) and because we 
can fix also a zero boundary condition on the face, there would be faces 
in the volume through which no flux of the curl of the interpolated 
field would flow, which is a condition we do not want at all.

You can find explanations about that in the paper

"A discrete sequence associated with mixed finite elements and its gauge 
condition for vector potentials", P. Dular, A. Nicolet, A. Genon, W. 
Legros, IEEE Transactions on Magnetics, Vol.31, No.3, (1995), pp. 1356-1359.

>c) By playing around we found that for faces where A has a symmetry boundary condition and no source terms (so b is normal to the surface), there should be no TreeStartingOn gauge condition defined. Is this correct?
>  
>
Yes. With the magnetic vector potential formulation, a symmetry 
condition with a normal b field is a natural boundary condition (not an 
essential one). This means that the constraint in b) is not necessary.

>d) How to correctly specify a boundary condition for a face where there is some current density normal to the face, and where  b should be considered only parallel to surface (the "physically gauged" A  then would be normal to the surface)? An additional Galerkin Term (CompZ[{a}],{a'}) In SampleSurfaceXY  with no gauge condition starting from this face seemed not to do work.
>  
>
This is a boundary condition (BC) of type n.b=0 (i.e., the normal 
component of b is zero). If you have non-connected surfaces of this type 
between which there is no flux of b, the BC on the magnetic vector 
potential A is nxA=0 (i.e., its tangential component is zero). This BC 
is thus directly defined in the Constraint in the FunctionSpace of A 
(fixing a zero value for all the edge degrees of freedom on the surface).

>e) Does the "tree-gauge-trick" also work on hexahedral and mixed (i.E. containing tetrahedral, hexahedral and "pyramidal") element meshes? Because I tried, and failed for some unknown reason: The resolution had two steps: Calculate el. potential from voltage boundary conditions and a spatial conductivity profile, all non-transient. Afterwards,  with the A-edge-tree gauge, calculate mag. vec. pot. A (using rect-shell transformations) from the current -cond*{d V}. I get nice results for a tetrahedral mesh, but when I go hexahedral or mixed with the _same_  .pro file, the calculated results become very counter-intuitive.
>  
>
It should work. The problem should be elsewhere in your .pro file.

>f) Is there a way to output the edge values in the post-processing, instead of some probably interpolated nodal values?
>  
>
There is no PostOperation for that. This information can nevertheless be 
found in the .res file to be put in parallel with the .pre file to 
localize the edge numbers and their position in the system of equations.

>g) Is there a good read on the "h-phi" formulation _apart_ from Bossavit? Finding some papers on the net is a pain, because of the unfortunate name of the formulation.
>  
>
Perhaps these papers can help:

* "A generalized source magnetic field calculation method for inductors 
of any shape", P. Dular, F. Henrotte, F. Robert, A. Genon, W. Legros, 
IEEE Transactions on Magnetics, Vol.33, No.2, (1997), pp. 1398-1401.
* "A natural method for coupling magnetodynamic h-formulations and 
circuit equations", P. Dular, C. Geuzaine, W. Legros, IEEE Transactions 
on Magnetics, Vol.35, No.3, (1999), pp. 1626-1629.
* "Dual magnetodynamic formulations and their source fields associated 
with massive and stranded inductors", P. Dular, P. Kuo-Peng, C. 
Geuzaine, N. Sadowski, J.P.A. Bastos, IEEE Transactions on Magnetics, 
Vol.36, No.4, (2000), pp. 1293-1299.

>h) Could one expect a big performance leap by exploiting the h-phi formulation compared to a "global A" formulation, if the source electric current for the magnetic field is only existing in one fith of the elements,  while the remaining cells are in a simple connected region?
>  
>
Yes.

>i) Are the following statements reasonable: 
>The matrix element M_ij originating from a Galerkin term (L_1 Dof{A},{L_2 A'}) is defined by 
>\int dv (L_1 l_i)(L_2 l_j)
>with the basis functions l of the field A. 
>The actual value is approximated using a numerical Gauss-Integration, by calculating a weighted sum of the values (L_1 l_i)*(L_2 l_j) at a finite number of well chosen positions wihin the shared support. So, if the operator L1 is changing polynomial with degree n in space, the number of points in the integration block of the .pro file should be larger than m, where m*2-1=n+2, since a N point Gauss quadrature is exact for polynomials with degree 2N-1, and the two bases functions are of order 1 \cite{Numerical Recipies in C, shortly before formula 4.5.9}.
>  
>
The number of Gauss integration points must indeed be fixed in order the 
integral of the total considered polynomial order is correct.

>j) Is more a proposal: in [1] there is the description of a loop gauge that (at least in the paper) seems to produce "nicer" vector-potentials than the edge-gauge one can use from getDP. Maybe this could be an interesting thing to add to the already great capabilites of getDP (e.g. by an additional Basis-Function-Type "BF_EdgesOfLoops").
>  
>
Why not...

>I feel very guilty for posing this bombardement of questions, but there is simply nobody else I know who I can ask. As well, this might be important for people who want to learn to work with getDP, or at least this is what I hope. So thank you in advance, in case your valueable time allows you to read/answer some of this (possibly stupid) questions.
>
>Nice greetings from Vienna,
>Bernhard Kubicek
>  
>
I hope these answers will help you.

Kind regards,

Patrick


>Some literature on this themes:
>
>[1] Ticar, I.; Biro, O.; Preis, K.: Vector potential expanded by edge basis functions associated with loops on finite element facets.  IEEE transactions on magnetics 38 (2002) , S. 437-440 
>
>[2] Vasile Catrinel Gradinaru: Whitney Elements on Sparse Grids, Dissertation, Tübingen 2002. 
>
>[3] J.C. Nedelec: Mixed Finite Elements in R3, Numerische Mathematik 1980. 
>
>[4] R. Albanese, G. Rubinacci: Solution of three dimensional eddy current Problems by integral and differential methods, IEEE Trans. Magnetics, 24(1), 1988. 
>
>[5] Biro, O.; Preis, K.; Richter, K. R.: On the use of the magnetic vector potential in the nodal and edge finite element analysis of 3D magnetostatic fields. IEEE transactions on magnetics 32 (1996) , S. 651-654
>
>_______________________________________________
>getdp mailing list
>getdp at geuz.org
>http://www.geuz.org/mailman/listinfo/getdp
>  
>
-- 
Patrick Dular, Dr. Ir., Research associate, F.N.R.S.
Department of Electrical Engineering and Computer Science
Unit of Applied Electricity
University of Liege - Montefiore Institute - B28 - Parking P32
B-4000 Liege - Belgium - Tel. +32-4 3663710 - Fax +32-4 3662910
E-mail: Patrick.Dular at ulg.ac.be