Construction of FE_Nedelec basis and ordering of DoFs for solving the Maxwell equation

[Ross Kynch]

Hi (again, second question of the day)..

I've been looking into preconditioning of the maxwell system (curl curl E + E = J type) and am wondering how the basis functions are constructed.

In [1], pdf here, they describe a basis set made up of lower order, faces (based on gradients of H1 functions and non-gradients), edges (grad/non-grad) and interior, elements (see page 13, eq (22)). I am wondering if the implementation in FE_Nedelec follows this or uses some order construction?

Also (assuming the answer to the above is yes) in the same paper (see page 16, eq (32)) they propose a block preconditioner which exploits the gradient-based parts of the basis on edges and faces. I'm interested in implementing this preconditioner so I am wondering if there is an easy way to access the specific blocks of the linear system in order to construct it. Would the function DoFRenumbering::block_wise be the right way to go about getting the degrees of freedom ordered for this?

Thanks for any thoughts!

Ross

[1] P. D. Ledger, S. Zaglmayr. hp-Finite Element Simulation of Three-Dimensional Eddy Current Problems on Multiply Connected Domains, NuMa-Report 2010/2,

[Ross Kynch]

Apologies for bumping this up but I'm still in need for an answer to this question. Even if I can only identify which degrees of freedom are associated with the low-order basis and which are associated with higher order it would be useful for preconditioning, although knowing which are associated with edges and faces too would be helpful.

[Wolfgang Bangerth]

> I've been looking into preconditioning of the maxwell system (curl curl E + E
> = J type) and am wondering how the basis functions are constructed.
>
> In [1], pdf here

> <http://www.numerik.math.tu-graz.ac.at/~zaglmayr/pub/Bericht0210.pdf>, they

> describe a basis set made up of lower order, faces (based on gradients of H1
> functions and non-gradients), edges (grad/non-grad) and interior, elements
> (see page 13, eq (22)). I am wondering if the implementation in FE_Nedelec
> follows this or uses some order construction?
>
> Also (assuming the answer to the above is yes) in the same paper (see page 16,
> eq (32)) they propose a block preconditioner which exploits the gradient-based
> parts of the basis on edges and faces. I'm interested in implementing this
> preconditioner so I am wondering if there is an easy way to access the
> specific blocks of the linear system in order to construct it. Would the
> function DoFRenumbering::block_wise be the right way to go about getting the
> degrees of freedom ordered for this?

I am somewhere over the Pacific northeast of Hawaii (where I'd also rather be,
for that matter), so I can't look up your paper. But what I think you're
suggesting is a way to preconditioning based on distinguishing lower and
higher order shape functions of the Nedelec element (or some other distinction
between the shape functions). If I'm right, then this is not an approach we
commonly take in deal.II. In particular, any of the DoFRenumbering functions
we currently have only renumber by the vector component or block a shape
function belongs to. (See the glossary for a lengthy discussion of what
components and blocks are.) In your case, for the Nedelec element, shape
functions correspond to multiple components and all belong to the same block.
So, if the only variable in your set of equations is the electric field, then
renumbering by blocks does nothing at all, because all shape functions belong
to block 0.

That said, the equation you are looking at lead to a positive definite system
matrix. These are typically efficiently preconditioned using multigrid or
algebraic multigrid methods. Why don't you try the ones in Trilinos or PETSc?

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

[Guido Kanschat]

Dear Ross,

the Schöberl/Zaglmayr approach for preconditioning high order methods is surely very neat and to my mind it will eventually be 'the way to go'. We have not considered it so far, though.

The good news first: on each cell, the degrees of freedom are ordered as follows:

1. degrees on vertices, by vertex

2. degrees on edges

3. degrees on faces

4. degrees on volumes

This order is currently reliably hard coded into deal.II and should help somewhat. If you want to implement their concept for sparse matrices, what is needed from my point of view is:

1. DoFRenumbering::by_dimension, which renumbers globally according to this scheme. That would give you the structure needed for (28) in the preprint cited in your email.

2. DoFTools::count_dofs_by_dimension, which would give you the sizes of the blocks in these matrices.

That should be it. One remark though: the Maxwell system is not uniformly (with respect to mesh size) positive definite, so the condition number will deteriorate even for multigrid methods if you do not use a solver for the mixed system (Arnold-Falk-Winther (AFW) 2000 or Hiptmair 1999, where I do not see how to generalize the latter efficiently for higher order). I assume you have a divergence condition like div E = 0? I have been using AFW for a while and it is very efficient. the BlockRelaxation methods in deal.II allow for a crude but working implementation.

Best,

Guido

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
--- You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--
Prof. Dr. Guido Kanschat
Interdisziplinäres Zentrum für Wissenschaftliches Rechnen
Universität Heidelberg
Im Neuenheimer Feld 368, 69120 Heidelberg

[Ross Kynch]

Thanks both for the replies. I'm currently making sure my complex-valued maxwell solver is working correctly and will be implementing the preconditioning once I am sure I can solve the eddy current problem.

Guido, when you say the DoFs are numbered in that order, and they also structured in a hierarchical ordering?

i.e. for FE_Nedelec<3>(p) the basis would be ordered as so:

1. edges at 0th order (12 of them)

2. edges at higher order (12*p of them)

3. faces at 1st order (12*p*(p+1) of them)

4. volumes at 1st order. (3*p*p*(p+1) of them)

With each of 2-4 being in groups, e.g. edges would be in groups of 12 (12 at 1st order, 12 at 2nd order and so on)?

I don't think we require a divE=0 condition as we're regularising the system by considering a small conductivity.

At the moment I'm just testing things out for a wave propagation problem ( E = p*exp(i k * x), with p, k being orthogonal vectors) and am getting some odd results when using Nedelec elements of order other than 0. I have a feeling the project_boundary_values_curl_conforming call may be going wrong on the faces (where it's n \times E for the dirichlet condition rather than n \dot E on the edges), I'll put together a small example once I am sure.

One additional issue is that the computational time for calculating the basis function for orders above 2 or 3 are very very high. I attempt to run a small example at order 4 and it was still generating the basis after an hour or so - could possibly include this in my example as well to see if it's just my install (currently up-to-date with the svn) or not.

Thanks again,

Ross

[Wolfgang Bangerth]

> Guido, when you say the DoFs are numbered in that order, and they also
> structured in a hierarchical ordering?
>
> i.e. for FE_Nedelec<3>(p) the basis would be ordered as so:
>
> 1. edges at 0th order (12 of them)
> 2. edges at higher order (12*p of them)
> 3. faces at 1st order (12*p*(p+1) of them)
> 4. volumes at 1st order. (3*p*p*(p+1) of them)
>
> With each of 2-4 being in groups, e.g. edges would be in groups of 12 (12 at
> 1st order, 12 at 2nd order and so on)?

No. The cell-local ordering first enumerates all shape functions associated
with vertex 0, then all vertex 1, ..., then all associated with edge 0, then
edge 1, ..., then all at face 0, face 1, ...

If you want information about the degree of each shape function, you will need
to ask the finite element for this. In fact, I'm not even sure that there are
0th order edge functions and higher order edge functions -- they may all be
high order (but, of course, spanning the same space). You will have to verify
this.

> At the moment I'm just testing things out for a wave propagation problem ( E =
> p*exp(i k * x), with p, k being orthogonal vectors) and am getting some odd
> results when using Nedelec elements of order other than 0. I have a feeling
> the project_boundary_values_curl_conforming call may be going wrong on the
> faces (where it's n \times E for the dirichlet condition rather than n \dot E
> on the edges), I'll put together a small example once I am sure.

Yes, testcases help.

> One additional issue is that the computational time for calculating the basis
> function for orders above 2 or 3 are very very high. I attempt to run a small
> example at order 4 and it was still generating the basis after an hour or so -
> could possibly include this in my example as well to see if it's just my
> install (currently up-to-date with the svn) or not.

When you say "generating the basis", do you mean "running the constructor of
the FE_Nedelec class"? Yes, again, testcases are the way to go. If you have
independent concerns (e.g., the two issues you mention), create independent
testcases for each and open bug reports (see the link to the bug database on
the website).

[Ross Kynch]

> No. The cell-local ordering first enumerates all shape functions associated 

> with vertex 0, then all vertex 1, ..., then all associated with edge 0, then 
> edge 1, ..., then all at face 0, face 1, ... 

Ok, so for Nedelec elements (no vertex shape functions in this case) it would be possible to reorder by putting the first function associated with each of the 12 edges in the first group of DoFs and then keeping the subsequent ordering after that. However:

> If you want information about the degree of each shape function, you will need 
> to ask the finite element for this. In fact, I'm not even sure that there are 
> 0th order edge functions and higher order edge functions -- they may all be 
> high order (but, of course, spanning the same space). You will have to verify 
> this. 

From what you've said here the basis set is not hierarchical (i.e. there are not shape functions of polynomial order 0, 1, 2, ... p but instead a set of functions all at order p), in which case my idea for preconditioning based on targeting the low-order edge elements is not going to work here. I'll look further into it as you suggest.

> When you say "generating the basis", do you mean "running the constructor of 

> the FE_Nedelec class"? Yes, again, testcases are the way to go. If you have 
> independent concerns (e.g., the two issues you mention), create independent 
> testcases for each and open bug reports (see the link to the bug database on 
> the website). 

Yes, the constructor takes 1000+ seconds for order 4, even in release mode. I've submitted a bug report.

Thanks again for the help.

Ross

[Wolfgang Bangerth]

> Ok, so for Nedelec elements (no vertex shape functions in this case) it would
> be possible to reorder by putting the first function associated with each of
> the 12 edges in the first group of DoFs and then keeping the subsequent
> ordering after that. However:

Yes, if the basis is hierarchical, then this will work.

> From what you've said here the basis set is not hierarchical (i.e. there are
> not shape functions of polynomial order 0, 1, 2, ... p but instead a set of
> functions all at order p), in which case my idea for preconditioning based on
> targeting the low-order edge elements is not going to work here. I'll look
> further into it as you suggest.

Yes. The answer to the first sentence of the paragraph is "I don't know". The
basis may be hierarchical or not, I simply don't know the answer. But you
should make sure it is what you expect, in any case.

> Yes, the constructor takes 1000+ seconds for order 4, even in release mode.
> I've submitted a bug report.

Thanks