FreeFOAM The Cross-Platform CFD Toolkit
directMappedPatchBase.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
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 
26 #include "directMappedPatchBase.H"
28 #include <OpenFOAM/ListListOps.H>
29 #include <meshTools/meshSearch.H>
30 #include <meshTools/meshTools.H>
31 #include <OpenFOAM/OFstream.H>
32 #include <OpenFOAM/Random.H>
33 #include <meshTools/treeDataFace.H>
35 #include <OpenFOAM/mapDistribute.H>
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41  defineTypeNameAndDebug(directMappedPatchBase, 0);
42 
43  template<>
45  {
46  "nearestCell",
47  "nearestPatchFace",
48  "nearestFace"
49  };
50 
52  directMappedPatchBase::sampleModeNames_;
53 }
54 
55 
56 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
57 
58 void Foam::directMappedPatchBase::collectSamples
59 (
61  labelList& patchFaceProcs,
62  labelList& patchFaces,
63  pointField& patchFc
64 ) const
65 {
66 
67  // Collect all sample points and the faces they come from.
68  List<pointField> globalFc(Pstream::nProcs());
69  List<pointField> globalSamples(Pstream::nProcs());
70  labelListList globalFaces(Pstream::nProcs());
71 
72  globalFc[Pstream::myProcNo()] = patch_.faceCentres();
73  globalSamples[Pstream::myProcNo()] = globalFc[Pstream::myProcNo()]+offsets_;
74  globalFaces[Pstream::myProcNo()] = identity(patch_.size());
75 
76  // Distribute to all processors
77  Pstream::gatherList(globalSamples);
78  Pstream::scatterList(globalSamples);
79  Pstream::gatherList(globalFaces);
80  Pstream::scatterList(globalFaces);
81  Pstream::gatherList(globalFc);
82  Pstream::scatterList(globalFc);
83 
84  // Rework into straight list
85  samples = ListListOps::combine<pointField>
86  (
87  globalSamples,
88  accessOp<pointField>()
89  );
90  patchFaces = ListListOps::combine<labelList>
91  (
92  globalFaces,
93  accessOp<labelList>()
94  );
95  patchFc = ListListOps::combine<pointField>
96  (
97  globalFc,
98  accessOp<pointField>()
99  );
100 
101  patchFaceProcs.setSize(patchFaces.size());
102  labelList nPerProc
103  (
105  (
106  globalFaces,
107  accessOp<labelList>()
108  )
109  );
110  label sampleI = 0;
111  forAll(nPerProc, procI)
112  {
113  for (label i = 0; i < nPerProc[procI]; i++)
114  {
115  patchFaceProcs[sampleI++] = procI;
116  }
117  }
118 }
119 
120 
121 // Find the processor/cell containing the samples. Does not account
122 // for samples being found in two processors.
123 void Foam::directMappedPatchBase::findSamples
124 (
125  const pointField& samples,
126  labelList& sampleProcs,
127  labelList& sampleIndices,
128  pointField& sampleLocations
129 ) const
130 {
131  // Lookup the correct region
132  const polyMesh& mesh = sampleMesh();
133 
134  // All the info for nearest. Construct to miss
135  List<nearInfo> nearest(samples.size());
136 
137  switch (mode_)
138  {
139  case NEARESTCELL:
140  {
141  if (samplePatch_.size() && samplePatch_ != "none")
142  {
144  (
145  "directMappedPatchBase::findSamples(const pointField&,"
146  " labelList&, labelList&, pointField&) const"
147  ) << "No need to supply a patch name when in "
148  << sampleModeNames_[mode_] << " mode." << exit(FatalError);
149  }
150 
151  // Octree based search engine
152  meshSearch meshSearchEngine(mesh, false);
153 
154  forAll(samples, sampleI)
155  {
156  const point& sample = samples[sampleI];
157 
158  label cellI = meshSearchEngine.findCell(sample);
159 
160  if (cellI == -1)
161  {
162  nearest[sampleI].second().first() = Foam::sqr(GREAT);
163  nearest[sampleI].second().second() = Pstream::myProcNo();
164  }
165  else
166  {
167  const point& cc = mesh.cellCentres()[cellI];
168 
169  nearest[sampleI].first() = pointIndexHit
170  (
171  true,
172  cc,
173  cellI
174  );
175  nearest[sampleI].second().first() = magSqr(cc-sample);
176  nearest[sampleI].second().second() = Pstream::myProcNo();
177  }
178  }
179  break;
180  }
181 
182  case NEARESTPATCHFACE:
183  {
184  Random rndGen(123456);
185 
186  const polyPatch& pp = samplePolyPatch();
187 
188  if (pp.empty())
189  {
190  forAll(samples, sampleI)
191  {
192  nearest[sampleI].second().first() = Foam::sqr(GREAT);
193  nearest[sampleI].second().second() = Pstream::myProcNo();
194  }
195  }
196  else
197  {
198  // patch faces
199  const labelList patchFaces(identity(pp.size()) + pp.start());
200 
201  treeBoundBox patchBb
202  (
203  treeBoundBox(pp.points(), pp.meshPoints()).extend
204  (
205  rndGen,
206  1E-4
207  )
208  );
209  patchBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
210  patchBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
211 
212  indexedOctree<treeDataFace> boundaryTree
213  (
214  treeDataFace // all information needed to search faces
215  (
216  false, // do not cache bb
217  mesh,
218  patchFaces // boundary faces only
219  ),
220  patchBb, // overall search domain
221  8, // maxLevel
222  10, // leafsize
223  3.0 // duplicity
224  );
225 
226  forAll(samples, sampleI)
227  {
228  const point& sample = samples[sampleI];
229 
230  pointIndexHit& nearInfo = nearest[sampleI].first();
231  nearInfo = boundaryTree.findNearest
232  (
233  sample,
234  magSqr(patchBb.span())
235  );
236 
237  if (!nearInfo.hit())
238  {
239  nearest[sampleI].second().first() = Foam::sqr(GREAT);
240  nearest[sampleI].second().second() =
242  }
243  else
244  {
245  point fc(pp[nearInfo.index()].centre(pp.points()));
246  nearInfo.setPoint(fc);
247  nearest[sampleI].second().first() = magSqr(fc-sample);
248  nearest[sampleI].second().second() =
250  }
251  }
252  }
253  break;
254  }
255 
256  case NEARESTFACE:
257  {
258  if (samplePatch_.size() && samplePatch_ != "none")
259  {
261  (
262  "directMappedPatchBase::findSamples(const pointField&,"
263  " labelList&, labelList&, pointField&) const"
264  ) << "No need to supply a patch name when in "
265  << sampleModeNames_[mode_] << " mode." << exit(FatalError);
266  }
267 
268  // Octree based search engine
269  meshSearch meshSearchEngine(mesh, false);
270 
271  forAll(samples, sampleI)
272  {
273  const point& sample = samples[sampleI];
274 
275  label faceI = meshSearchEngine.findNearestFace(sample);
276 
277  if (faceI == -1)
278  {
279  nearest[sampleI].second().first() = Foam::sqr(GREAT);
280  nearest[sampleI].second().second() = Pstream::myProcNo();
281  }
282  else
283  {
284  const point& fc = mesh.faceCentres()[faceI];
285 
286  nearest[sampleI].first() = pointIndexHit
287  (
288  true,
289  fc,
290  faceI
291  );
292  nearest[sampleI].second().first() = magSqr(fc-sample);
293  nearest[sampleI].second().second() = Pstream::myProcNo();
294  }
295  }
296  break;
297  }
298 
299  default:
300  {
301  FatalErrorIn("directMappedPatchBase::findSamples(..)")
302  << "problem." << abort(FatalError);
303  }
304  }
305 
306 
307  // Find nearest.
308  Pstream::listCombineGather(nearest, nearestEqOp());
310 
311  if (debug)
312  {
313  Info<< "directMappedPatchBase::findSamples on mesh " << sampleRegion_
314  << " : " << endl;
315  forAll(nearest, sampleI)
316  {
317  label procI = nearest[sampleI].second().second();
318  label localI = nearest[sampleI].first().index();
319 
320  Info<< " " << sampleI << " coord:"<< samples[sampleI]
321  << " found on processor:" << procI
322  << " in local cell/face:" << localI
323  << " with cc:" << nearest[sampleI].first().rawPoint() << endl;
324  }
325  }
326 
327  // Check for samples not being found
328  forAll(nearest, sampleI)
329  {
330  if (!nearest[sampleI].first().hit())
331  {
333  (
334  "directMappedPatchBase::findSamples"
335  "(const pointField&, labelList&"
336  ", labelList&, pointField&)"
337  ) << "Did not find sample " << samples[sampleI]
338  << " on any processor of region " << sampleRegion_
339  << exit(FatalError);
340  }
341  }
342 
343 
344  // Convert back into proc+local index
345  sampleProcs.setSize(samples.size());
346  sampleIndices.setSize(samples.size());
347  sampleLocations.setSize(samples.size());
348 
349  forAll(nearest, sampleI)
350  {
351  sampleProcs[sampleI] = nearest[sampleI].second().second();
352  sampleIndices[sampleI] = nearest[sampleI].first().index();
353  sampleLocations[sampleI] = nearest[sampleI].first().hitPoint();
354  }
355 }
356 
357 
358 void Foam::directMappedPatchBase::calcMapping() const
359 {
360  if (mapPtr_.valid())
361  {
362  FatalErrorIn("directMappedPatchBase::calcMapping() const")
363  << "Mapping already calculated" << exit(FatalError);
364  }
365 
366  if
367  (
368  gAverage(mag(offsets_)) <= ROOTVSMALL
369  && mode_ == NEARESTPATCHFACE
370  && sampleRegion_ == patch_.boundaryMesh().mesh().name()
371  && samplePatch_ == patch_.name()
372  )
373  {
374  WarningIn("directMappedPatchBase::calcMapping() const")
375  << "Invalid offset " << offsets_ << endl
376  << "Offset is the vector added to the patch face centres to"
377  << " find the patch face supplying the data." << endl
378  << "Setting it to " << offsets_
379  << " on the same patch, on the same region"
380  << " will find the faces themselves which does not make sense"
381  << " for anything but testing." << endl
382  << "patch_:" << patch_.name() << endl
383  << "sampleRegion_:" << sampleRegion_ << endl
384  << "mode_:" << sampleModeNames_[mode_] << endl
385  << "samplePatch_:" << samplePatch_ << endl
386  << "offsets_:" << offsets_ << endl;
387  }
388 
389 
390  // Get global list of all samples and the processor and face they come from.
392  labelList patchFaceProcs;
393  labelList patchFaces;
394  pointField patchFc;
395  collectSamples(samples, patchFaceProcs, patchFaces, patchFc);
396 
397  // Find processor and cell/face samples are in and actual location.
398  labelList sampleProcs;
399  labelList sampleIndices;
400  pointField sampleLocations;
401  findSamples(samples, sampleProcs, sampleIndices, sampleLocations);
402 
403 
404  // Now we have all the data we need:
405  // - where sample originates from (so destination when mapping):
406  // patchFaces, patchFaceProcs.
407  // - cell/face sample is in (so source when mapping)
408  // sampleIndices, sampleProcs.
409 
410  //forAll(samples, i)
411  //{
412  // Info<< i << " need data in region "
413  // << patch_.boundaryMesh().mesh().name()
414  // << " for proc:" << patchFaceProcs[i]
415  // << " face:" << patchFaces[i]
416  // << " at:" << patchFc[i] << endl
417  // << "Found data in region " << sampleRegion_
418  // << " at proc:" << sampleProcs[i]
419  // << " face:" << sampleIndices[i]
420  // << " at:" << sampleLocations[i]
421  // << nl << endl;
422  //}
423 
424 
425 
426  if (debug && Pstream::master())
427  {
428  OFstream str
429  (
430  patch_.boundaryMesh().mesh().time().path()
431  / patch_.name()
432  + "_directMapped.obj"
433  );
434  Pout<< "Dumping mapping as lines from patch faceCentres to"
435  << " sampled cell/faceCentres to file " << str.name() << endl;
436 
437  label vertI = 0;
438 
439  forAll(patchFc, i)
440  {
441  meshTools::writeOBJ(str, patchFc[i]);
442  vertI++;
443  meshTools::writeOBJ(str, sampleLocations[i]);
444  vertI++;
445  str << "l " << vertI-1 << ' ' << vertI << nl;
446  }
447  }
448 
449 
452  //if (Pstream::master())
453  //{
454  // const scalarField magOffset(mag(sampleLocations - patchFc));
455  // const scalar avgOffset(average(magOffset));
456  //
457  // forAll(magOffset, sampleI)
458  // {
459  // if
460  // (
461  // mag(magOffset[sampleI]-avgOffset)
462  // > max(SMALL, 0.001*avgOffset)
463  // )
464  // {
465  // WarningIn("directMappedPatchBase::calcMapping() const")
466  // << "The actual cell/face centres picked up using offset "
467  // << offsets_ << " are not" << endl
468  // << " on a single plane."
469  // << " This might give numerical problems." << endl
470  // << " At patchface " << patchFc[sampleI]
471  // << " the sampled cell/face " << sampleLocations[sampleI]
472  // << endl
473  // << " is not on a plane " << avgOffset
474  // << " offset from the patch." << endl
475  // << " You might want to shift your plane offset."
476  // << " Set the debug flag to get a dump of sampled cells."
477  // << endl;
478  // break;
479  // }
480  // }
481  //}
482 
483 
484  // Determine schedule.
485  mapPtr_.reset(new mapDistribute(sampleProcs, patchFaceProcs));
486 
487  // Rework the schedule from indices into samples to cell data to send,
488  // face data to receive.
489 
490  labelListList& subMap = mapPtr_().subMap();
491  labelListList& constructMap = mapPtr_().constructMap();
492 
493  forAll(subMap, procI)
494  {
495  subMap[procI] = UIndirectList<label>
496  (
497  sampleIndices,
498  subMap[procI]
499  );
500  constructMap[procI] = UIndirectList<label>
501  (
502  patchFaces,
503  constructMap[procI]
504  );
505 
506  if (debug)
507  {
508  Pout<< "To proc:" << procI << " sending values of cells/faces:"
509  << subMap[procI] << endl;
510  Pout<< "From proc:" << procI << " receiving values of patch faces:"
511  << constructMap[procI] << endl;
512  }
513  }
514 
515  // Redo constructSize
516  mapPtr_().constructSize() = patch_.size();
517 
518  if (debug)
519  {
520  // Check that all elements get a value.
521  PackedBoolList used(patch_.size());
522  forAll(constructMap, procI)
523  {
524  const labelList& map = constructMap[procI];
525 
526  forAll(map, i)
527  {
528  label faceI = map[i];
529 
530  if (used[faceI] == 0)
531  {
532  used[faceI] = 1;
533  }
534  else
535  {
536  FatalErrorIn("directMappedPatchBase::calcMapping() const")
537  << "On patch " << patch_.name()
538  << " patchface " << faceI
539  << " is assigned to more than once."
540  << abort(FatalError);
541  }
542  }
543  }
544  forAll(used, faceI)
545  {
546  if (used[faceI] == 0)
547  {
548  FatalErrorIn("directMappedPatchBase::calcMapping() const")
549  << "On patch " << patch_.name()
550  << " patchface " << faceI
551  << " is never assigned to."
552  << abort(FatalError);
553  }
554  }
555  }
556 }
557 
558 
559 // * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * * * //
560 
562 (
563  const polyPatch& pp
564 )
565 :
566  patch_(pp),
567  sampleRegion_(patch_.boundaryMesh().mesh().name()),
568  mode_(NEARESTPATCHFACE),
569  samplePatch_("none"),
570  uniformOffset_(true),
571  offset_(vector::zero),
572  offsets_(pp.size(), offset_),
573  sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
574  mapPtr_(NULL)
575 {}
576 
577 
579 (
580  const polyPatch& pp,
581  const word& sampleRegion,
582  const sampleMode mode,
583  const word& samplePatch,
584  const vectorField& offsets
585 )
586 :
587  patch_(pp),
588  sampleRegion_(sampleRegion),
589  mode_(mode),
590  samplePatch_(samplePatch),
591  uniformOffset_(false),
592  offsets_(offsets),
593  sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
594  mapPtr_(NULL)
595 {}
596 
597 
599 (
600  const polyPatch& pp,
601  const word& sampleRegion,
602  const sampleMode mode,
603  const word& samplePatch,
604  const vector& offset
605 )
606 :
607  patch_(pp),
608  sampleRegion_(sampleRegion),
609  mode_(mode),
610  samplePatch_(samplePatch),
611  uniformOffset_(true),
612  offset_(offset),
613  offsets_(pp.size(), offset_),
614  sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
615  mapPtr_(NULL)
616 {}
617 
618 
620 (
621  const polyPatch& pp,
622  const dictionary& dict
623 )
624 :
625  patch_(pp),
626  sampleRegion_
627  (
628  dict.lookupOrDefault
629  (
630  "sampleRegion",
631  patch_.boundaryMesh().mesh().name()
632  )
633  ),
634  mode_(sampleModeNames_.read(dict.lookup("sampleMode"))),
635  samplePatch_(dict.lookup("samplePatch")),
636  uniformOffset_(dict.found("offset")),
637  offset_
638  (
639  uniformOffset_
640  ? point(dict.lookup("offset"))
641  : vector::zero
642  ),
643  offsets_
644  (
645  uniformOffset_
646  ? pointField(pp.size(), offset_)
647  : dict.lookup("offsets")
648  ),
649  sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
650  mapPtr_(NULL)
651 {}
652 
653 
655 (
656  const polyPatch& pp,
657  const directMappedPatchBase& dmp
658 )
659 :
660  patch_(pp),
661  sampleRegion_(dmp.sampleRegion_),
662  mode_(dmp.mode_),
663  samplePatch_(dmp.samplePatch_),
664  uniformOffset_(dmp.uniformOffset_),
665  offset_(dmp.offset_),
666  offsets_(dmp.offsets_),
667  sameRegion_(dmp.sameRegion_),
668  mapPtr_(NULL)
669 {}
670 
671 
672 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
673 
675 {
676  clearOut();
677 }
678 
679 
681 {
682  mapPtr_.clear();
683 }
684 
685 
686 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
687 
689 {
690  return patch_.boundaryMesh().mesh().time().lookupObject<polyMesh>
691  (
692  sampleRegion_
693  );
694 }
695 
696 
698 {
699  const polyMesh& nbrMesh = sampleMesh();
700 
701  label patchI = nbrMesh.boundaryMesh().findPatchID(samplePatch_);
702 
703  if (patchI == -1)
704  {
705  FatalErrorIn("directMappedPatchBase::samplePolyPatch() ")
706  << "Cannot find patch " << samplePatch_
707  << " in region " << sampleRegion_ << endl
708  << "Valid patches are " << nbrMesh.boundaryMesh().names()
709  << exit(FatalError);
710  }
711 
712  return nbrMesh.boundaryMesh()[patchI];
713 }
714 
715 
717 {
718  os.writeKeyword("sampleMode") << sampleModeNames_[mode_]
719  << token::END_STATEMENT << nl;
720  os.writeKeyword("sampleRegion") << sampleRegion_
721  << token::END_STATEMENT << nl;
722  os.writeKeyword("samplePatch") << samplePatch_
723  << token::END_STATEMENT << nl;
724  if (uniformOffset_)
725  {
726  os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl;
727  }
728  else
729  {
730  os.writeKeyword("offsets") << offsets_ << token::END_STATEMENT << nl;
731  }
732 }
733 
734 
735 // ************************ vim: set sw=4 sts=4 et: ************************ //