54 "Coordinate dimension should be at least 3 for a 3D geometry.");
68 const int MaxIterations = 51;
79 ScaledTol *= Tol * Tol;
86 while (cnt++ < MaxIterations)
91 var[4] = Lcoords[0] * Lcoords[1];
92 var[5] = Lcoords[1] * Lcoords[2];
93 var[6] = Lcoords[0] * Lcoords[2];
94 var[7] = var[4] * Lcoords[2];
103 tmp[0] = coords[0] - xmap;
104 tmp[1] = coords[1] - ymap;
105 tmp[2] = coords[2] - zmap;
107 res = tmp[0] * tmp[0] + tmp[1] * tmp[1] + tmp[2] * tmp[2];
154 Lcoords[0] +=
Vmath::Dot(3, mat->GetPtr(), 1, tmp, 1);
155 Lcoords[1] +=
Vmath::Dot(3, mat->GetPtr() + 3, 1, tmp, 1);
156 Lcoords[2] +=
Vmath::Dot(3, mat->GetPtr() + 6, 1, tmp, 1);
158 if (!(std::isfinite(Lcoords[0]) && std::isfinite(Lcoords[1]) &&
159 std::isfinite(Lcoords[2])) ||
160 fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv ||
161 fabs(Lcoords[2]) > LcoordDiv)
163 std::ostringstream ss;
164 ss <<
"Iteration has diverged in NewtonIterationForLocCoord in "
172 if (cnt >= MaxIterations)
174 std::ostringstream ss;
176 ss <<
"Reached MaxIterations (" << MaxIterations
177 <<
") in Newton iteration ";
179 WARNINGL1(cnt < MaxIterations, ss.str());
213 NekDouble derx_1, derx_2, derx_3, dery_1, dery_2, dery_3, derz_1, derz_2,
217 NekDouble init0 = Lcoords[0], init1 = Lcoords[1], init2 = Lcoords[2];
230 m_xmap->PhysDeriv(ptsx, DxD1, DxD2, DxD3);
231 m_xmap->PhysDeriv(ptsy, DyD1, DyD2, DyD3);
232 m_xmap->PhysDeriv(ptsz, DzD1, DzD2, DzD3);
240 while (cnt++ < MaxIterations)
243 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
244 I[0] =
m_xmap->GetBasis(0)->GetI(eta);
245 I[1] =
m_xmap->GetBasis(1)->GetI(eta + 1);
246 I[2] =
m_xmap->GetBasis(2)->GetI(eta + 2);
249 xmap =
m_xmap->PhysEvaluate(I, ptsx);
250 ymap =
m_xmap->PhysEvaluate(I, ptsy);
251 zmap =
m_xmap->PhysEvaluate(I, ptsz);
253 F1 = coords[0] - xmap;
254 F2 = coords[1] - ymap;
255 F3 = coords[2] - zmap;
257 if (F1 * F1 + F2 * F2 + F3 * F3 < ScaledTol)
259 resid =
sqrt(F1 * F1 + F2 * F2 + F3 * F3);
264 derx_1 =
m_xmap->PhysEvaluate(I, DxD1);
265 derx_2 =
m_xmap->PhysEvaluate(I, DxD2);
266 derx_3 =
m_xmap->PhysEvaluate(I, DxD3);
267 dery_1 =
m_xmap->PhysEvaluate(I, DyD1);
268 dery_2 =
m_xmap->PhysEvaluate(I, DyD2);
269 dery_3 =
m_xmap->PhysEvaluate(I, DyD3);
270 derz_1 =
m_xmap->PhysEvaluate(I, DzD1);
271 derz_2 =
m_xmap->PhysEvaluate(I, DzD2);
272 derz_3 =
m_xmap->PhysEvaluate(I, DzD3);
274 jac = derx_1 * (dery_2 * derz_3 - dery_3 * derz_2) -
275 derx_2 * (dery_1 * derz_3 - dery_3 * derz_1) +
276 derx_3 * (dery_1 * derz_2 - dery_2 * derz_1);
282 ((dery_2 * derz_3 - dery_3 * derz_2) * (coords[0] - xmap) -
283 (derx_2 * derz_3 - derx_3 * derz_2) * (coords[1] - ymap) +
284 (derx_2 * dery_3 - derx_3 * dery_2) * (coords[2] - zmap)) /
289 ((dery_1 * derz_3 - dery_3 * derz_1) * (coords[0] - xmap) -
290 (derx_1 * derz_3 - derx_3 * derz_1) * (coords[1] - ymap) +
291 (derx_1 * dery_3 - derx_3 * dery_1) * (coords[2] - zmap)) /
296 ((dery_1 * derz_2 - dery_2 * derz_1) * (coords[0] - xmap) -
297 (derx_1 * derz_2 - derx_2 * derz_1) * (coords[1] - ymap) +
298 (derx_1 * dery_2 - derx_2 * dery_1) * (coords[2] - zmap)) /
301 if (!(std::isfinite(Lcoords[0]) && std::isfinite(Lcoords[1]) &&
302 std::isfinite(Lcoords[2])))
305 std::ostringstream ss;
306 ss <<
"nan or inf found in NewtonIterationForLocCoord in element "
311 if (fabs(Lcoords[0]) > LcoordDiv || fabs(Lcoords[1]) > LcoordDiv ||
312 fabs(Lcoords[2]) > LcoordDiv)
318 m_xmap->LocCoordToLocCollapsed(Lcoords, eta);
321 I[0] =
m_xmap->GetBasis(0)->GetI(eta);
322 I[1] =
m_xmap->GetBasis(1)->GetI(eta + 1);
323 I[2] =
m_xmap->GetBasis(2)->GetI(eta + 2);
325 xmap =
m_xmap->PhysEvaluate(I, ptsx);
326 ymap =
m_xmap->PhysEvaluate(I, ptsy);
327 zmap =
m_xmap->PhysEvaluate(I, ptsz);
328 F1 = coords[0] - xmap;
329 F2 = coords[1] - ymap;
330 F3 = coords[2] - zmap;
331 dist =
sqrt(F1 * F1 + F2 * F2 + F3 * F3);
338 if (cnt >= MaxIterations)
341 m_xmap->LocCoordToLocCollapsed(Lcoords, collCoords);
344 if ((collCoords[0] >= -1.0 && collCoords[0] <= 1.0) &&
345 (collCoords[1] >= -1.0 && collCoords[1] <= 1.0) &&
346 (collCoords[2] >= -1.0 && collCoords[2] <= 1.0))
348 std::ostringstream ss;
350 ss <<
"Reached MaxIterations (" << MaxIterations
351 <<
") in Newton iteration ";
352 ss <<
"Init value (" << std::setprecision(4) << init0 <<
","
353 << init1 <<
"," << init2 <<
") ";
354 ss <<
"Fin value (" << Lcoords[0] <<
"," << Lcoords[1] <<
","
355 << Lcoords[2] <<
") ";
356 ss <<
"Resid = " << resid
357 <<
" Tolerance = " <<
std::sqrt(ScaledTol);
359 WARNINGL1(cnt < MaxIterations, ss.str());
367 NekDouble dist = std::numeric_limits<double>::max();
369 tmpcoords[0] = coords[0];
370 tmpcoords[1] = coords[1];
371 tmpcoords[2] = coords[2];
396 int npts =
m_xmap->GetTotPoints();
419 int qa = za.size(), qb = zb.size();
422 eta[2] = zc[min_i / (qa * qb)];
423 min_i = min_i % (qa * qb);
424 eta[1] = zb[min_i / qa];
425 eta[0] = za[min_i % qa];
426 m_xmap->LocCollapsedToLocCoord(eta, Lcoords);
455 int nFaceCoeffs =
m_faces[i]->GetXmap()->GetNcoeffs();
462 m_xmap->GetTraceToElementMap(
469 m_xmap->GetTraceToElementMap(
480 for (k = 0; k < nFaceCoeffs; k++)
503 return m_xmap->PhysEvaluate(Lcoord, tmp);
538 "Edge ID must be between 0 and " +
539 std::to_string(
m_edges.size() - 1));
545 ASSERTL2((i >= 0) && (i <= 5),
"Edge id must be between 0 and 4");
552 "Edge ID must be between 0 and " +
553 std::to_string(
m_edges.size() - 1));
560 "Face ID must be between 0 and " +
561 std::to_string(
m_faces.size() - 1));
569 for (
int i = 0; i < 3; ++i)
571 for (
int j = 0; j < 3; ++j)
579 for (
int i = 0; i < 3; ++i)
582 for (
int j = 0; j < 3; ++j)
#define WARNINGL1(condition, msg)
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Geometry2DSharedPtr v_GetFace(int i) const override
Returns face i of this object.
int v_GetNumFaces() const override
Get the number of faces of this object.
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
int v_GetNumEdges() const override
Get the number of edges of this object.
std::vector< StdRegions::Orientation > m_forient
void v_CalculateInverseIsoParam() override
void v_FillGeom() override
Put all quadrature information into face/edge structure and backward transform.
NekDouble v_GetCoord(const int i, const Array< OneD, const NekDouble > &Lcoord) override
Given local collapsed coordinate Lcoord return the value of physical coordinate in direction i.
PointGeomSharedPtr v_GetVertex(int i) const override
void NewtonIterationForLocCoord(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &ptsx, const Array< OneD, const NekDouble > &ptsy, const Array< OneD, const NekDouble > &ptsz, Array< OneD, NekDouble > &Lcoords, NekDouble &dist)
Geometry1DSharedPtr v_GetEdge(int i) const override
Returns edge i of this object.
StdRegions::Orientation v_GetForient(const int i) const override
Returns the orientation of face i with respect to the ordering of faces in the standard element.
int v_GetShapeDim() const override
Get the object's shape dimension.
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.
std::vector< StdRegions::Orientation > m_eorient
int v_GetNumVerts() const override
Get the number of vertices of this object.
Base class for shape geometry information.
GeomState m_state
Enumeration to dictate whether coefficients are filled.
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.
static const unsigned int kNewtonIterations
@ eRegular
Geometry is straight-sided with constant geometric factors.
std::shared_ptr< PointGeom > PointGeomSharedPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
@ ePtsFilled
Geometric information has been generated.
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
T Dot(int n, const T *w, const T *x)
dot product
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > sqrt(scalarT< T > in)