56 "Coordinate dimension should be at least 2 for a 2D geometry");
74 for (
size_t i = 0; i <
m_verts.size(); ++i)
90 for (
size_t i = 0; i <
m_verts.size(); ++i)
112 int i0 = 0, i1 = 1, j1 = 0, j2 = 1;
132 Lcoords[j2] = -h / (0.5 +
sqrt(0.25 - epsilon * h));
133 Lcoords[j1] = -
beta * Lcoords[j2] + tty;
157 NekDouble derx_1, derx_2, dery_1, dery_2, jac;
160 NekDouble init0 = Lcoords[0], init1 = Lcoords[1];
168 m_xmap->PhysDeriv(ptsx, DxD1, DxD2);
169 m_xmap->PhysDeriv(ptsy, DyD1, DyD2);
177 while (cnt++ < MaxIterations)
180 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
181 I[0] =
m_xmap->GetBasis(0)->GetI(eta);
182 I[1] =
m_xmap->GetBasis(1)->GetI(eta + 1);
185 xmap =
m_xmap->PhysEvaluate(I, ptsx);
186 ymap =
m_xmap->PhysEvaluate(I, ptsy);
188 F1 = coords[0] - xmap;
189 F2 = coords[1] - ymap;
191 if (F1 * F1 + F2 * F2 < ScaledTol)
193 resid =
sqrt(F1 * F1 + F2 * F2);
198 derx_1 =
m_xmap->PhysEvaluate(I, DxD1);
199 derx_2 =
m_xmap->PhysEvaluate(I, DxD2);
200 dery_1 =
m_xmap->PhysEvaluate(I, DyD1);
201 dery_2 =
m_xmap->PhysEvaluate(I, DyD2);
203 jac = dery_2 * derx_1 - dery_1 * derx_2;
209 (dery_2 * (coords[0] - xmap) - derx_2 * (coords[1] - ymap)) / jac;
213 (-dery_1 * (coords[0] - xmap) + derx_1 * (coords[1] - ymap)) / jac;
215 if (!(std::isfinite(Lcoords[0]) && std::isfinite(Lcoords[1])))
218 std::ostringstream ss;
219 ss <<
"nan or inf found in NewtonIterationForLocCoord in element "
224 if (fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv)
230 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
233 I[0] =
m_xmap->GetBasis(0)->GetI(eta);
234 I[1] =
m_xmap->GetBasis(1)->GetI(eta + 1);
236 xmap =
m_xmap->PhysEvaluate(I, ptsx);
237 ymap =
m_xmap->PhysEvaluate(I, ptsy);
238 F1 = coords[0] - xmap;
239 F2 = coords[1] - ymap;
240 dist =
sqrt(F1 * F1 + F2 * F2);
247 if (cnt >= MaxIterations)
250 m_xmap->LocCoordToLocCollapsed(Lcoords, collCoords);
253 if ((collCoords[0] >= -1.0 && collCoords[0] <= 1.0) &&
254 (collCoords[1] >= -1.0 && collCoords[1] <= 1.0))
256 std::ostringstream ss;
258 ss <<
"Reached MaxIterations (" << MaxIterations
259 <<
") in Newton iteration ";
260 ss <<
"Init value (" << setprecision(4) << init0 <<
"," << init1
263 ss <<
"Fin value (" << Lcoords[0] <<
"," << Lcoords[1] <<
","
265 ss <<
"Resid = " << resid <<
" Tolerance = " <<
sqrt(ScaledTol);
267 WARNINGL1(cnt < MaxIterations, ss.str());
275 NekDouble dist = std::numeric_limits<double>::max();
296 int npts =
m_xmap->GetTotPoints();
305 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
308 m_xmap->LocCollapsedToLocCoord(eta, Lcoords);
364 return sqrt((xs[0] - gloCoord[0]) * (xs[0] - gloCoord[0]) +
365 (xs[1] - gloCoord[1]) * (xs[1] - gloCoord[1]) +
366 (xs[2] - gloCoord[2]) * (xs[2] - gloCoord[2]));
374 m_xmap->LocCollapsedToLocCoord(eta, xi);
380 int nq =
m_xmap->GetTotPoints();
388 zderxi1(nq, 0.0), xderxi2(nq, 0.0), yderxi2(nq, 0.0),
389 zderxi2(nq, 0.0), xderxi1xi1(nq, 0.0), yderxi1xi1(nq, 0.0),
390 zderxi1xi1(nq, 0.0), xderxi1xi2(nq, 0.0), yderxi1xi2(nq, 0.0),
391 zderxi1xi2(nq, 0.0), xderxi2xi1(nq, 0.0), yderxi2xi1(nq, 0.0),
392 zderxi2xi1(nq, 0.0), xderxi2xi2(nq, 0.0), yderxi2xi2(nq, 0.0),
396 std::array<NekDouble, 3> xc_derxi, yc_derxi, zc_derxi;
398 m_xmap->PhysDeriv(x, xderxi1, xderxi2);
399 m_xmap->PhysDeriv(y, yderxi1, yderxi2);
400 m_xmap->PhysDeriv(
z, zderxi1, zderxi2);
402 m_xmap->PhysDeriv(xderxi1, xderxi1xi1, xderxi1xi2);
403 m_xmap->PhysDeriv(yderxi1, yderxi1xi1, yderxi1xi2);
404 m_xmap->PhysDeriv(zderxi1, zderxi1xi1, zderxi1xi2);
406 m_xmap->PhysDeriv(yderxi2, yderxi2xi1, yderxi2xi2);
407 m_xmap->PhysDeriv(xderxi2, xderxi2xi1, xderxi2xi2);
408 m_xmap->PhysDeriv(zderxi2, zderxi2xi1, zderxi2xi2);
411 NekDouble fx_prev = std::numeric_limits<NekDouble>::max();
436 NekDouble fx = xdiff * xdiff + ydiff * ydiff + zdiff * zdiff;
438 NekDouble fx_derxi1 = 2.0 * xdiff * xc_derxi[0] +
439 2.0 * ydiff * yc_derxi[0] +
440 2.0 * zdiff * zc_derxi[0];
442 NekDouble fx_derxi2 = 2.0 * xdiff * xc_derxi[1] +
443 2.0 * ydiff * yc_derxi[1] +
444 2.0 * zdiff * zc_derxi[1];
447 2.0 * xdiff * xc_derxi1xi1 + 2.0 * xc_derxi[0] * xc_derxi[0] +
448 2.0 * ydiff * yc_derxi1xi1 + 2.0 * yc_derxi[0] * yc_derxi[0] +
449 2.0 * zdiff * zc_derxi1xi1 + 2.0 * zc_derxi[0] * zc_derxi[0];
452 2.0 * xdiff * xc_derxi1xi2 + 2.0 * xc_derxi[1] * xc_derxi[0] +
453 2.0 * ydiff * yc_derxi1xi2 + 2.0 * yc_derxi[1] * yc_derxi[0] +
454 2.0 * zdiff * zc_derxi1xi2 + 2.0 * zc_derxi[1] * zc_derxi[0];
457 2.0 * xdiff * xc_derxi2xi2 + 2.0 * xc_derxi[1] * xc_derxi[1] +
458 2.0 * ydiff * yc_derxi2xi2 + 2.0 * yc_derxi[1] * yc_derxi[1] +
459 2.0 * zdiff * zc_derxi2xi2 + 2.0 * zc_derxi[1] * zc_derxi[1];
470 1 / (fx_derxi1xi1 * fx_derxi2xi2 - fx_derxi1xi2 * fx_derxi1xi2);
471 hessInv[0][0] = det * fx_derxi2xi2;
472 hessInv[0][1] = det * -fx_derxi1xi2;
473 hessInv[1][0] = det * -fx_derxi1xi2;
474 hessInv[1][1] = det * fx_derxi1xi1;
477 if (
abs(fx - fx_prev) < 1e-12)
492 pk[0] = -(hessInv[0][0] * jac[0] + hessInv[1][0] * jac[1]);
493 pk[1] = -(hessInv[0][1] * jac[0] + hessInv[1][1] * jac[1]);
496 while (gamma > 1e-10)
499 xi_pk[0] = xi[0] + pk[0] * gamma;
500 xi_pk[1] = xi[1] + pk[1] * gamma;
503 m_xmap->LocCoordToLocCollapsed(xi_pk, eta_pk);
506 (-1 - std::numeric_limits<NekDouble>::epsilon()) ||
508 (1 + std::numeric_limits<NekDouble>::epsilon()) ||
510 (-1 - std::numeric_limits<NekDouble>::epsilon()) ||
511 eta_pk[1] > (1 + std::numeric_limits<NekDouble>::epsilon()))
517 std::array<NekDouble, 3> xc_pk_derxi, yc_pk_derxi, zc_pk_derxi;
527 NekDouble fx_pk = xc_pk_diff * xc_pk_diff +
528 yc_pk_diff * yc_pk_diff +
529 zc_pk_diff * zc_pk_diff;
531 NekDouble fx_pk_derxi1 = 2.0 * xc_pk_diff * xc_pk_derxi[0] +
532 2.0 * yc_pk_diff * yc_pk_derxi[0] +
533 2.0 * zc_pk_diff * zc_pk_derxi[0];
535 NekDouble fx_pk_derxi2 = 2.0 * xc_pk_diff * xc_pk_derxi[1] +
536 2.0 * yc_pk_diff * yc_pk_derxi[1] +
537 2.0 * zc_pk_diff * zc_pk_derxi[1];
541 NekDouble tmp = pk[0] * fx_derxi1 + pk[1] * fx_derxi2;
542 NekDouble tmp2 = pk[0] * fx_pk_derxi1 + pk[1] * fx_pk_derxi2;
543 if ((fx_pk - (fx + c1 * gamma * tmp)) <
544 std::numeric_limits<NekDouble>::epsilon() &&
545 (-tmp2 - (-c2 * tmp)) <
546 std::numeric_limits<NekDouble>::epsilon())
560 xi[0] += gamma * pk[0];
561 xi[1] += gamma * pk[1];
565 return sqrt(fx_prev);
569 ASSERTL0(
false,
"Geometry type unknown")
#define WARNINGL1(condition, msg)
#define ASSERTL0(condition, msg)
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
int v_GetShapeDim() const override
Get the object's shape dimension.
void NewtonIterationForLocCoord(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &ptsx, const Array< OneD, const NekDouble > &ptsy, Array< OneD, NekDouble > &Lcoords, NekDouble &dist)
Array< OneD, int > m_manifold
Array< OneD, Array< OneD, NekDouble > > m_edgeNormal
void SolveStraightEdgeQuad(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
PointGeomSharedPtr v_GetVertex(int i) const override
Geometry1DSharedPtr v_GetEdge(int i) const override
Returns edge i of this object.
std::vector< StdRegions::Orientation > m_eorient
NekDouble v_GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords) override
Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this ge...
int v_AllLeftCheck(const Array< OneD, const NekDouble > &gloCoord) override
void v_CalculateInverseIsoParam() override
StdRegions::Orientation v_GetEorient(const int i) const override
Returns the orientation of edge i with respect to the ordering of edges in the standard element.
NekDouble v_FindDistance(const Array< OneD, const NekDouble > &xs, Array< OneD, NekDouble > &xi) override
int v_GetNumEdges() const override
Get the number of edges of this object.
int v_GetNumVerts() const override
Get the number of vertices of this object.
Base class for shape geometry information.
NekDouble GetCoord(const int i, const Array< OneD, const NekDouble > &Lcoord)
Given local collapsed coordinate Lcoord, return the value of physical coordinate in direction i.
NekDouble GetLocCoords(const Array< OneD, const NekDouble > &coords, Array< OneD, NekDouble > &Lcoords)
Determine the local collapsed coordinates that correspond to a given Cartesian coordinate for this ge...
virtual void v_FillGeom()
Populate the coordinate mapping Geometry::m_coeffs information from any children geometry elements.
Array< OneD, Array< OneD, NekDouble > > m_isoParameter
Array< OneD, Array< OneD, NekDouble > > m_invIsoParam
bool ClampLocCoords(Array< OneD, NekDouble > &locCoord, NekDouble tol=std::numeric_limits< NekDouble >::epsilon())
Clamp local coords to be within standard regions [-1, 1]^dim.
int GetGlobalID(void) const
Get the ID of this object.
Array< OneD, Array< OneD, NekDouble > > m_coeffs
Array containing expansion coefficients of m_xmap.
GeomFactorsSharedPtr GetMetricInfo()
Get the geometric factors for this object.
StdRegions::StdExpansionSharedPtr m_xmap
mapping containing isoparametric transformation.
StdRegions::StdExpansionSharedPtr GetXmap() const
Return the mapping object Geometry::m_xmap that represents the coordinate transformation from standar...
GeomFactorsSharedPtr m_geomFactors
Geometric factors.
int m_coordim
Coordinate dimension of this geometry object.
@ beta
Gauss Radau pinned at x=-1,.
static const unsigned int kNewtonIterations
std::shared_ptr< Curve > CurveSharedPtr
@ eRegular
Geometry is straight-sided with constant geometric factors.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< PointGeom > PointGeomSharedPtr
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
std::vector< double > z(NPUPPER)
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
scalarT< T > abs(scalarT< T > in)
scalarT< T > sqrt(scalarT< T > in)