Re: [deal.II] interpolate FE_Nelelec

[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  

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/c2773841-4e03-42a1-9f1e-e1a41c0933c9n%40googlegroups.com.

Re: [deal.II] interpolate FE_Nelelec

[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

> On 21. Apr 2021, at 16:23, John Smith  wrote:
> 
> 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
> 
> 
> 
> On Tuesday, April 20, 2021 at 9:43:31 PM UTC+2 Konrad Simon wrote:
> 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
> 
> On Tuesday, April 20, 2021 at 9:23:23 PM UTC+2 Jean-Paul Pelteret wrote:
> 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
> 
> 
>> On 20. Apr 2021, at 15:47, John Smith > wrote:
>> 
> 
>> 
>> 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. 

Re: [deal.II] interpolate FE_Nelelec

[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:

> 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
>
> On 20. Apr 2021, at 15:47, John Smith  wrote:
>
>
> 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? 
>
>
>
> On Monday, April 19, 2021 at 8:47:29 AM UTC+2 Konrad Simon wrote:
>
>> 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
>>
>> On Monday, April 19, 2021 at 5:51:17 AM UTC+2 Wolfgang Bangerth wrote:
>>
>>> 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 

Re: [deal.II] interpolate FE_Nelelec

[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

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

> 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
>
> On 20. Apr 2021, at 15:47, John Smith  wrote:
>
>
> 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? 
>
>
>
> On Monday, April 19, 2021 at 8:47:29 AM UTC+2 Konrad Simon wrote:
>
>> 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 

Re: [deal.II] interpolate FE_Nelelec

[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

> On 20. Apr 2021, at 15:47, John Smith  wrote:
> 
> 
> 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? 
> 
> 
> 
> 
> On Monday, April 19, 2021 at 8:47:29 AM UTC+2 Konrad Simon wrote:
> 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
> 
> On Monday, April 19, 2021 at 5:51:17 AM UTC+2 Wolfgang Bangerth wrote:
> 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 

Re: [deal.II] interpolate FE_Nelelec

[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

> On 16. Apr 2021, at 22:26, Wolfgang Bangerth  wrote:
> 
> 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 
>> 
>> 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: bange...@colostate.edu
>   www: http://www.math.colostate.edu/~bangerth/
> 
> -- 
> 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+unsubscr...@googlegroups.com.
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/567f6646-df23-8c97-f210-87fd5a74e3ae%40colostate.edu.

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/DDD82757-089D-4391-A367-5464AA304081%40gmail.com.

Re: [deal.II] interpolate FE_Nelelec

[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  

On Friday, April 16, 2021 at 8:50:07 PM UTC+2 Jean-Paul Pelteret wrote:

> 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 
> void WireMagneticField::verify_source_terms () const
> {
> TimerOutput::Scope timer_scope (computing_timer, "Verify source terms");
>
> hp::FEFaceValues 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::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::faces_per_cell; ++face)
> {
> hp_fe_face_values.reinit(cell,face);
> const FEFaceValues _face_values = hp_fe_face_values.
> get_present_fe_values();
> const unsigned int _fq_points = fe_face_values.n_quadrature_points;
>
> std::vector > source_values (n_fq_points, Vector >(dim));
> function_free_current_density.vector_value_list (fe_face_values.
> get_quadrature_points(),
> source_values);
>
> for (unsigned int fq_point=0; 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>  = 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!"));
> }
> }
>
> On 16. Apr 2021, at 19:36, John Smith  wrote:
>
> 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
>
> 

Re: [deal.II] interpolate FE_Nelelec

[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 
void WireMagneticField::verify_source_terms () const
{
  TimerOutput::Scope timer_scope (computing_timer, "Verify source terms");

  hp::FEFaceValues 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::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::faces_per_cell; ++face)
{
  hp_fe_face_values.reinit(cell,face);
  const FEFaceValues _face_values = 
hp_fe_face_values.get_present_fe_values();
  const unsigned int  _fq_points = fe_face_values.n_quadrature_points;

  std::vector > source_values (n_fq_points, 
Vector(dim));
  function_free_current_density.vector_value_list 
(fe_face_values.get_quadrature_points(),
   source_values);

  for (unsigned int fq_point=0; fq_point 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>  = 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!"));
  }
}

> On 16. Apr 2021, at 19:36, John Smith  wrote:
> 
> 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
> 
> On Fri, Apr 16, 2021 at 5:28 PM simon...@gmail.com 
>   

[deal.II] 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 
 in function*

* void dealii::VectorTools::interpolate(const dealii::Mapping&,*

* const DoFHandlerType&, const dealii::Function&, VectorType&, const dealii::ComponentMask&) 
[with int *

*dim = 3; int spacedim = 3; VectorType = dealii::Vector; 
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

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/1d2f6b83-bf76-4762-978b-1102d2481568n%40googlegroups.com.