Problem with Dirichlet BCs for FE_Nedelec for 3D?

[Ross Kynch]

Hi all,

I'm solving a complex-valued vector-wave equation (the curl-curl formulation of the Maxwell equations):

curl(curl(E)) + kappa*E=0

with BCs given by

n x E = g on \Gamma_D (Dirichlet)

n x curl(E) = h on \Gamma_N (Neumann)

note, I've set mu =1 for this problem.

the wave propagation solution is given by:

 E = p*exp(i*k*x),  x, k and p in R^3

 

with p orthogonal to k, |k| < 1, |p|=1.

I've chosen k=(-0.1 -0.1, 0.2) and p=(1/sqrt(3))*[1 1 1].

I've implemented this for both Neumann and Dirichet BCs using Nedelec elements and splitting the system into real/imaginary parts. I am seeing the correct rates of convergence in  the L2 and H(curl) norms of the error for the Neumann version with both h and p refinement.

However, when using Dirichlet BCs, after p=0, something seems to go wrong. Am I enforcing the boundary conditions correctly for the project_boundary_values_curl_conforming routine? - I have a function with dim+dim components, with the first dim(3) being the real part and the latter dim being the imaginary part - is this the correct way to do this?

I've included an example with Dirichlet BCs, although if you uncomment the 666-678 in the code then the boundaries will be changed to Neumann. It takes the polynomial order as input, so run it with "./wave_propagation -p 2" for order 2 elements, etc.

For comparison, here are the results when I've run the code p=0,1,2 - using a direct solver for now, so going above p=2 causes memory problems, also runtime increases due to the slow generation of the Nedelec elements at higher orders.

Thanks,

Ross

Dirichlet:

./wave_propagation -p 0

cycle cells dofs    L2 Error    H(curl) Error  

    0     8  108 1.15904018e-01 1.24176508e-01 

    1    64  600 5.77896943e-02 6.19578648e-02 

    2   512 3888 2.88743606e-02 3.09624645e-02 

./wave_propagation -p 1

cycle cells dofs     L2 Error    H(curl) Error  

    0     8   600 8.35079734e-01 2.59980668e+00 

    1    64  3888 7.93792725e-01 3.36290438e+00 

    2   512 27744 7.77603612e-01 4.42928910e+00 

./wave_propagation -p 2

cycle cells dofs     L2 Error    H(curl) Error  

    0     8  1764 8.20717491e-01 2.50234579e+00 

    1    64 12168 7.87196814e-01 3.17447088e+00 

    2   512 90000 7.74601489e-01 4.11893779e+00 

Neumann:

./wave_propagation -p 0

    0     8  108 1.15471804e-01 1.23771671e-01 

    1    64  600 5.77352583e-02 6.19069021e-02 

    2   512 3888 2.88675425e-02 3.09560822e-02 

./wave_propagation -p 1

cycle cells dofs     L2 Error    H(curl) Error  

    0     8   600 2.58268149e-03 2.79282760e-03 

    1    64  3888 6.45549584e-04 6.98219104e-04 

    2   512 27744 1.61377755e-04 1.74553607e-04 

./wave_propagation -p 2

cycle cells dofs     L2 Error    H(curl) Error  

    0     8  1764 4.17828432e-05 4.53286783e-05 

    1    64 12168 5.22313050e-06 5.66700077e-06 

    2   512 90000 6.52897524e-07 7.08401532e-07 

[Guido Kanschat]

Dear Ross,

First question that comes to mind: did you check the error on the Dirichlet boundary? Is it zero or converging with the right order?

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+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Ross Kynch]

I've not explicitly calculated the error on the boundaries, but in all cases with Dirichlet conditions it does not look right and is getting worse with increasing p. This is just using paraview to visualise the solution on the surface of the cube and comparing it to how I expect the solution to look.

I've also noticed a small mistaken in my code - the off-diagonal entries (multipled by kappa_im) should be values[E_re]...*values[E_im]... and vice versa, not real*real and imaginary*imaginary as they are written in the code. This doesn't effect the solution in this case as the value of kappa_im is zero for this problem, but it's important if someone were to use the code as a basis for anything more complicated.

I'll look into calculating the surface errors, but I'm fairly sure they're going to confirm there is divergence from the solution on the boundary.

Ross

[Ross Kynch]

One other thing:

Can you confirm that the way I've implemented the function for the boundary values is correct? (i.e. that I've got 6 components instead of 3, where the first 3 are real and the second 3 are imaginary parts of the solution, E).

Thanks

Ross

[Guido Kanschat]

Dear Ross,

If you see and are sure that the tangential component of the solution is wrong, we may not need the error calculation.

Have you checked whether the function that computes the projected values behaves significant different depending on whether it receives one or two Nedelec elements? You may have found a bug here.

Best,
Guido

Best,

[Ross Kynch]

Do you mean in the constructor of FESystem?

i.e. using

fe (FE_Nedelec<dim>(order), 2)

or

fe (FE_Nedelec<dim>(order), 1, FE_Nedelec<dim>(order), 1)

If so - this makes no difference. If you mean reducing the complex-valued problem to a real valued problem so that only a single Nedelec element is required, then yes, I have tried a simple toy problem with non zero Dirichetl BCs - it too had issues with orders > 0. This was a while ago - unfortunately I got sidetracked with some other work at the time. I'll have to dig out the code and tidy it up if you'd like that example as well?

My colleague (supervisor on the EPRSC project I'm a PDRA on) has suggested that the handling of edges and faces in 3D is different to the 2D problem and he suspects there's a bug in the handling of the enforcement of the tangential components too. I'm unsure of where to look in the source code - my C++ knowledge isn't good enough to fully debug the problem.

Cheers

Ross

[Wolfgang Bangerth]

On 07/25/2014 05:10 AM, Ross Kynch wrote:
> Do you mean in the constructor of FESystem?
>
> i.e. using
> fe (FE_Nedelec<dim>(order), 2)
> or
>
> fe (FE_Nedelec<dim>(order), 1, FE_Nedelec<dim>(order), 1)

This should make no difference.

> If so - this makes no difference. If you mean reducing the complex-valued
> problem to a real valued problem so that only a single Nedelec element is
> required, then yes, I have tried a simple toy problem with non zero Dirichetl
> BCs - it too had issues with orders > 0. This was a while ago - unfortunately
> I got sidetracked with some other work at the time. I'll have to dig out the
> code and tidy it up if you'd like that example as well?
>
> My colleague (supervisor on the EPRSC project I'm a PDRA on) has suggested
> that the handling of edges and faces in 3D is different to the 2D problem and
> he suspects there's a bug in the handling of the enforcement of the tangential
> components too. I'm unsure of where to look in the source code - my C++
> knowledge isn't good enough to fully debug the problem.

It is possible that project_boundary_values_curl_conforming() has a bug for
p>0. To find out one way or the other, here is the way I usually attack these
issues:

First, make the problem as simple as possible. In your case, go to a problem
with only one, not two elements. If the boundary values are still wrong,
simplify the mesh -- maybe to a single cell. Are the boundary values still
wrong? Can you verify that this is truly the case and not just an effect of
visualization? For this, see for example

https://code.google.com/p/dealii/wiki/FrequentlyAskedQuestions#In_my_graphical_output,_the_solution_appears_discontinuous_at_ha

If the output definitely looks wrong, since you're on a single cell you can
probably compute the exact solution of your problem on a piece of paper. Does
your numerical solution differ? Is this true not just for the final solution,
but also for the constraints the project_boundary_values_curl_conforming()
function puts into the constraints object?

If still yes, then you know that there is a bug in the function that computes
the boundary values. Now that you know the exact solution and you have a small
testcase, the next step is to fix it. The point of all this simply is that
it's very difficult to find these bugs with a program of 726 lines. You want
to test the bug on a program that only calls the function in question on a
single cell.

Best
W.


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

[Ross Kynch]

Hi again,

Firstly, thanks to both of you for the advice on debugging - it proved useful in finding the problem here:

I've now done some testing and found that the problem only occurs when I use FESystem, even if only using a single FENedelec within it.

I've attached 2 codes which solve (in 3D):

curl(curl(E)) + E = Js, with E = (x^2, y^2, z^2) and Js=E.

As is obvious from the names: one uses FESystem which then contains and the other uses FENedelec on its own. The rest of the code should be identical between the two. So using the nedelec element on its own, all is fine and there's the proper convergence in L2 & H(curl) norms, with the exact solution found at p=2 as expected. Switch to the FESystem type and it gives the same error for p-0, but goes wrong for p>0. Also, if I switch the boundary conditions to Neumann then both work fine (obviously adding the appropriate part to the assemble routine) and I see the correct convergence rates. Run using:

./using_fe_xxx -p 1

for order 1 and so on.

I can only assume that the bug lies within FESystem somehow or that I've misunderstood how to use FESystem.

Cheers

Ross

[Alexander]

Hi Ross,

the bug is in the project_boundary_values_curl_conforming and related to the wrong ordering of the face DoFs (which is relevant for p > 0 only). 

I can try and fix it when I have time (probably next week) if nobody will do it before.

Best,

Alexander

[Ross Kynch]

Hi Alexander,

Thanks for finding it. Please let me know when it's fixed!

I'm also having some issues with non-cubic meshes (i.e. the elements are not bricks but general hex) for p>0, so am hoping this bug will fix that too!

Regards,

Ross

[Alexander]

Ross,

>I'm also having some issues with non-cubic meshes (i.e. the elements are not bricks but 

>general hex) for p>0, so am hoping this bug will fix that too!

have you read the main doc page of FE_Nedelec? It explains why FE_Nedelec does not support general hex elements. Your problem can be related to this issue rather. 

Alexander

[Ross Kynch]

Thanks for point that out, I assume this is the section you refer to?:

Todo:

Even if this element is implemented for two and three space dimensions, the definition of the node values relies on consistently oriented faces in 3D. Therefore, care should be taken on complicated meshes.

Restriction on transformations

In some sense, the implementation of this element is not complete, but you will rarely notice. Here is the fact: since the element is vector-valued already on the unit cell, the Jacobian matrix (or its inverse) is needed already to generate the values of the shape functions on the cells in real space. This is in contrast to most other elements, where you only need the Jacobian for the gradients. Thus, to generate the gradients of Nédélec shape functions, one would need to have the derivatives of the inverse of the Jacobian matrix.

Basically, the Nédélec shape functions can be understood as the gradients of scalar shape functions on the real cell. They are thus the inverse Jacobian matrix times the gradients of scalar shape functions on the unit cell. The gradient of Nédélec shape functions is then, by the product rule, the sum of first the derivative (with respect to true coordinates) of the inverse Jacobian times the gradient (in unit coordinates) of the scalar shape function, plus second the inverse Jacobian times the derivative (in true coordinates) of the gradient (in unit coordinates) of the scalar shape functions. Note that each of the derivatives in true coordinates can be expressed as inverse Jacobian times gradient in unit coordinates.

The problem is the derivative of the inverse Jacobian. This rank-3 tensor can actually be computed (and we did so in very early versions of the library), but is a large task and very time consuming, so we dropped it. Since it is not available, we simply drop this first term.

What this means for the present case: first the computation of gradients of Nédélec shape functions is wrong in general. Second, in the following two cases you will not notice this:

• If the cell is a parallelogram, then the usual bi-/trilinear mapping is in fact affine. In that case, the gradient of the Jacobian vanishes and the gradient of the shape functions is computed exactly, since the first term is zero.
• With the Nédélec elements, you will usually want to compute the curl, not the general derivative tensor. However, the curl of the Jacobian vanishes, so for the curl of shape functions the first term is irrelevant, and the curl will always be computed correctly even on cells that are not parallelograms.

However, my equations don't involve gradients. My aim is to solve curl mu^-1 curl E - kE = J where E and k are both complex valued. Surely I should be ok, even with a general cell?

The issue I'm having is with a mesh of a sphere inside a cube generated in trelis (previously cubit) - I'm attempting to compute the eddy currents within the sphere by imposing the analytical solution for the magnetic field (i.e. using neumann conditions) at a distance from the sphere. Is it likely that the face orientations are to blame (as stated in the TODO at the top)?

If that's the case - any suggestions on how to assure that the faces remain consistently orientated?

Thanks again for the help

Ross

[Alexander]

Ross,

you're likely facing the issue with non-standard faces. I do not know an easy way to overcome this problem. Adding support for this in the library is probably not trivial, especially for high-order elements, but you can try. 

Alexander

[Ross Kynch]

I'll discuss it with my colleague, who has a 3D time-harmonic Maxwell code using nedelec elements. His code is for tetrahedral elements, but I suppose the orientation issues are similar. I'm happy to contribute to the library and improve things for the nedelec element (e.g. I'd like to implement the hierarchal basis functions of Zaglmayr & Schöberl so that more efficient preconditioners can be used) but I'm unsure if my programming skills are good enough. My main problem is that I find some of the library's code quite hard to follow. My experience in C++ is fairly limited - I did all of my PhD work using F90 so I'm learning a lot of things as I go with deal.II. It'd be helpful if someone could chat through things on skype or google hangouts, although I'm not sure who is the best person to ask?

Ross

[Guido Kanschat]

Hi Ross,

first I must confess, that I am not fully aware of the details of face orientations or the implementation of the FE_Nedelec. But so much is for sure (3D):

1. Degrees of freedom on edges

  a) If the face of the other cell is flipped, then tangentials on two edges point in opposite directions.

  b) If a face is rotated, so are the degrees of freedom.

2. Degrees of freedom on faces (higher order): according to the documentation, they are defined as moments with respect to Nedelec polynomials of one order and one dimension lower. So, I guess we need to look at all those basis functions and see into which they transform under rotation and mirroring, including signs. Hopefully, there is some system behind this for higher order.

Comparison to other elements:

-FE_Q only has the reordering of degrees of freedom on the faces, since they are just scalars. But for higher order, we should be able to learn where exactly the reordering is done (I guess in get_dof_indices)

-FE_RaviartThomas also should have the change of sign, if for instance two faces with the same index on each cell, respectively, are adjacent. Is it verified that this works?

By the way, the Zaglmayr/Schöberl technique is definitely recommendable. Since our degrees of freedom are already ordered vertex/line/quad/hex, does that imply we are almost there?

Best,
Guido

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

[Ross Kynch]

With regard to the orientations:

One method for dealing with the orientations in tetrahedra is to characterise elements to be of two types - the way the basis functions are set up then depends on which type the element is, and then there a few combinations of how the elements can fit together. As you said, it seems this isn't so easy for hex and I can see why it's been avoided.

With regard to the basis functions:

It seems that your DoFs would fit nicely with the technique and should be fairly easy to implement. The only tricky part would be the restriction/prolongation operators - I think this is the part which takes all of the time to set up - something I mentioned previously about the high runtime for constructing a high-order nedelec element.

Anyway, I've been testing things with the current setup using the simple quadratic solution (real-valued using FE_Nedelec, not using FESystem), above, with p=2 (i.e. solution should be exact near to machine precision):

I've been experimenting with as few number of elements as possible and I'm finding that Dirichlet boundary conditions (implemented via project_boundary_values_curl_conforming()) do not work at all for elements which arn't cube-like (even fails on a parallelepiped). Neumann conditions (implemented by computing the integral of n \cross curl(E)) works fine for cubes, cuboids and a general hex with linear faces (tested by moving the vertices of the cube away). Neumann BCs even work fine on a cylinder with the correct cylindrical boundary, but does not work for a sphere (with a hyper ball boundary).

So some questions:

Is this difference only due to the bug which Alexander described?

On which shaped elements/boundaries generated within deal.II should I expect deal's Nedelec elements to work? - I need to be able to handle cylinders and spheres within my geometry.

How are the curls of the basis computed?:

Would I be right in assuming it's done as follows: curl N = 1/(detJ)*J*( curl_hat N_hat) ?

where curl and N are the physical curl and basis function and curl_hat and N_hat are the curl and basis function on the reference element.

I ask this because do_function_curls() in fe_values.cc doesn't seem to compute them this way, unless the function is computing the (curl_hat N_hat) and the jacobian is applied later. I can't see where this is done, if it is.

[Alexander]

Ross,

you're talking about many different things in your reply. I suggest we concentrate on boundaries here and issues related to face orientations or curl computations can be discussed separately. This increases chances to get answers. 

The issue with FESystem was addressed here: https://github.com/dealii/dealii/pull/82

Could you check that it works for you? Especially for the case of FESystem with two elements in it. 

From what I can see, the project_boundary_values_curl_conforming seems assuming that the faces on boundaries are rectangles, but I'm not sure actually. Do you have non-rectangular faces? And what it means it does not work? 

[Wolfgang Bangerth]

Ross,

> One method for dealing with the orientations in tetrahedra is to characterise
> elements to be of two types - the way the basis functions are set up then
> depends on which type the element is, and then there a few combinations of how
> the elements can fit together. As you said, it seems this isn't so easy for
> hex and I can see why it's been avoided.

The issue is more subtle than that. You can avoid these problems in 2d
altogether, see
http://www.dealii.org/developer/doxygen/deal.II/classGridReordering.html
but in 3d there is no way around the fact that one has to store orientations
on edges and faces of cells. I spent an enormous amount of time on this about
10 years ago, but I have since then found it incredibly hard to deal with
these issues every time they come up. My brain just doesn't seem to process
these geometric problems adequately (they only appear on nontrivial 3d
geometries) and this seems to be the case also for everyone else who tried.

So my guess is simply that these cases aren't correctly taken into account in
the implementation of the Nedelec element. The simplest case one could try
would be a torus of cubes that are twisted by 90 or 180 degrees, but even that
is nontrivial. One problem is that, just like for the geometry issues outlined
above, we as a project seem to have lost the expertise necessary to deal with
the higher order Nedelec elements. Alexander seems to be reading up on how it
works, and maybe Markus can still help, but other than that nobody knows
enough of the details of the element when using higher orders -- the perils of
a small project :-(

> On which shaped elements/boundaries generated within deal.II should I expect
> deal's Nedelec elements to work? - I need to be able to handle cylinders and
> spheres within my geometry.

Cylinders can be meshed without face or edge flips if you just extrude a 2d
mesh. But spheres can't.

In both cases, you will have cells that aren't parallelograms/parallel
epipeds, so the issues in the note you quoted do apply.

> How are the curls of the basis computed?:
> Would I be right in assuming it's done as follows: curl N = 1/(detJ)*J*(
> curl_hat N_hat) ?

No. The curl consists of some of the components of the gradient tensor and
that's how it is implemented: compute the gradient and pick individual elements.

[Ross Kynch]

Alexander:

Apologies, I am having a few (possibly connected) issues so am probably confusing things together. I agree that the issues seem to be boundaries at the moment.

> The issue with FESystem was addressed here: https://github.com/dealii/dealii/pull/82

> Could you check that it works for you? Especially for the case of FESystem with two elements in it. 

Thanks for this. I have tried out your patch and it appears to fix the problems with FESystem, at least for some simple test cases on cubes.

> From what I can see, the project_boundary_values_curl_conforming seems assuming that the faces on boundaries are rectangles, but I'm not sure actually. Do you have non-rectangular faces? And what it means it does not work?

When I say things don't work, I mean that I don't see the expected convergence in the L2 or the H(curl) norms. For the tests I've written I've used a quadratic solution, so I expect to recover the exact solution when using order 2 elements. 

I think you're probably right about the assumption of rectangular faces. For example, I have tried solving for both Dirichlet and Neumann conditions on a parallelepiped. Currently I see the expected convergence for Neumann (n x curl(E) type). However, if I switch to Dirichlet conditions via project_boundary_values_curl_conforming, then the convergence breaks down. I've attached an example similar to the previous ones (uncomment lines 463-476 for Neumann conditions) which shows that things break down on a parallelpiped (i.e non-rectangular faces) for Dirichlet conditions when using p>0. Note the error increase in the H(curl) norm.

Wolfgang:

I understand the issues with orientations in 3D - have the same problems myself.. everything is a lot easier to visualise in 2D! My background is in fluid dynamics and spectral element methods so I am fairly new to Nedelec elements and curl operators in the context of FEM, but I'm happy to help improve things as best I can.

My question about the curl computation was that I was thinking that the issue with needing the derivatives of the inverse Jacobian could be avoided by computing the curl in this manner - you no longer need to take the derivatives outside of the reference element. Maybe this would not help in either case.

Regards,

Ross

[Wolfgang Bangerth]

On 08/18/2014 11:22 AM, Ross Kynch wrote:
> My question about the curl computation was that I was thinking that the issue
> with needing the derivatives of the inverse Jacobian could be avoided by
> computing the curl in this manner - you no longer need to take the derivatives
> outside of the reference element. Maybe this would not help in either case.

It's possible, but this isn't how it's implemented right now. The code paths
to get at the curl are, at least currently, the same as those getting at the
gradients of shape functions. It would no doubt be a significant amount of
work to add a different way to compute them.

As for your test program: If you can verify that you can reproduce a linear
solution with p=1 elements but can't reproduce a quadratic solution with p=2
elements, then there must indeed be a bug somewhere. Please open a bug report
at the issues tracker in that case. I don't think we can promise to fix this
anytime soon, but at least the issue is recorded and documented.

Best
Wolfgang

[Alexander]

Wolfgang,

I would be interested actually adding the inverse Jacobian in FE_Nedelec? I would look closer at it if someone could give a couple of guidelines. Do other non-trivial FE implement this (RaviartThomas)? Can you foresee any caveats? 

Thanks,

Alexander

[Martin Kronbichler]

Alexander,

you mean the derivative of the inverse Jacobian? The computation of that rank-3 tensor, given the derivative of the Jacobian transformation (which you get by fe_values. jacobian_grad), is implemented in the evaluation routines of the matrix-free class. The computation of the inverse Jacobian is found in include/deal.II/matrix_free/mapping_info.templates.h, lines 450-497. Note that in that class one uses tensors of vectorized arrays instead of 'usual' tensors, but the computations actually work the same. The important thing to take care of is the order of indices (also when later applying the inverse Jacobian). Also, the algorithm needs modifications in case dim != spacedim.

Best,

Martin

[Wolfgang Bangerth]

Alexander,

> I would be interested actually adding the inverse Jacobian in
> FE_Nedelec? I would look closer at it if someone could give a couple of
> guidelines. Do other non-trivial FE implement this (RaviartThomas)? Can
> you foresee any caveats?

as you probably have noticed already, many of the issues surrounding
finite element implementations and the interface between FEValues and
the finite element class aren't completely familiar to anyone in the
project anymore. You may therefore have to dig around a bit to find
answers if you can't get them on the mailing list.

About your concrete question: the inverse of the Jacobian is something
that the mapping, not the finite element computes. If the element needs
it, it needs to make this clear to the mapping -- this happens at the
time an FEValues class is setup: it asks the finite element what it
needs from the mapping to compute what the update_flags say should be
computed. In this case, you need to make sure that the FE_Nedelec tells
the mapping that it needs the inverse of the Jacobian.

I suspect you actually meant the derivative of the inverse of the
Jacobian. Like Martin already pointed out, there is code that can do
this elsewhere already. I do not recall whether that is the current
status, but there used to be a time when we computed the derivative of
the Jacobian or its inverse via finite differences. This is because we
were young and stupid at the time because one can use the formula

grad (J^{-1}) = -J^{-1} (grad J) J^{-1}

(with an appropriate order of indices over which to contract). This
formula results from

grad (I) = grad (J J^{-1}) = 0

after applying the product rule on the middle expression. In other
words, if you need the derivative of the inverse of the Jacobian, all
you need is the derivative of the Jacobian and the inverse of the
Jacobian. Both are already available.

As I said, I don't recall whether that's actually implemented or if we
are still at using finite differences. If you want to go down this route
to extending FE_Nedelec, it would be useful to do it in stages and
submit a number of smaller, self-contained patches that can be
individually reviewed and applied.

[Jianan Zhang]

Dear Ross,

I am using dealii to solve the curl-curl equation:

curl(mu^(-1)curl(E)) + (-omega^2*epsilon+j*omega*sigma)*E=0,

with boundary conditions: n x E = n X F ,                          (Dirichlet)

                                    or  n X (curl(E)) = n X (curl(F))       (Neumann)

where F is a known exact solution.

Based on your code, I implemented my own exact solution and the assembly of system of equations. My exact solution (including real and imaginary parts) is as follows:

E_re = [0, sin(pi*x)/a*cos(k_z10_re*z)*exp(k_z10_im*z), 0],

E_im = [0, -sin(pi*x)/a*sin(k_z10_re*z)*exp(k_z10_im*z), 0]

And the exact solution for your case (if I take p=[0 1 0] and k=[0 0 0.1]) is given by:

E_re = [0, cos(0.1*z), 0],

E_im = [0, sin(0.1*z), 0] 

The domain I want to solve this problem is a 3D hyperrectangle  with width = 19.05e-3, height = 9.524e-3, and length = 15.24e-3. I've tried to solve the above equation with the two exact solutions. When I used the build-in mesh generator in dealii (i.e., GridGenerator::hyper_rectangle (triangulation, p1, p2);) to generate the mesh. The problem can be solved correctly for both cases:

my exact solution, built-in mesh generator:

your exact solution, built-in mesh generator:

Cycle 0:Number of degrees of freedom: 3888

 Done

cycle cells dofs    L2 Error    H(curl) Error  

    0    64 3888 9.73783895e-12 1.03298196e-11 

When I use the mesh generated from Gmsh software (after conversion into hexes), if I provide your exact solution, the problem can still be solved correctly, as shown bellow:

Cycle 0:number of cells:244

Number of degrees of freedom: 13176

 Done

cycle cells dofs     L2 Error    H(curl) Error  

    0   244 13176 1.36020912e-07 8.58355818e-05

However, when I switch to my exact solution, the results are totally wrong. The electric field obtained has totally wrong direction (seems to be arbitrary but the correct one should be along y-axis). The results are as follows

    

My question is like this, based on your experience of using FE_Nedelec, can we use FE_Nedelec  on a mesh imported from other mesh generators like Gmsh instead of the one generated from dealii itself? Does dealii of the current version has a good support for this kind of mesh? FYI, the following is the mesh generated from Gmsh (after transformation).

Any suggestions and help is much appreciated

Thanks in advance.

Jianan Zhang

在 2014年7月15日星期二 UTC-4下午12:58:46，Ross Kynch写道：

Hi all,

I'm solving a complex-valued vector-wave equation (the curl-curl formulation of the Maxwell equations):

curl(curl(E)) + kappa*E=0

with BCs given by

n x E = g on \Gamma_D (Dirichlet)

n x curl(E) = h on \Gamma_N (Neumann)

note, I've set mu =1 for this problem.

the wave propagation solution is given by:

 E = p*exp(i*k*x),  x, k and p in R^3

 

with p orthogonal to k, |k| < 1, |p|=1.

I've chosen k=(-0.1 -0.1, 0.2) and p=(1/sqrt(3))*[1 1 1].

I've implemented this for both Neumann and Dirichet BCs using Nedelec elements and splitting the system into real/imaginary parts. I am seeing the correct rates of convergence in  the L2 and H(curl) norms of the error for the Neumann version with both h and p refinement.

However, when using Dirichlet BCs, after p=0, something seems to go wrong. Am I enforcing the boundary conditions correctly for the project_boundary_values_curl_conforming routine? - I have a function with dim+dim components, with the first dim(3) being the real part and the latter dim being the imaginary part - is this the correct way to do this?

I've included an example with Dirichlet BCs, although if you uncomment the 666-678 in the code then the boundaries will be changed to Neumann. It takes the polynomial order as input, so run it with "./wave_propagation -p 2" for order 2 elements, etc.

For comparison, here are the results when I've run the code p=0,1,2 - using a direct solver for now, so going above p=2 causes memory problems, also runtime increases due to the slow generation of the Nedelec elements at higher orders.

Thanks,

Ross

Dirichlet:

./wave_propagation -p 0

cycle cells dofs    L2 Error    H(curl) Error  

    0     8  108 1.15904018e-01 1.24176508e-01 

    1    64  600 5.77896943e-02 6.19578648e-02 

    2   512 3888 2.88743606e-02 3.09624645e-02 

./wave_propagation -p 1

cycle cells dofs     L2 Error    H(curl) Error  

    0     8   600 8.35079734e-01 2.59980668e+00 

    1    64  3888 7.93792725e-01 3.36290438e+00 

    2   512 27744 7.77603612e-01 4.42928910e+00 

./wave_propagation -p 2

cycle cells dofs     L2 Error    H(curl) Error  

    0     8  1764 8.20717491e-01 2.50234579e+00 

    1    64 12168 7.87196814e-01 3.17447088e+00 

    2   512 90000 7.74601489e-01 4.11893779e+00 

Neumann:

./wave_propagation -p 0

    0     8  108 1.15471804e-01 1.23771671e-01 

    1    64  600 5.77352583e-02 6.19069021e-02 

    2   512 3888 2.88675425e-02 3.09560822e-02 

./wave_propagation -p 1

cycle cells dofs     L2 Error    H(curl) Error  

    0     8   600 2.58268149e-03 2.79282760e-03 

    1    64  3888 6.45549584e-04 6.98219104e-04 

    2   512 27744 1.61377755e-04 1.74553607e-04 

./wave_propagation -p 2

cycle cells dofs     L2 Error    H(curl) Error  

    0     8  1764 4.17828432e-05 4.53286783e-05 

    1    64 12168 5.22313050e-06 5.66700077e-06 

    2   512 90000 6.52897524e-07 7.08401532e-07