diff --git a/docs/html/mnt_regrid_edges_8h.html b/docs/html/mnt_regrid_edges_8h.html index c7aea798..3dd4d5c9 100644 --- a/docs/html/mnt_regrid_edges_8h.html +++ b/docs/html/mnt_regrid_edges_8h.html @@ -223,6 +223,7 @@

Note
the orientation of edges is left to right, bottom to top when using CELL_BY_CELL 3--->—2 | | ^ ^ | | 0--->—1 The line integrals should be compatible with the orientation of the edges.
Returns
error code (0 is OK)
diff --git a/docs/html/mnt_regrid_edges_8h_source.html b/docs/html/mnt_regrid_edges_8h_source.html index 48217b9d..e51f93a4 100644 --- a/docs/html/mnt_regrid_edges_8h_source.html +++ b/docs/html/mnt_regrid_edges_8h_source.html @@ -237,29 +237,29 @@
359LIBRARY_API
360int mnt_regridedges_getNumDstEdges(RegridEdges_t** self, std::size_t* n);
361
-
370LIBRARY_API
-
371int mnt_regridedges_apply(RegridEdges_t** self,
-
372 const double src_data[], double dst_data[], int placement);
-
373
-
384LIBRARY_API
-
385int mnt_regridedges_vectorApply(RegridEdges_t** self,
-
386 const double src_u[], const double src_v[],
-
387 double dst_u[], double dst_v[],
-
388 int fs);
-
389
-
399LIBRARY_API
-
400int mnt_regridedges_loadWeights(RegridEdges_t** self,
-
401 const char* fort_filename, int n);
-
402
-
411LIBRARY_API
-
412int mnt_regridedges_dumpWeights(RegridEdges_t** self,
-
413 const char* fort_filename, int n);
-
414
-
420LIBRARY_API
-
421int mnt_regridedges_print(RegridEdges_t** self);
-
422
-
423
-
424#endif // MNT_REGRID_EDGES
+
377LIBRARY_API
+
378int mnt_regridedges_apply(RegridEdges_t** self,
+
379 const double src_data[], double dst_data[], int placement);
+
380
+
391LIBRARY_API
+
392int mnt_regridedges_vectorApply(RegridEdges_t** self,
+
393 const double src_u[], const double src_v[],
+
394 double dst_u[], double dst_v[],
+
395 int fs);
+
396
+
406LIBRARY_API
+
407int mnt_regridedges_loadWeights(RegridEdges_t** self,
+
408 const char* fort_filename, int n);
+
409
+
418LIBRARY_API
+
419int mnt_regridedges_dumpWeights(RegridEdges_t** self,
+
420 const char* fort_filename, int n);
+
421
+
427LIBRARY_API
+
428int mnt_regridedges_print(RegridEdges_t** self);
+
429
+
430
+
431#endif // MNT_REGRID_EDGES
mnt_regridedges_del
LIBRARY_API int mnt_regridedges_del(RegridEdges_t **self)
Destructor.
mnt_regridedges_setDstGridFlags
LIBRARY_API int mnt_regridedges_setDstGridFlags(RegridEdges_t **self, int fixLonAcrossDateline, int averageLonAtPole)
Set the destination grid flags.
diff --git a/mint/tests/test_vector_interp_cubedsphere.py b/mint/tests/test_vector_interp_cubedsphere.py index 1c560414..cb02413d 100644 --- a/mint/tests/test_vector_interp_cubedsphere.py +++ b/mint/tests/test_vector_interp_cubedsphere.py @@ -74,7 +74,10 @@ import numpy import pytest -from mint import VectorInterp, CELL_BY_CELL_DATA +from mint import Grid, VectorInterp, CELL_BY_CELL_DATA +from tempfile import TemporaryDirectory + +DATA_DIR = Path(__file__).absolute().parent.parent.parent / Path('data') # scripts/generate_cubedsphere_grid.py is not part of the mint package; # reach it via sys.path, as in test_cubedsphere.py. @@ -236,6 +239,33 @@ def test_vector_interp_convergence(): f'({e_prev:.3e} -> {e_next:.3e}): the reconstruction is not converging') +def test_lonlat_vector_interp_fv3_grid(): + # Read in the FV3 cubed-sphere grid. + src_grid = Grid() + src_grid.setFlags(fixLonAcrossDateline=1, averageLonAtPole=1, degrees=True) + filename = str(DATA_DIR / Path('C24_SCRIP_desc.181018.nc')) + src_grid.loadFromFV32DFile(filename) + + # Construct target points that are regularly spaced + # in lon-lat coordinates + lon_reg_A = numpy.linspace(1, 359, 180) + lat_reg_A = numpy.linspace(-89, 89, 90) + + LAT_REG_A, LON_REG_A = numpy.meshgrid(lat_reg_A, lon_reg_A) + ntarget = LON_REG_A.shape[0] * LAT_REG_A.shape[1] + target_points_A = numpy.zeros((ntarget, 3), numpy.float64) + target_points_A[:, 0] = LON_REG_A.flat + target_points_A[:, 1] = LAT_REG_A.flat + + # Check the points are valid + vi = VectorInterp() + vi.setGrid(src_grid) + vi.buildLocator(numCellsPerBucket=128, periodX=360., enableFolding=True) + numbad = vi.findPoints(target_points_A, tol2=1.e-12) + + assert numbad == 0, f'{numbad} bad cells found' + + if __name__ == '__main__': test_exact_gradient_matches_finite_difference() for M in (4, 8, 16, 32): diff --git a/src/MvMat_double.cpp b/src/MvMat_double.cpp index 29d2a583..d0643481 100644 --- a/src/MvMat_double.cpp +++ b/src/MvMat_double.cpp @@ -29,41 +29,37 @@ ColMat eye(std::size_t n){ ***********************************************************************/ -ColMat load(const std::string& cfile) { +ColMat load(const std::string& cfile) { - // copy to temp1 + // strip lines starting with % (comments), in memory -- no temp files, no shelling + // out to cp/sed/rm (which also meant cfile was being pasted, unescaped, into a shell + // command) + ifstream file(cfile.c_str()); + std::ostringstream stripped; + std::string line; + while (std::getline(file, line)) { + if (line.empty() || line[0] != '%') { + stripped << line << '\n'; + } + } - std::string temp1 = "temp1"; - std::string s1 = "cp " + cfile + " " + temp1 + '\n'; - system(s1.c_str()); + std::istringstream in(stripped.str()); - // get rid of lines starting with % (comments) + unsigned int nr, nc; + in >> nr; + in >> nc; - std::string temp2 = "temp2"; - std::string s2 = "sed 's;^%.*$;;g' " + temp1 + " > " + temp2 + '\n'; - system(s2.c_str()); - - const char *c_temp2 = temp2.c_str(); - ifstream file(c_temp2); - unsigned int nr, nc; - file >> nr; - file >> nc; + ColMat a(nr, nc, 0.); + for(unsigned int i = 0; i < nr; ++i){ + for(unsigned int j = 0; j < nc; ++j){ + in >> a(i,j); + } + } - ColMat a(nr, nc, 0.); - - for(unsigned int i = 0; i < nr; ++i){ - for(unsigned int j = 0; j < nc; ++j){ - file >> a(i,j); - } - } - - std::string s3 = "rm " + temp1 + ' ' + temp2 + '\n'; - system(s3.c_str()); - - return a; -} + return a; +} // real matrix dot complex vector Vec_cmplx dot(const Mat &a, const Vec_cmplx &b) diff --git a/src/mntPolylineIntegral.cpp b/src/mntPolylineIntegral.cpp index a0843ae8..c2628d26 100644 --- a/src/mntPolylineIntegral.cpp +++ b/src/mntPolylineIntegral.cpp @@ -61,6 +61,14 @@ int mnt_polylineintegral_buildLocator(PolylineIntegral_t** self, if (enableFolding == 1) { (*self)->loc->enableFolding(); } + // NOTE: deliberately not calling setCubedSphere here. containsPoint's spherical + // treatment (see vmtCellLocator.h) is only consistent with FindCell/VectorInterp, + // whose pcoords/weights come from the very same spherical model. PolylineIntegral's + // line/cell-edge intersections (collectIntersectionPoints) use a separate, still + // flat-(lon,lat)-straight-line algorithm (LineLineIntersector) that setCubedSphere + // would put out of sync with containsPoint near a pole -- fixing that would mean + // reworking the intersection math itself to use great-circle arcs, not just this + // locator, so it's left as flat/unchanged here rather than partially fixed. (*self)->loc->BuildLocator(); return ier; diff --git a/src/mntRegridEdges.cpp b/src/mntRegridEdges.cpp index 01a8c5ce..ef467d5d 100644 --- a/src/mntRegridEdges.cpp +++ b/src/mntRegridEdges.cpp @@ -562,6 +562,14 @@ int mnt_regridedges_buildLocator(RegridEdges_t** self, int numCellsPerBucket, if (enableFolding == 1) { (*self)->srcLoc->enableFolding(); } + // NOTE: deliberately not calling setCubedSphere here. containsPoint's spherical + // treatment (see vmtCellLocator.h) is only consistent with FindCell/VectorInterp, + // whose pcoords/weights come from the very same spherical model. RegridEdges' + // line/cell-edge intersections (collectIntersectionPoints) use a separate, still + // flat-(lon,lat)-straight-line algorithm (LineLineIntersector) that setCubedSphere + // would put out of sync with containsPoint near a pole -- fixing that would mean + // reworking the intersection math itself to use great-circle arcs, not just this + // locator, so it's left as flat/unchanged here rather than partially fixed. (*self)->srcLoc->BuildLocator(); return 0; diff --git a/src/mntVectorInterp.cpp b/src/mntVectorInterp.cpp index d0de357c..c14dce21 100644 --- a/src/mntVectorInterp.cpp +++ b/src/mntVectorInterp.cpp @@ -56,6 +56,11 @@ int mnt_vectorinterp_buildLocator(VectorInterp_t** self, int numCellsPerBucket, if (enableFolding == 1) { (*self)->locator->enableFolding(); } + // fixLonAcrossDateline and averageLonAtPole are only ever both set for a gnomonic + // cubed-sphere grid (see Grid_t/mnt_grid_setFlags) -- a plain (possibly rotated) + // lon-lat grid sets neither, even when it reaches a pole + (*self)->locator->setCubedSphere((*self)->grid->fixLonAcrossDateline && + (*self)->grid->averageLonAtPole); (*self)->locator->BuildLocator(); return 0; diff --git a/src/vmtCellLocator.cpp b/src/vmtCellLocator.cpp index 6ba30eb6..03bec6e5 100644 --- a/src/vmtCellLocator.cpp +++ b/src/vmtCellLocator.cpp @@ -50,6 +50,13 @@ vmtCellLocator::vmtCellLocator() { this->grid = NULL; + // 0 (not periodic) until setPeriodicityLengthX says otherwise -- containsPoint + // uses this to decide whether (lon, lat)-in-degrees semantics apply at all + this->periodX = 0.0; + + // false until setCubedSphere says otherwise + this->isCubedSphere = false; + // will be determined in SetDataSet this->numBucketsX = 1; this->numBucketsY = 1; @@ -115,6 +122,13 @@ vmtCellLocator::BuildLocator() { // assign each face to one or more buckets depending on where the face's nodes fall // WARNING: this could fail if the buckets are much smaller than some cells! vtkIdType numFaces = this->grid->GetNumberOfCells(); + + if (this->periodX > 0 && this->isCubedSphere) { + this->cubedSphereVerts.resize(numFaces*4); + this->cubedSphereCentroid.resize(numFaces); + this->cubedSphereRadius2.resize(numFaces); + } + for (vtkIdType faceId = 0; faceId < numFaces; ++faceId) { std::vector nodes = getFacePoints(faceId); for (const Vec3& p : nodes) { @@ -123,6 +137,24 @@ vmtCellLocator::BuildLocator() { int bucketId = this->getBucketId(&p[0]); this->bucket2Faces[bucketId].insert(faceId); } + + if (this->periodX > 0 && this->isCubedSphere) { + // precompute this face's corners in Cartesian XYZ, plus its centroid and + // radius for containsPointCubedSphere's cheap pre-filter, once here instead + // of on every containment check against this face + Vec3* verts = &this->cubedSphereVerts[faceId*4]; + for (int i = 0; i < 4; ++i) { + verts[i] = this->lonLatDegToXYZ(&nodes[i][0]); + } + Vec3 centroid = (verts[0] + verts[1] + verts[2] + verts[3]) / 4.; + double radius2 = 0.0; + for (int i = 0; i < 4; ++i) { + Vec3 d = verts[i] - centroid; + radius2 = std::max(radius2, dot(d, d)); + } + this->cubedSphereCentroid[faceId] = centroid; + this->cubedSphereRadius2[faceId] = radius2; + } } } @@ -154,16 +186,66 @@ vmtCellLocator::enableFolding() { this->kFolding[1] = 1; } -bool +bool vmtCellLocator::containsPoint(vtkIdType faceId, const double point[3], double tol) const { tol = std::abs(tol); + + if (this->periodX > 0 && this->isCubedSphere) { + // a genuine gnomonic cubed sphere (see setCubedSphere) -- use the exact + // spherical bilinear patch rather than a straight (lon, lat) chord; pcoords and + // weights are not needed here, only the containment decision + double pcoords[3], weights[8]; + return this->containsPointCubedSphere(faceId, point, tol, pcoords, weights); + } + std::vector nodes = this->getFacePoints(faceId); Vec3 targetPoint(point); + return isPointInQuad(targetPoint, nodes, tol); +} - bool res = isPointInQuad(targetPoint, nodes, tol); - return res; +bool +vmtCellLocator::invertSphericalBilinearPatch(const Vec3& target, const Vec3 verts[4], + double& xsi, double& eta) const { + + const int maxIter = 30; + const double newtonTol2 = 1.e-24; // (1.e-12)^2 + const double h = 1.e-6; + + xsi = 0.5; + eta = 0.5; + + for (int iter = 0; iter < maxIter; ++iter) { + + Vec3 res = this->sphericalBilinearMap(xsi, eta, verts) - target; + + // finite-difference Jacobian columns, d(map)/dXsi and d(map)/dEta + Vec3 rXsi = (this->sphericalBilinearMap(xsi + h, eta, verts) - + this->sphericalBilinearMap(xsi - h, eta, verts)) / (2.*h); + Vec3 rEta = (this->sphericalBilinearMap(xsi, eta + h, verts) - + this->sphericalBilinearMap(xsi, eta - h, verts)) / (2.*h); + + // Gauss-Newton step: solve the 2x2 normal equations (J^T J) delta = -J^T res, + // where J = [rXsi, rEta] is the map's (3 x 2) Jacobian + double a11 = dot(rXsi, rXsi), a12 = dot(rXsi, rEta), a22 = dot(rEta, rEta); + double b1 = -dot(rXsi, res), b2 = -dot(rEta, res); + double det = a11*a22 - a12*a12; + if (std::abs(det) < 1.e-30) { + // degenerate Jacobian -- should not happen for a non-degenerate quad + return false; + } + double dXsi = (a22*b1 - a12*b2) / det; + double dEta = (a11*b2 - a12*b1) / det; + + xsi += dXsi; + eta += dEta; + + if (dXsi*dXsi + dEta*dEta < newtonTol2) { + return true; + } + } + return false; } @@ -203,18 +285,55 @@ vmtCellLocator::containsPointMultiValued(vtkIdType faceId, const double point[3] vtkIdType vmtCellLocator::FindCell(const double point[3], double tol, vtkGenericCell *notUsed, double pcoords[3], double *weights) { - int bucketId = this->getBucketId(point); double closestPoint[3]; int subId; double dist2; - const std::set& faces = this->bucket2Faces.find(bucketId)->second; + // Try the point as given, then folded across the pole and/or shifted by + // +-periodX, so a point near a periodic seam or pole singularity is + // looked up in the bucket(s) where its true cell actually lives -- + // mirrors the shift-before-bucket-lookup pattern already used by + // findIntersectionsWithLine for line integrals. + for (const int& kFold : this->kFolding) { + + Vec3 p(point); + + if (kFold == 1) { + if (std::abs(p[1]) <= 90) { + // folding only makes sense for points that already fell + // outside the +-90 deg latitude range + continue; + } + this->foldAtPole(&p[0]); + } + + for (const double& modPx : this->modPeriodX) { + + p[0] += modPx; + + int bucketId = this->getBucketId(&p[0]); + const std::set& faces = this->bucket2Faces.find(bucketId)->second; + + for (const vtkIdType& cId : faces) { + if (this->periodX > 0 && this->isCubedSphere) { + // get pcoords/weights from the same spherical model used for the + // containment decision itself (see containsPointCubedSphere) -- + // computing them together, here, avoids solving twice and, more + // importantly, avoids handing this point to vtkQuad::EvaluatePosition's + // unrelated flat, straight-(lon,lat)-chord model below, which can + // extrapolate (xsi, eta) wildly for a point close to a pole + if (this->containsPointCubedSphere(cId, &p[0], tol, pcoords, weights)) { + return cId; + } + } + else if (this->containsPoint(cId, &p[0], tol)) { + vtkCell* quad = this->grid->GetCell(cId); + quad->EvaluatePosition(&p[0], closestPoint, subId, pcoords, dist2, weights); + return cId; + } + } - for (const vtkIdType& cId : faces) { - if (this->containsPoint(cId, point, tol)) { - vtkCell* quad = this->grid->GetCell(cId); - quad->EvaluatePosition((double*) point, closestPoint, subId, pcoords, dist2, weights); - return cId; + p[0] -= modPx; } } diff --git a/src/vmtCellLocator.h b/src/vmtCellLocator.h index ffb8913e..f8d5794c 100644 --- a/src/vmtCellLocator.h +++ b/src/vmtCellLocator.h @@ -111,6 +111,25 @@ class vmtCellLocator { */ void enableFolding(); + /** + * Declare whether the grid is a gnomonic (central-projection) cubed sphere, e.g. FV3's + * @param isCubedSphere true if the grid is a cubed sphere + * @note only meaningful together with a periodic-longitude-in-degrees grid (periodX + * set > 0); for such a grid, a straight (lon, lat) chord between two cell + * corners is an excellent approximation of that cube-face edge everywhere + * except close to a pole -- where, since gnomonic projection maps straight + * lines on a cube face to great circles, the actual edge is a great circle that + * a straight (lon, lat) chord can badly misrepresent (see + * invertSphericalBilinearPatch). A plain (possibly rotated) lon-lat grid has no + * such cube faces -- its cell edges (lines of constant longitude or latitude) + * are exactly what a straight (lon, lat) chord already represents (a meridian + * is itself a great circle; a circle of latitude is not, but is still not a + * cubed-sphere edge either) -- so leave this false (the default) for those. + */ + void setCubedSphere(bool isCubedSphere) { + this->isCubedSphere = isCubedSphere; + } + /** * Find all intersection points between line and the grid * @param pBeg start point of the line @@ -179,6 +198,22 @@ class vmtCellLocator { // periodicity in x double periodX; + // true if the grid is a gnomonic cubed sphere (see setCubedSphere) + bool isCubedSphere; + + // for a cubed-sphere grid, every face's 4 corners as unit vectors in Cartesian XYZ + // (flat array, faceId*4 + corner), computed once in BuildLocator instead of on every + // containsPointCubedSphere call -- the same handful of candidate faces per bucket get + // tested against many different target points, so this easily pays for itself + std::vector cubedSphereVerts; + + // for a cubed-sphere grid, every face's centroid (mean of its 4 corners' unit + // vectors, not itself a unit vector) and squared radius (largest squared distance + // from that centroid to any of the 4 corners) -- the cheap pre-filter in + // containsPointCubedSphere compares against these instead of recomputing them + std::vector cubedSphereCentroid; + std::vector cubedSphereRadius2; + /** * Adjust the longitude and latitude to account for the folding at the pole * @param point lon, lat in input and transformed lon, lat on output @@ -190,6 +225,137 @@ class vmtCellLocator { point[1] = sgnTheta*180 - point[1]; } + /** + * Convert a (lon, lat) point, in degrees, to a unit vector in 3D Cartesian space + * @param lonLatDeg pointer to lon, lat, in degrees + * @return unit vector on the sphere + */ + inline Vec3 lonLatDegToXYZ(const double lonLatDeg[2]) const { + const double deg2rad = M_PI / 180.0; + double lam = lonLatDeg[0] * deg2rad; + double the = lonLatDeg[1] * deg2rad; + double cosThe = std::cos(the); + Vec3 xyz; + xyz[0] = cosThe * std::cos(lam); + xyz[1] = cosThe * std::sin(lam); + xyz[2] = std::sin(the); + return xyz; + } + + /** + * Spherical linear interpolation between two unit vectors + * @param u start unit vector + * @param w end unit vector + * @param t interpolation parameter, in [0, 1] + * @return unit vector, u at t=0, w at t=1, along the great-circle arc between them + */ + inline Vec3 slerp(const Vec3& u, const Vec3& w, double t) const { + double cosOmega = std::max(-1.0, std::min(1.0, dot(u, w))); + double omega = std::acos(cosOmega); + if (omega < 1.e-9) { + // u and w (nearly) coincide -- avoid the 0/0 below, any t gives the same point + return u; + } + double s = std::sin(omega); + return (std::sin((1. - t)*omega)/s)*u + (std::sin(t*omega)/s)*w; + } + + /** + * Map parametric coordinates to a point on the exact spherical quad spanned by a + * cubed-sphere face's 4 corners (a "spherical bilinear patch": interpolate along the + * 0->1 and 3->2 edges at xsi, then between those two points at eta -- all by great + * circle, via slerp) + * @param xsi xsi parametric coordinate, in [0, 1] inside the face + * @param eta eta parametric coordinate, in [0, 1] inside the face + * @param verts the face's 4 corners, unit vectors in Cartesian XYZ, in the same + * 0->1->2->3 corner order used everywhere else in this class + * @return unit vector on the sphere + */ + inline Vec3 sphericalBilinearMap(double xsi, double eta, const Vec3 verts[4]) const { + Vec3 a = this->slerp(verts[0], verts[1], xsi); + Vec3 b = this->slerp(verts[3], verts[2], xsi); + Vec3 r = this->slerp(a, b, eta); + return r / std::sqrt(dot(r, r)); + } + + /** + * Find the parametric coordinates (xsi, eta) at which the spherical bilinear patch + * spanned by a cubed-sphere face's 4 corners passes through a target point, by + * Gauss-Newton iteration (finite-difference Jacobian) + * @param target target point, unit vector in Cartesian XYZ + * @param verts the face's 4 corners, unit vectors in Cartesian XYZ + * @param xsi xsi parametric coordinate (output), whether or not the iteration converges + * @param eta eta parametric coordinate (output), whether or not the iteration converges + * @return true if the iteration converged + * @note convergence to (xsi, eta) outside [0, 1] is expected and correct for a target + * point outside the face -- the values stay bounded and well behaved (verified + * against faces spanning up to a full pole-adjacent cubed-sphere panel corner), + * they are just not a valid location inside this particular face + */ + bool invertSphericalBilinearPatch(const Vec3& target, const Vec3 verts[4], + double& xsi, double& eta) const; + + /** + * Check if a point is inside a cubed-sphere face, using the exact spherical bilinear + * patch spanned by its 4 corners for both the containment decision and the + * parametric coordinates -- see invertSphericalBilinearPatch + * @param faceId face/cell Id + * @param point point, (lon, lat) in degrees + * @param tol tolerance on xsi, eta + * @param pcoords parametric coordinates (output): pcoords[0] = xsi, pcoords[1] = eta, + * pcoords[2] = 0; valid whether or not the point turns out to be inside + * @param weights bilinear interpolation weights from (xsi, eta) (output), one per + * corner, in the same 0->1->2->3 order + * @return true if inside + * @note deriving both the containment decision and the parametric coordinates from + * the same spherical model, instead of pairing an accurate containment test + * with vtkQuad's flat, straight-(lon,lat)-chord EvaluatePosition, is the point: + * the two can never disagree with each other close to a pole, where they used + * to -- a point the flat model sees as just outside its own corners can get + * assigned wildly extrapolated (xsi, eta), and from there arbitrarily large + * interpolation weights. + */ + inline bool containsPointCubedSphere(vtkIdType faceId, const double point[3], double tol, + double pcoords[3], double weights[8]) const { + // corners, centroid and radius are all precomputed once, in BuildLocator -- see + // cubedSphereVerts + const Vec3* verts = &this->cubedSphereVerts[faceId*4]; + Vec3 target = this->lonLatDegToXYZ(point); + + // cheap pre-filter: the locator's buckets are sized for a handful-of-flops + // containment test, so a bucket can hold many candidate cells for every one + // that actually contains a given point -- reject the (typically large) majority + // that are nowhere close before paying for a multi-iteration Newton solve. + // Compare squared Cartesian distance to this cell's own centroid against its own + // "radius" (the farthest corner from that centroid), inflated by a generous + // safety factor: this can only ever over-accept (falling through to the exact + // solve below), never wrongly reject a point that's actually inside -- every + // point of the spherical bilinear patch stays within the spherical convex hull + // of its 4 corners, which in turn stays within their centroid's own radius, so + // safetyFactor > 1 is already conservative; the extra margin is just insurance. + const double safetyFactor = 2.0; + Vec3 dt = target - this->cubedSphereCentroid[faceId]; + if (dot(dt, dt) > safetyFactor*safetyFactor*this->cubedSphereRadius2[faceId]) { + pcoords[0] = pcoords[1] = pcoords[2] = -1.0; // clearly outside, not solved + weights[0] = weights[1] = weights[2] = weights[3] = 0.0; + return false; + } + + double xsi = 0.5, eta = 0.5; + bool converged = this->invertSphericalBilinearPatch(target, verts, xsi, eta); + + pcoords[0] = xsi; + pcoords[1] = eta; + pcoords[2] = 0.0; + + weights[0] = (1. - xsi)*(1. - eta); + weights[1] = xsi*(1. - eta); + weights[2] = xsi*eta; + weights[3] = (1. - xsi)*eta; + + return converged && xsi >= -tol && xsi <= 1. + tol && eta >= -tol && eta <= 1. + tol; + } + /** * Get the flat array index of a bucket containing a given point diff --git a/tests/testLineLineIntersector.cxx b/tests/testLineLineIntersector.cxx index 5cc577a1..4fdcc74c 100644 --- a/tests/testLineLineIntersector.cxx +++ b/tests/testLineLineIntersector.cxx @@ -16,8 +16,8 @@ void test1() { lli.setPoints(2, p0, p1, q0, q1); Vec2 xi = lli.getSolution(); std::cout << "test1: xi = " << xi << '\n'; - assert(abs(xi[0] - 1./2.) < tol); - assert(abs(xi[1] - 1./3.) < tol); + assert(std::abs(xi[0] - 1./2.) < tol); + assert(std::abs(xi[1] - 1./3.) < tol); } void test1_3d() { @@ -30,8 +30,8 @@ void test1_3d() { LineLineIntersector lli; lli.setPoints(3, p0, p1, q0, q1); Vec2 xi = lli.getSolution(); - assert(abs(xi[0] - 1./2.) < tol); - assert(abs(xi[1] - 1./3.) < tol); + assert(std::abs(xi[0] - 1./2.) < tol); + assert(std::abs(xi[1] - 1./3.) < tol); } void test1_3DOffset() { @@ -44,8 +44,8 @@ void test1_3DOffset() { LineLineIntersector lli; lli.setPoints(3, p0, p1, q0, q1); Vec2 xi = lli.getSolution(); - assert(abs(xi[0] - 1./2.) < tol); - assert(abs(xi[1] - 1./3.) < tol); + assert(std::abs(xi[0] - 1./2.) < tol); + assert(std::abs(xi[1] - 1./3.) < tol); } @@ -59,7 +59,7 @@ void test2() { LineLineIntersector lli; lli.setPoints(2, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); assert(lli.hasSolution(tol)); } @@ -73,7 +73,7 @@ void test2_3DOffset() { LineLineIntersector lli; lli.setPoints(3, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); assert(! lli.hasSolution(tol)); } @@ -87,7 +87,7 @@ void test3() { LineLineIntersector lli; lli.setPoints(2, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); std::cout << "test3: det = " << det << '\n'; assert(! lli.hasSolution(tol)); } @@ -102,7 +102,7 @@ void test3_3DOffset() { LineLineIntersector lli; lli.setPoints(3, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); // no solution as the two lines run parallel without ever touch each other assert(! lli.hasSolution(tol)); } @@ -118,7 +118,7 @@ void test3_crossAt0() { lli.setPoints(3, p0, p1, q0, q1); double det = lli.getDet(); std::cerr << "test3_crossAt0: det = " << det << '\n'; - assert(abs(det) > tol); + assert(std::abs(det) > tol); assert(lli.hasSolution(tol)); } @@ -133,7 +133,7 @@ void test3_touchMiddle() { lli.setPoints(3, p0, p1, q0, q1); double det = lli.getDet(); std::cerr << "test3_touchMiddle: det = " << det << '\n'; - assert(abs(det) > tol); + assert(std::abs(det) > tol); assert(lli.hasSolution(tol)); } @@ -148,7 +148,7 @@ void test3_noTouchMiddle() { lli.setPoints(3, p0, p1, q0, q1); double det = lli.getDet(); std::cerr << "test3_noTouchMiddle: det = " << det << '\n'; - assert(abs(det) > tol); + assert(std::abs(det) > tol); assert(lli.hasSolution(tol)); } @@ -162,7 +162,7 @@ void testNoOverlap() { LineLineIntersector lli; lli.setPoints(2, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); assert(! lli.hasSolution(tol)); } @@ -176,7 +176,7 @@ void testNoOverlap_3DOffset() { LineLineIntersector lli; lli.setPoints(3, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); assert(! lli.hasSolution(tol)); } @@ -190,7 +190,7 @@ void testNoOverlap2() { LineLineIntersector lli; lli.setPoints(2, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); assert(! lli.hasSolution(tol)); } @@ -204,7 +204,7 @@ void testNoOverlap2_3DOffset() { LineLineIntersector lli; lli.setPoints(3, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); assert(! lli.hasSolution(tol)); } @@ -217,7 +217,7 @@ void testPartialOverlap() { LineLineIntersector lli; lli.setPoints(2, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); assert(lli.hasSolution(tol)); std::pair< double, double > p = lli.getBegEndParamCoords(); double lamA = p.first; @@ -250,7 +250,7 @@ void testPartialOverlap_3DOffset() { LineLineIntersector lli; lli.setPoints(3, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); // lines don't touch assert(! lli.hasSolution(tol)); std::pair< double, double > p = lli.getBegEndParamCoords(); @@ -305,8 +305,8 @@ void testPartialOverlap2() { } u /= sqrt(dot(dp10, dp10)); - assert(abs(dot(dpaq0, u)) < tol); - assert(abs(dot(dpbp1, u)) < tol); + assert(std::abs(dot(dpaq0, u)) < tol); + assert(std::abs(dot(dpbp1, u)) < tol); } void testPartialOverlap2_3DOffset() { @@ -339,8 +339,8 @@ void testPartialOverlap2_3DOffset() { } u /= sqrt(dot(dp10, dp10)); - assert(abs(dot(dpaq0, u)) < tol); - assert(abs(dot(dpbp1, u)) < tol); + assert(std::abs(dot(dpaq0, u)) < tol); + assert(std::abs(dot(dpbp1, u)) < tol); } void testPartialOverlap3() { @@ -352,7 +352,7 @@ void testPartialOverlap3() { LineLineIntersector lli; lli.setPoints(2, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); std::cerr << "testPartialOverlap3: det = " << det << '\n'; assert(lli.hasSolution(tol)); std::pair< double, double > p = lli.getBegEndParamCoords(); @@ -372,8 +372,8 @@ void testPartialOverlap3() { dp10[i] = p1[i] - p0[i]; } u /= sqrt(dot(dp10, dp10)); - assert(abs(dot(dpap0, u)) < tol); - assert(abs(dot(dpbq0, u)) < tol); + assert(std::abs(dot(dpap0, u)) < tol); + assert(std::abs(dot(dpbq0, u)) < tol); } void testPartialOverlap3_3DOffset() { @@ -386,7 +386,7 @@ void testPartialOverlap3_3DOffset() { LineLineIntersector lli; lli.setPoints(3, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); // lines don't touch assert(! lli.hasSolution(tol)); std::pair< double, double > p = lli.getBegEndParamCoords(); @@ -406,8 +406,8 @@ void testPartialOverlap3_3DOffset() { dp10[i] = p1[i] - p0[i]; } u /= sqrt(dot(dp10, dp10)); - assert(abs(dot(dpap0, u)) < tol); - assert(abs(dot(dpbq0, u)) < tol); + assert(std::abs(dot(dpap0, u)) < tol); + assert(std::abs(dot(dpbq0, u)) < tol); } @@ -421,7 +421,7 @@ void testQInsideP() { LineLineIntersector lli; lli.setPoints(2, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); assert(lli.hasSolution(tol)); std::pair< double, double > p = lli.getBegEndParamCoords(); double lamA = p.first; @@ -440,8 +440,8 @@ void testQInsideP() { dpbq1[i] = pb[i] - q1[i]; } u /= sqrt(dot(dp10, dp10)); - assert(abs(dot(dpaq0, u)) < tol); - assert(abs(dot(dpbq1, u)) < tol); + assert(std::abs(dot(dpaq0, u)) < tol); + assert(std::abs(dot(dpbq1, u)) < tol); } void testQInsideP_3DOffset() { @@ -454,7 +454,7 @@ void testQInsideP_3DOffset() { LineLineIntersector lli; lli.setPoints(3, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); // lines don't touch assert(! lli.hasSolution(tol)); std::pair< double, double > p = lli.getBegEndParamCoords(); @@ -474,8 +474,8 @@ void testQInsideP_3DOffset() { dpbq1[i] = pb[i] - q1[i]; } u /= sqrt(dot(dp10, dp10)); - assert(abs(dot(dpaq0, u)) < tol); - assert(abs(dot(dpbq1, u)) < tol); + assert(std::abs(dot(dpaq0, u)) < tol); + assert(std::abs(dot(dpbq1, u)) < tol); } void testPInsideQ() { @@ -488,7 +488,7 @@ void testPInsideQ() { LineLineIntersector lli; lli.setPoints(2, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); assert(lli.hasSolution(tol)); std::pair< double, double > p = lli.getBegEndParamCoords(); double lamA = p.first; @@ -507,8 +507,8 @@ void testPInsideQ() { dpbp1[i] = pb[i] - p1[i]; } u /= sqrt(dot(dp10, dp10)); - assert(abs(dot(dpap0, u)) < tol); - assert(abs(dot(dpbp1, u)) < tol); + assert(std::abs(dot(dpap0, u)) < tol); + assert(std::abs(dot(dpbp1, u)) < tol); } void testPInsideQ_3DOffset() { @@ -521,7 +521,7 @@ void testPInsideQ_3DOffset() { LineLineIntersector lli; lli.setPoints(3, p0, p1, q0, q1); double det = lli.getDet(); - assert(abs(det) < tol); + assert(std::abs(det) < tol); // lines don't touch assert(! lli.hasSolution(tol)); std::pair< double, double > p = lli.getBegEndParamCoords(); @@ -541,8 +541,8 @@ void testPInsideQ_3DOffset() { dpbp1[i] = pb[i] - p1[i]; } u /= sqrt(dot(dp10, dp10)); - assert(abs(dot(dpap0, u)) < tol); - assert(abs(dot(dpbp1, u)) < tol); + assert(std::abs(dot(dpap0, u)) < tol); + assert(std::abs(dot(dpbp1, u)) < tol); } diff --git a/version.txt b/version.txt index 7aa332e4..2b17ffd5 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -1.33.0 +1.34.0