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 meshToMesh0& meshToMesh0Interp)
84 {
85  // Determine which particles are in meshTarget
86  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
87 
88  // target to source cell map
89  const labelList& cellAddressing = meshToMesh0Interp.cellAddressing();
90 
91  // Invert celladdressing to get source to target(s).
92  // Note: could use sparse addressing but that is too storage inefficient
93  // (Map<labelList>)
94  labelListList sourceToTargets
95  (
96  invertOneToMany(meshToMesh0Interp.fromMesh().nCells(), cellAddressing)
97  );
98 
99  const fvMesh& meshSource = meshToMesh0Interp.fromMesh();
100  const fvMesh& meshTarget = meshToMesh0Interp.toMesh();
101 
102 
103  fileNameList cloudDirs
104  (
105  readDir
106  (
107  meshSource.time().timePath()/cloud::prefix,
109  )
110  );
111 
112  forAll(cloudDirs, cloudI)
113  {
114  // Search for list of lagrangian objects for this time
115  IOobjectList objects
116  (
117  meshSource,
118  meshSource.time().timeName(),
119  cloud::prefix/cloudDirs[cloudI]
120  );
121 
122  IOobject* positionsPtr = objects.lookup("positions");
123 
124  if (positionsPtr)
125  {
126  Info<< nl << " processing cloud " << cloudDirs[cloudI] << endl;
127 
128  // Read positions & cell
129  passiveParticleCloud sourceParcels
130  (
131  meshSource,
132  cloudDirs[cloudI],
133  false
134  );
135  Info<< " read " << sourceParcels.size()
136  << " parcels from source mesh." << endl;
137 
138  // Construct empty target cloud
139  passiveParticleCloud targetParcels
140  (
141  meshTarget,
142  cloudDirs[cloudI],
143  IDLList<passiveParticle>()
144  );
145 
146  passiveParticle::trackingData td(targetParcels);
147 
148  label sourceParticleI = 0;
149 
150  // Indices of source particles that get added to targetParcels
151  DynamicList<label> addParticles(sourceParcels.size());
152 
153  // Unmapped particles
154  labelHashSet unmappedSource(sourceParcels.size());
155 
156 
157  // Initial: track from fine-mesh cell centre to particle position
158  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
159  // This requires there to be no boundary in the way.
160 
161 
162  forAllConstIter(Cloud<passiveParticle>, sourceParcels, iter)
163  {
164  bool foundCell = false;
165 
166  // Assume that cell from read parcel is the correct one...
167  if (iter().cell() >= 0)
168  {
169  const labelList& targetCells =
170  sourceToTargets[iter().cell()];
171 
172  // Particle probably in one of the targetcells. Try
173  // all by tracking from their cell centre to the parcel
174  // position.
175 
176  forAll(targetCells, i)
177  {
178  // Track from its cellcentre to position to make sure.
179  autoPtr<passiveParticle> newPtr
180  (
181  new passiveParticle
182  (
183  meshTarget,
184  barycentric(1, 0, 0, 0),
185  targetCells[i],
186  meshTarget.cells()[targetCells[i]][0],
187  1
188  )
189  );
190  passiveParticle& newP = newPtr();
191 
192  newP.track(iter().position() - newP.position(), 0);
193 
194  if (!newP.onFace())
195  {
196  // Hit position.
197  foundCell = true;
198  addParticles.append(sourceParticleI);
199  targetParcels.addParticle(newPtr.ptr());
200  break;
201  }
202  }
203  }
204 
205  if (!foundCell)
206  {
207  // Store for closer analysis
208  unmappedSource.insert(sourceParticleI);
209  }
210 
211  sourceParticleI++;
212  }
213 
214  Info<< " after meshToMesh0 addressing found "
215  << targetParcels.size()
216  << " parcels in target mesh." << endl;
217 
218 
219  // Do closer inspection for unmapped particles
220  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
221 
222  if (unmappedSource.size())
223  {
224  sourceParticleI = 0;
225 
226  forAllIter(Cloud<passiveParticle>, sourceParcels, iter)
227  {
228  if (unmappedSource.found(sourceParticleI))
229  {
230  label targetCell =
231  findCell(targetParcels, iter().position());
232 
233  if (targetCell >= 0)
234  {
235  unmappedSource.erase(sourceParticleI);
236  addParticles.append(sourceParticleI);
237  targetParcels.addParticle
238  (
239  new passiveParticle
240  (
241  meshTarget,
242  iter().position(),
243  targetCell
244  )
245  );
246  sourceParcels.remove(&iter());
247  }
248  }
249  sourceParticleI++;
250  }
251  }
252  addParticles.shrink();
253 
254  Info<< " after additional mesh searching found "
255  << targetParcels.size() << " parcels in target mesh." << endl;
256 
257  if (addParticles.size())
258  {
259  IOPosition<passiveParticleCloud>(targetParcels).write();
260 
261  // addParticles now contains the indices of the sourceMesh
262  // particles that were appended to the target mesh.
263 
264  // Map lagrangian fields
265  // ~~~~~~~~~~~~~~~~~~~~~
266 
267  MapLagrangianFields<label>
268  (cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
269  MapLagrangianFields<scalar>
270  (cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
271  MapLagrangianFields<vector>
272  (cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
273  MapLagrangianFields<sphericalTensor>
274  (cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
275  MapLagrangianFields<symmTensor>
276  (cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
277  MapLagrangianFields<tensor>
278  (cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
279  }
280  }
281  }
282 }
283 
284 
285 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
286 
287 } // End namespace Foam
288 
289 // ************************************************************************* //
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
labelListList invertOneToMany(const label len, const labelUList &)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:67
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.