Skip to content

Commit f2e45ef

Browse files
author
Lutz Gross
committed
Merge branch 'oxley-finley-export' into master
oxley stops assembling PDEs and hands the work to finley instead: the forest is converted to a conforming simplex mesh - Tri3 in 2D, Tet4 in 3D - and finley owns the overlap, the solve and the MPI. The 2:1 seams are resolved by the triangulation, so the exported space is an ordinary P1 space with no hanging degree of freedom in it, and data moves between the two domains through explicit transfer operators for all four function spaces. Both dimensions are complete, conforming and graded, on any number of ranks. The mesh carries its element, face and node tags, its tag names and its Dirac points. The split's central invariant - that two octants cut their shared face into the same triangles - is checked combinatorially on every export and across ranks, which is the only thing that can see a mismatched diagonal: a linear field lies in both triangulations of a planar quad, so no physics test would.
2 parents 0e6323c + 4c51a62 commit f2e45ef

86 files changed

Lines changed: 11368 additions & 6578 deletions

File tree

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎.gitignore‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,11 @@ test_*.log
4949
test_*.txt
5050
scons_py_tests*.log
5151
uniform_mesh_ae.silo/removed_code/
52+
# ... and what the oxley validators and the mesh IO tests write. Named rather
53+
# than a blanket *.vtu, because doc/examples ships .vtu files that ARE source.
54+
oxley_a3_*.vtu
55+
oxley_hanging_*.vtu
56+
_meshio_*
5257

5358
# design working notes, kept out of the repo
5459
oxley-designs/DESIGN.md

‎escriptcore/py_src/testing.py‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -133,6 +133,26 @@ def run_tests(modulename, classes = [], exit_on_failure = False):
133133
else:
134134
for test_class in classes:
135135
suite.addTest(unittest.TestLoader().loadTestsFromTestCase(test_class))
136+
137+
# A suite that collects nothing is not a passing suite. unittest calls an
138+
# empty result successful, so without this the script exits 0, scons writes
139+
# its .passed file and the suite is indistinguishable from one whose tests
140+
# all passed - which is how several suites went years without running.
141+
# Note this is about COLLECTING nothing: a skipped test is still collected,
142+
# reported as skipped, and does not trip this.
143+
if suite.countTestCases() == 0:
144+
if rank == 0:
145+
sys.stderr.write(
146+
"ERROR: %s collected no tests. Test classes must be named "
147+
"Test* to be discovered; if the suite is deliberately "
148+
"disabled, leave one Test class marked "
149+
"@unittest.skip(reason) so it reports as skipped.\n"
150+
% modulename)
151+
sys.stderr.flush()
152+
if exit_on_failure:
153+
MPIBarrierWorld()
154+
sys.exit(1)
155+
136156
s=unittest.TextTestRunner(stream=stream,verbosity=verb).run(suite)
137157
if exit_on_failure and not s.wasSuccessful():
138158
sys.stderr.flush()

‎escriptcore/py_src/util.py‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2816,6 +2816,22 @@ def grad(arg,where=None):
28162816
:type where: ``None`` or `escript.FunctionSpace`
28172817
:return: gradient of ``arg``
28182818
:rtype: `escript.Data` or `Symbol`
2819+
:note: A gradient on ``FunctionOnBoundary`` (or its reduced form) needs
2820+
face elements that carry the nodes of the element BEHIND the face,
2821+
not only the nodes lying on it: the derivative normal to the face
2822+
cannot be recovered from values on the face alone. Where the mesh
2823+
has only the face's own nodes the result is the TANGENTIAL part of
2824+
the gradient, with the normal component zero, and **no exception is
2825+
raised** - so a wrong answer here looks like a plausible one.
2826+
2827+
Which meshes provide it: the structured `finley` factories take
2828+
``useElementsOnFace=True``, which is their default. Meshes read from
2829+
gmsh do not - ``ReadGmsh`` has no such option, and the boundary
2830+
elements come from the file as written, carrying their own nodes
2831+
only. Building the parent-shaped faces afterwards means finding each
2832+
face's element and reordering its nodes, which is costly in 3D, so
2833+
for such a mesh this is a limitation to work around rather than a
2834+
setting to change.
28192835
"""
28202836
if isinstance(arg,sym.Symbol):
28212837
if where is None:

‎finley/src/FinleyDomain.h‎

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -85,6 +85,50 @@ enum SystemMatrixType {
8585
SMT_UNROLL = 1<<17
8686
};
8787

88+
/**
89+
\brief
90+
A mesh handed over as flat arrays, for FinleyDomain::createFromArrays().
91+
92+
Node ids, and the node ids appearing in the element tables, are GLOBAL. Each
93+
rank describes only the part of the mesh it holds: it need not supply every
94+
node its own elements refer to, and two ranks may supply the same node.
95+
createFromArrays() resolves the references, distributes the mesh and builds
96+
the parallel overlap, so a caller needs no MPI communication of its own.
97+
*/
98+
struct MeshArrays
99+
{
100+
/// spatial dimension, 2 or 3
101+
int numDim = 0;
102+
/// type of the volume elements, e.g. Tri3, Rec4, Tet4, Hex8
103+
ElementTypeId elementType = NoRef;
104+
/// type of the face elements, e.g. Line2 in 2D, Tri3 or Rec4 in 3D
105+
ElementTypeId faceElementType = NoRef;
106+
107+
/// global id of each supplied node, length numNodes
108+
std::vector<index_t> nodeId;
109+
/// coordinate d of node i at [i*numDim+d], length numNodes*numDim
110+
std::vector<double> nodeCoords;
111+
/// tag of each node, length numNodes, or empty for untagged
112+
std::vector<int> nodeTag;
113+
114+
/// global node ids of each element, length numElements*<nodes per element>
115+
std::vector<index_t> elementNodes;
116+
/// global id of each element, length numElements, or empty to number them
117+
std::vector<index_t> elementId;
118+
/// tag of each element, length numElements, or empty for untagged
119+
std::vector<int> elementTag;
120+
121+
/// global node ids of each face element
122+
std::vector<index_t> faceNodes;
123+
/// global id of each face element, or empty to number them
124+
std::vector<index_t> faceId;
125+
/// tag of each face element, or empty for untagged
126+
std::vector<int> faceTag;
127+
128+
/// tag name to tag value, copied onto the finished domain
129+
TagMap tagMap;
130+
};
131+
88132
/**
89133
\brief
90134
FinleyDomain implements the AbstractContinuousDomain interface for the
@@ -249,6 +293,30 @@ class FINLEY_DLL_API FinleyDomain : public escript::AbstractContinuousDomain
249293
bool useMacroElements, bool optimize,
250294
escript::JMPI jmpi);
251295

296+
/**
297+
\brief
298+
Creates a domain from a mesh supplied as flat arrays.
299+
300+
Each rank passes the part of the mesh it holds, with nodes identified by
301+
global id; nodes referred to by an element need not be present on the rank
302+
that supplies the element. The node and element tables are resolved and
303+
distributed here, so the caller performs no communication itself. This is
304+
the entry point used by mesh generators that live outside finley.
305+
306+
\param arrays Input - the mesh, see MeshArrays
307+
\param name Input - a descriptive name for the domain
308+
\param order Input - integration order (1 or 2)
309+
\param reducedOrder Input - reduced integration order (1 or 2)
310+
\param optimize Input - whether to optimize node/DOF labelling. Note that
311+
this repartitions with ParMETIS, so a caller that
312+
wants to keep its own partition should pass false.
313+
\param jmpi Input - shared pointer to MPI information to be used
314+
*/
315+
static escript::Domain_ptr createFromArrays(const MeshArrays& arrays,
316+
const std::string& name,
317+
int order, int reducedOrder,
318+
bool optimize, escript::JMPI jmpi);
319+
252320
/**
253321
\brief
254322
Constructor for FinleyDomain
@@ -382,6 +450,16 @@ class FINLEY_DLL_API FinleyDomain : public escript::AbstractContinuousDomain
382450
*/
383451
virtual std::string getDescription() const;
384452

453+
/**
454+
\brief
455+
returns the name this mesh was created with.
456+
457+
getDescription() answers "FinleyMesh" for every mesh, so it cannot tell
458+
two meshes apart; the name can, which matters to whoever built the mesh
459+
and later has to recognise it.
460+
*/
461+
const std::string& getName() const { return m_name; }
462+
385463
/**
386464
\brief
387465
Return a description for the given function space type code

‎finley/src/Mesh_fromArrays.cpp‎

Lines changed: 200 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,200 @@
1+
2+
/*****************************************************************************
3+
*
4+
* Copyright (c) 2003-2026 by the esys.escript Group
5+
* https://github.com/LutzGross/esys-escript.github.io
6+
*
7+
* Primary Business: Queensland, Australia
8+
* Licensed under the Apache License, version 2.0
9+
* http://www.apache.org/licenses/LICENSE-2.0
10+
*
11+
* See CREDITS file for contributors and development history
12+
**
13+
*****************************************************************************/
14+
15+
#include "FinleyDomain.h"
16+
17+
#include <escript/index.h>
18+
19+
#include <sstream>
20+
21+
namespace finley {
22+
23+
namespace {
24+
25+
/// the contact element type that goes with a given face element type
26+
ElementTypeId contactTypeFor(ElementTypeId faceType)
27+
{
28+
switch (faceType) {
29+
case Point1: return Point1_Contact;
30+
case Line2: return Line2_Contact;
31+
case Line3: return Line3_Contact;
32+
case Tri3: return Tri3_Contact;
33+
case Tri6: return Tri6_Contact;
34+
case Rec4: return Rec4_Contact;
35+
case Rec8: return Rec8_Contact;
36+
case Rec9: return Rec9_Contact;
37+
default: break;
38+
}
39+
std::stringstream ss;
40+
ss << "createFromArrays: no contact element type is defined for face "
41+
"element type " << faceType;
42+
throw escript::ValueError(ss.str());
43+
}
44+
45+
/// fills an element table from flat arrays. `nodes` holds global node ids and
46+
/// its length fixes the element count. Ids are made globally unique when the
47+
/// caller does not supply them.
48+
void fillElementTable(ElementFile* ef, const std::vector<index_t>& nodes,
49+
const std::vector<index_t>& ids,
50+
const std::vector<int>& tags,
51+
const char* what, escript::JMPI mpiInfo)
52+
{
53+
const int NN = ef->numNodes;
54+
if (NN < 1) {
55+
std::stringstream ss;
56+
ss << "createFromArrays: " << what << " reference element has no nodes";
57+
throw escript::ValueError(ss.str());
58+
}
59+
if (nodes.size() % NN != 0) {
60+
std::stringstream ss;
61+
ss << "createFromArrays: the " << what << " node table has "
62+
<< nodes.size() << " entries, which is not a multiple of the " << NN
63+
<< " nodes per element";
64+
throw escript::ValueError(ss.str());
65+
}
66+
const dim_t numElements = nodes.size() / NN;
67+
if (!ids.empty() && (dim_t)ids.size() != numElements) {
68+
std::stringstream ss;
69+
ss << "createFromArrays: " << what << " has " << numElements
70+
<< " elements but " << ids.size() << " ids";
71+
throw escript::ValueError(ss.str());
72+
}
73+
if (!tags.empty() && (dim_t)tags.size() != numElements) {
74+
std::stringstream ss;
75+
ss << "createFromArrays: " << what << " has " << numElements
76+
<< " elements but " << tags.size() << " tags";
77+
throw escript::ValueError(ss.str());
78+
}
79+
80+
// When ids are not supplied, number the elements consecutively across ranks
81+
// so that they stay unique once the mesh is distributed.
82+
//
83+
// The scan runs UNCONDITIONALLY, and only its result is conditional. It is a
84+
// collective, and `ids.empty()` is a per-rank test: a rank that simply has
85+
// none of this kind of element - no boundary faces, say, because it owns
86+
// only interior cells - supplies an empty id list too, and cannot be told
87+
// apart from a caller that omitted them. Guarding the collective with that
88+
// test let such a rank enter the scan alone while the others went on, and
89+
// the run deadlocked here with the ranks in different collectives.
90+
index_t idOffset = 0;
91+
#ifdef ESYS_MPI
92+
if (mpiInfo->size > 1) {
93+
index_t local = numElements;
94+
index_t scan = 0;
95+
MPI_Exscan(&local, &scan, 1, MPI_DIM_T, MPI_SUM, mpiInfo->comm);
96+
if (ids.empty() && mpiInfo->rank != 0)
97+
idOffset = scan;
98+
}
99+
#endif
100+
101+
ef->allocTable(numElements);
102+
ef->minColor = 0;
103+
ef->maxColor = numElements > 0 ? numElements - 1 : -1;
104+
105+
#pragma omp parallel for
106+
for (index_t e = 0; e < numElements; e++) {
107+
ef->Id[e] = ids.empty() ? (idOffset + e) : ids[e];
108+
ef->Tag[e] = tags.empty() ? 0 : tags[e];
109+
ef->Owner[e] = mpiInfo->rank;
110+
ef->Color[e] = e;
111+
for (int k = 0; k < NN; k++)
112+
ef->Nodes[INDEX2(k, e, NN)] = nodes[e * NN + k];
113+
}
114+
}
115+
116+
} // anonymous namespace
117+
118+
escript::Domain_ptr FinleyDomain::createFromArrays(const MeshArrays& in,
119+
const std::string& name,
120+
int order, int reducedOrder,
121+
bool optimize,
122+
escript::JMPI mpiInfo)
123+
{
124+
if (in.numDim != 2 && in.numDim != 3) {
125+
std::stringstream ss;
126+
ss << "createFromArrays: numDim is " << in.numDim << ", must be 2 or 3";
127+
throw escript::ValueError(ss.str());
128+
}
129+
if (in.elementType == NoRef || in.faceElementType == NoRef)
130+
throw escript::ValueError("createFromArrays: element type and face "
131+
"element type must both be set");
132+
133+
const dim_t numNodes = in.nodeId.size();
134+
if (in.nodeCoords.size() != (size_t)numNodes * in.numDim) {
135+
std::stringstream ss;
136+
ss << "createFromArrays: " << numNodes << " nodes in " << in.numDim
137+
<< " dimensions need " << (size_t)numNodes * in.numDim
138+
<< " coordinates, got " << in.nodeCoords.size();
139+
throw escript::ValueError(ss.str());
140+
}
141+
if (!in.nodeTag.empty() && (dim_t)in.nodeTag.size() != numNodes) {
142+
std::stringstream ss;
143+
ss << "createFromArrays: " << numNodes << " nodes but "
144+
<< in.nodeTag.size() << " node tags";
145+
throw escript::ValueError(ss.str());
146+
}
147+
148+
FinleyDomain* out = new FinleyDomain(name, in.numDim, mpiInfo);
149+
150+
const_ReferenceElementSet_ptr refElements(
151+
new ReferenceElementSet(in.elementType, order, reducedOrder));
152+
const_ReferenceElementSet_ptr refFaceElements(
153+
new ReferenceElementSet(in.faceElementType, order, reducedOrder));
154+
const_ReferenceElementSet_ptr refContactElements(
155+
new ReferenceElementSet(contactTypeFor(in.faceElementType), order,
156+
reducedOrder));
157+
const_ReferenceElementSet_ptr refPoints(
158+
new ReferenceElementSet(Point1, order, reducedOrder));
159+
160+
ElementFile* elements = new ElementFile(refElements, mpiInfo);
161+
out->setElements(elements);
162+
ElementFile* faces = new ElementFile(refFaceElements, mpiInfo);
163+
out->setFaceElements(faces);
164+
out->setContactElements(new ElementFile(refContactElements, mpiInfo));
165+
out->setPoints(new ElementFile(refPoints, mpiInfo));
166+
167+
// node table. The global id doubles as the degree of freedom: unlike the
168+
// structured generators there is no periodicity to fold away here.
169+
NodeFile* nodes = out->getNodes();
170+
nodes->allocTable(numNodes);
171+
#pragma omp parallel for
172+
for (index_t i = 0; i < numNodes; i++) {
173+
nodes->Id[i] = in.nodeId[i];
174+
nodes->Tag[i] = in.nodeTag.empty() ? 0 : in.nodeTag[i];
175+
nodes->globalDegreesOfFreedom[i] = in.nodeId[i];
176+
for (int d = 0; d < in.numDim; d++)
177+
nodes->Coordinates[INDEX2(d, i, in.numDim)] =
178+
in.nodeCoords[(size_t)i * in.numDim + d];
179+
}
180+
181+
fillElementTable(elements, in.elementNodes, in.elementId, in.elementTag,
182+
"element", mpiInfo);
183+
fillElementTable(faces, in.faceNodes, in.faceId, in.faceTag,
184+
"face element", mpiInfo);
185+
out->getContactElements()->allocTable(0);
186+
out->getPoints()->allocTable(0);
187+
188+
for (TagMap::const_iterator it = in.tagMap.begin();
189+
it != in.tagMap.end(); ++it) {
190+
out->setTagMap(it->first, it->second);
191+
}
192+
193+
// resolve the global node references, then distribute and build the
194+
// overlap, the mappings and the element colouring
195+
out->resolveNodeIds();
196+
out->prepare(optimize);
197+
return out->getPtr();
198+
}
199+
200+
} // namespace finley

‎finley/src/SConscript‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,7 @@ sources = """
4949
IndexList.cpp
5050
Mesh_addPoints.cpp
5151
Mesh_findMatchingFaces.cpp
52+
Mesh_fromArrays.cpp
5253
Mesh_getPasoPattern.cpp
5354
Mesh_getTrilinosGraph.cpp
5455
Mesh_glueFaces.cpp

0 commit comments

Comments
 (0)