Non-nodal FEs (e.g., Raviart-Thomas and Nedelec)

[santiagobadia]

In the near future, we should provide non-nodal FE spaces as Raviart-Thomas for H(div) approximation and Nedelec for H(curl) approximation, by simply defining the right prebasis and DoFs in the ConformingFESpace

[santiagobadia]

@fverdugo No pressure but we still need the anisotropic monomial dev in this issue to start implementing this functionality.

[santiagobadia]

@fverdugo I have worked on the polytope and now we can generate nonuniform nodes arrays. I have one question. I feel that in some parts of the code I am re-implementing the things we have at CellValues. For instance, to provide the DOFs at the closure of a NFace, I would use:

function closurenfacedofs(this::LagrangianRefFE{D,T} where {D,T})
  cv1 = CellValueFromArray(this.polytope.nf_nfs) # cell to index
  cv2 = CellValueFromArray(this.nfacedofs)
  return CellVectorByComposition(cv1,cv2)
end

I think that it is reasonable to use such structures at the RefFE / Polytope level. What do you think? Any performance issue? I would say not a big point.

[fverdugo]

@fverdugo No pressure but we still need the anisotropic monomial dev in this issue to start implementing this functionality.

Arbitrary anisotropic Q-monomials are not implemented, but I have implemented two particular cases: See GradMonomialBasis and CurlGradMonomialBasis in Gridap. I think this is actually what you need right? It is really necessary to have arbitrarily anisoropic monomials, or with the two particular cases is enough for our purposes?

[fverdugo]

@fverdugo I have worked on the polytope and now we can generate nonuniform nodes arrays. I have one question. I feel that in some parts of the code I am re-implementing the things we have at CellValues. For instance, to provide the DOFs at the closure of a NFace, I would use:

function closurenfacedofs(this::LagrangianRefFE{D,T} where {D,T})
  cv1 = CellValueFromArray(this.polytope.nf_nfs) # cell to index
  cv2 = CellValueFromArray(this.nfacedofs)
  return CellVectorByComposition(cv1,cv2)
end

I think that it is reasonable to use such structures at the RefFE / Polytope level. What do you think? Any performance issue? I would say not a big point.

With which granularity this function is to be called? If you call it a lot, perhaps it is better to pre-compute it in the constructor of the struct.

By the way, if you want to work with native julia arrays, you can always collect the values, e.g.: return collect(CellVectorByComposition(cv1,cv2)). This will return a Vector{Vector{Int}}. In that case, it is better to precompute the result.

For very small data like this one, perhaps it is better to precompute it as native julia arrays. (Here we do not really need lazyness, and native arrays will simplify things)

[santiagobadia]

In this case it is not a problem because I think that this method will be called few times. But I wanted to know whether you think it is reasonable to use these cell arrays for fine grain operations.

[santiagobadia]

With regard to anisotropic Lagrangian methods I have a constructor that allows one to provide a vector of orders per direction. For n-cubes all the polytope machinery is working for anisotropic order, the only issue is the generation of anisotropic monomials. So, I would include this option in the future. For the moment I am enforcing equal order.

You are right. With the ones you already have I can implement RaviartThomas. I will create a constructor especific for such RefFEs.

[fverdugo]

In this case it is not a problem because I think that this method will be called few times. But I wanted to know whether you think it is reasonable to use these cell arrays for fine grain operations.

The rule of thumb is: construct the object only once and use it as much as you want. They should be type stable and efficient (unless there is a bug in the code). However, for such small data, I will definitively collect the values in order to have native julia arrays (this simplifies things, speeds up compilation times, etc, and for small data lazyiness is not really needed).

[fverdugo]

@santiagobadia just another important thing to be aware of w.r.t CellValue usage.

See issue #91

[santiagobadia]

Raviart-Thomas FE methods working for arbitrary order!!!

[santiagobadia]

I would like to re-write some parts to make it easy to extend it for other methods like Nedelec, BDM, etc.

[fverdugo]

great @santiagobadia !

In the rewiring process, it would be nice to add some extra tests for the RT elements. I have seen that the coverage has dropped about 3%.

[santiagobadia]

I forgot to include the RT test in the runtests. Now it has increased a 3%.

[santiagobadia]

Pending tasks:

• [x] To define a set of methods for n-face moments to be used by the RefFE generator. This way, we can define a new RefFE by just providing these methods.

• [ ] To define a new method local_to_global_dof_value and global_to_local_dof_value to be able to deal with more general meshes

• [x] Piola and contra-Piola maps for general meshes

• [x] Nedelec spaces for quads

• [ ] RT spaces for tets

[fverdugo]

@fverdugo No pressure but we still need the anisotropic monomial dev in this issue to start implementing this functionality.

Done!

[santiagobadia]

@fverdugo I have worked on the Darcy tutorial in Tutorials repos. I have checked that the method provides the exact solution but when plotted using writevtk the result in Paraview has no sense. I guess I have to look at how it is written, probably there is something to change there. Not sure I will have time before leaving. I won't report flux plots.

[fverdugo]

This is something that we need to look at.

The machinery to interpolate the solution for visualization is the same used to compute the values of the solution at integration points. Thus, if the error norms are OK, visualizing should be OK as well (up to a dummy error somewhere)...

[fverdugo]

I can take a look in what is going on. And include the plot of the flux if I fix the issue.

[santiagobadia]

Not, the problem is somewhere else, the visualization is OK. The test in Gridap works perfect but the driver does not work as expected. I have to investigate it.

[santiagobadia]

I am trying to implement non-trivial field maps from the reference basis to the physical basis. See, e.g., the beginning of p 9 here.

For the moment, we are assuming an identity operator, as for Lagrangian FEM, i.e., do-nothing.

I want to implement the Piola map. See Eq. (17) in the same reference.

I was trying to express the physical basis given a physical point x already mapped in the reference space hatx, which is all we need. I have written the following code

using Gridap
import Gridap: ∇

import LinearAlgebra: det

model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,2))
trian = Triangulation(model)
basis = CellBasis(trian)
phi = CellGeomap(trian)
jac = gradient(phi)
detjac = det(jac)
# I would call it field_geomap
piola_map = (1.0/detjac)*jac

V = FESpace( reffe=:RaviartThomas, conformity=:HDiv, order=1, model=model, diritags = [5,6])
basis = V.cellbasis

quad = CellQuadrature(trian,degree=1)
z = coordinates(quad)

pv = evaluate(piola_map,z)

bv = evaluate(basis,z)

phys_basis = piola_map*basis

dbv = evaluate(phys_basis,z)

It seems to work.

So, now we should decide were to put it in our FE integration machinery. Probably just to change the constructor for div-conforming FE spaces, changing the current cell basis for the one with the Piola map. Not sure we need the untouched one at all.

[fverdugo]

If I have understood correctly, you end up with a basis which is evaluated at the quadrature points in the reference space, but the value of the basis refers to quantities in the physical space, right?

For Lagrangian elems we have something similar, we map the gradients to the physical space (but we still have an object evaluated at the reference space). So, we can accommodate both things in the same framework.

It seems that each type of reference element needs to map the basis to the physical space in a different way. For Lagrangian only modify the gradient, for RT use piola map etc. Thus, we need a new function with mutltiple dispatch.

Something like:

map_to_physical_space(reffe::RefFE,basis::CellBasis,geomap::CellGeomal) -> GellBasis

For Lagrangian elems it will be implemented like this

function map_to_physical_space(reffe::LagrangianRefFE,basis::CellBasis,geomap::CellGeomal)
  attachgeomap(basis,geomap)
end

And for RT, it will be implemented as you have detailed.

Then, this new function has the be called in the ConformingFESpace constructor instead of calling attachgeomap as it is done now.

We can start with something like this. But, it can be improved, since this interface assumes a unique reference elem in the entire mesh. Perhaps, we can change reffe::RefFE by reffes::CellValue{RefFE} to fix this.

[santiagobadia]

Hi Francesc

The attachgeomap is also needed here, since it is needed to compute derivatives of the shape functions. In fact, for non-affine maps, we would also need to compute the derivatives of the Piola map (again, we need the geomap here) and use the chain rule. I think it is very complicating and not sure we need it now because RT or Nedelec on non-affine maps have problems.\

[santiagobadia]

I have the following issue.

For grad-conforming spaces, the a shape function in the physical space g(x) is expressed as a reference function g'(x) and a geometrical map m (from ref to physical space) as g = g(inv(m(x))). Thus, using differentiation rules, we can check that its derivative is grad(g) = grad(g')grad(inv(m)), i.e., grad(g')inv(Jac). This is what we have implemented.

Now the situation is more involved because g= P(g'(m)), where P is the Piola map. So, its derivative is the sum of P(grad(g')inv(Jac)) and one term that involves the derivative of P. For non-affine maps, these methods have problems, and for affine maps the term that involves derivatives of P vanish and we simply have P grad(g')inv(Jac). Thus, we only need to implement the application of the Piola map over the gradient in the physical space obtained as usual.

Problem 1 is that we have not implemented the gradient of CellMap as composition (it is not that simple). So, if I want to avoid this implementation (using the fact that the expression can be simplified for affine maps) and consider only the case of affine geometrical maps. I think that I need: 1) Define a cell basis concrete version PiolaCellMap <: CellMap with attributes P and the basis with attached geomap already in Gridap representing g'(m) and a method gradient for this type that simply computes the gradient of g'(m). Do we agree on this?

Not sure it is clear enough.