interpolate FE_Nelelec

[John Smith]

Hello,

It seems I am unable to find a function similar to VectorTools::interpolate for FE_Nedelec finite elements. The existing implementation of this functions gives the following run-time error:

An error occurred in line <556> of file </home/john/dealii-9.1.1/include/deal.II

/numerics/vector_tools.templates.h> in function

void dealii::VectorTools::interpolate(const dealii::Mapping<dim, spacedim>&,

const DoFHandlerType<dim, spacedim>&, const dealii::Function<spacedim, typename

VectorType::value_type>&, VectorType&, const dealii::ComponentMask&) [with int

dim = 3; int spacedim = 3; VectorType = dealii::Vector<double>; DoFHandlerType =

dealii::DoFHandler; typename VectorType::value_type = double]

The violated condition was:

dof_handler.get_fe().n_components() == function.n_components

Additional information:

Dimension 3 not equal to 1.

Looks like VectorTools::interpolate does not like the vector-valued finite elements. Could you please point me in the right direction?

Best,

John

[simon...@gmail.com]

Hi,
I don't have any experience with FE_Nedelec, but the error message states that the function you are trying to interpolate has 1 component, while the element has dim components.

If you implemented your own Function, you need to make sure it has dim-components by calling the constructor of Function with dim as inargument:

template <int dim>
class CustomFunction : public Function<dim>
{
public:
CustomFunction()
: Function<dim>(dim) // The function has dim components
{}

double
value(const Point<dim> & point,
const unsigned int component = 0) const override
{
// Return value will depend on component here
return 0;
}
};


Best,
Simon

[John Smith]

Dear Simon,

Thank you for your reply. Your suggestion definitely works. I used functions that output one component when I was interpolating scalar potentials. It works well.

However, the FE_Nedelec is different. It implements edge elements. They are vector-based. That is, functions are represented by a superposition of vector-valued shape functions:

\vec{A} = \sum u_i vec{N}_i .

Therefore, the output of the "value" method in “class CustomFunction : public Function” must be vector-valued. Three components in a three-dimensional space.  Otherwise, there is no point in interpolation.

This kind of vectorial approximations is a bread-and-butter topic in magnetics. See, for example, equation (7) in:

https://ieeexplore.ieee.org/document/497322

In short, the source, i.e the current vector potential, must be projected on the space spanned by the vector-valued shape functions. Otherwise, the simulation is numerically unstable. The last, from my experience, is definitely true.

Best,

John

--
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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/wG3j4yD1rpw/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/6d051c9e-8840-4e5d-8021-f8007a1dc543n%40googlegroups.com.

[Jean-Paul Pelteret]

Dear John,

I’m not sure that there is an interpolation function that would work in that way for Nedelec elements (there is VectorTools::project_boundary_values_curl_conforming_l2() for Neumann boundary conditions, not that this is particularly useful for you here). As always we’d be happy to accept a contribution to the library for the extension that you require!

You don’t mention specifically what you’re trying to model, or where this interpolated function is to be used, so I’m just going to throw out a couple of things for you to think about and decide if they’re valid or of use to your specific application. What I’ve done in the past is to use a dim-component Function for the free current source, and I’d written a check (which I’ve pasted below) to verify that I’d not set up this source function incorrectly. Admittedly, I had a relatively straight-forward geometry so it wasn’t challenging to write such a source function. If you’ve got a more complex setup with the current source derived from some secondary finite element problem then maybe such an approach is not appropriate. But if it so happens that you can define or solve for a scalar potential whose gradient is the source current, then this is divergence free if the scalar potential is approximated using Lagrange polynomials (i.e. FE_Q . You can then use an FEFieldFunction to compute this gradient at each point within the current source, to act as a source term for the magnetic vector potential problem.

Best,

Jean-Paul

template <int dim>

void WireMagneticField<dim>::verify_source_terms () const

{

TimerOutput::Scope timer_scope (computing_timer, "Verify source terms");

hp::FEFaceValues<dim> hp_fe_face_values (mapping_collection,

fe_collection,

qf_collection_face,

update_quadrature_points |

update_normal_vectors |

update_JxW_values);

// To verify that the source terms are divergence free,

// we make use of the divergence theorem:

// 0 = \int_{\Omega_{e}} div (J_f)

// = \int_{\partial\Omega_{e}} J_f . n

typename hp::DoFHandler<dim>::active_cell_iterator

cell = hp_dof_handler.begin_active(),

endc = hp_dof_handler.end();

for (; cell!=endc; ++cell)

{

if (cell->subdomain_id() != this_mpi_process) continue;

double div_J_f = 0.0;

for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)

{

hp_fe_face_values.reinit(cell,face);

const FEFaceValues<dim> &fe_face_values = hp_fe_face_values.get_present_fe_values();

const unsigned int &n_fq_points = fe_face_values.n_quadrature_points;

std::vector<Vector<double> > source_values (n_fq_points, Vector<double>(dim));

function_free_current_density.vector_value_list (fe_face_values.get_quadrature_points(),

source_values);

for (unsigned int fq_point=0; fq_point<n_fq_points; ++fq_point)

{

const Tensor<1,dim> J_f ({source_values[fq_point][0],

source_values[fq_point][1],

source_values[fq_point][2]}); // Note: J_f must be divergence free!

const double JxW = fe_face_values.JxW(fq_point);

const Tensor<1,dim> &N = fe_face_values.normal_vector(fq_point);

div_J_f += (J_f*N) * JxW;

}

}

AssertThrow(std::abs(div_J_f) < 1e-9, ExcMessage("Source term is not divergence free!"));

}

}

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.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/CAH3L9AVtN%3Dtaqi2YiTjfoyfhxo-%2B9XnNzpBqNXSHqk5wX4iUtQ%40mail.gmail.com.

[John Smith]

Ah! Now I got it. The function you suggested indeed vector-valued. It works just fine now. Thanks a lot! 

I was trying to output the function values via “vector_value” method like this:

class CustomFunction : public Function<3>

{
public:
virtual void vector_value( const Point<3> & p,
            Vector< double > & values ) const override;
};

That was a bad idea ...

Thanks! Now it works! Problem solved!

Best,
John.

[John Smith]

Dear Jean-Paul,

Thank you for your reply. I am solving a very standard curl(1/mu(curl(A)) = J problem as described in the paper of Oszkar Biro.

https://ieeexplore.ieee.org/document/497322

No gauge. I need interpolation to project the vector current potential on the shape functions.

I just tried as a test to project the magnetic field of a magnetic dipole on a cuboidal domain using the dim- component Function following example of Simon and VectorTools::interpolate. It looks ok … But I will be extra careful checking the projections made by this interpolation function.

Thanks!
John 

[Wolfgang Bangerth]

On 4/16/21 11:36 AM, John Smith wrote:
> However, the FE_Nedelec is different. It implements edge elements. They are
> vector-based. That is, functions are represented by a superposition of
> vector-valued shape functions:
>
> \vec{A} = \sum u_i vec{N}_i .
>
> Therefore, the output of the "value" method in “class CustomFunction : public
> Function” must be vector-valued. Three components in a three-dimensional
> space.  Otherwise, there is no point in interpolation.
>
> This kind of vectorial approximations is a bread-and-butter topic in
> magnetics. See, for example, equation (7) in:
>
> https://ieeexplore.ieee.org/document/497322

> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fieeexplore.ieee.org%2Fdocument%2F497322&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cf7f578c8c6fc44890c3308d900fe30ef%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C1%7C637541914012798794%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=tJVhiTO272UFKlS0mWzTCrVp1VzPV3zbr2I16KxHkFI%3D&reserved=0>

>
> In short, the source, i.e the current vector potential, must be projected on
> the space spanned by the vector-valued shape functions. Otherwise, the
> simulation is numerically unstable. The last, from my experience, is
> definitely true.

I think you already found your solution, but just for clarity: When we use the
term "interpolation", we typically ask for (scalar) coefficients U_i so that

u_h(x_j) = sum U_i \phi_i(x_j) = g(x_j)

where g is given and x_j are the node points. The problem is that for the
Nedelec element, this is not always possible: Not all possible vectors g(x_j)
can be represented. This makes sense because if you have N node points in 3d,
then you have N scalar coefficients U_i, but you have 3N components of the
values g(x_j).

One way to deal with this is to associate a vector, let's say y_j with every
node point x_j, and require that

y_j \cdot u_h(x_j) = sum U_i y_j \cdot \phi_i(x_j) = y_j \cdot g(x_j)

These y_j could, for example, be the tangential direction associated with the
shape function phi_j. One *could* call this a variation of the term
"interpolation", but it is not what VectorTools::interpolate() implements. It
would probably not be terribly difficult to implement this kind of function,
however!

Best
W.



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

[Jean-Paul Pelteret]

To add to Wolfgang’s comment, you can have a look at the compute_edge_projection_l2() function that forms a part of the function that computes curl-conforming boundary conditions. You might be able to gain some inspiration from there? The associated paper that explains the method that’s being implemented is 

// See, for example, section 4.2 of:
// Electromagnetic scattering simulation using an H (curl) conforming hp-
// finite element method in three dimensions, PD Ledger, K Morgan, O
// Hassan, Int. J. Num. Meth. Fluids, Volume 53, Issue 8, pages
// 1267-1296, 20 March 2007:
// http://onlinelibrary.wiley.com/doi/10.1002/fld.1223/abstract

Best,

Jean-Paul

--
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.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/567f6646-df23-8c97-f210-87fd5a74e3ae%40colostate.edu.

[John Smith]

Dear Wolfgang,

Thank you for your reply. It is a bit difficult to read formulas in text. So I have put a few questions I have in a pdf file. Formulas are better there. It is attached to this message. May I ask you to have a look at it?

Best,

John

[Wolfgang Bangerth]

On 4/18/21 12:24 PM, John Smith wrote:
>
> Thank you for your reply. It is a bit difficult to read formulas in text. So I
> have put a few questions I have in a pdf file. Formulas are better there. It
> is attached to this message. May I ask you to have a look at it?

John:

If I understand your first question right, then you are given a vector field J
and you are looking for a vector field T so that
curl T = J
I don't really have anything to offer to this kind of question. There are many
vector fields T that satisfy this because of the large null space of curl. You
have a number of condition you would like to "approximate" but there are many
ways to "approximate" something. In essence, you have two goals: To satisfy
the equation above and to approximate something. You have to come up with a
way to weigh these two goals. For example, you could look to minimize
F(T) = ||curl T - J||^2 + alpha ||T-T_desired||^2
where T_desired is the right hand side of (0.0.5) and you have to determine
what alpha should be.

As for your other questions: (0.0.9) is indeed called "interpolation"

The issue with the Nedelec element (as opposed to a Q(k)^dim) field is that
for the Nedelec element,
\vec phi(x_j)
is not independent of
\vec phi(x_k)
and so you can't choose the matrix in (0.0.8) as the identity matrix. You
realize this in (0.0.13). The point I was trying to make is that (0.0.13)
cannot be exactly satisfied if T(x_j) is not a vector parallel to \hat e_j,
which in general it will not be. You have to come up with a different notion
of what it means to solve (0.0.13) because you cannot expect the left and
right hand side to be equal.

Best
W>

[Konrad Simon]

Dear John,

I had a look at the pdf you sent. I noticed some conceptual inconsistencies that are important for the discretization. Let me start like this: The curl-curl-problem that you are trying to solve is - according to your description of the gauge - actually a curl-curl + grad-div problem and hence a Laplace-problem that describes a Helmholtz-decomposition. In other words in order to get a well-posed simulation problem you need more boundary conditions. You already noticed that since your curl-curl-matrix is singular (the kernel is quite complicated and contains all longitudinal waves). 

You have of course the option to solve (0.0.2) with a Krylov-solver but then you need to make sure that the right-hand side is in the orthogonal complement of this kernel before each iteration which is quite difficult. I would not recommend that option.

Another point is that if you have a solution for your field A like in (0.0.4) you can not have a similar representation for T. This is mathematically not possible.

Since you not care about the gauge let me tell you how I would tackle this: The reason  is  - to come back to my original point - you are missing a boundary condition that makes you gauge unique. Since you apply natural boundary conditions (BCs) on A you must do so as well to determine div(A). This second BC for A is applied to a different part of the system that is usually neglected in the literature and A is partly determined from this part (this part then describes the missing transversal waves). The conditions you mention on curl(A) are some what contradictory to the space in which A lives. Nedelec elements which you must use (you can not use FE_Q to enforce the conditions) cannot generate your desired T_0.

There are some principles when discretizing these problems (which are not obvious) that you MUST adhere to (choice of finite elements, boundary conditions, what is the exact system etc) if you want to get a stable solution and these are only understood recently. I am solving very similar problems (with Deal.II) in a fluid context and will be happy to share my experiences with you. Just email me: konrad...@uni-hamburg.de

Best,
Konrad

[John Smith]

Dear Konrad,

Thank you for your help! I have recently modified the step-3 tutorial program to save FE_Nedelec<2> fe(0) shape function into *.csv files. It is attached to this message. The result of this program is not what I would expect it to be. I have also attached a pdf file which compares my expectations (marked as “Expectations” in the pdf file) with the results of the modified step-3 program (marked as “Reality” in the pdf file). My expectations are based on the literature I have. My question is: why the shape functions of FE_Nedelec<2> fe(0) differ from my expectations? I guess my expectations are all wrong. I would also appreciate if you will share with me a reference to a paper or a book, so I can extend my modest library.

Thanks!

John

P.S. I have sent you an e-mail as you have requested. Could you please check your spam folder?

[Jean-Paul Pelteret]

Hi John,

Unfortunately, you’ve fallen into the trap of confusing what the entries in the solution vector mean for the different element types. For Nedelec elements, they really are coefficients of a polynomial function, so you can’t simply set each coefficient to 1 to visualise the shape function (like one might be inclined to try for Lagrange type polynomials). If you want a little piece of code that will plot the Nedelec shape functions for you, then take a look over here:

https://github.com/dealii/dealii/issues/8634#issuecomment-595715677

BTW. In case you’re not aware, a newer (and probably more robust) implementation of a Nedelec element is the FE_NedelecSZ that uses hierarchical shape functions to form the higher order bases. The PhD thesis of S. Zaglmayr that describes the element is mentioned in the class documentation.  

Best,

Jean-Paul

-- 
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.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/348bd8f9-0891-42e5-b150-c3723eacb84bn%40googlegroups.com.
<FE_Nedelec.pdf><step-3.cc>

[Konrad Simon]

Dear John,

As Jean-Paul pointed out the entries of the solution vector for a Nedelec element have a different meaning. The entries are weights that result from edge moments (and in case of higher order also face moments and volume moments), i.e, integrals on edges. Nedelec elements have a deep geometric meaning: they discretize quantities that you can integrate (a physicist would say "measure") along 1D-lines in 3D.

Btw, if you are using Nedelec in 2d be aware that it has some difficulties on complicated meshes. In 3D the lowest order is fine. NedelecSZ (for 2D lowest order and 3D all orders), as Jean-Paul pointed out, are another option. They are meant to by-pass certain mesh orientation issues on complicated geometries (I can tell you a bit about that since currently I am chasing some problems there).

Best,
Konrad

[John Smith]

Dear Jean-Paul,

Thanks! It is indeed a trap.

John

On Tuesday, April 20, 2021 at 9:23:23 PM UTC+2 Jean-Paul Pelteret wrote:

[John Smith]

Dear Konrad,

I do not understand when a geometry gets complicated. Is a toroid inside a sphere, both centred at the origin, a complicated geometry? To start with, I will use a relatively simple mesh. I will not use local refinement, only global. I can do without refinement as well, i.e. making the mesh in gmsh.

Yes, I would like to know more about orientations issues on complicated geometries. Could please tell me about it? I would very much prefer to use FE_Nedelec without hierarchical shape functions. The first order 3D FE_Nedelec will do nicely provided the orientation issues will be irrelevant to the mesh.

John

[Jean-Paul Pelteret]

Dear John,

The class documentation states

> Several aspects of the implementation are experimental. For the moment, it is safe to use the element on globally refined meshes with consistent orientation of faces.

My interpretation would be that the mesh should be Cartesian aligned, or at least “grid-like". My limited experience with the standard FE_Nedelec element on a mesh of similar geometric “complexity” to yours was not satisfactory. Of course, I wouldn’t discount that I made a mistake somewhere, so it would be good if others who have spent more time using the FE_Nedelec element would share their experiences. Konrad is actively working on improving the current situation (see https://github.com/dealii/dealii/pull/12016, thanks Konrad!), but naturally it's an ongoing project. I believe that the FE_NedelecSZ does not suffer from the face/edge orientation issues, which is why I recommended that you try it (especially if you’re going to stick to the lowest order elements that, all things being equal, should render identical results to the canonical Nedelec element).

Best,

Jean-Paul

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/b82d4595-8508-4201-98d7-c554019f6574n%40googlegroups.com.

[John Smith]

Dear Jean-Paul,

Thank you for sharing this information with me!

Could you please check if I understand it correctly:

If I use the FE_NedelecSZ<3> fe(0) on a mesh imported from gmsh without any refinement in deal.ii, then the FE_NedelecSZ<3> fe(0) finite elements will not be hierarchical (as there is no refinement) . This will allow me to escape the face/edge orientation issues. In this case, the FE_NedelecSZ<3> fe(0) will be as FE_Nedelec<0> fe(0) minus the issues with the orientation of the faces and edges.

Is this correct? If it is, then I switch to FE_NedelecSZ<3> fe(0).

I have attached an example of a mesh I am thinking about. There are three figures: inner part of the mesh, outer sphere (infinity), and an adaptor (it just keeps the mesh cells nicely aliened) . I think it is “edge-elements friendly”…

John

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/A86FD959-B364-460E-B77C-CA24133C4F38%40gmail.com.

[Konrad Simon]

Hi John,

Maybe I can give some hints of how I would approach this problem (these are just some quick thoughts):

Write this problem as a vector Laplace problem

curl(nu*curl(A)) - grad(div(A)) = J

which you then have to write as a system of PDEs. Note, this can only be the same as

curl(nu*curl(A)) = J

if J is a curl itself.

If you do not want to impose conditions on the divergence then you have two options:

1.) You use sigma = div(A) as an auxiliary variable and impose natural boundary conditions (i.e., the ones that you get through integration by parts) on A*n and nu*curl(A)xn, where n is the outward unit normal. This means that sigma is a 0-form and A is a 1-form -> sigma is then discretized (for example) by FE_Q(n+1) and A by FE_Nedelec(n)

System:

(A1) sigma + div(A) = 0
(B1) grad(sigma) + curl(nu*curl(A)) = J

For the weak form you multiply (A1) with a FE_Q test function and integrate the second summand by parts. Secondly, you multiply (B1) with a Nedelec test function and also integrate the second summand by parts. You will see your natural boundary conditions pop up.

2.) You use sigma = nu*curl(A) as an auxiliary variable and impose essential boundary conditions on A*n and nu*curl(A)xn, where n is the outward unit normal. This means that sigma is interpreted as a 1-form and A is a 2-form -> sigma is then discretized (for example) by FE_Nedelec(n) and A by FE_RaviartThomas(n)

System:

(A2) (1/nu)*sigma - curl(A) = 0
(B2) curl(sigma) -  grad(div(A)) = J

For the weak form you multiply (A2) with a Nedelec test function and integrate the second summand by parts. Secondly, you multiply (B2)
with a Raviart-Thomas test function and also integrate the second summand by parts. You will NOT see your natural boundary conditions pop up since conditions on A*n and nu*curl(A)xn are essential in this case. You need to enforce them on the function spaces directly. In deal.ii you do so by using project_boundary_values_curl_conforming_l2 or project_boundary_values_div_conforming

Note that the additional boundary conditions make the system invertible and if J is a curl it will turn out that div(A)=0.

There is a field in mathematics that is called finite element exterior calculus (FEEC) that answers the question of stability when solving such problems. See this book. I guess Chapters 4, 5 and 6 are most interesting for you.

I do not understand when a geometry gets complicated. Is a toroid inside a sphere, both centred at the origin, a complicated geometry? To start with, I will use a relatively simple mesh. I will not use local refinement, only global. I can do without refinement as well, i.e. making the mesh in gmsh.

Yes, I would like to know more about orientations issues on complicated geometries. Could please tell me about it? I would very much prefer to use FE_Nedelec without hierarchical shape functions. The first order 3D FE_Nedelec will do nicely provided the orientation issues will be irrelevant to the mesh.

 Lowest order Nedelec and all other elements I mentioned are fine on the meshes you need I guess.

Hope that helps.

Best,

Konrad