Implementing Crouzeix-Raviart (non-conforming) interpolation

[bhaveshshrimali]

This follows #363 (comment)

Just reminding myself to pick on this thread up soon :) @termi-official

Thanks

[termi-official]

I also looked into how to introduce non-conforming elements into Ferrite (especially DG) - it should be possible with some moderate modifications to __close! . The problem here is that the used algorithm assumes 1:1 correspondences between shared entities and the dofs. Technically we could be fine by exchanging Dict{Int,Int} with something like Dict{Int,Array{Int}} for NC element formulations on conforming geometries and add some additional routines to determine whether a dof is continuous. I am not super confident that this will also do the correct thing for NC geometries (in combination with NC elements).

Edit: It should also be noted that this may degrade performance, because I do not see a way to reuse ht_keyindex2! Nvm. Derped here.

[koehlerson]

Maybe the *Index structs should be incorporated there as well to be more consistent?

[termi-official]

While I agree conceptually it does not work. The *Index structs refer to local entities as for example half-faces (e.g. cell id-local side pairs). However, when distributing dofs we need global entities (e.g. for faces the sorted vertex indices).

[termi-official]

Probably this struct should be renamed (do not have a good name here, maybe Local*Index and introduce something like a Global*Index?) and better documented. However, Global*Index is then again ambiguous in the case of NC geometry...

[koehlerson]

ah nvm it's not a good idea :D

[bhaveshshrimali]

@koehlerson @termi-official
Now that #363 is merged, is there a way to, say, assemble the jump of a field over the internal facets of a mesh? Something like dot(jump(u), jump(v)) and dot(jump(du), jump(v)) where u, v and du are the function, test function and trial function respectively.

I will look at the code in some detail tomorrow and try to construct a working example.

[koehlerson]

A sketch is given in https://github.com/Ferrite-FEM/Ferrite.jl/blob/master/test/test_grid_dofhandler_vtk.jl#L362-L393 - note that this code currently just works if either the gradienf field is constant or if you have 1 quadrature point per face, because the quadrature points do not match. Some adjustments have to be made to match these points.

I think, in most cases, the adjustment for higher order fields should be rather
easy. If the field is not constant the faces will be flipped, i.e. FaceIndex(2,2) and FaceIndex(3,4) are neighbors and the first face in global representation could be (3,7) and the latter (7,3). So, it should be sufficient to loop over the gauss points in the cell by 1:nquadpoints and zip another index variable with nquadpoints:-1:1 for the neighbor face gauss point stuff.

[termi-official]

Note that Maxi's statement is, in general, only true for 2D. In 3D we have also to account for the rotation of the faces. I.e. the triangle with vertices 1,2,3 can be described by the sequences:

1-2-3
2-3-1
3-1-2
1-3-2
3-2-1
2-1-3

However, although I do not recommend this approach, you should be able to "unrotate" the faces in preprocessing step to overcome this issue. I have an idea for a cleaner solution to this problem, but won't find any time soon to implement it anytime in the next weeks.

Edit: In case someone wants to fix this - my idea for a clean solution is to add a helper routine which maps the integration points to the global face representation, which in turn defined by the sortface function.

[koehlerson]

yes thats true, but I think grid_generator grids will be fine. Depends really on what you are planning to do @bhaveshshrimali

[bhaveshshrimali]

Thanks Maxi and Dennis!

I can try things with a grid_generator generated grid for now in 2D. I will take a look at 3D later since a 3D-CR interpolation needs to be added to Ferrite anyway right? The goal right now is to reimplement something like the incompressible hyperelasticity example but with these CR interpolations instead of mixed (disp-press) interpolation, say a trimmed down version of https://github.com/bhaveshshrimali/FEniCS_Shrimali_Ghosh_Kumar_Lopez-Pamies

Once I get down to the details, I will update this discussion..

[termi-official]

I am not aware that anyone has tried implementing 3D CR elements yet. Adding these will be a valuable contribution, because they can also be a template for deriving an interesting class of error estimators and a discontinuous Galerkin example! Happy to hear updates on this.