Implementing Discontinuous Lagrange Interpolations

[bhaveshshrimali]

Hi,

It would be helpful if discontinuous interpolations on simplices could be supported. I am taking a look at the codebase. If someone could guide me on the relevant modules that need to be changed, I could take an attempt at implementing those. Thanks!

[kimauth]

Hi Bhavesh,
you probably found the interpolations.jl file already. You could just define a new subtype of Interpolation{dim,shape,order}. Depending on your interpolation and use-case, you might then want to have a look at the CellValues, which store pre-computed shape function values and derivatives. But if your interpolation is similar to the existing Lagrange interpolations, the CellValues might work out of the box. (I use a custom subtype of Interpolation for cohesive zone modelling and it didn't need much change.)
Eventually you also want to check out the DofHandler or MixedDofHandler if you need a non-standard dof-distribution. In that case it is useful to notice that the InterpolationInfo is used for the dof-distribution.

[koehlerson]

You can get some inspiration from https://github.com/Paulms/HDiscontinuousGalerkin.jl/, especially https://github.com/Paulms/HDiscontinuousGalerkin.jl/blob/master/src/basis.jl. Further you can check the modifications in the dofhandler and what additional information are stored in the FESpace which is a related entity to CellValues

Edit: As far as I can tell HDiscontinuousGalerkin.jl tries to do more spectral element stuff. If you go for a nodal representation DG scheme, you can just wrap the lagrangean interpolations. However, what you need to change is how the degrees of freedom are counted in MixedDofHandler or DofHandler, because you need to double them. Since you interconnect later those duplicate DOFs through some suitable flux at the boundary of each element I'd suggest to create a new CellValues subtype. Probably you need to store face and cell information inside them.

[bhaveshshrimali]

Hi Bhavesh,
you probably found the interpolations.jl file already. You could just define a new subtype of Interpolation{dim,shape,order}. Depending on your interpolation and use-case, you might then want to have a look at the CellValues, which store pre-computed shape function values and derivatives. But if your interpolation is similar to the existing Lagrange interpolations, the CellValues might work out of the box. (I use a custom subtype of Interpolation for cohesive zone modelling and it didn't need much change.)
Eventually you also want to check out the DofHandler or MixedDofHandler if you need a non-standard dof-distribution. In that case it is useful to notice that the InterpolationInfo is used for the dof-distribution.

Thanks Kim. I did look at interpolations.jl and in fact much of what I need is already in Lagrange bases. Will start by looking at DofHandler/MixedDofHandler to explore further.

[bhaveshshrimali]

You can get some inspiration from Paulms/HDiscontinuousGalerkin.jl, especially Paulms/HDiscontinuousGalerkin.jl@master/src/basis.jl. Further you can check the modifications in the dofhandler and what additional information are stored in the FESpace which is a related entity to CellValues

Edit: As far as I can tell HDiscontinuousGalerkin.jl tries to do more spectral element stuff. If you go for a nodal representation DG scheme, you can just wrap the lagrangean interpolations. However, what you need to change is how the degrees of freedom are counted in MixedDofHandler or DofHandler, because you need to double them. Since you interconnect later those duplicate DOFs through some suitable flux at the boundary of each element I'd suggest to create a new CellValues subtype. Probably you need to store face and cell information inside them.

Thanks for the pointers @koehlerson These are indeed very helpful. I wasn't aware of the HDiscontinuousGalerkin package. But you guessed it right, I am looking for a nodal discontinuous interpolation. I will take a look at DofHandler to start with. Let's see how it goes... Thanks again for these pointers!

[koehlerson]

To extend this a bit: Let's consider for simplicity just dofs at the vertices.
https://github.com/KristofferC/JuAFEM.jl/blob/master/src/Dofs/DofHandler.jl#L174 this statement checks if the dof is already in the vertexdict. The vertexdict maps a given vertex to a degree of freedom, which makes sense for continuous Galerkin methods. Now what you want to do is to double each dof. So, the most primitive way to do so is to remove that if statement. That would complete the first thing you need for DG.

Secondly, you need a new cellvalue type that gives you normals of the cell just as FaceValues. I don't know what's the best way to include it, maybe the most native thing to do in the beginning is a composition of CellValues and FaceValue and from this call just the functions that are already implemented.

Thirdly, the most challenging part probably, is to get the neighboring information from a given cell. They are not stored in the grid as it is right now, however, you can recompute the neighbors based on the connectivity that is stored in the grid.cells vector. If you have all of those three points you should be ready to go to assemble some conservation system with some flux element connection on the rhs.

[bhaveshshrimali]

Secondly, you need a new cellvalue type that gives you normals of the cell just as FaceValues. I don't know what's the best way to include it, maybe the most native thing to do in the beginning is a composition of CellValues and FaceValue and from this call just the functions that are already implemented.

Thirdly, the most challenging part probably, is to get the neighboring information from a given cell. They are not stored in the grid as it is right now, however, you can recompute the neighbors based on the connectivity that is stored in the grid.cells vector. If you have all of those three points you should be ready to go to assemble some conservation system with some flux element connection on the rhs.

I see, thanks again @koehlerson. I want to first implement (and test) just a discontinuous lagrange approximation. This can be useful in some mixed-formulations in incompressible elasticity/stokes. I wouldn't need facet normal information or inter-element flux for that. So I guess it is just modifying the vertexdict map I suppose.

Once I have that working, I will proceed with points (2) and (3) to implement a discontinuous galerkin method (which is also useful in incompressible elasticity -- without any mixed-formulation).

Does that sound correct? I guess, I would fully know only when I start implementing...

[koehlerson]

I'm not so familiar with incompressible elasticity, but checking https://www.sciencedirect.com/science/article/pii/S0045782501003589 or https://www.sciencedirect.com/science/article/pii/S0045782505002641 show both some interelement "flux" term in the weak form. I'm not quite sure how you would obtain a well posed solution at the "duplicate" dofs otherwise. Nonetheless, if there is some formulation that doesn't require (2) or (3), you still can go on with this roadmap. As I wrote before, you can test a simple linear case by removing the if statement and check if everything else is working down the line (which is what I expect). If so, you can get your head around how to make the dof counting "agnostic" for discontinuous/continuous interpolations. Most likely the part where you want to declare it is when you push a field in a DofHandler

[bhaveshshrimali]

I'm not so familiar with incompressible elasticity, but checking sciencedirect.com/science/article/pii/S0045782501003589 or sciencedirect.com/science/article/pii/S0045782505002641 show both some interelement "flux" term in the weak form

Absolutely! In a Discontinuous Galerking formulation there are "flux" terms. What I was referring to is a mixed formulation (displacement-pressure, say) for incompressible elasticity. Something similar to a Taylor-Hood approximation for Stokes, but instead of "continuous" pressure approximation there a discontinuous one in this case. With displacement enriched with bubbles and pressure being "discontinuous"-P1, the element satisfies the inf-sup condition. The space of pressure is a "discontinuous"-P1. It's also sometimes called the Crouzeix-Raviart (conforming) element after it was presented in their 1973 paper, see 4.13 therein. Since the basis functions are known over the reference simplex, I will start with first writing these and their gradients and then start looking at DofHandler.

Thanks

[koehlerson]

Ok, if you want Crouzeix-Raviart elements, then you actually have to do something slightly different. If we consider the linear case, where you have at each face in the middle sitting a DOF instead of the vertex position, then you only have to create a new interpolation with the correct dispatches on https://github.com/Ferrite-FEM/Ferrite.jl/blob/master/src/interpolations.jl#L93-L96 https://github.com/Ferrite-FEM/Ferrite.jl/blob/master/src/interpolations.jl#L192-L219 and so on. Then, everything should probably work ootb

[bhaveshshrimali]

where you have at each face in the middle sitting a DOF instead of the vertex position

Thanks @koehlerson I believe this is for the P1 non-conforming Crouzeix-Raviart element where the dofs are (only) at the centroid of each facet. The 1973 paper presents both the conforming and non-conforming elements for Stokes (which also work for incompressible elasticity). I will take a look this too. What I was referring to is it's conforming counterpart.

The terminology isn't super popular, and Crouzeix-Raviart almost always ends up being used for the non-conforming P1 element introduced in that paper.

For the conforming Crouziex-Raviart the dofs are at the vertices and edge-midpoints for a 2D case and vertices, edge-midpoints, facet-midpoints and at the element centroid for a 3D case. The basis functions for 3D can be found here and also in scikit-fem where I implemented these as a test case to run small problems here.