All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mapLagrangian.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-2018 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 
26 #include "MapLagrangianFields.H"
27 #include "passiveParticleCloud.H"
28 #include "meshSearch.H"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 static const scalar perturbFactor = 1e-6;
36 
37 
38 // Special version of findCell that generates a cell guaranteed to be
39 // compatible with tracking.
40 static label findCell(const Cloud<passiveParticle>& cloud, const point& pt)
41 {
42  label celli = -1;
43  label tetFacei = -1;
44  label tetPtI = -1;
45 
46  const polyMesh& mesh = cloud.pMesh();
47 
48  mesh.findCellFacePt(pt, celli, tetFacei, tetPtI);
49 
50  if (celli >= 0)
51  {
52  return celli;
53  }
54  else
55  {
56  // See if particle on face by finding nearest face and shifting
57  // particle.
58 
59  meshSearch meshSearcher
60  (
61  mesh,
62  polyMesh::FACE_PLANES // no decomposition needed
63  );
64 
65  label facei = meshSearcher.findNearestBoundaryFace(pt);
66 
67  if (facei >= 0)
68  {
69  const point& cc = mesh.cellCentres()[mesh.faceOwner()[facei]];
70 
71  const point perturbPt = (1-perturbFactor)*pt+perturbFactor*cc;
72 
73  mesh.findCellFacePt(perturbPt, celli, tetFacei, tetPtI);
74 
75  return celli;
76  }
77  }
78 
79  return -1;
80 }
81 
82 
83 void mapLagrangian(const meshToMesh& interp)
84 {
85  // Determine which particles are in meshTarget
86  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
87 
88  const polyMesh& meshSource = interp.srcRegion();
89  const polyMesh& meshTarget = interp.tgtRegion();
90  const labelListList& sourceToTarget = interp.srcToTgtCellAddr();
91 
92  fileNameList cloudDirs
93  (
94  readDir
95  (
96  meshSource.time().timePath()/cloud::prefix,
98  )
99  );
100 
101  forAll(cloudDirs, cloudI)
102  {
103  // Search for list of lagrangian objects for this time
104  IOobjectList objects
105  (
106  meshSource,
107  meshSource.time().timeName(),
108  cloud::prefix/cloudDirs[cloudI]
109  );
110 
111  IOobject* positionsPtr = objects.lookup(word("positions"));
112 
113  if (positionsPtr)
114  {
115  Info<< nl << " processing cloud " << cloudDirs[cloudI] << endl;
116 
117  // Read positions & cell
118  passiveParticleCloud sourceParcels
119  (
120  meshSource,
121  cloudDirs[cloudI],
122  false
123  );
124  Info<< " read " << sourceParcels.size()
125  << " parcels from source mesh." << endl;
126 
127  // Construct empty target cloud
128  passiveParticleCloud targetParcels
129  (
130  meshTarget,
131  cloudDirs[cloudI],
132  IDLList<passiveParticle>()
133  );
134 
135  passiveParticle::trackingData td(targetParcels);
136 
137  label sourceParticleI = 0;
138 
139  // Indices of source particles that get added to targetParcels
140  DynamicList<label> addParticles(sourceParcels.size());
141 
142  // Unmapped particles
143  labelHashSet unmappedSource(sourceParcels.size());
144 
145 
146  // Initial: track from fine-mesh cell centre to particle position
147  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
148  // This requires there to be no boundary in the way.
149 
150 
151  forAllConstIter(Cloud<passiveParticle>, sourceParcels, iter)
152  {
153  bool foundCell = false;
154 
155  // Assume that cell from read parcel is the correct one...
156  if (iter().cell() >= 0)
157  {
158  const labelList& targetCells =
159  sourceToTarget[iter().cell()];
160 
161  // Particle probably in one of the targetcells. Try
162  // all by tracking from their cell centre to the parcel
163  // position.
164 
165  forAll(targetCells, i)
166  {
167  // Track from its cellcentre to position to make sure.
168  autoPtr<passiveParticle> newPtr
169  (
170  new passiveParticle
171  (
172  meshTarget,
173  barycentric(1, 0, 0, 0),
174  targetCells[i],
175  meshTarget.cells()[targetCells[i]][0],
176  1
177  )
178  );
179  passiveParticle& newP = newPtr();
180 
181  newP.track(iter().position() - newP.position(), 0);
182 
183  if (!newP.onFace())
184  {
185  // Hit position.
186  foundCell = true;
187  addParticles.append(sourceParticleI);
188  targetParcels.addParticle(newPtr.ptr());
189  break;
190  }
191  }
192  }
193 
194  if (!foundCell)
195  {
196  // Store for closer analysis
197  unmappedSource.insert(sourceParticleI);
198  }
199 
200  sourceParticleI++;
201  }
202 
203  Info<< " after meshToMesh addressing found "
204  << targetParcels.size()
205  << " parcels in target mesh." << endl;
206 
207 
208  // Do closer inspection for unmapped particles
209  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
210 
211  if (unmappedSource.size())
212  {
213  sourceParticleI = 0;
214 
215  forAllIter(Cloud<passiveParticle>, sourceParcels, iter)
216  {
217  if (unmappedSource.found(sourceParticleI))
218  {
219  label targetCell =
220  findCell(targetParcels, iter().position());
221 
222  if (targetCell >= 0)
223  {
224  unmappedSource.erase(sourceParticleI);
225  addParticles.append(sourceParticleI);
226  targetParcels.addParticle
227  (
228  new passiveParticle
229  (
230  meshTarget,
231  iter().position(),
232  targetCell
233  )
234  );
235  sourceParcels.remove(&iter());
236  }
237  }
238  sourceParticleI++;
239  }
240  }
241  addParticles.shrink();
242 
243  Info<< " after additional mesh searching found "
244  << targetParcels.size() << " parcels in target mesh." << endl;
245 
246  if (addParticles.size())
247  {
248  IOPosition<passiveParticleCloud>(targetParcels).write();
249 
250  // addParticles now contains the indices of the sourceMesh
251  // particles that were appended to the target mesh.
252 
253  // Map lagrangian fields
254  // ~~~~~~~~~~~~~~~~~~~~~
255 
256  MapLagrangianFields<label>
257  (
258  cloudDirs[cloudI],
259  objects,
260  meshTarget,
261  addParticles
262  );
263  MapLagrangianFields<scalar>
264  (
265  cloudDirs[cloudI],
266  objects,
267  meshTarget,
268  addParticles
269  );
270  MapLagrangianFields<vector>
271  (
272  cloudDirs[cloudI],
273  objects,
274  meshTarget,
275  addParticles
276  );
277  MapLagrangianFields<sphericalTensor>
278  (
279  cloudDirs[cloudI],
280  objects,
281  meshTarget,
282  addParticles
283  );
284  MapLagrangianFields<symmTensor>
285  (
286  cloudDirs[cloudI],
287  objects,
288  meshTarget,
289  addParticles
290  );
291  MapLagrangianFields<tensor>
292  (
293  cloudDirs[cloudI],
294  objects,
295  meshTarget,
296  addParticles
297  );
298  }
299  }
300  }
301 }
302 
303 
304 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
305 
306 } // End namespace Foam
307 
308 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:459
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Barycentric< scalar > barycentric
A scalar version of the templated Barycentric.
Definition: barycentric.H:45
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:211
void write(std::ostream &os, const bool binary, List< floatScalar > &fField)
Write floats ascii or binary.
List< label > labelList
A List of labels.
Definition: labelList.H:56
fileNameList readDir(const fileName &, const fileType=fileType::file, const bool filterVariants=true, const bool followLink=true)
Read a directory and return the entries as a string list.
Definition: POSIX.C:662
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
static const char nl
Definition: Ostream.H:260
objects
static const word prefix
The prefix to local: lagrangian.
Definition: cloud.H:62
void mapLagrangian(const meshToMesh0 &meshToMesh0Interp)
Maps lagrangian positions and fields.
vector point
Point is a vector.
Definition: point.H:41
messageStream Info
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
List< fileName > fileNameList
A List of fileNames.
Definition: fileNameList.H:50
Namespace for OpenFOAM.