Core support of quadratic elements

[KennethEJansen]

In the early mid 90s to mid 2000's SCOREC was ground zero for the development of pde solvers using the hierarchic basis (HB). Quite a bit of software was written that used them but it looks like most of that did not come into core. Below is a summary of email dialogs on what SCOREC/core does and does not support along with a wish list for HB.

Jansen: Does core have any support yet for higher order fields?

Specifically, for the hierarchical basis, attaching solution coefficients to edges as was done in NSpre?

Here is the wish list and maybe you can tell me what the prospects are:

1) chef prepares hierarchical basis runs for phasta which requires: i) expanding the notion of connectivity to include global mode numbers for all edges ii) expanding the notion of connectivity to have the ien array know the global edge mode number for each edge of the element (ien(1:numel, 1:10) for tets or ien(1:numel, 1:20) for hexes). iii) adding boundary conditions for these global edge modes iv) adding initial conditions for these global edges.
v) obviously write it all out Ideally we can rob from NSpre the code that does all of these things (note iii) and iv) are solved as local problems on the entity so it is not all that much harder to do than what is currently done for nodes from the Simmetrix smd file expressions).

2) in the hay day of HB work we also had code that would subsample and HB basis restart file + geom.sms to a given refinement level so that we could visualize the higher order fields in ParaView. Has anybody done anything like that yet with SCOREC/core? Again, I think I can find the old Simmetrix-based code and I suppose if we bring it back to serial or low process counts (it probably only worked in serial before) then either update this code OR we can port it into SCOREC/core to work directly with the mds file (preferred because the stream approach avoids Simmetrix entirely)

[KennethEJansen]

Morteza's response:

I don't quite understand what is needed here but here is the summary of what can actually be done in CORE in regards to higher order fields at the moment.

1st and 2nd order Lagrange fields on a mesh can be defined using Field createField (Mesh m, const char name, int valueType, FieldShape shape)

FieldShape* apf::getLagrange ( int order )   Using order 2 would enable storing one value associated with each edge.

For going to orders higher than 2 one needs to use the bezier shapes. Creating a bezier field is done in the same way as above except that for the shape parameter in createField one has to use apf::FieldShape* crv::getBezier(int order )

where order could be anyting upto 6.

As for the visualization in paraview, I have never tried visualizeing any higher order fields except for the coordinate field, and as you mentioned that is done by subsampling points on the actual curved entities and basically creaiting a very fine linear mesh that paraview can understand.

I should mention that paraview can visualize quadratic meshes. In other words, if you have a mds mesh with quadratic Lagrange shapes you sould be able to just visualize it in paraview by wrting is to VTK using this api void apf::writeVtkFiles ( const char prefix,   Mesh m,   int cellDim = -1   ) This api woulld also write all the fields to the VTK files given they are "complete". The "complete" in this contex means that for quadratic fields there must be a value associated with each vertex and also a value associated with each edge in the mesh.

[KennethEJansen]

Jansen Reponse: Thanks for the explanation and the pointers.

You mention for quadratic, Lagrange is used and one dof is added per edge. From my experience, Lagrange adds an interpolatory NODE to the edge AND it changes all the node functions to be quadratic and interpolatory — that is, for a tet, the 4 original, linear nodal functions are discarded and replaced by 4 new quadratic functions that together with the 6 new “edge” functions are 1 at their entity and zero at all others. Hierarchical basis functions preserve the four linear nodal functions, add 6 quadratic functions that are not interpolatory. It is obviously possible to map between the solution coefficients since they are the same function space but the 6 edge coefficient values will be different so it is important for me to be sure I understand what you mean by Lagrange fields. Is my understanding that getLagrange elevates a linear mesh to a globally quadratic mesh that is equivalent in some sense to adding node/modes to each edge (I guess you would assume it is added at the midpoint) of the mesh AND elevating the nodal basis to quadratic?

Can you comment on whether the functions you are referring to support element shapes beyond tetrahedra (e.g., whether they could be used for hexes)? For hexes, adding a mode to each edge does make it a complete cubic but in this case is is what is sometimes referred to as a serendipity element since it lacks the by quadratic mode that would otherwise appear on each hex face and the tri-quadratic mode that would appear on each cell center from a standard, tensor product, Lagrange basis. Again, hierarchical basis are well defined on hexes with either just the edge modes or adding more if you like with the same additive construction that leaves the vertex modes as linear.

I suppose at this point, Cameron is going to say we should have made a ticket out of this which is of course correct. I will do that Cameron but can you invite John Evans and Corey Nelson to the SCORE/core ticket process?

[KennethEJansen]

Morteza Response: You are right. For Lagrange shapes, going from order 1 to 2 will replace all the linear shape functions at the vertices with quadratic ones, and they are also interpolatory as you mentioned.

In general, all the non-bezier shape functions that we support in CORE are defined in apf/apfShape.cc We have "Linear" shapes for edges, tris, quads, tets, prisms, pyramids, and hexes https://github.com/SCOREC/core/blob/0f41f64649f82fd5daf3c57cfd2d413903ac70ee/apf/apfShape.cc#L71 We have "Quadratic" shapes for edges/tris and tets https://github.com/SCOREC/core/blob/0f41f64649f82fd5daf3c57cfd2d413903ac70ee/apf/apfShape.cc#L356 We have "LagrangeQuadratic" shapes for quads with 9 nodes https://github.com/SCOREC/core/blob/0f41f64649f82fd5daf3c57cfd2d413903ac70ee/apf/apfShape.cc#L477 We have "SerendipityQuadratic" for quads with 8 nodes https://github.com/SCOREC/core/blob/0f41f64649f82fd5daf3c57cfd2d413903ac70ee/apf/apfShape.cc#L576 We have "LagrangeCubic" (I was not aware of this until now!) for edges, tris, and tets https://github.com/SCOREC/core/blob/0f41f64649f82fd5daf3c57cfd2d413903ac70ee/apf/apfShape.cc#L646

So to answer your second question, based on the above list it seems that we do not support anything beyond first order for hexes.

[KennethEJansen]

Jansen (last time I do this because now author implies who posted): Thanks for digging out this index of what is there. If you are familiar with the process of adding these elements, can you comment on whether the infrastructure is in place to support the following "elements": 1) Serendipity hex (8 nodes on vertices and 12 "nodes" on edges) This is the quadratic hex element supported by ParaView by the way. 2) HB-p2-Hex (same as 1 but in the hierarchic basis where we would be leaving the 8 vertex modes as they are in class Linear for hex and the quadratic modes are added to each edge following Dey's thesis which is SCOREC report 1997-1 3) HB-p2_TP-Hex (same as 2) but adding bi-quadratic mode on each face of the hex and tri-quadratic function that is non-zero only on the volume)

I am looking for advice as to whether, after we add these to field shape we would then be able to, with modest effort, accomplish the objectives listed in the opening comment regarding pre-processing HB cases for PHASTA AND develop code that would do a series of mesh subdivisions (and yes we are talking about hexes here so I know there would probably either need to be some development here OR at least make some code that tetrahedronizes the hex mesh and then does subdivisions of that tetrahedronized mesh). Note this is all mute if we conclude that adding those three higher order hex elements to core is a major challenge.

Given that we only want to do this on structured meshes, one might imagine that it is not a very difficult task to do all of the chef-type prep based on linear hexes and then do all the higher order mode construction in PHASTA, possibly with the aid of a few helper data structures (e.g., global face number connectivity on each element, global edge number connectivity on each element).

[mortezah]

I am not very familiar with those elements, but in general, as long as one has the expressions for all the shape functions (and their derivatives) in terms of the parametric coordinates of the (parent) entity, then it should be relatively straight-forward to add those entity shapes to the code.

For example, the LagrangeCubic edges, tris, and tets seem to be the entity shapes implemented most recently in the code. See the following commits:

https://github.com/SCOREC/core/commit/7d1bdcd02a117c1b72fb59ea59f6791912c3e11b https://github.com/SCOREC/core/commit/20f5295cd01a6fe47d4e13bacc904ddb44b30e9b

All that is done, in these commits is to add the expressions for the shape functions and their derivatives.

You should be able to follow the steps in the commits mentioned above to implement the desired entity shapes.

[KennethEJansen]

Thanks for that feedback. Cameron, can you comment on any known obstacles beyond adding these element types and the appropriate shape functions and derivatives (which we already have coded by the way within the /phasta)?

I realize this is a loaded/ambiguous question but I am just trying to get a sense of what additional obstacles we would hit with bringing higher order hexes through the chef process outlined in the first message. I think expanding FieldShape will go very quickly for Corey as he is a very competent C++ programmer and he has already worked with the HB functions. He is much less familiar with the rest of the steps chef does (again, I think all the change "deltas" are described in the issue-opening comment) and so I am by this line of questioning trying to map out what major obstacles experienced folks might know about or get a forecast for smooth sailing if that is what we expect.

I can free up a half hour to walk the code if that helps.

Again, we have "working" NSpre (built off of current Simmetrix libraries) that does the prep work described in opening comment 1). We have old and moldy stuff that does the post stuff in comment 2).

[cwsmith]

@KennethEJansen I'm not very familiar with the fields portion of the code and, at best, could only speculate on possible issues. It sounds like if the examples @mortezah has pointed out map to your use things should be OK.

Do you need to adapt meshes with these HB fields?

[KennethEJansen]

No plan for h adapt (near term at least). The HB is great for p-adapt though.

[cwsmith]

I'm going to close this as it seems the discussion has settled. Please reopen as needed.