Supplying element-neighbor pairs for (periodic) interface calculations

[ttruster]

I am planning some implementations regarding periodic boundary conditions. Though there are several existing methods for them in MOOSE, I feel that adding the weak enforcement using surface integrals has some benefits not yet realized.

To start out, I studied the behavior of input files having BreakMeshByBlockGenerator. If the Exodus mesh is written out after calling this command, and then a new input file reads this mesh, then the computed result is different. As mentioned in the discussion here, #18240 (reply in thread), that is because a 'trick' is used in BreakMeshByBlockGenerator whereby the nodes are duplicated and the nodes-to-element connectivity is updated, but the element-neighbor information is not updated (i.e. is retained) internally. So it makes sense that if the modified mesh is written, though it contains a Side-Set, that set is one-sided and lacks the neighbor-element pairing required to recreate the internal state of the memory of the previous analysis.

My question is, is the only information I would be required to read in through a text file, the element-neighbor pairing (element ID and faces)? Or is there some other orientation information needed to update the Libmesh quadrature point projections? @arovinelli
suggested that only this information would be needed, and gave 4 lines of code for it. If this sounds like the approach, then I'll try hard-coding those lines within a new Mesh-Generator object, and see if that allows the read-in mesh to provide the same result as the original BreakMeshByBlockGenerator result.

Alternatively, I saw a discussion here, #18580 (reply in thread), about Libmesh periodic features. Are those periodic listings node-based (like multi-point conditions) or face-based? From this discussion #20232, it sounds like they are face-based. I can't easily tell from the source code https://mooseframework.inl.gov/docs/doxygen/libmesh/periodic__boundaries_8C_source.html#l00038 what kind of periodic treatment Libmesh is doing (nodal, penalty, mortar, Nitsche); so suggested reading would help here as well.

[ttruster]

Question at the bottom about inverse_map.

So I was able to verify that reassigning the element-neighbor pair through a Mesh-Generator class was able to restore the connection needed for the cohesive simulation.
Elem * current_elem = mesh->elem_ptr(elem_id);
Elem * neighbor_elem = mesh->elem_ptr(neighbor_id);
current_elem->set_neighbor(elem_side, neighbor_elem);
neighbor_elem->set_neighbor(neighbor_side, current_elem);
Thus, maybe on day in the future, we could add a feature to BreakMeshByBlockGenerator to export a CSV list of this data while running, and having another mesh class that reads it back in.

On the second topic for periodic BC, I spent an hour looking through the code for how neighbor_data is initialized, and I finally found within the Assembly class the reinitElemAndNeighbor function and within it the Libmesh FEInterface::inverse_map call. That's what I expected to find, because somehow the neighbor quadrature point needs to know which element quadrature point to link up with. The face alone isn't enough; there needs to be some projection/lookup/reorientation that is checked.

Is anyone familiar with the source code of inverse_map who can discuss how I would add the calculation of a translate/offset during the projection to handle periodic pairs that are parallel but translated (i.e. not adjacent like in usual meshes with Interface Kernels)? I can imagine doing this from the average coordinate of the corner nodes or of the quadrature point coordinates, such that it is internal and doesn't need any change to the argument list.

[ttruster]

The PeriodicBoundary may be associated with a set of variables

That's even better, because for cohesive zone (CZM) stuff, they would be the same for disp_x, y and z; so giving it once would be convenient.

I don't know the implementation of these details, only the interfaces to them. But supposing we could add an enum STRONG_ENFORCEMENT/WEAK_ENFORCEMENT member to PeriodicBoundary and I made a copy of AddPeriodicBCAction which sets that member to weak and doesn't make a constraint matrix for the variable (or maybe add the option for week as a parameter in that AddPeriodicBCAction class instead of a copy of it). Then in the onSeparatedInterface loop I could do a proof of concept if you let me know what the interface is to look up a variable and its PBC from the System. Namely, the proof of concept would be to recover the result of a two-element CZM analyses from a restart, which currently isn't possible.

[lindsayad]

or maybe add the option for week as a parameter in that AddPeriodicBCAction class instead of a copy of it

I think the parameter would be nicer than having a copy

[lindsayad]

You can get periodic boundaries from the DofMap. You can get the DofMap from the libMesh System using get_dof_map(). You can see what variables a given PeriodicBoundary(Base) is associated with using its get_variables() API (see periodic_boundary_base.h).

[ttruster]

@roystgnr and @lindsayad: I finally got some time to work on this within the past week. The source code is here https://github.com/Truster-Research-Group/moose/tree/topological-neighbor-PBC and https://github.com/ttruster/libmesh/tree/topological-neighbor-PBC. Perhaps we can open an issue to continue the discussion of this? I had also heard from @reverendbedford that the discussion of stateful mortar constraints was continuing, so I noticed a new issue was started here for that #24030.

Summary: the implementation is crude (no error checks) with the following components along the ideas mentioned in posts.
Libmesh: I added an enforcement type for STRONG_ENFORCEMENT/WEAK_ENFORCEMENT to the periodic boundary class and the AddPeriodicBCAction class in MOOSE. There is a set/get method pair for it. Then in fe_base.C I added a conditional to not create a DOF constraint matrix for weak PBC. Lastly in elem.C, I added a copy of topological_neighbor which also sets an integer pointer for the facet of the neighbor on the PBC, which is passed up from the neighbor function in the periodic_boundaries object.
MOOSE: Duplicated the method onInterface into onPeriodicBoundary within the ThreadedElementLoopBase.h class to be called for facets that don't have a neighbor. Implemented the method similarly in both NonlinearThread and ``ComputeMaterialsObjectThreadbut instead callingtopological_neighbor` (once) to determine the neighbor and the facet, and then using that data for the reinitialize methods. Modified the `reinitializeNeighbor` methods in `Assembly, FEProblemBase, DisplacedProblem`, and `Subproblem` to use that data and to compute an offset vector from element to neighbor side which allows for a coordinate offset between those two sides before calling `inverse_map`.

Both branches have been pushed with the updates. There is an example file bilinear_topo.i in the CZM test folder of tensor mechanics module that produces identical behavior to bilinear_split.i but uses a .e file for the mesh input and a Periodic BC for letting the two sides be connected. A refined mesh with 4 elements along the interface also matches (though the Newton algorithm fails at step 4 or later, in both cases). The --mesh-only option needs to be called on the original input file to get the _in file referenced in the test input file. The key to determine if the result is correct is to check the transverse displacement disp_x. If that field is zero, then the two sides are not "talking to each other", meaning no physics is coming from the Interface Kernel.

Some of the discussions to occur in the issue include

• Error checks to verify which variables are periodic and that they are attached to the interface kernel
• A test of whether any PBC exist in the model, to test at the top of onPeriodicBoundary
• Devise a caching scheme in MOOSE (seems better than libmesh) to list the order that the PBC pairs will be integrated, so that the results from topological_neighbor can be stored instead of recomputed at each nonlinear iteration. Perhaps a sorted table based on element and side could be used, so that a binary search could be used?
• any other concerns with my crude implementation

Reminder about the benefits for this (option 2 of the 3 tier hierarchy, where stateful mortar is the complicated one for nonconforming meshes): only involves changes to 5 libmesh files while retaining default behavior if no flags are set, allows existing Interface Kernels to be used in MOOSE for offset boundary problems with no change to the input file, and works within the AddPeriodicBCAction. The drawbacks are: the neighbor data is currently recalculated at every iteration for no other reason than we are accessing the topology of the mesh in a roundabout way through field variables and DOF map tables; and relatedly that the DOF map (maybe, unless you tell me where) doesn't exist for aux variables, or at least there is not a _nl member for User Objects, and so the Interface user objects can't involve this offset option.

[lindsayad]

There is a DofMap for every system, so we can work on that. Yes I think it would be good to create an issue or even better probably to get the ball rolling would be to open a PR with your libMesh changes.

[lindsayad]

Do you not want to use mortar?

[ttruster]

Good question! So, I had forgotten my notes from the summer when talking with @reverendbedford that the Mortar objects don't allow stateful materials. This is an issue for applying PBC using the Nitsche/DG method for Tensor Mechanics using the Mortar approach. So I was developing for the heat conduction/diffusion case using the mortar objects up to now (which you're familiar with that), and then we tried porting to Tensor Mechanics and hit the stateful error message.

I can add more details later. What the two paths are for getting to our formulation goal is:

1. Add stateful materials to mortar objects in the case when quadrature points are fixed for all time (the case for mesh tying such as PBC; NOT the case for contact mechanics). Then add the ability to have a vector offset for the search/projection direction between primary and secondary sides; this is needed for non-rectangular RVE domains.
2. Use interface kernels for the PBC Nitsche method and Lagrange multiplier methods, which work for stateful materials. This requires the reading in of element-neighbor pairs (demonstrated above) and the ability for the pair to not be physically adjacent. Obviously the mortar method functions already do that for quadrature points (using projection along the unit outward normal) since the surfaces can be separated. So I just need to access the inverse_map in a similar fashion.

It seemed to me that, option 2 would require less implementation (as I have existing codes that generates all the needed element-neighbor pairs; just need a MeshGenerator that reads a CSV file of those and updates the mesh). And I can use the diffusion mortar PBC implementations to check the new implementation.

[hugary1995]

I had a very similar discussion with Alex a while ago actually. 1 is tough, quoting Alex. Though it might be easier if we assume the mortar mesh never changes. Plus it doesn't require libmesh level changes, unlike option 2 which requires modifying inverse_map().

I can see how 2 might work as well. Since implementing a mesh generator which enables restarting a CZM simulation also helps many of our simulations, I will also be supportive of option 2.

[lindsayad]

Let's try to respond within relevant threads to keep the discussion thread clear. For cases where the mesh can change, I've often felt that the current way we do stateful material properties is a little disappointing. Yes we have an implementation for projection of material properties, but this would be very natural if the "stateful properties" were finite element approximations. Then we have incredibly generic projections (the libMesh GenericProjector for example) immediately at our disposal. Additionally within the mortar method, even if the mortar segments change, the ability to evaluate finite element bases connected to the parent mesh is flawless. So aux variables for instance work "flawlessly" with mortar, and you have access to old and older states (up to arbitrarily old states). I quote "flawlessly" because it's not quite flawless when you want to do AD as we do not yet have first class capability for storing aux variables with AD information

[lindsayad]

For problems in which you do not have mesh adaptivity or something like mortar, the way we do stateful material properties is great... you do not introduce any error by projecting your material property evaluations into a finite element basis approximation.

[ttruster]

I started discussing the mortar case sense you brought it up :) but my other problem is that I missed the "Write a reply" box.

The mortar system was helpful for testing out the weak imposition of periodic conditions, but it is more general than our typical use case which will be conforming mesh tying across the external surfaces. Mortar would allow for non-conforming PBC, but is a bit more expensive due to interface segments and also has a different parallelization compared to Interface Kernels. Also, the topic of stateful materials for constitutive models in tensor mechanics is a lengthy discussion which I agree is off topic here.

So, I will up-vote my response above about the element-neighbor pairing and comment on the inverse_map below, which can be moved to a separate discussion and issue.

[ttruster]

I'm not sure how inverse_map() works, but if it has access to either the coordinates of 1. nodes of the face or 2. the quadrature points of the face or else 3. the centroid of the face, then the easy way to compute the offset is just to average them and subtract them. That can be done at EVERY interface, since for adjacent ones the offset will be (0,0,0). That way, no flags have to get passed in, at the expense of some vector additions/divisions for each element pair.

[hugary1995]

The inverse map solves a normal equation to locate a point in the parametric space. So 1. Yes. 2. It could be any point, not limited to qps. 3. I'm sure that information is indirectly available somehow.

But inverse_map() is a performance critical low level routine, you'll probably have to overload it, not modifying the existing implementation which adds any computational cost.

[ttruster]

I found the source code with the source code in FEMap::inverse_map() that solves the normal equation you mentioned. So if I look up the data structure of element_data and neighbor_data, I will look for a low-impact way to incorporate the change, which will likely involve a flag and might get handled either in inverse_map or even nicer if in Assembly::reinitElemAndNeighbor.