Chef makes boundary elements on periodic/matched faces

[KennethEJansen]

Perhaps this has always been there for chef but I don't think it was there in the predecessor that chef was built to replace (phNSpre). Periodic faces are not (should not be) part of the boundary that the FEM should evaluate boundary integrals arising from volume terms that were integrated by parts. We match the mesh on periodic faces and think of elements that contribute to nodes on that boundary are connected (assembled) to the nodes on the matched face. Computing the boundary flux from a periodic solution (e.g., one that has been matched as the periodic BC does maintain) will cancel under exact arithmetic but this un-necessary work and un-necessary round off error. Cutting from phoOutput.cc

247   int boundaryDim = m->getDimension() - 1;
 248   apf::MeshEntity* f;
 249   apf::MeshIterator* it = m->begin(boundaryDim);
 250   while ((f = m->iterate(it))) {
 251     apf::ModelEntity* me = m->toModel(f);
 252     if (m->getModelType(me) != boundaryDim)
 253       continue;
 254     if (getBCValue(m->getModel(), bcs.fields["DG interface"], (gmi_ent*) me) != 0){
 255       apf::DgCopies dgCopies;
 256       m->getDgCopies(f, dgCopies);
 257       if (dgCopies.getSize() == 1) // This prevents adding interface elements...
 258         continue;
 259     }
 260     if (m->countUpward(f)>1)   // don't want interior region boundaries here...
 261       continue;
 262     gmi_ent* gf = (gmi_ent*)me;
 263     apf::MeshEntity* e = m->getUpward(f, 0);

We see that interior faces are skipped on line 260 (as they must be), and on 254-259 DG interfaces are skipped (as they must be) but periodic faces are not skipped (as they should be).

I am looking over the code in core/phasta to try to see if there are attached flags or other data structures that can be checked to escape the boundary element creation in a similar way. This looks promising but I don't understand it well enough to be sure

266 
267 static SavedMatches* savedVertexMatches = 0;
268 static SavedMatches* savedFaceMatches = 0;
269 
270 void enterFilteredMatching(apf::Mesh2* m, Input& in, BCs& bcs)
271 {
272   if (!in.filterMatches)
273     return;
274   savedVertexMatches = new SavedMatches();
275   saveMatches(m, 0, *savedVertexMatches);
276   if (in.formElementGraph) {
277     savedFaceMatches = new SavedMatches();
278     saveMatches(m, 2, *savedFaceMatches);
279   }
280   ModelMatching mm;
281   getFullAttributeMatching(m->getModel(), bcs, mm);
282   filterMatching(m, mm, 0);
283   if (in.formElementGraph)
284     filterMatching(m, mm, 2);
285 }

from phFilterMatching.cc . Thinking about this a bit more, I recall the following history. It used to be that we would request periodicity as an analysis attribute. I see the code still exists in core to do that. Then, suddenly, some body (have not searched git blame) decided it should ALWAYS be applied when users generated a mesh with the meshing attribute matched mesh. I protested saying that sometimes, we run multiple analyses on the same mesh and we might compare periodicity to symmetry (on both faces that were previously linked with periodicity). I lost the battle to revert back to having it be an analysis attribute but, to enable us to still do non-periodic cases with matched meshes someone created the adapt.inp input formElementGraph. I have no idea what that means but I think it causes the matching attribute to be filtered as we see in this code.

If somebody sees/understands enough to help me identify faces that have matching to escape the loop I would appreciate help. You can't be too explicit because I am not trained in C.

[cwsmith]

I think these are the queries you are looking for:

core/apf/apfMesh.h

Lines 351 to 354 in ba7c15c

	/** \brief return true if the mesh has matched entities */
	virtual bool hasMatching() = 0;
	/** \brief get the periodic copies of an entity */
	virtual void getMatches(MeshEntity* e, Matches& m) = 0;

core/apf/apfMigrate.cc

Lines 63 to 69 in ba7c15c

	if (m->hasMatching())
	{
	Matches matches;
	m->getMatches(adjacent[i],matches);
	for (size_t j=0; j < matches.getSize(); ++j)
	PCU_COMM_PACK(matches[j].peer,matches[j].entity);
	}

[KennethEJansen]

If I understand these correctly, hasMatching is true if there are any matched entities. I guess I can see how a query on getMatches for the current face in that loop will give me its matches and if not zero I will want to skip that face (exempt it from the boundary element creation).

This raises another question. Are matches bi-diretional? that is if face 27 answers getMatches with a non-zero match and say that is face 289, does the same query on face 289 answer getMatches with a non-zero match and list face 27?

If yes, I can exempt simply by finding a non-zero match. If not, I guess I will need to traverse the faces once to generate a second list of pointed to matches to check for exemption as well since I need the boundary element skipped for both.

I guess I also need to be sure that matching is not limited to nodes but is also known to edges and faces. Nodes can have multiple matches as can edges (when multiple directions are matched) but faces can have only one but obviously, the above plan relies on faces answering that query correctly when e is a mesh face entity.

[KennethEJansen]

That change compiled fine with the following code insertion

diff --git a/phasta/phOutput.cc b/phasta/phOutput.cc
index e27e95bc..28d52f4e 100644
--- a/phasta/phOutput.cc
+++ b/phasta/phOutput.cc
@@ -257,6 +257,13 @@ static void getBoundary(Output& o, BCs& bcs, apf::Numbering* n)
       if (dgCopies.getSize() == 1) // This prevents adding interface elements...
         continue;
     }
+    if (m->hasMatching()) 
+    { 
+      apf::Matches matches; 
+      m->getMatches(f,matches); 
+      if(matches.getSize()>0)  // periodic boundary faces should not be in boundary element list
+        continue;
+    } 
     if (m->countUpward(f)>1)   // don't want interior region boundaries here...
       continue;
     gmi_ent* gf = (gmi_ent*)me;
diff --git a/pumi-meshes b/pumi-meshes
index 73950fdf..635fff26 160000
--- a/pumi-meshes
+++ b/pumi-meshes
@@ -1 +1 @@
-Subproject commit 73950fdfee8547925a60dc3f456767d52a4e80ad

But I am not optimistic about this working for the workflow that motivated it -- matchedNodeElmReader. Recall that workflow does not involve Simmetrtix but instead streams in coordinates, connectivity, and NODAL matches. Unless there is code to compute matches of edges and faces inferred from nodal matches (together with the topological model info), I am expecting to not find correct answers to a getMatches query on a face entity.

I suppose a work around to this is to get the vertices of the current face and see if they ALL have matches. I say ALL because mesh faces that contact a model EDGE will have two candidate boundary element faces. WE do not want to exempt the face that has 2 nodes with matches and one that does not because that face only has an edge on the matched boundary and the other node(s) complete a face on boundary that needs a boundary element (non-matched).

Having said all of this there could be a problem with tri-faces (from tets) in cases with two matched boundaries. The right way to do this is to ask if all the nodes of the face are classified on the classified on the closure of the model face that mesh face is classified on.