54 "Coordinate dimension should be at least 2 for a 2D geometry");
67 int i0 = 1, i1 = 0, direction = 1;
68 for (
size_t i = 0; i <
m_verts.size(); ++i)
83 for (
size_t i = 0; i <
m_verts.size(); ++i)
100 for (
size_t i = 0; i <
m_verts.size(); ++i)
112 for (
size_t i = 0; i <
m_verts.size(); ++i)
114 int i1 = (i + 1) %
m_verts.size();
121 m_verts[i1]->GetCoords(vertex);
146 int i0 = 0, i1 = 1, j1 = 0, j2 = 1;
166 Lcoords[j2] = -h / (0.5 +
sqrt(0.25 - epsilon * h));
167 Lcoords[j1] = -
beta * Lcoords[j2] + tty;
191 NekDouble derx_1, derx_2, dery_1, dery_2, jac;
194 NekDouble init0 = Lcoords[0], init1 = Lcoords[1];
202 m_xmap->PhysDeriv(ptsx, DxD1, DxD2);
203 m_xmap->PhysDeriv(ptsy, DyD1, DyD2);
211 while (cnt++ < MaxIterations)
214 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
215 I[0] =
m_xmap->GetBasis(0)->GetI(eta);
216 I[1] =
m_xmap->GetBasis(1)->GetI(eta + 1);
219 xmap =
m_xmap->PhysEvaluate(I, ptsx);
220 ymap =
m_xmap->PhysEvaluate(I, ptsy);
222 F1 = coords[0] - xmap;
223 F2 = coords[1] - ymap;
225 if (F1 * F1 + F2 * F2 < ScaledTol)
227 resid =
sqrt(F1 * F1 + F2 * F2);
232 derx_1 =
m_xmap->PhysEvaluate(I, DxD1);
233 derx_2 =
m_xmap->PhysEvaluate(I, DxD2);
234 dery_1 =
m_xmap->PhysEvaluate(I, DyD1);
235 dery_2 =
m_xmap->PhysEvaluate(I, DyD2);
237 jac = dery_2 * derx_1 - dery_1 * derx_2;
243 (dery_2 * (coords[0] - xmap) - derx_2 * (coords[1] - ymap)) / jac;
247 (-dery_1 * (coords[0] - xmap) + derx_1 * (coords[1] - ymap)) / jac;
249 if (!(std::isfinite(Lcoords[0]) && std::isfinite(Lcoords[1])))
252 std::ostringstream ss;
253 ss <<
"nan or inf found in NewtonIterationForLocCoord in element "
258 if (fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv)
264 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
267 I[0] =
m_xmap->GetBasis(0)->GetI(eta);
268 I[1] =
m_xmap->GetBasis(1)->GetI(eta + 1);
270 xmap =
m_xmap->PhysEvaluate(I, ptsx);
271 ymap =
m_xmap->PhysEvaluate(I, ptsy);
272 F1 = coords[0] - xmap;
273 F2 = coords[1] - ymap;
274 dist =
sqrt(F1 * F1 + F2 * F2);
281 if (cnt >= MaxIterations)
284 m_xmap->LocCoordToLocCollapsed(Lcoords, collCoords);
287 if ((collCoords[0] >= -1.0 && collCoords[0] <= 1.0) &&
288 (collCoords[1] >= -1.0 && collCoords[1] <= 1.0))
290 std::ostringstream ss;
292 ss <<
"Reached MaxIterations (" << MaxIterations
293 <<
") in Newton iteration ";
294 ss <<
"Init value (" << std::setprecision(4) << init0 <<
","
297 ss <<
"Fin value (" << Lcoords[0] <<
"," << Lcoords[1] <<
","
299 ss <<
"Resid = " << resid
300 <<
" Tolerance = " <<
std::sqrt(ScaledTol);
302 WARNINGL1(cnt < MaxIterations, ss.str());
310 NekDouble dist = std::numeric_limits<double>::max();
331 int npts =
m_xmap->GetTotPoints();
340 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
343 m_xmap->LocCollapsedToLocCoord(eta, Lcoords);
351 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
353 m_xmap->LocCollapsedToLocCoord(eta, xi);
354 int npts =
m_xmap->GetTotPoints();
360 dist =
sqrt(
z *
z + dist * dist);
418 return sqrt((xs[0] - gloCoord[0]) * (xs[0] - gloCoord[0]) +
419 (xs[1] - gloCoord[1]) * (xs[1] - gloCoord[1]) +
420 (xs[2] - gloCoord[2]) * (xs[2] - gloCoord[2]));
428 m_xmap->LocCollapsedToLocCoord(eta, xi);
434 int nq =
m_xmap->GetTotPoints();
442 zderxi1(nq, 0.0), xderxi2(nq, 0.0), yderxi2(nq, 0.0),
443 zderxi2(nq, 0.0), xderxi1xi1(nq, 0.0), yderxi1xi1(nq, 0.0),
444 zderxi1xi1(nq, 0.0), xderxi1xi2(nq, 0.0), yderxi1xi2(nq, 0.0),
445 zderxi1xi2(nq, 0.0), xderxi2xi1(nq, 0.0), yderxi2xi1(nq, 0.0),
446 zderxi2xi1(nq, 0.0), xderxi2xi2(nq, 0.0), yderxi2xi2(nq, 0.0),
450 std::array<NekDouble, 3> xc_derxi, yc_derxi, zc_derxi;
452 m_xmap->PhysDeriv(x, xderxi1, xderxi2);
453 m_xmap->PhysDeriv(y, yderxi1, yderxi2);
454 m_xmap->PhysDeriv(
z, zderxi1, zderxi2);
456 m_xmap->PhysDeriv(xderxi1, xderxi1xi1, xderxi1xi2);
457 m_xmap->PhysDeriv(yderxi1, yderxi1xi1, yderxi1xi2);
458 m_xmap->PhysDeriv(zderxi1, zderxi1xi1, zderxi1xi2);
460 m_xmap->PhysDeriv(yderxi2, yderxi2xi1, yderxi2xi2);
461 m_xmap->PhysDeriv(xderxi2, xderxi2xi1, xderxi2xi2);
462 m_xmap->PhysDeriv(zderxi2, zderxi2xi1, zderxi2xi2);
465 NekDouble fx_prev = std::numeric_limits<NekDouble>::max();
490 NekDouble fx = xdiff * xdiff + ydiff * ydiff + zdiff * zdiff;
492 NekDouble fx_derxi1 = 2.0 * xdiff * xc_derxi[0] +
493 2.0 * ydiff * yc_derxi[0] +
494 2.0 * zdiff * zc_derxi[0];
496 NekDouble fx_derxi2 = 2.0 * xdiff * xc_derxi[1] +
497 2.0 * ydiff * yc_derxi[1] +
498 2.0 * zdiff * zc_derxi[1];
501 2.0 * xdiff * xc_derxi1xi1 + 2.0 * xc_derxi[0] * xc_derxi[0] +
502 2.0 * ydiff * yc_derxi1xi1 + 2.0 * yc_derxi[0] * yc_derxi[0] +
503 2.0 * zdiff * zc_derxi1xi1 + 2.0 * zc_derxi[0] * zc_derxi[0];
506 2.0 * xdiff * xc_derxi1xi2 + 2.0 * xc_derxi[1] * xc_derxi[0] +
507 2.0 * ydiff * yc_derxi1xi2 + 2.0 * yc_derxi[1] * yc_derxi[0] +
508 2.0 * zdiff * zc_derxi1xi2 + 2.0 * zc_derxi[1] * zc_derxi[0];
511 2.0 * xdiff * xc_derxi2xi2 + 2.0 * xc_derxi[1] * xc_derxi[1] +
512 2.0 * ydiff * yc_derxi2xi2 + 2.0 * yc_derxi[1] * yc_derxi[1] +
513 2.0 * zdiff * zc_derxi2xi2 + 2.0 * zc_derxi[1] * zc_derxi[1];
524 1 / (fx_derxi1xi1 * fx_derxi2xi2 - fx_derxi1xi2 * fx_derxi1xi2);
525 hessInv[0][0] = det * fx_derxi2xi2;
526 hessInv[0][1] = det * -fx_derxi1xi2;
527 hessInv[1][0] = det * -fx_derxi1xi2;
528 hessInv[1][1] = det * fx_derxi1xi1;
531 if (
abs(fx - fx_prev) < 1e-12)
546 pk[0] = -(hessInv[0][0] * jac[0] + hessInv[1][0] * jac[1]);
547 pk[1] = -(hessInv[0][1] * jac[0] + hessInv[1][1] * jac[1]);
550 while (gamma > 1e-10)
553 xi_pk[0] = xi[0] + pk[0] * gamma;
554 xi_pk[1] = xi[1] + pk[1] * gamma;
557 m_xmap->LocCoordToLocCollapsed(xi_pk, eta_pk);
560 (-1 - std::numeric_limits<NekDouble>::epsilon()) ||
562 (1 + std::numeric_limits<NekDouble>::epsilon()) ||
564 (-1 - std::numeric_limits<NekDouble>::epsilon()) ||
565 eta_pk[1] > (1 + std::numeric_limits<NekDouble>::epsilon()))
571 std::array<NekDouble, 3> xc_pk_derxi, yc_pk_derxi, zc_pk_derxi;
581 NekDouble fx_pk = xc_pk_diff * xc_pk_diff +
582 yc_pk_diff * yc_pk_diff +
583 zc_pk_diff * zc_pk_diff;
585 NekDouble fx_pk_derxi1 = 2.0 * xc_pk_diff * xc_pk_derxi[0] +
586 2.0 * yc_pk_diff * yc_pk_derxi[0] +
587 2.0 * zc_pk_diff * zc_pk_derxi[0];
589 NekDouble fx_pk_derxi2 = 2.0 * xc_pk_diff * xc_pk_derxi[1] +
590 2.0 * yc_pk_diff * yc_pk_derxi[1] +
591 2.0 * zc_pk_diff * zc_pk_derxi[1];
595 NekDouble tmp = pk[0] * fx_derxi1 + pk[1] * fx_derxi2;
596 NekDouble tmp2 = pk[0] * fx_pk_derxi1 + pk[1] * fx_pk_derxi2;
597 if ((fx_pk - (fx + c1 * gamma * tmp)) <
598 std::numeric_limits<NekDouble>::epsilon() &&
599 (-tmp2 - (-c2 * tmp)) <
600 std::numeric_limits<NekDouble>::epsilon())
614 xi[0] += gamma * pk[0];
615 xi[1] += gamma * pk[1];
619 return sqrt(fx_prev);
623 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.
int GetVid(int i) const
Get the ID of vertex i of this object.
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)