FreeFOAM The Cross-Platform CFD Toolkit
mapLagrangian.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 "MapLagrangianFields.H"
27 #include <lagrangian/Cloud.H>
29 #include <meshTools/meshSearch.H>
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 static const scalar perturbFactor = 1E-6;
37 
38 
39 // Special version of findCell that generates a cell guaranteed to be
40 // compatible with tracking.
41 static label findCell(const meshSearch& meshSearcher, const point& pt)
42 {
43  const polyMesh& mesh = meshSearcher.mesh();
44 
45  // Use tracking to find cell containing pt
46  label cellI = meshSearcher.findCell(pt);
47 
48  if (cellI >= 0)
49  {
50  return cellI;
51  }
52  else
53  {
54  // See if particle on face by finding nearest face and shifting
55  // particle.
56 
57  label faceI = meshSearcher.findNearestBoundaryFace(pt);
58 
59  if (faceI >= 0)
60  {
61  const point& cc = mesh.cellCentres()[mesh.faceOwner()[faceI]];
62  const point perturbPt = (1-perturbFactor)*pt+perturbFactor*cc;
63 
64  return meshSearcher.findCell(perturbPt);
65  }
66  }
67  return -1;
68 }
69 
70 
71 void mapLagrangian(const meshToMesh& meshToMeshInterp)
72 {
73  // Determine which particles are in meshTarget
74  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
75 
76  // target to source cell map
77  const labelList& cellAddressing = meshToMeshInterp.cellAddressing();
78 
79  // Invert celladdressing to get source to target(s).
80  // Note: could use sparse addressing but that is too storage inefficient
81  // (Map<labelList>)
82  labelListList sourceToTargets
83  (
84  invertOneToMany(meshToMeshInterp.fromMesh().nCells(), cellAddressing)
85  );
86 
87  const fvMesh& meshSource = meshToMeshInterp.fromMesh();
88  const fvMesh& meshTarget = meshToMeshInterp.toMesh();
89  const pointField& targetCc = meshTarget.cellCentres();
90 
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("positions");
112 
113  if (positionsPtr)
114  {
115  Info<< nl << " processing cloud " << cloudDirs[cloudI] << endl;
116 
117  // Read positions & cell
118  Cloud<passiveParticle> 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  Cloud<passiveParticle> targetParcels
129  (
130  meshTarget,
131  cloudDirs[cloudI],
132  IDLList<passiveParticle>()
133  );
134 
135  label sourceParticleI = 0;
136 
137  // Indices of source particles that get added to targetParcels
138  DynamicList<label> addParticles(sourceParcels.size());
139 
140  // Unmapped particles
141  labelHashSet unmappedSource(sourceParcels.size());
142 
143 
144  // Initial: track from fine-mesh cell centre to particle position
145  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
146  // This requires there to be no boundary in the way.
147 
148 
149  forAllConstIter(Cloud<passiveParticle>, sourceParcels, iter)
150  {
151  bool foundCell = false;
152 
153  // Assume that cell from read parcel is the correct one...
154  if (iter().cell() >= 0)
155  {
156  const labelList& targetCells =
157  sourceToTargets[iter().cell()];
158 
159  // Particle probably in one of the targetcells. Try
160  // all by tracking from their cell centre to the parcel
161  // position.
162 
163  forAll(targetCells, i)
164  {
165  // Track from its cellcentre to position to make sure.
166  autoPtr<passiveParticle> newPtr
167  (
168  new passiveParticle
169  (
170  targetParcels,
171  targetCc[targetCells[i]],
172  targetCells[i]
173  )
174  );
175  passiveParticle& newP = newPtr();
176 
177  scalar fraction = 0;
178  label faceI = newP.track(iter().position(), fraction);
179 
180  if (faceI < 0 && newP.cell() >= 0)
181  {
182  // Hit position.
183  foundCell = true;
184  addParticles.append(sourceParticleI);
185  targetParcels.addParticle(newPtr.ptr());
186  break;
187  }
188  }
189  }
190 
191  if (!foundCell)
192  {
193  // Store for closer analysis
194  unmappedSource.insert(sourceParticleI);
195  }
196 
197  sourceParticleI++;
198  }
199 
200  Info<< " after meshToMesh addressing found "
201  << targetParcels.size()
202  << " parcels in target mesh." << endl;
203 
204 
205  // Do closer inspection for unmapped particles
206  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
207 
208  if (unmappedSource.size())
209  {
210  meshSearch targetSearcher(meshTarget, false);
211 
212  sourceParticleI = 0;
213 
214  forAllIter(Cloud<passiveParticle>, sourceParcels, iter)
215  {
216  if (unmappedSource.found(sourceParticleI))
217  {
218  label targetCell =
219  findCell(targetSearcher, iter().position());
220 
221  if (targetCell >= 0)
222  {
223  unmappedSource.erase(sourceParticleI);
224  addParticles.append(sourceParticleI);
225  iter().cell()=targetCell;
226  targetParcels.addParticle
227  (
228  sourceParcels.remove(&iter())
229  );
230  }
231  }
232  sourceParticleI++;
233  }
234  }
235  addParticles.shrink();
236 
237  Info<< " after additional mesh searching found "
238  << targetParcels.size() << " parcels in target mesh." << endl;
239 
240  if (addParticles.size())
241  {
242  IOPosition<passiveParticle>(targetParcels).write();
243 
244  // addParticles now contains the indices of the sourceMesh
245  // particles that were appended to the target mesh.
246 
247  // Map lagrangian fields
248  // ~~~~~~~~~~~~~~~~~~~~~
249 
250  MapLagrangianFields<label>
251  (cloudDirs[cloudI], objects, meshToMeshInterp, addParticles);
252  MapLagrangianFields<scalar>
253  (cloudDirs[cloudI], objects, meshToMeshInterp, addParticles);
254  MapLagrangianFields<vector>
255  (cloudDirs[cloudI], objects, meshToMeshInterp, addParticles);
256  MapLagrangianFields<sphericalTensor>
257  (cloudDirs[cloudI], objects, meshToMeshInterp, addParticles);
258  MapLagrangianFields<symmTensor>
259  (cloudDirs[cloudI], objects, meshToMeshInterp, addParticles);
260  MapLagrangianFields<tensor>
261  (cloudDirs[cloudI], objects, meshToMeshInterp, addParticles);
262  }
263  }
264  }
265 }
266 
267 
268 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
269 
270 } // End namespace Foam
271 
272 // ************************ vim: set sw=4 sts=4 et: ************************ //