Nektar++
ProcessWallNormalData.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessWallNormalData.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Get the wall-normal data at a given origin.
32 //
33 ////////////////////////////////////////////////////////////////////////////////
34 
35 #include <iostream>
36 #include <string>
37 
41 #include <MultiRegions/ExpList.h>
42 
43 #include "ProcessWallNormalData.h"
44 
45 using namespace std;
46 
47 namespace Nektar
48 {
49 namespace FieldUtils
50 {
51 
52 ModuleKey ProcessWallNormalData::className =
54  ModuleKey(eProcessModule, "wallNormalData"),
55  ProcessWallNormalData::create,
56  "Export data in wall-normal direction from a single point on the "
57  "wall.");
58 
59 ProcessWallNormalData::ProcessWallNormalData(FieldSharedPtr f)
61 {
62  f->m_writeBndFld = false; // turned on in upstream ProcessBoundaryExtract
63 
64  m_config["xorig"] =
65  ConfigOption(false, "0.5,0.0,0.0",
66  "The point to be projected onto the wall to get the \
67  sampling origin. default=[0.5,0,0]");
68  m_config["projDir"] =
69  ConfigOption(false, "0.0,1.0,0.0",
70  "The direction to project the point onto the wall to \
71  get the sampling origin. default=[0,1,0]");
72  m_config["maxDist"] = ConfigOption(
73  false, "1.0", "Distance to limit projection distance to find the \
74  desired sampling origin. defalut=1.0");
75  m_config["distH"] = ConfigOption(
76  false, "0.01", "Sampling distance along the wall normal at the \
77  sampling origin. default=0.1");
78  m_config["nptsH"] = ConfigOption(
79  false, "5", "Number of sampling points along the wall normal. \
80  default=5");
81  m_config["d"] = ConfigOption(
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");
86 }
87 
89 {
90 }
91 
92 /*Note
93  * This module is used to get field data in the wall-normal direction.
94  * The input cases can be 2D, 2.5D and 3D.
95  * The data will be exported with .pts extension.
96  *
97  * The user defined parameters are: bnd, xorig, searchDir, maxDist, distH,
98  * nptsH, and d.
99  * - bnd=0 is the boundary id. This boundary should contain the desired origin.
100  * - xorig="x,y,z" are the coordinates of the input sampling origin. They are
101  * used as the references to get exact origin on the wall.
102  * - projDir="0,1,0" is the projection direction to find the point on the wall
103  * - maxDist=1.0 is the distance that constrains the projection on the wall.
104  * The projection distance should be smaller than this distance. This
105  * parameter is set in case the boundary is wavey so that multiple
106  * projections exist at the same time.
107  * - distH=0.01 is the sampling depth in the wall-normal direction.
108  * - nptsH=11 is the number of sampling points along wall-normal direction.
109  * - d=0.1 is a destribution control parameter of the sampling points. It
110  * should be in the range (0,inf). d in range (0,0.95] for controlled array.
111  * d>0.95 gives evenly spaced array.
112  */
113 void ProcessWallNormalData::Process(po::variables_map &vm)
114 {
116 
117  // Get dim to store data
118  const int nfields = m_f->m_variables.size();
119  const int nCoordDim = m_f->m_exp[0]->GetCoordim(0);
120  m_spacedim = nCoordDim + m_f->m_numHomogeneousDir;
121  const int nBndLcoordDim = nCoordDim - 1;
122  const int totVars = m_spacedim + m_f->m_variables.size();
123 
124  // Initialize the sampling parameters
125  std::vector<NekDouble> xorig, searchDir;
126  ASSERTL0(ParseUtils::GenerateVector(m_config["xorig"].as<string>(), xorig),
127  "Failed to interpret origin coordinates string");
128  ASSERTL0(
129  ParseUtils::GenerateVector(m_config["projDir"].as<string>(), searchDir),
130  "Failed to interpret search direction string");
131  const NekDouble maxDist = m_config["maxDist"].as<NekDouble>();
132  const NekDouble distH = m_config["distH"].as<NekDouble>();
133  const int nptsH = m_config["nptsH"].as<int>();
134  const NekDouble delta = m_config["d"].as<NekDouble>();
135 
136  Array<OneD, NekDouble> orig(3), projDir(3); // gloCoord of the origin
137  for (int i = 0; i < 3; ++i)
138  {
139  orig[i] = (xorig.size() >= (i + 1)) ? xorig[i] : 0.0;
140  projDir[i] = (searchDir.size() >= (i + 1)) ? searchDir[i] : 0.0;
141  }
142  if (nCoordDim == 2)
143  {
144  projDir[2] = 0.0;
145  }
146  Vmath::Smul(3, 1.0 / sqrt(Vmath::Dot(3, projDir, 1, projDir, 1)), projDir,
147  1, projDir, 1);
148 
149  // Update z according to the closest plane in 2.5D cases
150  if (m_f->m_numHomogeneousDir == 1)
151  {
152  int nPlanes = m_f->m_exp[0]->GetHomogeneousBasis()->GetZ().size();
153  NekDouble lHom = m_f->m_exp[0]->GetHomoLen();
154  if (orig[2] < 0.0 || orig[2] > lHom)
155  {
156  orig[2] = 0.0;
157  }
158  else
159  {
160  NekDouble dZ = lHom / nPlanes;
161  NekDouble zTmp, zCur = 0.0, distTmp, distCur = 999.0;
162  for (int i = 0; i <= nPlanes; ++i)
163  {
164  zTmp = dZ * i;
165  distTmp = fabs(orig[2] - zTmp);
166  if (distTmp < distCur)
167  {
168  distCur = distTmp;
169  zCur = zTmp;
170  if (i == nPlanes)
171  {
172  zCur = 0.0;
173  }
174  }
175  }
176  orig[2] = zCur;
177  }
178  }
179 
180  // Get the bnd id
182  m_f->m_exp[0]->GetGraph());
184  bcs.GetBoundaryRegions();
185  map<int, int> BndRegionMap;
186  int cnt = 0;
187  for (auto &breg_it : bregions)
188  {
189  BndRegionMap[breg_it.first] = cnt++;
190  }
191  int bnd = BndRegionMap[m_f->m_bndRegionsToWrite[0]];
192 
193  // Get expansion list for boundary
195  for (int i = 0; i < nfields; ++i)
196  {
197  BndExp[i] = m_f->m_exp[i]->UpdateBndCondExpansion(bnd);
198  }
199 
200  //-------------------------------------------------------------------------
201  // Find the element that contains the origin, and give the precise origin
203  Array<OneD, NekDouble> locCoord(nBndLcoordDim, -999.0);
204  NekDouble projDist;
205  int elmtid;
206  bool isInside = false;
207 
208  // Search and get precise locCoord
209  const int nElmts = BndExp[0]->GetNumElmts();
210  for (elmtid = 0; elmtid < nElmts; ++elmtid) // nElmts
211  {
212  bndGeom = BndExp[0]->GetExp(elmtid)->GetGeom();
213  isInside = BndElmtContainsPoint(bndGeom, orig, projDir, locCoord,
214  projDist, maxDist);
215  if (isInside)
216  {
217  break;
218  }
219  }
220  ASSERTL0(isInside, "Failed to find the sampling origin on the boundary.");
221 
222  // Then Update the precise sampling position
223  StdRegions::StdExpansionSharedPtr bndXmap = bndGeom->GetXmap();
224  const int npts = bndXmap->GetTotPoints();
225  Array<OneD, Array<OneD, NekDouble>> pts(nCoordDim);
226  Array<OneD, Array<OneD, const NekDouble>> bndCoeffs(nCoordDim);
227 
228  for (int i = 0; i < nCoordDim; ++i)
229  {
230  pts[i] = Array<OneD, NekDouble>(npts);
231  bndCoeffs[i] = bndGeom->GetCoeffs(i); // 0/1/2 for x/y/z
232  bndXmap->BwdTrans(bndCoeffs[i], pts[i]);
233  orig[i] = bndXmap->PhysEvaluate(locCoord, pts[i]); // Update wall point
234  }
235 
236  // Get outward-pointing normal vectors for all quadrature points on bnd
237  // Use these normals as diretion references
239  m_f->m_exp[0]->GetBoundaryNormals(bnd, normalsQ);
240 
241  // Get the normals in the element
242  // Set point key to get point id
243  Array<OneD, LibUtilities::PointsKey> from_key(nBndLcoordDim);
244  int from_nPtsPerElmt = 1;
245  for (int i = 0; i < nBndLcoordDim; ++i)
246  {
247  from_key[i] = BndExp[0]->GetExp(elmtid)->GetBasis(i)->GetPointsKey();
248  from_nPtsPerElmt *= from_key[i].GetNumPoints();
249  }
250 
251  // Get ref normal
252  Array<OneD, NekDouble> normalsRef(3, 0.0);
253  int refId = elmtid * from_nPtsPerElmt + 0; // the 1st normal in the bnd elmt
254  for (int i = 0; i < m_spacedim; ++i)
255  {
256  normalsRef[i] = normalsQ[i][refId];
257  }
258 
259  // Get the precise normals and correct the direction to be inward
260  Array<OneD, NekDouble> normals(3, 0.0);
261  GetNormals(bndGeom, locCoord, normals);
262 
263  if (Vmath::Dot(3, normals, normalsRef) > 0.0)
264  {
265  Vmath::Neg(3, normals, 1);
266  }
267 
268  // Output the info if -v is set
269  if (m_f->m_verbose)
270  {
271  cout << "------ wallNormalData module ------\n";
272  cout << "Input point:\n";
273  cout << " - [Px,Py,Pz] = [" << xorig[0] << ", " << xorig[1];
274  if (xorig.size() >= 3)
275  {
276  cout << ", " << xorig[2] << "]\n";
277  }
278  else
279  {
280  cout << ", " << 0.0 << "]\n";
281  }
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] << ", "
287  << orig[2] << "]\n";
288  cout << "Normals at the origin:\n";
289  cout << " - [nx,ny,nz] = [" << normals[0] << ", " << normals[1] << ", "
290  << normals[2] << "]\n";
291  cout
292  << "Ref normals (at quadrature points in the projected element):\n";
293  for (int i = 0; i < from_nPtsPerElmt; ++i)
294  {
295  cout << " - " << i << ": [nx,ny,nz] = ["
296  << -normalsQ[0][elmtid * from_nPtsPerElmt + i] << ", "
297  << -normalsQ[1][elmtid * from_nPtsPerElmt + i];
298  if (m_spacedim == 3)
299  {
300  cout << ", " << -normalsQ[2][elmtid * from_nPtsPerElmt + i]
301  << "]\n";
302  }
303  else
304  {
305  cout << "]\n";
306  }
307  }
308  cout << endl;
309  }
310 
311  //-------------------------------------------------------------------------
312  // Set the depth of sampling array
313  Array<OneD, NekDouble> h(nptsH, 0.0);
314  if (delta > 0 && delta <= 0.95)
315  {
316  // Use expression in Agrawal's paper:
317  // h = 1-tanh((1-ksi)*atanh(sqrt(1-delta)))/sqrt(1-delta), ksi in [0,1]
318  NekDouble tmp1;
319  const NekDouble tmp2 = 1.0 / (static_cast<NekDouble>(nptsH) - 1.0);
320  const NekDouble tmp3 = sqrt(1.0 - delta);
321  const NekDouble tmp4 = atanh(tmp3);
322  const NekDouble tmp5 = 1.0 / tmp3;
323  for (int i = 1; i < nptsH; ++i)
324  {
325  tmp1 = 1.0 - i * tmp2; // tmp1 = 1-ksi
326  h[i] = 1 - tanh(tmp1 * tmp4) * tmp5;
327  h[i] *= distH; // physical distance in normal direction
328  }
329  }
330  else if (delta > 0.95)
331  {
332  // Evenly spaced array
333  const NekDouble tmp1 = 1.0 / (static_cast<NekDouble>(nptsH) - 1.0);
334  for (int i = 1; i < nptsH; ++i)
335  {
336  h[i] = i * tmp1;
337  h[i] *= distH; // physical distance in normal direction
338  }
339  }
340  else
341  {
342  ASSERTL0(false, "Input error. Delta needs to be greater than 0.0.");
343  }
344 
345  // Set pts coordinates and interpoate the data
346  Array<OneD, Array<OneD, NekDouble>> ptsH(totVars);
347  for (int i = 0; i < totVars; ++i)
348  {
349  ptsH[i] = Array<OneD, NekDouble>(nptsH, 0.0);
350  }
351 
352  for (int i = 0; i < m_spacedim; ++i)
353  {
354  for (int j = 0; j < nptsH; ++j)
355  {
356  ptsH[i][j] = orig[i] + h[j] * normals[i]; // x0+dist*nx
357  }
358  }
359 
361  m_spacedim, m_f->m_variables, ptsH, LibUtilities::NullPtsInfoMap);
362 
363  Interpolator interp;
364  interp.Interpolate(m_f->m_exp, m_f->m_fieldPts,
366 }
367 
368 /**
369  * @brief Project a single point along the given direction to a plane
370  * @param gloCoord Global coordinate of the point. size=3.
371  * @param projDir Projection direction, which is also the normal vector of
372  * the target plane. size=3, norm=1.
373  * @param distToOrig The distance from the origin (0,0,0) to the target plane
374  * @param projGloCoord The global coordinate of the projecion result.
375  */
377  const Array<OneD, const NekDouble> &gloCoord,
378  const Array<OneD, const NekDouble> &projDir, const NekDouble distToOrig,
379  Array<OneD, NekDouble> &projGloCoord)
380 {
381  NekDouble tmp1;
382  Array<OneD, NekDouble> tmp2(3);
383 
384  tmp1 = Vmath::Dot(3, gloCoord, 1, projDir, 1); // |xn| = x dot n
385  Vmath::Smul(3, distToOrig - tmp1, projDir, 1, tmp2,
386  1); // d_xn = (dist-|xn|) n
387  Vmath::Vadd(3, gloCoord, 1, tmp2, 1, projGloCoord, 1); // x' = x + d_xn
388 }
389 
390 /**
391  * @brief Project the vertices of the elmt along the given direction to a plane
392  * @param pts Global coordinate of the vertices of the elmt. size=2/3.
393  * @param projDir Projection direction, which is also the normal vector of
394  * the target plane. size=3, norm=1.
395  * @param distToOrig The distance from the origin (0,0,0) to the target plane.
396  * @param projPts The global coordinate of the projecion result.
397  */
399  const Array<OneD, const Array<OneD, NekDouble>> &pts,
400  const Array<OneD, const NekDouble> &projDir, const NekDouble distToOrig,
401  Array<OneD, Array<OneD, NekDouble>> &projPts)
402 {
403  const int nCoordDim = pts.size(); // 2 for 2.5D cases
404  const int npts = pts[0].size();
405 
406  NekDouble tmp1;
407  Array<OneD, NekDouble> singlePnt(nCoordDim), tmp2(nCoordDim);
408 
409  for (int i = 0; i < npts; ++i)
410  {
411  // Get a point
412  for (int j = 0; j < nCoordDim; ++j)
413  {
414  singlePnt[j] = pts[j][i];
415  }
416 
417  // Projection
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);
421 
422  // Save to projPts
423  for (int j = 0; j < nCoordDim; ++j)
424  {
425  projPts[j][i] = singlePnt[j];
426  }
427  }
428 }
429 
430 /**
431  * @brief Determine if the projected point is inside the projected element.
432  * @param projGloCoord The global coordinate of the projected single point.
433  * @param projPts The global coordinate of the projected vertices,size=2/3
434  * @param projDir Projection direction, which is used as the reference
435  * direction in the 3D routine. size=3, norm=1.
436  * @param paralTol Tolerence to check if two vectors are parallel.
437  * @param angleTol Tolerence to check if the total angle is 2*pi.
438  * @return Inside (true) or not (false)
439  */
441  const Array<OneD, const NekDouble> &projGloCoord,
442  const Array<OneD, const Array<OneD, NekDouble>> &projPts,
443  const NekDouble paralTol)
444 {
445  const NekDouble npts = projPts[0].size();
446 
447  Array<OneD, NekDouble> vec1(2, 0.0), vec2(2, 0.0);
448 
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];
453 
454  Vmath::Smul(2, 1.0 / sqrt(Vmath::Dot(2, vec1, 1, vec1, 1)), vec1, 1, vec1,
455  1);
456  Vmath::Smul(2, 1.0 / sqrt(Vmath::Dot(2, vec2, 1, vec2, 1)), vec2, 1, vec2,
457  1);
458 
459  if ((fabs(vec1[0] + vec2[0]) + fabs(vec1[1] + vec2[1])) < paralTol)
460  {
461  // On the line segement, true
462  return true;
463  }
464  else
465  {
466  return false;
467  }
468 }
469 
471  const Array<OneD, const NekDouble> &projGloCoord,
472  const Array<OneD, const Array<OneD, NekDouble>> &projPts,
473  const Array<OneD, const NekDouble> &projDir, const NekDouble paralTol,
474  const NekDouble angleTol)
475 {
476 
477  // Generate a polygen (quad), re-order the points on the edge
478  const NekDouble npts = projPts[0].size();
479  const int nptsEdge =
480  sqrt(npts); // num of points on edge for geo representation
481  const int nptsPolygon = 4 * nptsEdge - 4;
482  Array<OneD, Array<OneD, NekDouble>> ptsPolygon(3);
483  for (int i = 0; i < 3; ++i)
484  {
485  ptsPolygon[i] = Array<OneD, NekDouble>(nptsPolygon);
486 
487  for (int j = 0; j < nptsEdge; ++j)
488  {
489  ptsPolygon[i][j] = projPts[i][j];
490  }
491  for (int j = 0; j < nptsEdge - 2; ++j)
492  {
493  ptsPolygon[i][nptsEdge + j] = projPts[i][(j + 2) * nptsEdge - 1];
494  }
495  for (int j = 0; j < nptsEdge; ++j)
496  {
497  ptsPolygon[i][2 * nptsEdge - 2 + j] = projPts[i][npts - 1 - j];
498  }
499  for (int j = 0; j < nptsEdge - 2; ++j)
500  {
501  ptsPolygon[i][3 * nptsEdge - 2 + j] =
502  projPts[i][nptsEdge * (nptsEdge - j - 2)];
503  }
504  }
505 
506  // Determine relation using angle method (sum=2*pi)
507  NekDouble angleCos, angleAbs, angleSign, angleSum = 0.0;
508  Array<OneD, NekDouble> vec1(3, 0.0), vec2(3, 0.0), vec3(3, 0.0);
509  int id1, id2;
510 
511  for (int i = 0; i < nptsPolygon; ++i)
512  {
513  id1 = i;
514  id2 = (id1 == (nptsPolygon - 1)) ? 0 : (id1 + 1);
515 
516  for (int j = 0; j < 3; ++j)
517  {
518  vec1[j] = ptsPolygon[j][id1] - projGloCoord[j];
519  vec2[j] = ptsPolygon[j][id2] - projGloCoord[j];
520  }
521 
522  Vmath::Smul(3, 1.0 / sqrt(Vmath::Dot(3, vec1, 1, vec1, 1)), vec1, 1,
523  vec1, 1);
524  Vmath::Smul(3, 1.0 / sqrt(Vmath::Dot(3, vec2, 1, vec2, 1)), vec2, 1,
525  vec2, 1);
526 
527  if ((fabs(vec1[0] + vec2[0]) + fabs(vec1[1] + vec2[1]) +
528  fabs(vec1[2] + vec2[2])) < paralTol)
529  {
530  // On the line segement, true
531  // This branch is used to aviod angle=pi but angleSign=-1 (not 1)
532  return true;
533  }
534  else
535  {
536  // Off the line segement, compute angle
537  // vec3 = vec1 x vec2, not scaled
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];
541 
542  // Use the projDir as reference direction (positive)
543  angleSign = Vmath::Dot(3, vec3, 1, projDir, 1) > 0.0 ? 1.0 : -1.0;
544  angleCos = Vmath::Dot(3, vec1, 1, vec2, 1);
545 
546  // Avoid error in using acos()
547  if (angleCos > 1.0)
548  {
549  angleCos = 1.0;
550  }
551  else if (angleCos < -1.0)
552  {
553  angleCos = -1.0;
554  }
555 
556  angleAbs = acos(angleCos);
557  angleSum += angleSign * angleAbs;
558  }
559  }
560 
561  // Check
562  angleSum = fabs(angleSum);
563  if (fabs(angleSum - 2.0 * M_PI) < angleTol)
564  {
565  return true;
566  }
567  else
568  {
569  return false;
570  }
571 }
572 
573 /**
574  * @brief Use iteration to get the locCoord. This routine should be used after
575  * we have checked the projected point is inside the projected element.
576  * @param bndGeom Geometry to get the xmap.
577  * @param gloCoord Global coordinate of the point. size=3.
578  * @param pts Global coordinate of the vertices of the elmt. size=2/3.
579  * @param dieUse The main direction(s) used to compute local coordinate
580  * @param locCoord Iteration results for local coordinate(s)
581  * @param dist Returned distance in physical space if the collapsed
582  * locCoord is out of range [-1,1].
583  * @param iterTol Tolerence for iteration.
584  * @param iterMax Maximum iteration steps
585  * @return Converged (true) or not (false)
586  */
589  const Array<OneD, const NekDouble> &gloCoord,
590  const Array<OneD, const Array<OneD, NekDouble>> &pts,
591  const Array<OneD, const int> &dirUse, Array<OneD, NekDouble> &locCoord,
592  const NekDouble iterTol, const int iterMax)
593 {
594  // Initial settings
595  Array<OneD, NekDouble> etaLR(2); // range [-1,1]
596  etaLR[0] = -1.0; // left
597  etaLR[1] = 1.0; // right
598  NekDouble tmpL, tmpR; // tmp values for L/R
599 
600  StdRegions::StdExpansionSharedPtr bndXmap = bndGeom->GetXmap();
601  locCoord[0] = -2.0;
602  int cnt = 0;
603  bool isConverge = false;
604  while (cnt < iterMax)
605  {
606  tmpL = bndXmap->PhysEvaluate(etaLR, pts[dirUse[0]]);
607  tmpR = bndXmap->PhysEvaluate(etaLR + 1, pts[dirUse[0]]);
608 
609  if (fabs(gloCoord[dirUse[0]] - tmpL) >=
610  fabs(gloCoord[dirUse[0]] - tmpR))
611  {
612  etaLR[0] = 0.5 * (etaLR[0] + etaLR[1]);
613  }
614  else
615  {
616  etaLR[1] = 0.5 * (etaLR[0] + etaLR[1]);
617  }
618 
619  if ((etaLR[1] - etaLR[0]) < iterTol)
620  {
621  locCoord[0] = 0.5 * (etaLR[0] + etaLR[1]);
622  isConverge = true;
623  break;
624  }
625 
626  ++cnt;
627  }
628 
629  // Warning if failed
630  if (cnt >= iterMax)
631  {
632  WARNINGL1(false, "Bisection iteration is not converged");
633  }
634 
635  return isConverge;
636 }
637 
640  const Array<OneD, const NekDouble> &gloCoord,
641  const Array<OneD, const Array<OneD, NekDouble>> &pts,
642  const Array<OneD, const int> &dirUse, Array<OneD, NekDouble> &locCoord,
643  NekDouble &dist, const NekDouble iterTol, const int iterMax)
644 {
645 
646  const NekDouble LcoordDiv = 15.0;
647 
648  StdRegions::StdExpansionSharedPtr bndXmap = bndGeom->GetXmap();
649 
651  bndGeom->GetMetricInfo()->GetJac(bndXmap->GetPointsKeys());
652  NekDouble scaledTol =
653  Vmath::Vsum(Jac.size(), Jac, 1) / ((NekDouble)Jac.size());
654  scaledTol *= iterTol;
655 
656  // Set the gloCoord used to compute locCoord
657  const int dir1 = dirUse[0];
658  const int dir2 = dirUse[1];
659 
660  Array<OneD, NekDouble> Dx1D1(pts[dir1].size());
661  Array<OneD, NekDouble> Dx1D2(pts[dir1].size());
662  Array<OneD, NekDouble> Dx2D1(pts[dir1].size());
663  Array<OneD, NekDouble> Dx2D2(pts[dir1].size());
664 
665  // Ideally this will be stored in m_geomfactors
666  bndXmap->PhysDeriv(pts[dir1], Dx1D1, Dx1D2);
667  bndXmap->PhysDeriv(pts[dir2], Dx2D1, Dx2D2);
668 
669  // Initial the locCoord, in [-1,1]
670  locCoord[0] = 0.0;
671  locCoord[1] = 0.0;
672 
673  NekDouble x1map, x2map, F1, F2;
674  NekDouble derx1_1, derx1_2, derx2_1, derx2_2, jac;
675  NekDouble resid;
676  int cnt = 0;
677  bool isConverge = false;
678 
679  F1 = F2 = 2000; // Starting value of Function
680  while (cnt++ < iterMax)
681  {
682  x1map = bndXmap->PhysEvaluate(locCoord, pts[dir1]);
683  x2map = bndXmap->PhysEvaluate(locCoord, pts[dir2]);
684 
685  F1 = gloCoord[dir1] - x1map;
686  F2 = gloCoord[dir2] - x2map;
687 
688  if (F1 * F1 + F2 * F2 < scaledTol)
689  {
690  resid = sqrt(F1 * F1 + F2 * F2);
691  isConverge = true;
692  break;
693  }
694 
695  // Interpolate derivative metric at locCoord
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);
700 
701  jac = derx2_2 * derx1_1 - derx2_1 * derx1_2;
702 
703  // Use analytical inverse of derivitives which are
704  // also similar to those of metric factors.
705  locCoord[0] = locCoord[0] + (derx2_2 * (gloCoord[dir1] - x1map) -
706  derx1_2 * (gloCoord[dir2] - x2map)) /
707  jac;
708 
709  locCoord[1] = locCoord[1] + (-derx2_1 * (gloCoord[dir1] - x1map) +
710  derx1_1 * (gloCoord[dir2] - x2map)) /
711  jac;
712 
713  // locCoord have diverged so stop iteration
714  if (!(std::isfinite(locCoord[0]) && std::isfinite(locCoord[1])))
715  {
716  dist = 1e16;
717  std::ostringstream ss;
718  ss << "nan or inf found in NewtonIterForLocCoordOnProjBndElmt in "
719  "element "
720  << bndGeom->GetGlobalID();
721  WARNINGL1(false, ss.str());
722  return false;
723  }
724  if (fabs(locCoord[0]) > LcoordDiv || fabs(locCoord[1]) > LcoordDiv)
725  {
726  break;
727  }
728  }
729 
730  // Check distance for collapsed coordinate
731  Array<OneD, NekDouble> eta(2);
732  bndXmap->LocCoordToLocCollapsed(locCoord, eta);
733 
734  if (bndGeom->ClampLocCoords(eta, 0.0))
735  {
736  // calculate the global point corresponding to locCoord
737  x1map = bndXmap->PhysEvaluate(eta, pts[dir1]);
738  x2map = bndXmap->PhysEvaluate(eta, pts[dir2]);
739 
740  F1 = gloCoord[dir1] - x1map;
741  F2 = gloCoord[dir2] - x2map;
742 
743  dist = sqrt(F1 * F1 + F2 * F2);
744  }
745  else
746  {
747  dist = 0.0;
748  }
749 
750  // Warning if failed
751  if (cnt >= iterMax)
752  {
753  Array<OneD, NekDouble> collCoords(2);
754  bndXmap->LocCoordToLocCollapsed(locCoord, collCoords);
755 
756  // if coordinate is inside element dump error!
757  if ((collCoords[0] >= -1.0 && collCoords[0] <= 1.0) &&
758  (collCoords[1] >= -1.0 && collCoords[1] <= 1.0))
759  {
760  std::ostringstream ss;
761 
762  ss << "Reached MaxIterations (" << iterMax
763  << ") in Newton iteration ";
764  ss << "Init value (" << setprecision(4) << 0 << "," << 0 << ","
765  << ") ";
766  ss << "Fin value (" << locCoord[0] << "," << locCoord[1] << ","
767  << ") ";
768  ss << "Resid = " << resid << " Tolerance = " << sqrt(scaledTol);
769 
770  WARNINGL1(cnt < iterMax, ss.str());
771  }
772  }
773 
774  return isConverge;
775 }
776 
777 /**
778  * @brief Check if a point can be projected onto an oundary element in a given
779  * direction. If yes, give the local coordinates of the projected point.
780  * we have checked the projected point is inside the projected element.
781  * @param bndGeom Pointer to the geometry of the boundary element.
782  * @param gloCoord Global coordinate of the point. size=3.
783  * @param projDir Projection direction, which is used as the reference
784  * direction in the 3D routine. size=3, norm=1.
785  * @param locCoord Iteration results for local coordinates (if inside).
786  * @param projDist Projection distance betweem the point to the wall point.
787  * @param maxDist Disntance to check if the wall point is desired.
788  * @param iterTol Tolerence for iteration.
789  * @return Inside (true) or not (false)
790  */
793  const Array<OneD, const NekDouble> &gloCoord,
794  const Array<OneD, const NekDouble> &projDir,
795  Array<OneD, NekDouble> &locCoord, NekDouble &projDist,
796  const NekDouble maxDist, const NekDouble iterTol)
797 {
798  // Get variables
799  StdRegions::StdExpansionSharedPtr bndXmap = bndGeom->GetXmap();
800  const int npts = bndXmap->GetTotPoints();
801  const int nCoordDim = m_f->m_exp[0]->GetCoordim(0); // =2 for 2.5D cases
802 
803  Array<OneD, Array<OneD, const NekDouble>> bndCoeffs(nCoordDim);
804  Array<OneD, Array<OneD, NekDouble>> pts(nCoordDim);
805  Array<OneD, Array<OneD, NekDouble>> projPts(nCoordDim);
806  Array<OneD, NekDouble> projGloCoord(3, 0.0);
807 
808  for (int i = 0; i < nCoordDim; ++i)
809  {
810  pts[i] = Array<OneD, NekDouble>(npts);
811  projPts[i] = Array<OneD, NekDouble>(npts);
812  bndCoeffs[i] = bndGeom->GetCoeffs(i); // 0/1/2 for x/y/z
813  bndXmap->BwdTrans(bndCoeffs[i], pts[i]);
814  }
815 
816  // Project the point and vertices of the element in the input direction
817  ProjectPoint(gloCoord, projDir, 0.0, projGloCoord);
818  ProjectVertices(pts, projDir, 0.0, projPts);
819 
820  // Set the main direction(s) and the minor direction
821  // The gloCoord for minor direction will not be used for locCoord iteration
822  // dirUse[0] is the main dir for 2D/2.5D cases, dirUse[1]/[2] is the minor
823  // one dirUse[0] and dirUse[1] are the main dir for 3D cases, dirUse[2] is
824  // hte minor one
825  int dirMaxId = 0; // id to get the dir with largest projDir component
826  for (int i = 1; i < nCoordDim; ++i)
827  {
828  if (fabs(projDir[i]) > fabs(projDir[dirMaxId]))
829  {
830  dirMaxId = i;
831  }
832  }
833 
834  Array<OneD, int> dirUse(3, 0);
835  if (nCoordDim == 2)
836  {
837  // 2D or 2.5D cases
838  if (dirMaxId == 0)
839  {
840  dirUse[0] = 1;
841  dirUse[1] = 0;
842  dirUse[2] = 2;
843  }
844  else
845  {
846  dirUse[0] = 0;
847  dirUse[1] = 1;
848  dirUse[2] = 2;
849  }
850  }
851  else
852  {
853  // 3D cases
854  if (dirMaxId == 0)
855  {
856  dirUse[0] = 1;
857  dirUse[1] = 2;
858  dirUse[2] = 0;
859  }
860  else if (dirMaxId == 1)
861  {
862  dirUse[0] = 2;
863  dirUse[1] = 0;
864  dirUse[2] = 1;
865  }
866  else
867  {
868  dirUse[0] = 0;
869  dirUse[1] = 1;
870  dirUse[2] = 2;
871  }
872  }
873 
874  // Check if the projected point is in the projected elmt
875  // If yes, then compute the locCoord and check if desired point is found
876  if (nCoordDim == 2)
877  {
878  if (isInProjectedArea2D(projGloCoord, projPts, 1.0e-12))
879  {
880  bool isConverge, isDesired;
881 
882  isConverge = BisectionForLocCoordOnBndElmt(
883  bndGeom, projGloCoord, projPts, dirUse, locCoord, iterTol);
884 
885  Array<OneD, NekDouble> tmp(2, 0.0);
886  tmp[0] = bndXmap->PhysEvaluate(locCoord, pts[0]) - gloCoord[0];
887  tmp[1] = bndXmap->PhysEvaluate(locCoord, pts[1]) - gloCoord[1];
888  projDist = Vmath::Dot(2, tmp, 1, projDir, 1); // can be negative
889 
890  isDesired = (projDist > 0.0) && (projDist < maxDist);
891 
892  return isConverge && isDesired;
893  }
894  else
895  {
896  return false;
897  }
898  }
899  else
900  {
901  if (isInProjectedArea3D(projGloCoord, projPts, projDir, 1.0e-12,
902  1.0e-6))
903  {
904  NekDouble dist;
905  bool isConverge, isDesired;
906 
907  isConverge =
908  NewtonIterForLocCoordOnBndElmt(bndGeom, projGloCoord, projPts,
909  dirUse, locCoord, dist, iterTol);
910 
911  if (dist > iterTol)
912  {
913  std::ostringstream ss;
914  ss << "Collapsed locCoord out of range.\n"
915  << "Newton iteration gives the distance: " << dist;
916  WARNINGL1(false, ss.str());
917  }
918 
919  Array<OneD, NekDouble> tmp(3, 0.0);
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];
923  projDist = Vmath::Dot(3, tmp, 1, projDir, 1); // can be negative
924 
925  isDesired = (projDist > 0.0) && (projDist < maxDist);
926 
927  return isConverge && isDesired;
928  }
929  else
930  {
931  return false;
932  }
933  }
934 }
935 
936 /**
937  * @brief Get the normals for a given locCoord
938  * @param bndGeom Pointer to the geometry of the boundary element.
939  * @param locCoord Iteration results for local coordinates (if inside).
940  * @param normals Wall normal as the result
941  */
944  const Array<OneD, const NekDouble> &locCoord,
945  Array<OneD, NekDouble> &normals)
946 {
947  const int nCoordDim = m_f->m_exp[0]->GetCoordim(0); // =2 for 2.5D cases
948  StdRegions::StdExpansionSharedPtr bndXmap = bndGeom->GetXmap();
951  int npts = bndXmap->GetTotPoints();
952 
953  // Get pts in the element
954  for (int i = 0; i < nCoordDim; ++i)
955  {
956  pts[i] = Array<OneD, NekDouble>(npts);
957  bndCoeffs[i] = bndGeom->GetCoeffs(i); // 0/1/2 for x/y/z
958  bndXmap->BwdTrans(bndCoeffs[i], pts[i]);
959  }
960 
961  // Get the outward-pointing normals according to the given locCoord
962  if (nCoordDim == 2)
963  {
964  Array<OneD, NekDouble> DxD1(pts[0].size());
965  Array<OneD, NekDouble> DyD1(pts[0].size());
966 
967  bndXmap->PhysDeriv(pts[0], DxD1);
968  bndXmap->PhysDeriv(pts[1], DyD1);
969 
970  NekDouble dxd1, dyd1, norm;
971  dxd1 = bndXmap->PhysEvaluate(locCoord, DxD1);
972  dyd1 = bndXmap->PhysEvaluate(locCoord, DyD1);
973  norm = sqrt(dxd1 * dxd1 + dyd1 * dyd1);
974 
975  normals[0] = dyd1 / norm;
976  normals[1] = -dxd1 / norm;
977  normals[2] = 0.0;
978  }
979  else
980  {
981  Array<OneD, NekDouble> DxD1(pts[0].size());
982  Array<OneD, NekDouble> DxD2(pts[0].size());
983  Array<OneD, NekDouble> DyD1(pts[0].size());
984  Array<OneD, NekDouble> DyD2(pts[0].size());
985  Array<OneD, NekDouble> DzD1(pts[0].size());
986  Array<OneD, NekDouble> DzD2(pts[0].size());
987 
988  bndXmap->PhysDeriv(pts[0], DxD1, DxD2);
989  bndXmap->PhysDeriv(pts[1], DyD1, DyD2);
990  bndXmap->PhysDeriv(pts[2], DzD1, DzD2);
991 
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);
999 
1000  NekDouble n1, n2, n3, norm;
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);
1005 
1006  normals[0] = n1 / norm;
1007  normals[1] = n2 / norm;
1008  normals[2] = n3 / norm;
1009  }
1010 }
1011 
1012 } // namespace FieldUtils
1013 } // namespace Nektar
#define WARNINGL1(condition, msg)
Definition: ErrorUtil.hpp:250
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
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.
Definition: Module.h:225
std::map< std::string, ConfigOption > m_config
List of configuration values.
Definition: Module.h:228
This processing module sets up for the boundary field to be extracted.
virtual void Process(po::variables_map &vm)
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...
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.
Definition: NekFactory.hpp:198
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.
Definition: ParseUtils.cpp:129
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:234
std::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:989
std::pair< ModuleType, std::string > ModuleKey
Definition: Module.h:285
ModuleFactory & GetModuleFactory()
Definition: Module.cpp:49
static std::map< PtsInfo, int > NullPtsInfoMap
Definition: PtsField.h:70
static const NekDouble kNekUnsetDouble
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:210
std::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
double NekDouble
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:895
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:1100
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.
Definition: Vmath.cpp:359
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291
Represents a command-line configuration option.
Definition: Module.h:131