SCOREC
Chef makes boundary elements on periodic/matched faces
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.
I think these are the queries you are looking for:
https://github.com/SCOREC/core/blob/ba7c15c42de21efe6fac82043db4e51a5a3d9c24/apf/apfMesh.h#L351-L354
https://github.com/SCOREC/core/blob/ba7c15c42de21efe6fac82043db4e51a5a3d9c24/apf/apfMigrate.cc#L63-L69
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.
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.