mappedInternalPatchBase.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
27 #include "SubField.H"
28 #include "Time.H"
29 #include "triPointRef.H"
30 #include "treeDataCell.H"
31 #include "indexedOctree.H"
32 #include "globalIndex.H"
33 #include "RemoteData.H"
34 #include "OBJstream.H"
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
42 
43  template<>
44  const char* Foam::NamedEnum
45  <
47  2
48  >::names[] =
49  {
50  "normal",
51  "direction"
52  };
53 }
54 
55 
58 
59 
60 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
61 
64 (
65  const dictionary& dict
66 ) const
67 {
68  if (dict.found("offsetMode"))
69  {
70  return offsetModeNames_.read(dict.lookup("offsetMode"));
71  }
72  else
73  {
74  const bool haveDistance = dict.found("distance");
75  const bool haveOffset = dict.found("offset");
76 
77  if (haveDistance == haveOffset)
78  {
79  // Error. Demand "offsetMode" setting to disambiguate.
80  return offsetModeNames_.read(dict.lookup("offsetMode"));
81  }
82  else if (haveDistance)
83  {
84  return NORMAL;
85  }
86  else
87  {
88  return DIRECTION;
89  }
90  }
91 }
92 
93 
95 {
96  if (mapPtr_.valid())
97  {
99  << "Mapping already calculated" << exit(FatalError);
100  }
101 
102  const polyMesh& nbrMesh = this->nbrMesh();
103 
104  const globalIndex patchGlobalIndex(patch_.size());
105 
106  // Find processor and cell/face indices of samples
107  labelList sampleGlobalPatchFaces, sampleIndices;
108  {
109  // Gather the sample points into a single globally indexed list
110  List<point> allPoints(patchGlobalIndex.size());
111  {
112  List<pointField> procSamplePoints(Pstream::nProcs());
113  procSamplePoints[Pstream::myProcNo()] = this->samplePoints();
114  Pstream::gatherList(procSamplePoints);
115  Pstream::scatterList(procSamplePoints);
116 
117  forAll(procSamplePoints, proci)
118  {
119  forAll(procSamplePoints[proci], procSamplei)
120  {
121  allPoints
122  [
123  patchGlobalIndex.toGlobal(proci, procSamplei)
124  ] = procSamplePoints[proci][procSamplei];
125  }
126  }
127  }
128 
129  // List of possibly remote sampling cells
130  List<RemoteData<scalar>> allNearest(patchGlobalIndex.size());
131 
132  // Find containing cell for every sampling point
133  const indexedOctree<Foam::treeDataCell>& tree = nbrMesh.cellTree();
134  forAll(allPoints, alli)
135  {
136  const point& p = allPoints[alli];
137 
138  const label celli = tree.findInside(p);
139 
140  if (celli != -1)
141  {
142  const point& cc = nbrMesh.cellCentres()[celli];
143 
144  allNearest[alli].proci = Pstream::myProcNo();
145  allNearest[alli].elementi = celli;
146  allNearest[alli].data = magSqr(cc - p);
147  }
148  }
149 
150  // Find nearest. Combine on master.
152  (
153  allNearest,
155  );
156  Pstream::listCombineScatter(allNearest);
157 
158  // Determine the number of missing samples (no reduction necessary as
159  // this is computed from synchronised data)
160  label nNotFound = 0;
161  forAll(allPoints, alli)
162  {
163  if (allNearest[alli].proci == -1)
164  {
165  nNotFound ++;
166  }
167  }
168 
169  // If any points were not found within cells then re-search for them
170  // using a nearest test, which should not fail. Warn that this is
171  // happening. If any points were not found for some other method, then
172  // fail.
173  if (nNotFound)
174  {
176  << "Did not find a containing cell for " << nNotFound
177  << " out of " << returnReduce(patch_.size(), sumOp<label>())
178  << " total faces. Using nearest cell for these faces instead."
179  << endl << "On patch " << patch_.name() << " on region "
180  << nbrRegionName() << " with offset mode "
181  << offsetModeNames_[offsetMode_] << "." << endl;
182 
183  const indexedOctree<Foam::treeDataCell>& tree = nbrMesh.cellTree();
184  forAll(allPoints, alli)
185  {
186  const point& p = allPoints[alli];
187 
188  if (allNearest[alli].proci == -1)
189  {
190  const pointIndexHit pih = tree.findNearest(p, sqr(great));
191 
192  allNearest[alli].proci = Pstream::myProcNo();
193  allNearest[alli].elementi = pih.index();
194  allNearest[alli].data = magSqr(pih.hitPoint() - p);
195  }
196  }
197 
199  (
200  allNearest,
202  );
203  Pstream::listCombineScatter(allNearest);
204  }
205 
206  // Build lists of samples
207  DynamicList<label> samplePatchGlobalFacesDyn;
208  DynamicList<label> sampleIndicesDyn;
209  forAll(allNearest, alli)
210  {
211  if (allNearest[alli].proci == Pstream::myProcNo())
212  {
213  samplePatchGlobalFacesDyn.append(alli);
214  sampleIndicesDyn.append(allNearest[alli].elementi);
215  }
216  }
217  sampleGlobalPatchFaces.transfer(samplePatchGlobalFacesDyn);
218  sampleIndices.transfer(sampleIndicesDyn);
219  }
220 
221  // Construct distribution schedule
222  List<Map<label>> compactMap;
223  mapPtr_.reset
224  (
225  new distributionMap
226  (
227  patchGlobalIndex,
228  sampleGlobalPatchFaces,
229  compactMap
230  )
231  );
232  const labelList oldToNew(move(sampleGlobalPatchFaces));
233  const labelList oldSampleIndices(move(sampleIndices));
234 
235  // Construct input mapping for data to be distributed
236  cellIndices_ = labelList(mapPtr_->constructSize(), -1);
237  UIndirectList<label>(cellIndices_, oldToNew) = oldSampleIndices;
238 
239  // Reverse the map. This means the map is "forward" when going from cells
240  // to this patch, which is logical.
241  mapPtr_.reset
242  (
243  new distributionMap
244  (
245  patch_.size(),
246  move(mapPtr_->constructMap()),
247  move(mapPtr_->subMap())
248  )
249  );
250 
251  // Dump connecting lines for debugging
252  if (debug)
253  {
254  OBJstream obj
255  (
256  patch_.name()
257  + "_processor"
259  + ".obj"
260  );
261 
262  // Cells -> Patch
263  pointField ccs(nbrMesh.cellCentres(), cellIndices_);
264  mapPtr_->distribute(ccs);
265  forAll(patch_, patchFacei)
266  {
267  const point& fc = patch_.faceCentres()[patchFacei];
268  const point mid = 0.51*fc + 0.49*ccs[patchFacei];
269  obj.write(linePointRef(fc, mid));
270  }
271 
272  // Patch -> Cells
273  pointField fcs(patch_.faceCentres());
274  mapPtr_->reverseDistribute(cellIndices_.size(), fcs);
275  forAll(cellIndices_, i)
276  {
277  const label celli = cellIndices_[i];
278  if (celli == -1) continue;
279  const point& cc = nbrMesh.cellCentres()[celli];
280  const point mid = 0.51*cc + 0.49*fcs[i];
281  obj.write(linePointRef(cc, mid));
282  }
283  }
284 }
285 
286 
287 // * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * * * //
288 
290 (
291  const polyPatch& pp
292 )
293 :
294  patch_(pp),
295  nbrRegionName_(patch_.boundaryMesh().mesh().name()),
296  offsetMode_(NORMAL),
297  distance_(NaN),
298  offset_(vector::uniform(NaN)),
299  mapPtr_(nullptr),
300  cellIndices_()
301 {}
302 
303 
305 (
306  const polyPatch& pp,
307  const word& neighbourRegion
308 )
309 :
310  patch_(pp),
311  nbrRegionName_(neighbourRegion),
312  offsetMode_(NORMAL),
313  distance_(NaN),
314  offset_(vector::uniform(NaN)),
315  mapPtr_(nullptr),
316  cellIndices_()
317 {}
318 
319 
321 (
322  const polyPatch& pp,
323  const dictionary& dict
324 )
325 :
326  patch_(pp),
327  nbrRegionName_
328  (
329  dict.lookupOrDefaultBackwardsCompatible<word>
330  (
331  {"neighbourRegion", "sampleRegion"},
332  patch_.boundaryMesh().mesh().name()
333  )
334  ),
335  offsetMode_(readOffsetMode(dict)),
336  distance_
337  (
338  offsetMode_ == NORMAL
339  ? dict.lookup<scalar>("distance")
340  : NaN
341  ),
342  offset_
343  (
344  offsetMode_ == DIRECTION
345  ? dict.lookup<vector>("offset")
346  : vector::uniform(NaN)
347  ),
348  mapPtr_(nullptr),
349  cellIndices_()
350 {}
351 
352 
354 (
355  const polyPatch& pp,
356  const mappedInternalPatchBase& mpb
357 )
358 :
359  patch_(pp),
360  nbrRegionName_(mpb.nbrRegionName_),
361  offsetMode_(mpb.offsetMode_),
362  distance_(mpb.distance_),
363  offset_(mpb.offset_),
364  mapPtr_(nullptr),
365  cellIndices_()
366 {}
367 
368 
369 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
370 
372 {}
373 
374 
375 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
376 
378 {
379  return patch_.boundaryMesh().mesh().time().lookupObject<polyMesh>
380  (
381  nbrRegionName()
382  );
383 }
384 
385 
387 {
388  const polyMesh& mesh = patch_.boundaryMesh().mesh();
389 
390  // Force construction of min-tet decomp
391  (void)mesh.tetBasePtIs();
392 
393  // Allocate the result
394  tmp<pointField> tresult(new pointField(patch_.size()));
395  pointField& result = tresult.ref();
396 
397  // Compute the face points. Assume the face centre to begin with, and then
398  // try and find a better match on complex faces by doing a ray intersection
399  // with a triangulation of the face. This triangulation matches the
400  // tetrahedralisation used for cell searches. This maximises the chances
401  // that the point will be successfully found in a cell.
402  forAll(patch_, patchFacei)
403  {
404  const pointField& ps = mesh.points();
405 
406  const face& f = patch_[patchFacei];
407  const point& fc = patch_.faceCentres()[patchFacei];
408 
409  result[patchFacei] = fc;
410 
411  if (f.size() > 3)
412  {
413  const label facei = patch_.start() + patchFacei;
414  const label celli = patch_.faceCells()[patchFacei];
415 
416  const point& cc = mesh.cellCentres()[celli];
417  const vector d = fc - cc;
418 
419  const label faceBasePtI = mesh.tetBasePtIs()[facei];
420 
421  for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI ++)
422  {
423  const label facePtI = (tetPtI + faceBasePtI) % f.size();
424  const label otherFacePtI = f.fcIndex(facePtI);
425 
426  const triPointRef tri
427  (
428  ps[f[faceBasePtI]],
429  ps[f[facePtI]],
430  ps[f[otherFacePtI]]
431  );
432 
433  const pointHit hitInfo =
435 
436  if (hitInfo.hit() && hitInfo.distance() > 0)
437  {
438  result[patchFacei] = hitInfo.hitPoint();
439  break;
440  }
441  }
442  }
443  }
444 
445  // Apply offset to get sample points
446  switch (offsetMode_)
447  {
448  case NORMAL:
449  result += distance_*patch_.faceNormals();
450  break;
451  case DIRECTION:
452  result += offset_;
453  break;
454  }
455 
456  return tresult;
457 }
458 
459 
461 {
462  mapPtr_.clear();
463  cellIndices_.clear();
464 }
465 
466 
468 {
469  return
470  dict.found("neighbourRegion")
471  || dict.found("offsetMode")
472  || dict.found("distance")
473  || dict.found("offset");
474 }
475 
476 
478 {
480  (
481  os,
482  "neighbourRegion",
483  patch_.boundaryMesh().mesh().name(),
484  nbrRegionName_
485  );
486 
487  writeEntry(os, "offsetMode", offsetModeNames_[offsetMode_]);
488 
489  switch (offsetMode_)
490  {
491  case NORMAL:
492  writeEntry(os, "distance", distance_);
493  break;
494  case DIRECTION:
495  writeEntry(os, "offset", offset_);
496  break;
497  }
498 }
499 
500 
501 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Macros for easy insertion into run-time selection tables.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
virtual const fileName & name() const
Return the name of the stream.
Definition: IOstream.H:294
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:54
Enum read(Istream &) const
Read a word from Istream and return the corresponding.
Definition: NamedEnum.C:61
OFstream which keeps track of vertices.
Definition: OBJstream.H:56
virtual Ostream & write(const char)
Write character.
Definition: OBJstream.C:81
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:57
const Point & hitPoint() const
Return hit point.
Definition: PointHit.H:126
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
bool hit() const
Is there a hit.
Definition: PointHit.H:120
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:54
const Point & hitPoint() const
Return hit point.
label index() const
Return index.
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
A List with indirect addressing.
Definition: UIndirectList.H:60
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
T * data()
Return a pointer to the first data element,.
Definition: UListI.H:149
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
static Form uniform(const Cmpt &s)
Return a VectorSpace with all elements = s.
Definition: VectorSpaceI.H:168
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:710
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:509
Class containing processor-to-processor mapping information.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:76
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:64
label size() const
Global sum of localSizes.
Definition: globalIndexI.H:66
label toGlobal(const label i) const
From local to global.
Definition: globalIndexI.H:82
label findInside(const point &) const
Find shape containing point. Only implemented for certain.
Engine which provides mapping from cells to patch faces.
virtual void write(Ostream &) const
Write as a dictionary.
const polyMesh & nbrMesh() const
Get the region mesh.
tmp< pointField > samplePoints() const
Get the sample points.
virtual ~mappedInternalPatchBase()
Destructor.
static const NamedEnum< offsetMode, 2 > offsetModeNames_
offsetMode
How to project face centres.
static bool specified(const dictionary &dict)
Return whether or not the given dictionary contains a.
void calcMapping() const
Calculate mapping.
mappedInternalPatchBase(const polyPatch &)
Construct from patch.
offsetMode readOffsetMode(const dictionary &dict) const
Read the offset mode from a dictionary.
void clearOut()
Clear out data on mesh change.
const Time & time() const
Return time.
const Type & lookupObject(const word &name) const
Lookup and return the object of the given Type and name.
const polyMesh & mesh() const
Return the mesh reference.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
const labelIOList & tetBasePtIs() const
Return the tetBasePtIs.
Definition: polyMesh.C:1023
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:404
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1313
const indexedOctree< treeDataCell > & cellTree() const
Return the cell search tree.
Definition: polyMesh.C:1057
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:70
const vectorField & cellCentres() const
A class for managing temporary objects.
Definition: tmp.H:55
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
A triangle primitive used to calculate face areas and swept volumes.
Definition: triangle.H:80
pointHit intersection(const point &p, const vector &q, const intersection::algorithm alg, const scalar tol=0.0) const
Fast intersection with a ray.
Definition: triangleI.H:409
A class for handling words, derived from string.
Definition: word.H:62
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
void writeEntryIfDifferent(Ostream &os, const word &entryName, const EntryType &value1, const EntryType &value2)
Helper function to write the keyword and entry only if the.
List< label > labelList
A List of labels.
Definition: labelList.H:56
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
void writeEntry(Ostream &os, const HashTable< T, Key, Hash > &ht)
Definition: HashTableIO.C:96
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
defineTypeNameAndDebug(combustionModel, 0)
error FatalError
dimensioned< scalar > magSqr(const dimensioned< Type > &)
labelList f(nPoints)
dictionary dict
volScalarField & p
Operator to take smallest valid value.
Definition: RemoteData.H:94