Vector Finite Elements

[lindsayad]

Rationale

For solving magnetostatic problems for example, we require continuity of the tangential component of the vector potential but not continuity of the normal component.

Description

A good element choice for the problem type above yields degrees of freedom corresponding to the tangential projection of the vector potential onto element edges. These edge elements are typically Raviart-Thomas elements for H(div) problems and Nedelec elements for H(curl) problems. LibMesh at least has an example with Nedelec elements. Implementing vector elements will take some good design thought, perhaps borrowing some concepts from #6881.

Impact

Increase our capability for solving electromagnetic problems.

[friedmud]

This is definitely difficult - but not because it's technically difficult... just that there are many design decisions... and the need to get it in there without effecting the performance of everything else.

I definitely think the ArrayKernel stuff shows some promise as to how to compute vectors of residuals... but there were still many open questions on the design even as it was.

One of the very first issues that you pretty much need to change computeQpResidual() to a "set" syntax instead of "return". Definitely take a look at what I did in ArrayKernel where I used Eigen's capabilities to remap memory behind an EigenVector to give us something that looks as "MOOSE-like" as possible.

[lindsayad]

So, the i'th residual for a vector fe problem is still a scalar, so in my mind our kernel design for vector finite elements really shouldn't look any different. However, the huge difference lies in our assembly and FEProblem etc. because our shape functions are now vectors instead of scalars. The nested class FEShapeData for example in Assembly assumes scalar shape functions; same thing for our zero interface in FEProblemBase. LibMesh handles this by making FEGenericBase templated. @idaholab/moose-team would we be open to templating our FEProblem? One consequence of this which I was hoping to avoid is that AFAIK this would make it impossible to do vector FE and scalar FE calculations in a single app; coupling would have to be done through multi-apps. Someone could correct me on that though.

[friedmud]

Let's not think about templating yet... did you look at what I did in the VectorKernel stuff? Basically... make new MooseVariable types... that way we can mix things together more seamlessly.

As for your statement about scalar values... it's not really a great way to look at it. To really take advantage of the vector valued FE stuff... you really want Kernels to be able to compute vector values all at one time (not one component per call). This is isn't insurmountable... you can see what I did in the VectorKernel stuff for instance.

On Tue, Oct 17, 2017 at 10:12 AM Alex Lindsay notifications@github.com wrote:

So, the i'th residual for a vector fe problem is still a scalar, so in my mind our kernel design for vector finite elements really shouldn't look any different. However, the huge difference lies in our assembly and FEProblem etc. because our shape functions are now vectors instead of scalars. The nested class FEShapeData for example in Assembly assumes scalar shape functions; same thing for our zero interface in FEProblemBase. LibMesh handles this by making FEGenericBase templated. @idaholab/moose-team https://github.com/orgs/idaholab/teams/moose-team would we be open to templating our FEProblem? One consequence of this which I was hoping to avoid is that AFAIK this would make it impossible to do vector FE and scalar FE calculations in a single app; coupling would have to be done through multi-apps. Someone could correct me on that though.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/idaholab/moose/issues/10049#issuecomment-337244524, or mute the thread https://github.com/notifications/unsubscribe-auth/AA1JMR5Z6aBre1UX8K_xV1vUOQUJU4vjks5stLXQgaJpZM4PzGgm .

[lindsayad]

We might not be on the same page. I'm saying that the weak form for a vector FE problem yields one scalar equation, not three component equations. For example, the weak form of:

is:

Thus, I wouldn't be calling a kernel three separate times for three different moose variables corresponding to components. A kernel would be called one single time for a single VectorMooseVariable whose values at quadrature points are vectors instead of scalars. This VectorMooseVariable is based off your ArrayMooseVariable :-)

[jwpeterson]

I think @lindsayad is on the right track here. In his example, we have to think of $\vec{E}$ as a "single" variable even though we are internally solving for its 3 components. computeQpResidual and computeQpJacobian should still return Real...

[YaqiWang]

computeQpResidual is used to fill in the global residual vector of 3 components. So It cannot return a scalar value. The three components of v is arbitrary and independent and do not appear in the computeQpResidual function. When making the design, please also consider the size of the vector variable could be arbitrary (not necessarily space dimension). It will be important for radiation transport.

[lindsayad]

Yes, $\vec{E}$ is truly a single variable. We're just moving from:

to:

[lindsayad]

What do you mean they're arbitrary? The shape functions for vector finite elements are explicitly defined. Here's the code in libmesh for first order Nedelec elements on 2D triangles:

              switch(i)
                {
                case 0:
                  {
                    if (elem->point(0) > elem->point(1))
                      return RealGradient( -1.0+eta, -xi );
                    else
                      return RealGradient( 1.0-eta, xi );
                  }
                case 1:
                  {
                    if (elem->point(1) > elem->point(2))
                      return RealGradient( eta, -xi );
                    else
                      return RealGradient( -eta, xi );
                  }

                case 2:
                  {
                    if (elem->point(2) > elem->point(0))
                      return RealGradient( eta, -xi+1.0 );
                    else
                      return RealGradient( -eta, xi-1.0 );
                  }

As you can see the shape functions themselves are vectors

[lindsayad]

Here's a libmesh code snippet assembling the residual and Jacobian for vector FEM:

  // Loop over quadrature points
  for (unsigned int qp=0; qp != n_qpoints; qp++)
    {
      Gradient u;
      Gradient curl_u;

      c.interior_value(u_var, qp, u);

      c.interior_curl(u_var, qp, curl_u);

      // Value of the forcing function at this quadrature point
      RealGradient f = this->forcing(qpoint[qp]);

      // First, an i-loop over the degrees of freedom.
      for (unsigned int i=0; i != n_u_dofs; i++)
        {
          Fu(i) += (curl_u*curl_phi[i][qp] + u*phi[i][qp] - f*phi[i][qp])*JxW[qp];

          // Matrix contributions for the uu and vv couplings.
          if (request_jacobian)
            for (unsigned int j=0; j != n_u_dofs; j++)
              Kuu(i,j) += (curl_phi[j][qp]*curl_phi[i][qp] +
                           phi[j][qp]*phi[i][qp])*JxW[qp];
        }
    } // end of the quadrature point qp-loop

• curl_u is a vector
• curl_phi[i][qp] is a vector.
• u is a vector
• phi[i][q] is a vector
• curl_u*curl_phi[i][qp] : scalar
• u*phi[i][qp] : scalar
• Fu(i) : scalar

[jwpeterson]

@YaqiWang there may be a terminology issue. The vector-valued FEs (Nedelec, Raviart-Thomas etc) are not arbitrarily long. They correspond to the spatial dimension in which you are solving the problem...

[lindsayad]

@YaqiWang I think we have to make a distinction between vector and perhaps array. To me vector FEM is a fundamentally spatial concept where the dimension should be no greater than our physical dimension of 3. The shape functions in vector FEM are fundamentally that: vectors.

[YaqiWang]

Oh, I have a misunderstanding. I thought the DoFs of E contains three components. Yes, vector and array are then two different things.

[andrsd]

I do not know if this was discussed, but how's this idea (just an idea, it might have issues):

1. Refactor Assembly, so you have one for scalar- and one for vector-valued FEs.
2. The current kernels will use the scalar-values one.
3. For vector-valued FE's, there would be something like VectorFEKernelBase that will get _phis from vector-valued Assembly and _us and _curl_u and what not from VectorMooseVariable (or whatever name it gets). Similarly for BCs.
4. Then, VectorMooseVariable computing _u, _curl_u, etc.

If 1. happens, we would probably need to allocate the right Assembly type when the corresponding FE variable is added. Not sure if 1. should be part of this. It can be probably done later.

I do not think I want to template FEProblemBase.

I only did simple electromagnetic eqns in 2009, so it has been a while, but the set of what we call kernels was different. That's why I am thinking of a parallel branch in Kernel hierarchy. It might accommodate Hdiv spaces as well eventually.

[lindsayad]

I personally like including 1.

[andrsd]

Ok. Then, for 1) notice that some data will be common between scalar- and vector-valued Assembly - things like _q_points, _JxW and probably others. So maybe have that in a base class and then inherit the scalar-values and vector-valued Assembly from that.

Then another thing comes to my mind: things like reinitFEs will need to be virtual, because of the different FETypes. But, you probably thought of this already...

I am also assuming that you want to include 1) because of the FEShapeData caching, right?

Last thing, since nobody really commented on it. I am fine with calling one system Vector and the other Array. I think Vector and VectorFE would be just too prone to misunderstandings. Unless somebody has a better suggestion.

[lindsayad]

Yes, I think 1) will just cut down on a lot of code duplication.

On Thu, Oct 19, 2017 at 10:14 AM, David Andrs notifications@github.com wrote:

Ok. Then, for 1) notice that some data will be common between scalar- and vector-valued Assembly - things like _q_points, _JxW and probably others. So maybe have that in a base class and then inherit the scalar-values and vector-valued Assembly from that.

Then another thing comes to my mind: things like reinitFEs will need to be virtual, because of the different FETypes. But, you probably thought of this already...

I am also assuming that you want to include 1) because of the FEShapeData caching, right?

Last thing, since nobody really commented on it. I am fine with calling one system Vector and the other Array. I think Vector and VectorFE would be just too prone to misunderstandings. Unless somebody has a better suggestion.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/idaholab/moose/issues/10049#issuecomment-337958812, or mute the thread https://github.com/notifications/unsubscribe-auth/AJxgcNR3lQj7tLt3Icz4dRUWnXzF8Aheks5st3VLgaJpZM4PzGgm .

[friedmud]

As long as we’re refactoring in there: I don’t see any reason to keep the FEShapeData caching stuff. I don’t know of anyone that actually uses that and it adds tons of complication.

It was a cool idea that I implemented a long time ago: but I’m just not sure it’s worth the code complexity.

Also: Vector and Array are good with me. On Thu, Oct 19, 2017 at 3:14 PM Alex Lindsay notifications@github.com wrote:

Yes, I think 1) will just cut down on a lot of code duplication.

On Thu, Oct 19, 2017 at 10:14 AM, David Andrs notifications@github.com wrote:

Ok. Then, for 1) notice that some data will be common between scalar- and vector-valued Assembly - things like _q_points, _JxW and probably others. So maybe have that in a base class and then inherit the scalar-values and vector-valued Assembly from that.

Then another thing comes to my mind: things like reinitFEs will need to be virtual, because of the different FETypes. But, you probably thought of this already...

I am also assuming that you want to include 1) because of the FEShapeData caching, right?

Last thing, since nobody really commented on it. I am fine with calling one system Vector and the other Array. I think Vector and VectorFE would be just too prone to misunderstandings. Unless somebody has a better suggestion.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/idaholab/moose/issues/10049#issuecomment-337958812, or mute the thread < https://github.com/notifications/unsubscribe-auth/AJxgcNR3lQj7tLt3Icz4dRUWnXzF8Aheks5st3VLgaJpZM4PzGgm

.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/idaholab/moose/issues/10049#issuecomment-338008709, or mute the thread https://github.com/notifications/unsubscribe-auth/AA1JMYaQ-3WKot7g3whk6HIfZjmqEpAwks5st5-egaJpZM4PzGgm .

[lindsayad]

Alright cool; yea that reduces a lot of the complexity

On Thu, Oct 19, 2017 at 1:35 PM, Derek Gaston notifications@github.com wrote:

As long as we’re refactoring in there: I don’t see any reason to keep the FEShapeData caching stuff. I don’t know of anyone that actually uses that and it adds tons of complication.

It was a cool idea that I implemented a long time ago: but I’m just not sure it’s worth the code complexity.

Also: Vector and Array are good with me. On Thu, Oct 19, 2017 at 3:14 PM Alex Lindsay notifications@github.com wrote:

Yes, I think 1) will just cut down on a lot of code duplication.

On Thu, Oct 19, 2017 at 10:14 AM, David Andrs notifications@github.com wrote:

Ok. Then, for 1) notice that some data will be common between scalar- and vector-valued Assembly - things like _q_points, _JxW and probably others. So maybe have that in a base class and then inherit the scalar-values and vector-valued Assembly from that.

Then another thing comes to my mind: things like reinitFEs will need to be virtual, because of the different FETypes. But, you probably thought of this already...

I am also assuming that you want to include 1) because of the FEShapeData caching, right?

Last thing, since nobody really commented on it. I am fine with calling one system Vector and the other Array. I think Vector and VectorFE would be just too prone to misunderstandings. Unless somebody has a better suggestion.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <https://github.com/idaholab/moose/issues/10049#issuecomment-337958812 , or mute the thread < https://github.com/notifications/unsubscribe-auth/ AJxgcNR3lQj7tLt3Icz4dRUWnXzF8Aheks5st3VLgaJpZM4PzGgm

.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/idaholab/moose/issues/10049#issuecomment-338008709, or mute the thread https://github.com/notifications/unsubscribe-auth/AA1JMYaQ- 3WKot7g3whk6HIfZjmqEpAwks5st5-egaJpZM4PzGgm

.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/idaholab/moose/issues/10049#issuecomment-338013748, or mute the thread https://github.com/notifications/unsubscribe-auth/AJxgcKIya4JG_baRBvHM2sjiG91IyiK2ks5st6RwgaJpZM4PzGgm .

[YaqiWang]

You may want to set up benchmarks to ensure refactoring does not degrade performance. I saw some fundamental changes in MOOSE recently. I think the sooner we set up benchmarks the better.

[mangerij]

Just to motivate you, I can think of some very high impact projects if the Maxwell problem can be solved when coupled to the output from some existing modules/applications.

I understand that this is a difficult task, but most of the FDTD solvers available are sad black boxes.

[gridley]

More motivation: this would make adding P1 neutronics to Moltres a lot more elegant!

[lindsayad]

Closed by #10238