From 7c14c53025443696c9d8ae1c251a73e94c70785c Mon Sep 17 00:00:00 2001 From: Marc Day Date: Mon, 3 Aug 2026 10:53:15 +0200 Subject: [PATCH 1/4] Fix constructor declarators in SDFGen vec.h --- Tools/SDFGen/vec.h | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/Tools/SDFGen/vec.h b/Tools/SDFGen/vec.h index d846af68..455d3cce 100644 --- a/Tools/SDFGen/vec.h +++ b/Tools/SDFGen/vec.h @@ -23,45 +23,45 @@ struct Vec { T v[N]; - Vec(void) + Vec(void) {} - explicit Vec(T value_for_all) + explicit Vec(T value_for_all) { for(unsigned int i=0; i - explicit Vec(const S *source) + explicit Vec(const S *source) { for(unsigned int i=0; i - explicit Vec(const Vec& source) + explicit Vec(const Vec& source) { for(unsigned int i=0; i(T v0, T v1) + Vec(T v0, T v1) { assert(N==2); v[0]=v0; v[1]=v1; } - Vec(T v0, T v1, T v2) + Vec(T v0, T v1, T v2) { assert(N==3); v[0]=v0; v[1]=v1; v[2]=v2; } - Vec(T v0, T v1, T v2, T v3) + Vec(T v0, T v1, T v2, T v3) { assert(N==4); v[0]=v0; v[1]=v1; v[2]=v2; v[3]=v3; } - Vec(T v0, T v1, T v2, T v3, T v4) + Vec(T v0, T v1, T v2, T v3, T v4) { assert(N==5); v[0]=v0; v[1]=v1; v[2]=v2; v[3]=v3; v[4]=v4; } - Vec(T v0, T v1, T v2, T v3, T v4, T v5) + Vec(T v0, T v1, T v2, T v3, T v4, T v5) { assert(N==6); v[0]=v0; v[1]=v1; v[2]=v2; v[3]=v3; v[4]=v4; v[5]=v5; From efae514beae6db121dc12d6a0607794eb71af262 Mon Sep 17 00:00:00 2001 From: Marc Day Date: Mon, 3 Aug 2026 10:53:33 +0200 Subject: [PATCH 2/4] isosurface: fix computeArea reporting 0, and support 2D --- Src/isosurface.cpp | 105 ++++++++++++++++++++++++++++++++------------- 1 file changed, 76 insertions(+), 29 deletions(-) diff --git a/Src/isosurface.cpp b/Src/isosurface.cpp index 8083698f..1fb3bfcc 100644 --- a/Src/isosurface.cpp +++ b/Src/isosurface.cpp @@ -1889,6 +1889,82 @@ main (int argc, const Real uniq_time = end_time_uniq - strt_time_uniq; Print() << "Uniquify time: " << uniq_time << '\n'; + // Measure of the extracted iso-level: area in 3D (triangles), arc length in + // 2D (2-node segments). + // + // Disjoint sections need no special handling in either dimension: the + // measure is additive over elements and independent of how they connect, and + // eltSet is deduplicated so nothing is counted twice. So several separate + // contour lines in 2D, or several surface sheets in 3D, all just add up. + // Only a per-section breakdown would need the connectivity, which would be a + // separate feature. + // + // This must run BEFORE the surface-output block below. That block releases + // eltSet (and, with surface_is_large=1, sortedNodes) to reclaim memory for + // large surfaces, so computing the measure afterwards silently reported 0 + // for any run that also wrote a surface. + bool computeArea = false; + pp.query("computeArea",computeArea); + if (computeArea) { + Real measure = 0; + int nSkipped = 0; + for (std::set::const_iterator it = eltSet.begin(); it != eltSet.end(); ++it) { + const Element& elt = *it; + if (elt.size() != nodesPerElt) { + nSkipped++; + continue; + } + + bool nodeOutOfRange = false; + for (int i = 0; i < nodesPerElt; ++i) { + if (elt[i] >= sortedNodes.size()) { + std::cerr << "Accessing node past end: element node " << i + << " = " << elt[i] << " of " << sortedNodes.size() << std::endl; + nodeOutOfRange = true; + } + } + if (nodeOutOfRange) { + nSkipped++; + continue; + } + +#if AMREX_SPACEDIM==2 + // Segment length. + const Real* p0 = sortedNodes[elt[0]]->m_vec; + const Real* p1 = sortedNodes[elt[1]]->m_vec; + + measure += std::sqrt( (p1[0] - p0[0])*(p1[0] - p0[0]) + + (p1[1] - p0[1])*(p1[1] - p0[1]) ); +#else + // Triangle area, as half the magnitude of the edge cross product. + const Real* p0 = sortedNodes[elt[0]]->m_vec; + const Real* p1 = sortedNodes[elt[1]]->m_vec; + const Real* p2 = sortedNodes[elt[2]]->m_vec; + + measure += 0.5*std::sqrt( + std::pow(( p1[1] - p0[1])*(p2[2]-p0[2]) + -(p1[2] - p0[2])*(p2[1]-p0[1]), 2) + + + std::pow(( p1[2] - p0[2])*(p2[0]-p0[0]) + -(p1[0] - p0[0])*(p2[2]-p0[2]), 2) + + + std::pow(( p1[0] - p0[0])*(p2[1]-p0[1]) + -(p1[1] - p0[1])*(p2[0]-p0[0]), 2) ); +#endif + } + + if (nSkipped > 0) { + std::cerr << "computeArea: skipped " << nSkipped << " of " << eltSet.size() + << " elements (wrong node count or node index out of range)" + << std::endl; + } +#if AMREX_SPACEDIM==2 + Print() << "Total length = " << measure << '\n'; +#else + Print() << "Total area = " << measure << '\n'; +#endif + } + const Real strt_time_sout = ParallelDescriptor::second(); bool writeSurf = true; pp.query("writeSurf",writeSurf); @@ -2233,35 +2309,6 @@ main (int argc, const Real sout_time = end_time_sout - strt_time_sout; std::cout << "Surface output time: " << sout_time << '\n'; - // Compute area of isosurface - bool computeArea = false; - pp.query("computeArea",computeArea); - if (computeArea && (AMREX_SPACEDIM==3)) { - Real Area = 0; - for (std::set::const_iterator it = eltSet.begin(); it != eltSet.end(); ++it) { - const Element& elt = *it; - if (elt.size()==3) { - if (elt[0]>=sortedNodes.size() || elt[1]>=sortedNodes.size() || elt[2]>=sortedNodes.size()) { - std::cerr << "Accessing node past end: " << elt[0] << ", " << elt[1] << ", " << elt[2] << std::endl; - } - - const Real* p0 = sortedNodes[elt[0]]->m_vec; - const Real* p1 = sortedNodes[elt[1]]->m_vec; - const Real* p2 = sortedNodes[elt[2]]->m_vec; - - Area += 0.5*sqrt( - pow(( p1[1] - p0[1])*(p2[2]-p0[2]) - -(p1[2] - p0[2])*(p2[1]-p0[1]), 2) - - + pow(( p1[2] - p0[2])*(p2[0]-p0[0]) - -(p1[0] - p0[0])*(p2[2]-p0[2]), 2) - - + pow(( p1[0] - p0[0])*(p2[1]-p0[1]) - -(p1[1] - p0[1])*(p2[0]-p0[0]), 2) ); - } - } - Print() << "Total area = " << Area << '\n'; - } } // IOProc } amrex::Finalize(); From df510b1acff4593e2e1f3c6105fb7fb203c53f02 Mon Sep 17 00:00:00 2001 From: Marc Day Date: Mon, 3 Aug 2026 10:53:53 +0200 Subject: [PATCH 3/4] binMEF: support 2D contour surfaces --- Src/binMEF.cpp | 302 +++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 268 insertions(+), 34 deletions(-) diff --git a/Src/binMEF.cpp b/Src/binMEF.cpp index 67eb847b..c2a45f9c 100755 --- a/Src/binMEF.cpp +++ b/Src/binMEF.cpp @@ -49,23 +49,44 @@ Real triangleArea(const vector& p0, { // Note: assumes (x,y,z) are first 3 components return 0.5*sqrt( - pow( ( p1[1] - p0[1])*(p2[2]-p0[2]) + pow( ( p1[1] - p0[1])*(p2[2]-p0[2]) -(p1[2] - p0[2])*(p2[1]-p0[1]), 2) - - + pow(( p1[2] - p0[2])*(p2[0]-p0[0]) + + + pow(( p1[2] - p0[2])*(p2[0]-p0[0]) -(p1[0] - p0[0])*(p2[2]-p0[2]), 2) - - + pow(( p1[0] - p0[0])*(p2[1]-p0[1]) + + + pow(( p1[0] - p0[0])*(p2[1]-p0[1]) -(p1[1] - p0[1])*(p2[0]-p0[0]), 2)); } +// Length of a 2-node element. The 2D counterpart of triangleArea: for a contour +// written by isosurface in 2D, elements are segments and the measure that +// partitions among bins is arc length rather than area. +// +// Uses the leading AMREX_SPACEDIM components as the coordinates, which is how +// isosurface lays out a MEF node ("X Y [Z] "). Note triangleArea above +// hard-codes components 0,1,2, so it must not be used on a 2D MEF -- component 2 +// there is the first interpolated field, not z. +static +Real segmentLength(const vector& p0, + const vector& p1) +{ + Real sumSq = 0; + for (int d=0; d& A, vector& Abin, vector& B, vector& Bbin, vector& C, vector& Cbin, int binID) { - // Order big to small. + // Order big to small. if (Bbin[binID] > Abin[binID]) { vector t = A; @@ -89,6 +110,62 @@ void orderNodes(vector& A, vector& Abin, } } +// Order a 2-node element so that A is in the higher bin for coordinate binID. +// The 2D counterpart of orderNodes. +static +void orderNodesSegment(vector& A, vector& Abin, + vector& B, vector& Bbin, int binID) +{ + if (Bbin[binID] > Abin[binID]) + { + vector t = A; + vector tbin = Abin; + A = B; Abin = Bbin; + B = t; Bbin = tbin; + } +} + +// Find the single point D where segment AB crosses the bin boundary below A, and +// interpolate all node states to it. +// +// The 2D counterpart of findFG. A triangle straddling a boundary needs two cut +// points (and yields three sub-triangles); a segment needs exactly one, and +// yields two sub-segments, which is why this is so much simpler than the +// triangle case. +static +void findD(const vector& A, + const vector& B, + vector& D, + const vector& binLO, + Real binMax, + int abin, + int comp) +{ + // Assumes A[comp] > B[comp] (i.e. orderNodesSegment has been applied) and + // that A and B share a component layout. + Real fAB; + if (abin >= binLO.size()) // A above the upper bin bound + { + fAB = (A[comp] - binMax)/(A[comp] - B[comp]); + } + else + { + // A is the higher node, so its bin index cannot be below the range here: + // that would require B's to be lower still, and getBin() floors at -1, so + // both would be -1 and the caller would have taken the "same bin" branch. + AMREX_ALWAYS_ASSERT(abin >= 0); + fAB = (A[comp] - binLO[abin])/(A[comp] - B[comp]); + } + + AMREX_ALWAYS_ASSERT(fAB>=0 && fAB<=1); + + // Interpolate all states to the interface + for (int i=0; i& A, const vector& B, @@ -191,7 +268,7 @@ getBin (const vector& val, std::upper_bound(binLO[j].begin(), binLO[j].end(), val[binComps[j]]); --it; - + retVal[j] = it - binLO[j].begin(); } } @@ -225,7 +302,36 @@ bool satisfyCondition(const vector& A, const vector& B, const vector return false; } +// The 2D counterpart of satisfyCondition: keep a segment only if BOTH its nodes +// satisfy the condition. +static +bool satisfyConditionSegment(const vector& A, const vector& B, + int condComp, Real condVal, int condSgn) +{ + if (condSgn > 0) + { + if (A[condComp] > condVal && B[condComp] > condVal) + return true; + } + else if (condSgn < 0) + { + if (A[condComp] < condVal && B[condComp] < condVal) + return true; + } + else + { + if (A[condComp] == condVal && B[condComp] == condVal) + return true; + } + + return false; +} + static long NmyTriangles = 0; +static long NmySegments = 0; + +// Elements this rank actually binned: triangles in 3D, segments in 2D. +static long NmyElements() { return NmyTriangles + NmySegments; } static void processTriangle(const vector& Ai, const vector& AbinI, @@ -268,7 +374,7 @@ void processTriangle(const vector& Ai, const vector& AbinI, return; } else if ( (AbinI[binID]==BbinI[binID]) && (BbinI[binID]==CbinI[binID]) ) - { + { processTriangle(Ai,AbinI,Bi,BbinI,Ci,CbinI,bins,binLO,binMax,areaEps,binComps, condComp,condVal,condSgn,condApply,binID+1); } @@ -327,7 +433,100 @@ void processTriangle(const vector& Ai, const vector& AbinI, processTriangle(F,Fbin,C,Cbin,G,Gbin,bins,binLO,binMax,areaEps,binComps, condComp,condVal,condSgn,condApply,binID); } - } + } +} + +// Distribute a 2-node element's arc length into bins, clipping it at every bin +// boundary. The 2D counterpart of processTriangle, and it follows the same shape: +// +// * recurse over the binning coordinates via binID; +// * when both nodes share a bin for the current coordinate, move to the next; +// * otherwise cut at the boundary and process each fragment. +// +// The cutting step is where it simplifies. A triangle straddling a boundary has +// two topologies (one node separated, or two) needing two cut points and yielding +// three sub-triangles; a segment has only one topology, one cut point D, and two +// sub-segments AD and DB. So there is no orderNodes-driven case split beyond +// putting A on the high side. +// +// As in the triangle version, the fragment on the far side of the cut is +// re-processed at the SAME binID with its bin index decremented, so that a +// segment spanning many bins is clipped repeatedly until each piece lies in one +// bin. Length is therefore partitioned exactly, not apportioned by an +// extent-overlap approximation. +static +void processSegment(const vector& Ai, const vector& AbinI, + const vector& Bi, const vector& BbinI, + map,Real >& bins, + const vector >& binLO, + const vector& binMax, + Real lengthEps, + const vector& binComps, + int condComp, Real condVal, int condSgn, bool condApply, + int binID=0) +{ + const Real length = segmentLength(Ai,Bi); + + if (length < lengthEps) { + return; + } + + if (binID >= AbinI.size()) + { + bool in_range = true; + for (int i=0; i=binLO[i].size()) + in_range = false; + + if (in_range) + { + NmySegments++; + + if ( !(condApply) || satisfyConditionSegment(Ai,Bi,condComp,condVal,condSgn)) + { + bins[AbinI] += length; + } + else + { + areaOutsideCondition += length; + } + } + return; + } + else if (AbinI[binID]==BbinI[binID]) + { + processSegment(Ai,AbinI,Bi,BbinI,bins,binLO,binMax,lengthEps,binComps, + condComp,condVal,condSgn,condApply,binID+1); + } + else + { + vector Abin = AbinI; + vector Bbin = BbinI; + vector A = Ai; + vector B = Bi; + orderNodesSegment(A,Abin,B,Bbin,binID); + + int nComp = A.size(); + vector D(nComp); + findD(A,B,D,binLO[binID],binMax[binID],Abin[binID],binComps[binID]); + vector Dbin = getBin(D,binComps,binLO,binMax); + for (int i=0; i<=binID; ++i) + { + Dbin[i] = Abin[i]; + } + + // A-side fragment: lies wholly inside A's bin for this coordinate, so + // move on to the next binning coordinate. + processSegment(A,Abin,D,Dbin,bins,binLO,binMax,lengthEps,binComps, + condComp,condVal,condSgn,condApply,binID+1); + + // Far-side fragment: D sits exactly on the boundary, so getBin would put + // it back in A's bin. Decrement it and re-process at this same binID, in + // case D..B crosses further boundaries. + Dbin[binID] = Abin[binID] - 1; + processSegment(D,Dbin,B,Bbin,bins,binLO,binMax,lengthEps,binComps, + condComp,condVal,condSgn,condApply,binID); + } } int @@ -373,6 +572,19 @@ main (int argc, int MYLEN; ifs >> nElts; ifs >> MYLEN; + + // 2 nodes per element is a 2D contour (arc length), 3 is a 3D surface + // (area). isosurface writes nodesPerElt = AMREX_SPACEDIM, so a mismatch means + // the MEF was produced by a build of a different dimensionality -- in which + // case the coordinate columns would also be miscounted, silently, since the + // leading AMREX_SPACEDIM components are taken as (x,y[,z]). + if (MYLEN != 2 && MYLEN != 3) + Abort("binMEF supports 2 nodes per element (2D contours) or 3 (3D surfaces)"); + if (MYLEN != AMREX_SPACEDIM) + Abort("This MEF has " + std::to_string(MYLEN) + " nodes per element but " + "binMEF was built for DIM=" + std::to_string(AMREX_SPACEDIM) + + ". Rebuild with DIM=" + std::to_string(MYLEN) + "."); + if (ParallelDescriptor::IOProcessor()) cerr << "...finished reading data header" << endl; @@ -416,13 +628,13 @@ main (int argc, vector binComps; int nc; - if (nc = pp.countval("binComps")) + if ((nc = pp.countval("binComps"))) { binComps.resize(nc); pp.getarr("binComps",binComps,0,nc); for (int i=0; i=nComp) - Abort("At least one element in binComps out of range"); + Abort("At least one element in binComps out of range"); } else Abort("Need to specify binComps array"); @@ -517,7 +729,7 @@ main (int argc, Real area = 0; map,Real > bins; - + int idx = 0; for (int i=0; i& E = eltVec[i]; const vector& A = nodeVec[E[0]-1]; const vector& B = nodeVec[E[1]-1]; - const vector& C = nodeVec[E[2]-1]; vector Abin = getBin(A,binComps,binLO,binMax); vector Bbin = getBin(B,binComps,binLO,binMax); - vector Cbin = getBin(C,binComps,binLO,binMax); - - area += triangleArea(A,B,C); - processTriangle(A,Abin,B,Bbin,C,Cbin,bins,binLO,binMax,areaEps,binComps, - condComp,condVal,condSgn,condApply); + if (MYLEN == 2) + { + // 2D contour: elements are segments and the measure is arc length. + area += segmentLength(A,B); + + processSegment(A,Abin,B,Bbin,bins,binLO,binMax,areaEps,binComps, + condComp,condVal,condSgn,condApply); + } + else + { + const vector& C = nodeVec[E[2]-1]; + vector Cbin = getBin(C,binComps,binLO,binMax); + + area += triangleArea(A,B,C); + + processTriangle(A,Abin,B,Bbin,C,Cbin,bins,binLO,binMax,areaEps,binComps, + condComp,condVal,condSgn,condApply); + } } } // @@ -543,7 +767,7 @@ main (int argc, // //std::cerr << ParallelDescriptor::MyProc() << ": finished processing triangles." << std::endl; ParallelDescriptor::ReduceRealSum(area, IOProc); - + vector binIdx(bins.size()*nc); vector binDat(bins.size()); int icnt = 0; @@ -559,7 +783,7 @@ main (int argc, // Now get data to IOProc. // std::vector pSizes = ParallelDescriptor::Gather(int(bins.size()),IOProc); - + for (int i=0; i,Real >::const_iterator it=bins.begin(); it!=bins.end(); ++it) binSum += it->second; - + // // Dump binned data. @@ -617,7 +841,7 @@ main (int argc, { box = Box(IntVect::TheZeroVector(),IntVect(AMREX_D_DECL(nBins[0]-1,nBins[1]-1,0))); } - + FArrayBox outFab(box,1); outFab.setVal(0.0); for (std::map,Real >::const_iterator it=bins.begin(); it!=bins.end(); ++it) @@ -662,29 +886,39 @@ main (int argc, cout << it->second << endl; } } - - cerr << "Total area of this surface: " << area << " (sum of bins: " << binSum << ")" << endl; + + // "area" for a 3D surface, "length" for a 2D contour. The surrounding + // format is deliberately unchanged, since it is what downstream scripts + // parse to check that the bins account for the whole surface. +#if AMREX_SPACEDIM==2 + const char* measureLabel = "length"; +#else + const char* measureLabel = "area"; +#endif + cerr << "Total " << measureLabel << " of this surface: " << area + << " (sum of bins: " << binSum << ")" << endl; if (condApply) - cerr << " area outside condition: " << areaOutsideCondition + cerr << " " << measureLabel << " outside condition: " << areaOutsideCondition << " (total: " << areaOutsideCondition + binSum << ")" << endl; } if (ParallelDescriptor::NProcs()>1) { if (ParallelDescriptor::IOProcessor()) std::cerr << "Load balance: " << std::endl; - + for (int i=0; i Date: Tue, 4 Aug 2026 12:55:16 +0200 Subject: [PATCH 4/4] isosurface: check element edge lengths and periodic face images Two cheap checks on the assembled surface, both added because a defect there is silent: it corrupts the reported measure by tens of percent while every element still looks individually plausible. Edge lengths: marching squares/cubes cuts one cell, so no edge can exceed the finest cell diagonal. Checked for the single segment in 2D and all three triangle edges in 3D, since a triangle can be badly shaped while its area looks unremarkable. Reports the count, the worst edge in cell diagonals and its endpoints, and states that the measure is unreliable. Tolerance via edge_length_tol (default 1.5). Periodic face images: after a correct periodic clip, each crossing of a periodic face leaves one node on each face at the same transverse position, so the two face node sets are images of each other. A dropped sliver, a duplicated piece, an inconsistently resolved saddle or a resolution mismatch across the seam all break that. Inactive until the clip exists -- no node lies on a face without it -- and it says so rather than passing silently. Motivated by a periodic 1600^2 case where the crossing between the last and first cell centres was interpolated from raw coordinates rather than the minimum image, placing a node mid-domain: segments of 65.8 mm and 14.1 mm where a cell is 0.05 mm, inflating the arc length by 45%, and 109% in the worst plotfile. It is intermittent across regrids, so a passing case proves nothing -- hence a run-time check rather than a test. Disable with check_surface=0. --- Src/isosurface.cpp | 184 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 184 insertions(+) diff --git a/Src/isosurface.cpp b/Src/isosurface.cpp index 1fb3bfcc..779e35cc 100644 --- a/Src/isosurface.cpp +++ b/Src/isosurface.cpp @@ -2,6 +2,8 @@ #include #include #include +#include +#include #include #include @@ -1889,6 +1891,188 @@ main (int argc, const Real uniq_time = end_time_uniq - strt_time_uniq; Print() << "Uniquify time: " << uniq_time << '\n'; + // --------------------------------------------------------------------- + // Surface sanity checks. + // + // Both are cheap, run on the assembled surface, and exist because a + // defect here is otherwise silent: it corrupts the reported measure by + // tens of percent while every element still looks individually plausible. + // + // Enable/disable with check_surface (default on). + // --------------------------------------------------------------------- + bool check_surface = true; + pp.query("check_surface",check_surface); + + if (check_surface && !sortedNodes.empty()) { + + // Finest cell size, per direction, by the same route used elsewhere in + // this file. + Array dxf; + for (int d=0; d worstA(AMREX_SPACEDIM,0), worstB(AMREX_SPACEDIM,0); + + for (std::set::const_iterator it=eltSet.begin(); it!=eltSet.end(); ++it) { + const Element& elt = *it; + if (elt.size() != nodesPerElt) continue; // reported separately below + bool inRange = true; + for (int i=0; i= long(sortedNodes.size())) inRange = false; + } + if (!inRange) continue; + + // 2D: the one segment. 3D: all three triangle edges. + const int nEdges = (nodesPerElt == 2 ? 1 : nodesPerElt); + for (int e=0; em_vec; + const Real* pb = sortedNodes[elt[(e+1) % nodesPerElt]]->m_vec; + Real len2 = 0; + for (int d=0; d edgeMax) { + nBadEdges++; + if (len > worstEdge) { + worstEdge = len; + for (int d=0; d 0) { + std::cerr << "WARNING: surface check: " << nBadEdges + << " element edge(s) exceed " << edge_length_tol + << " x the finest cell diagonal (" << cellDiag << ").\n" + << " Longest is " << worstEdge << " (" + << worstEdge/cellDiag << " cell diagonals), from ("; + for (int d=0; d> onLo, onHi; + for (long i=0; im_vec; + std::vector transverse; + for (int e=0; e matchTol) { + std::cerr << "WARNING: surface check: direction " << d + << " is periodic but its face node sets are not images of" + << " each other: worst partner separation " + << worstMismatch << " against a tolerance of " + << matchTol << " (" << onLo.size() << " nodes per face)." + << "\n"; + } else { + Print() << "Surface check: direction " << d << " periodic faces" + << " match, " << onLo.size() << " node(s) per face, worst" + << " separation " << worstMismatch << "." << std::endl; + } + } + } + // Measure of the extracted iso-level: area in 3D (triangles), arc length in // 2D (2-node segments). //