52 ModuleKey ProcessWallNormalData::className =
55 ProcessWallNormalData::create,
56 "Export data in wall-normal direction from a single point on the "
62 f->m_writeBndFld =
false;
66 "The point to be projected onto the wall to get the \
67 sampling origin. default=[0.5,0,0]");
70 "The direction to project the point onto the wall to \
71 get the sampling origin. default=[0,1,0]");
73 false,
"1.0",
"Distance to limit projection distance to find the \
74 desired sampling origin. defalut=1.0");
76 false,
"0.01",
"Sampling distance along the wall normal at the \
77 sampling origin. default=0.1");
79 false,
"5",
"Number of sampling points along the wall normal. \
82 false,
"0.1",
"The parameter that controls the sampling points' \
83 distribution. d should be in the range (0,inf). d \
84 in (0,0.95] gives controled points; d in (0.95,inf) \
85 gives evenly spaced points");
118 const int nfields =
m_f->m_variables.size();
119 const int nCoordDim =
m_f->m_exp[0]->GetCoordim(0);
121 const int nBndLcoordDim = nCoordDim - 1;
125 std::vector<NekDouble> xorig, searchDir;
127 "Failed to interpret origin coordinates string");
130 "Failed to interpret search direction string");
133 const int nptsH =
m_config[
"nptsH"].as<
int>();
137 for (
int i = 0; i < 3; ++i)
139 orig[i] = (xorig.size() >= (i + 1)) ? xorig[i] : 0.0;
140 projDir[i] = (searchDir.size() >= (i + 1)) ? searchDir[i] : 0.0;
150 if (
m_f->m_numHomogeneousDir == 1)
152 int nPlanes =
m_f->m_exp[0]->GetHomogeneousBasis()->GetZ().size();
154 if (orig[2] < 0.0 || orig[2] > lHom)
161 NekDouble zTmp, zCur = 0.0, distTmp, distCur = 999.0;
162 for (
int i = 0; i <= nPlanes; ++i)
165 distTmp = fabs(orig[2] - zTmp);
166 if (distTmp < distCur)
182 m_f->m_exp[0]->GetGraph());
185 map<int, int> BndRegionMap;
187 for (
auto &breg_it : bregions)
189 BndRegionMap[breg_it.first] = cnt++;
191 int bnd = BndRegionMap[
m_f->m_bndRegionsToWrite[0]];
195 for (
int i = 0; i < nfields; ++i)
197 BndExp[i] =
m_f->m_exp[i]->UpdateBndCondExpansion(bnd);
206 bool isInside =
false;
209 const int nElmts = BndExp[0]->GetNumElmts();
210 for (elmtid = 0; elmtid < nElmts; ++elmtid)
212 bndGeom = BndExp[0]->GetExp(elmtid)->GetGeom();
220 ASSERTL0(isInside,
"Failed to find the sampling origin on the boundary.");
224 const int npts = bndXmap->GetTotPoints();
228 for (
int i = 0; i < nCoordDim; ++i)
231 bndCoeffs[i] = bndGeom->GetCoeffs(i);
232 bndXmap->BwdTrans(bndCoeffs[i], pts[i]);
233 orig[i] = bndXmap->PhysEvaluate(locCoord, pts[i]);
239 m_f->m_exp[0]->GetBoundaryNormals(bnd, normalsQ);
244 int from_nPtsPerElmt = 1;
245 for (
int i = 0; i < nBndLcoordDim; ++i)
247 from_key[i] = BndExp[0]->GetExp(elmtid)->GetBasis(i)->GetPointsKey();
248 from_nPtsPerElmt *= from_key[i].GetNumPoints();
253 int refId = elmtid * from_nPtsPerElmt + 0;
256 normalsRef[i] = normalsQ[i][refId];
271 cout <<
"------ wallNormalData module ------\n";
272 cout <<
"Input point:\n";
273 cout <<
" - [Px,Py,Pz] = [" << xorig[0] <<
", " << xorig[1];
274 if (xorig.size() >= 3)
276 cout <<
", " << xorig[2] <<
"]\n";
280 cout <<
", " << 0.0 <<
"]\n";
282 cout <<
"Projection direction:\n";
283 cout <<
" - [vx,vy,vz] = [" << projDir[0] <<
", " << projDir[1] <<
", "
284 << projDir[2] <<
"]\n";
285 cout <<
"Sampling origin on the wall:\n";
286 cout <<
" - [Ox,Oy,Oz] = [" << orig[0] <<
", " << orig[1] <<
", "
288 cout <<
"Normals at the origin:\n";
289 cout <<
" - [nx,ny,nz] = [" << normals[0] <<
", " << normals[1] <<
", "
290 << normals[2] <<
"]\n";
292 <<
"Ref normals (at quadrature points in the projected element):\n";
293 for (
int i = 0; i < from_nPtsPerElmt; ++i)
295 cout <<
" - " << i <<
": [nx,ny,nz] = ["
296 << -normalsQ[0][elmtid * from_nPtsPerElmt + i] <<
", "
297 << -normalsQ[1][elmtid * from_nPtsPerElmt + i];
300 cout <<
", " << -normalsQ[2][elmtid * from_nPtsPerElmt + i]
314 if (delta > 0 && delta <= 0.95)
323 for (
int i = 1; i < nptsH; ++i)
325 tmp1 = 1.0 - i * tmp2;
326 h[i] = 1 - tanh(tmp1 * tmp4) * tmp5;
330 else if (delta > 0.95)
334 for (
int i = 1; i < nptsH; ++i)
342 ASSERTL0(
false,
"Input error. Delta needs to be greater than 0.0.");
347 for (
int i = 0; i < totVars; ++i)
354 for (
int j = 0; j < nptsH; ++j)
356 ptsH[i][j] = orig[i] + h[j] * normals[i];
384 tmp1 =
Vmath::Dot(3, gloCoord, 1, projDir, 1);
385 Vmath::Smul(3, distToOrig - tmp1, projDir, 1, tmp2,
387 Vmath::Vadd(3, gloCoord, 1, tmp2, 1, projGloCoord, 1);
403 const int nCoordDim = pts.size();
404 const int npts = pts[0].size();
409 for (
int i = 0; i < npts; ++i)
412 for (
int j = 0; j < nCoordDim; ++j)
414 singlePnt[j] = pts[j][i];
418 tmp1 =
Vmath::Dot(nCoordDim, singlePnt, 1, projDir, 1);
419 Vmath::Smul(nCoordDim, distToOrig - tmp1, projDir, 1, tmp2, 1);
420 Vmath::Vadd(nCoordDim, singlePnt, 1, tmp2, 1, singlePnt, 1);
423 for (
int j = 0; j < nCoordDim; ++j)
425 projPts[j][i] = singlePnt[j];
445 const NekDouble npts = projPts[0].size();
449 vec1[0] = projPts[0][0] - projGloCoord[0];
450 vec1[1] = projPts[1][0] - projGloCoord[1];
451 vec2[0] = projPts[0][npts - 1] - projGloCoord[0];
452 vec2[1] = projPts[1][npts - 1] - projGloCoord[1];
459 if ((fabs(vec1[0] + vec2[0]) + fabs(vec1[1] + vec2[1])) < paralTol)
478 const NekDouble npts = projPts[0].size();
481 const int nptsPolygon = 4 * nptsEdge - 4;
483 for (
int i = 0; i < 3; ++i)
487 for (
int j = 0; j < nptsEdge; ++j)
489 ptsPolygon[i][j] = projPts[i][j];
491 for (
int j = 0; j < nptsEdge - 2; ++j)
493 ptsPolygon[i][nptsEdge + j] = projPts[i][(j + 2) * nptsEdge - 1];
495 for (
int j = 0; j < nptsEdge; ++j)
497 ptsPolygon[i][2 * nptsEdge - 2 + j] = projPts[i][npts - 1 - j];
499 for (
int j = 0; j < nptsEdge - 2; ++j)
501 ptsPolygon[i][3 * nptsEdge - 2 + j] =
502 projPts[i][nptsEdge * (nptsEdge - j - 2)];
507 NekDouble angleCos, angleAbs, angleSign, angleSum = 0.0;
511 for (
int i = 0; i < nptsPolygon; ++i)
514 id2 = (id1 == (nptsPolygon - 1)) ? 0 : (id1 + 1);
516 for (
int j = 0; j < 3; ++j)
518 vec1[j] = ptsPolygon[j][id1] - projGloCoord[j];
519 vec2[j] = ptsPolygon[j][id2] - projGloCoord[j];
527 if ((fabs(vec1[0] + vec2[0]) + fabs(vec1[1] + vec2[1]) +
528 fabs(vec1[2] + vec2[2])) < paralTol)
538 vec3[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
539 vec3[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
540 vec3[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
543 angleSign =
Vmath::Dot(3, vec3, 1, projDir, 1) > 0.0 ? 1.0 : -1.0;
551 else if (angleCos < -1.0)
556 angleAbs = acos(angleCos);
557 angleSum += angleSign * angleAbs;
562 angleSum = fabs(angleSum);
563 if (fabs(angleSum - 2.0 * M_PI) < angleTol)
592 const NekDouble iterTol,
const int iterMax)
603 bool isConverge =
false;
604 while (cnt < iterMax)
606 tmpL = bndXmap->PhysEvaluate(etaLR, pts[dirUse[0]]);
607 tmpR = bndXmap->PhysEvaluate(etaLR + 1, pts[dirUse[0]]);
609 if (fabs(gloCoord[dirUse[0]] - tmpL) >=
610 fabs(gloCoord[dirUse[0]] - tmpR))
612 etaLR[0] = 0.5 * (etaLR[0] + etaLR[1]);
616 etaLR[1] = 0.5 * (etaLR[0] + etaLR[1]);
619 if ((etaLR[1] - etaLR[0]) < iterTol)
621 locCoord[0] = 0.5 * (etaLR[0] + etaLR[1]);
632 WARNINGL1(
false,
"Bisection iteration is not converged");
651 bndGeom->GetMetricInfo()->GetJac(bndXmap->GetPointsKeys());
654 scaledTol *= iterTol;
657 const int dir1 = dirUse[0];
658 const int dir2 = dirUse[1];
666 bndXmap->PhysDeriv(pts[dir1], Dx1D1, Dx1D2);
667 bndXmap->PhysDeriv(pts[dir2], Dx2D1, Dx2D2);
674 NekDouble derx1_1, derx1_2, derx2_1, derx2_2, jac;
677 bool isConverge =
false;
680 while (cnt++ < iterMax)
682 x1map = bndXmap->PhysEvaluate(locCoord, pts[dir1]);
683 x2map = bndXmap->PhysEvaluate(locCoord, pts[dir2]);
685 F1 = gloCoord[dir1] - x1map;
686 F2 = gloCoord[dir2] - x2map;
688 if (F1 * F1 + F2 * F2 < scaledTol)
690 resid =
sqrt(F1 * F1 + F2 * F2);
696 derx1_1 = bndXmap->PhysEvaluate(locCoord, Dx1D1);
697 derx1_2 = bndXmap->PhysEvaluate(locCoord, Dx1D2);
698 derx2_1 = bndXmap->PhysEvaluate(locCoord, Dx2D1);
699 derx2_2 = bndXmap->PhysEvaluate(locCoord, Dx2D2);
701 jac = derx2_2 * derx1_1 - derx2_1 * derx1_2;
705 locCoord[0] = locCoord[0] + (derx2_2 * (gloCoord[dir1] - x1map) -
706 derx1_2 * (gloCoord[dir2] - x2map)) /
709 locCoord[1] = locCoord[1] + (-derx2_1 * (gloCoord[dir1] - x1map) +
710 derx1_1 * (gloCoord[dir2] - x2map)) /
714 if (!(std::isfinite(locCoord[0]) && std::isfinite(locCoord[1])))
717 std::ostringstream ss;
718 ss <<
"nan or inf found in NewtonIterForLocCoordOnProjBndElmt in "
720 << bndGeom->GetGlobalID();
724 if (fabs(locCoord[0]) > LcoordDiv || fabs(locCoord[1]) > LcoordDiv)
732 bndXmap->LocCoordToLocCollapsed(locCoord, eta);
734 if (bndGeom->ClampLocCoords(eta, 0.0))
737 x1map = bndXmap->PhysEvaluate(eta, pts[dir1]);
738 x2map = bndXmap->PhysEvaluate(eta, pts[dir2]);
740 F1 = gloCoord[dir1] - x1map;
741 F2 = gloCoord[dir2] - x2map;
743 dist =
sqrt(F1 * F1 + F2 * F2);
754 bndXmap->LocCoordToLocCollapsed(locCoord, collCoords);
757 if ((collCoords[0] >= -1.0 && collCoords[0] <= 1.0) &&
758 (collCoords[1] >= -1.0 && collCoords[1] <= 1.0))
760 std::ostringstream ss;
762 ss <<
"Reached MaxIterations (" << iterMax
763 <<
") in Newton iteration ";
764 ss <<
"Init value (" << setprecision(4) << 0 <<
"," << 0 <<
","
766 ss <<
"Fin value (" << locCoord[0] <<
"," << locCoord[1] <<
","
768 ss <<
"Resid = " << resid <<
" Tolerance = " <<
sqrt(scaledTol);
800 const int npts = bndXmap->GetTotPoints();
801 const int nCoordDim =
m_f->m_exp[0]->GetCoordim(0);
808 for (
int i = 0; i < nCoordDim; ++i)
812 bndCoeffs[i] = bndGeom->GetCoeffs(i);
813 bndXmap->BwdTrans(bndCoeffs[i], pts[i]);
826 for (
int i = 1; i < nCoordDim; ++i)
828 if (fabs(projDir[i]) > fabs(projDir[dirMaxId]))
860 else if (dirMaxId == 1)
880 bool isConverge, isDesired;
883 bndGeom, projGloCoord, projPts, dirUse, locCoord, iterTol);
886 tmp[0] = bndXmap->PhysEvaluate(locCoord, pts[0]) - gloCoord[0];
887 tmp[1] = bndXmap->PhysEvaluate(locCoord, pts[1]) - gloCoord[1];
890 isDesired = (projDist > 0.0) && (projDist < maxDist);
892 return isConverge && isDesired;
905 bool isConverge, isDesired;
909 dirUse, locCoord, dist, iterTol);
913 std::ostringstream ss;
914 ss <<
"Collapsed locCoord out of range.\n"
915 <<
"Newton iteration gives the distance: " << dist;
920 tmp[0] = bndXmap->PhysEvaluate(locCoord, pts[0]) - gloCoord[0];
921 tmp[1] = bndXmap->PhysEvaluate(locCoord, pts[1]) - gloCoord[1];
922 tmp[2] = bndXmap->PhysEvaluate(locCoord, pts[2]) - gloCoord[2];
925 isDesired = (projDist > 0.0) && (projDist < maxDist);
927 return isConverge && isDesired;
947 const int nCoordDim =
m_f->m_exp[0]->GetCoordim(0);
951 int npts = bndXmap->GetTotPoints();
954 for (
int i = 0; i < nCoordDim; ++i)
957 bndCoeffs[i] = bndGeom->GetCoeffs(i);
958 bndXmap->BwdTrans(bndCoeffs[i], pts[i]);
967 bndXmap->PhysDeriv(pts[0], DxD1);
968 bndXmap->PhysDeriv(pts[1], DyD1);
971 dxd1 = bndXmap->PhysEvaluate(locCoord, DxD1);
972 dyd1 = bndXmap->PhysEvaluate(locCoord, DyD1);
973 norm =
sqrt(dxd1 * dxd1 + dyd1 * dyd1);
975 normals[0] = dyd1 / norm;
976 normals[1] = -dxd1 / norm;
988 bndXmap->PhysDeriv(pts[0], DxD1, DxD2);
989 bndXmap->PhysDeriv(pts[1], DyD1, DyD2);
990 bndXmap->PhysDeriv(pts[2], DzD1, DzD2);
992 NekDouble dxd1, dyd1, dzd1, dxd2, dyd2, dzd2;
993 dxd1 = bndXmap->PhysEvaluate(locCoord, DxD1);
994 dxd2 = bndXmap->PhysEvaluate(locCoord, DxD2);
995 dyd1 = bndXmap->PhysEvaluate(locCoord, DyD1);
996 dyd2 = bndXmap->PhysEvaluate(locCoord, DyD2);
997 dzd1 = bndXmap->PhysEvaluate(locCoord, DzD1);
998 dzd2 = bndXmap->PhysEvaluate(locCoord, DzD2);
1001 n1 = dyd1 * dzd2 - dyd2 * dzd1;
1002 n2 = dzd1 * dxd2 - dzd2 * dxd1;
1003 n3 = dxd1 * dyd2 - dxd2 * dyd1;
1004 norm =
sqrt(n1 * n1 + n2 * n2 + n3 * n3);
1006 normals[0] = n1 / norm;
1007 normals[1] = n2 / norm;
1008 normals[2] = n3 / norm;
#define WARNINGL1(condition, msg)
#define ASSERTL0(condition, msg)
A class that contains algorithms for interpolation between pts fields, expansions and different meshe...
FIELD_UTILS_EXPORT void Interpolate(const std::vector< MultiRegions::ExpListSharedPtr > expInField, std::vector< MultiRegions::ExpListSharedPtr > &expOutField, NekDouble def_value=0.0)
Interpolate from an expansion to an expansion.
FieldSharedPtr m_f
Field object.
std::map< std::string, ConfigOption > m_config
List of configuration values.
void ProjectPoint(const Array< OneD, const NekDouble > &gloCoord, const Array< OneD, const NekDouble > &projDir, const NekDouble distToOrig, Array< OneD, NekDouble > &projGloCoord)
Project a single point along the given direction to a plane.
virtual void Process(po::variables_map &vm)
Write mesh to output file.
bool BndElmtContainsPoint(SpatialDomains::GeometrySharedPtr bndGeom, const Array< OneD, const NekDouble > &gloCoord, const Array< OneD, const NekDouble > &projDir, Array< OneD, NekDouble > &locCoord, NekDouble &projDist, const NekDouble maxDist=1.0, const NekDouble iterTol=1.0e-8)
Check if a point can be projected onto an oundary element in a given direction. If yes,...
bool NewtonIterForLocCoordOnBndElmt(SpatialDomains::GeometrySharedPtr bndGeom, const Array< OneD, const NekDouble > &gloCoord, const Array< OneD, const Array< OneD, NekDouble >> &pts, const Array< OneD, const int > &dirUse, Array< OneD, NekDouble > &locCoord, NekDouble &dist, const NekDouble iterTol=1.0e-8, const int iterMax=51)
void ProjectVertices(const Array< OneD, const Array< OneD, NekDouble >> &pts, const Array< OneD, const NekDouble > &projDir, const NekDouble distToOrig, Array< OneD, Array< OneD, NekDouble >> &projPts)
Project a single point along the given direction to a plane.
void GetNormals(SpatialDomains::GeometrySharedPtr bndGeom, const Array< OneD, const NekDouble > &locCoord, Array< OneD, NekDouble > &normals)
Get the normals for a given locCoord.
bool BisectionForLocCoordOnBndElmt(SpatialDomains::GeometrySharedPtr bndGeom, const Array< OneD, const NekDouble > &gloCoord, const Array< OneD, const Array< OneD, NekDouble >> &pts, const Array< OneD, const int > &dirUse, Array< OneD, NekDouble > &locCoord, const NekDouble iterTol=1.0e-8, const int iterMax=51)
Use iteration to get the locCoord. This routine should be used after we have checked the projected po...
virtual ~ProcessWallNormalData()
bool isInProjectedArea3D(const Array< OneD, const NekDouble > &projGloCoord, const Array< OneD, const Array< OneD, NekDouble >> &projPts, const Array< OneD, const NekDouble > &projDir, const NekDouble paralTol=1.0e-12, const NekDouble angleTol=1.0e-6)
bool isInProjectedArea2D(const Array< OneD, const NekDouble > &projGloCoord, const Array< OneD, const Array< OneD, NekDouble >> &projPts, const NekDouble paralTol=1.0e-12)
Determine if the projected point is inside the projected element.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
static bool GenerateVector(const std::string &str, std::vector< T > &out)
Takes a comma-separated string and converts it to entries in a vector.
const BoundaryRegionCollection & GetBoundaryRegions(void) const
std::shared_ptr< Field > FieldSharedPtr
std::pair< ModuleType, std::string > ModuleKey
ModuleFactory & GetModuleFactory()
static std::map< PtsInfo, int > NullPtsInfoMap
static const NekDouble kNekUnsetDouble
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
std::shared_ptr< Geometry > GeometrySharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
The above copyright notice and this permission notice shall be included.
void Neg(int n, T *x, const int incx)
Negate x = -x.
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
scalarT< T > sqrt(scalarT< T > in)
Represents a command-line configuration option.