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 #include #include +#include +#include #include #include @@ -1889,6 +1891,264 @@ 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). + // + // 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 +2493,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(); 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;