Description
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...
2
57D9
61 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.