mirrorFvMesh.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 
26 #include "mirrorFvMesh.H"
27 #include "Time.H"
28 #include "plane.H"
29 
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
31 
32 Foam::mirrorFvMesh::mirrorFvMesh(const IOobject& io, const IOobject& dictIO)
33 :
34  fvMesh(io, false),
35  mirrorMeshDict_(dictIO)
36 {
37  plane mirrorPlane(mirrorMeshDict_);
38 
39  scalar planeTolerance
40  (
41  mirrorMeshDict_.lookup<scalar>("planeTolerance")
42  );
43 
44  const pointField& oldPoints = points();
45  const faceList& oldFaces = faces();
46  const cellList& oldCells = cells();
47  const label nOldInternalFaces = nInternalFaces();
48  const polyPatchList& oldPatches = boundaryMesh();
49 
50  // Mirror the points
51  Info<< "Mirroring points. Old points: " << oldPoints.size();
52 
53  pointField newPoints(2*oldPoints.size());
54  label nNewPoints = 0;
55 
56  labelList mirrorPointLookup(oldPoints.size(), -1);
57 
58  // Grab the old points
59  forAll(oldPoints, pointi)
60  {
61  newPoints[nNewPoints] = oldPoints[pointi];
62  nNewPoints++;
63  }
64 
65  forAll(oldPoints, pointi)
66  {
67  scalar alpha =
68  mirrorPlane.normalIntersect
69  (
70  oldPoints[pointi],
71  mirrorPlane.normal()
72  );
73 
74  // Check plane on tolerance
75  if (mag(alpha) > planeTolerance)
76  {
77  // The point gets mirrored
78  newPoints[nNewPoints] =
79  oldPoints[pointi] + 2.0*alpha*mirrorPlane.normal();
80 
81  // remember the point correspondence
82  mirrorPointLookup[pointi] = nNewPoints;
83  nNewPoints++;
84  }
85  else
86  {
87  // The point is on the plane and does not get mirrored
88  // Adjust plane location
89  newPoints[nNewPoints] =
90  oldPoints[pointi] + alpha*mirrorPlane.normal();
91 
92  mirrorPointLookup[pointi] = pointi;
93  }
94  }
95 
96  // Reset the size of the point list
97  Info<< " New points: " << nNewPoints << endl;
98  newPoints.setSize(nNewPoints);
99 
100  // Construct new to old map
101  pointMapPtr_.reset(new labelList(newPoints.size()));
102  labelList& pointMap = pointMapPtr_();
103  // Insert old points
104  forAll(oldPoints, oldPointi)
105  {
106  pointMap[oldPointi] = oldPointi;
107  }
108  forAll(mirrorPointLookup, oldPointi)
109  {
110  pointMap[mirrorPointLookup[oldPointi]] = oldPointi;
111  }
112 
113 
114 
115  Info<< "Mirroring faces. Old faces: " << oldFaces.size();
116 
117  // Algorithm:
118  // During mirroring, the faces that were previously boundary faces
119  // in the mirror plane may become ineternal faces. In order to
120  // deal with the ordering of the faces, the algorithm is split
121  // into two parts. For original faces, the internal faces are
122  // distributed to their owner cells. Once all internal faces are
123  // distributed, the boundary faces are visited and if they are in
124  // the mirror plane they are added to the master cells (the future
125  // boundary faces are not touched). After the first phase, the
126  // internal faces are collected in the cell order and numbering
127  // information is added. Then, the internal faces are mirrored
128  // and the face numbering data is stored for the mirrored section.
129  // Once all the internal faces are mirrored, the boundary faces
130  // are added by mirroring the faces patch by patch.
131 
132  // Distribute internal faces
133  labelListList newCellFaces(oldCells.size());
134 
135  const labelUList& oldOwnerStart = lduAddr().ownerStartAddr();
136 
137  forAll(newCellFaces, celli)
138  {
139  labelList& curFaces = newCellFaces[celli];
140 
141  const label s = oldOwnerStart[celli];
142  const label e = oldOwnerStart[celli + 1];
143 
144  curFaces.setSize(e - s);
145 
146  forAll(curFaces, i)
147  {
148  curFaces[i] = s + i;
149  }
150  }
151 
152  // Distribute boundary faces. Remember the faces that have been inserted
153  // as internal
154  boolListList insertedBouFace(oldPatches.size());
155 
156  forAll(oldPatches, patchi)
157  {
158  const polyPatch& curPatch = oldPatches[patchi];
159 
160  if (curPatch.coupled())
161  {
163  << "Found coupled patch " << curPatch.name() << endl
164  << " Mirroring faces on coupled patches destroys"
165  << " the ordering. This might be fixed by running a dummy"
166  << " createPatch afterwards." << endl;
167  }
168 
169  boolList& curInsBouFace = insertedBouFace[patchi];
170 
171  curInsBouFace.setSize(curPatch.size());
172  curInsBouFace = false;
173 
174  // Get faceCells for face insertion
175  const labelUList& curFaceCells = curPatch.faceCells();
176  const label curStart = curPatch.start();
177 
178  forAll(curPatch, facei)
179  {
180  // Find out if the mirrored face is identical to the
181  // original. If so, the face needs to become internal and
182  // added to its owner cell
183  const face& origFace = curPatch[facei];
184 
185  face mirrorFace(origFace.size());
186  forAll(mirrorFace, pointi)
187  {
188  mirrorFace[pointi] = mirrorPointLookup[origFace[pointi]];
189  }
190 
191  if (origFace == mirrorFace)
192  {
193  // The mirror is identical to current face. This will
194  // become an internal face
195  const label oldSize = newCellFaces[curFaceCells[facei]].size();
196 
197  newCellFaces[curFaceCells[facei]].setSize(oldSize + 1);
198  newCellFaces[curFaceCells[facei]][oldSize] = curStart + facei;
199 
200  curInsBouFace[facei] = true;
201  }
202  }
203  }
204 
205  // Construct the new list of faces. Boundary faces are added
206  // last, cush that each patch is mirrored separately. The
207  // addressing is stored in two separate arrays: first for the
208  // original cells (face order has changed) and then for the
209  // mirrored cells.
210  labelList masterFaceLookup(oldFaces.size(), -1);
211  labelList mirrorFaceLookup(oldFaces.size(), -1);
212 
213  faceList newFaces(2*oldFaces.size());
214  label nNewFaces = 0;
215 
216  // Insert original (internal) faces
217  forAll(newCellFaces, celli)
218  {
219  const labelList& curCellFaces = newCellFaces[celli];
220 
221  forAll(curCellFaces, cfI)
222  {
223  newFaces[nNewFaces] = oldFaces[curCellFaces[cfI]];
224  masterFaceLookup[curCellFaces[cfI]] = nNewFaces;
225 
226  nNewFaces++;
227  }
228  }
229 
230  // Mirror internal faces
231  for (label facei = 0; facei < nOldInternalFaces; facei++)
232  {
233  const face& oldFace = oldFaces[facei];
234  face& nf = newFaces[nNewFaces];
235  nf.setSize(oldFace.size());
236 
237  nf[0] = mirrorPointLookup[oldFace[0]];
238 
239  for (label i = 1; i < oldFace.size(); i++)
240  {
241  nf[i] = mirrorPointLookup[oldFace[oldFace.size() - i]];
242  }
243 
244  mirrorFaceLookup[facei] = nNewFaces;
245  nNewFaces++;
246  }
247 
248  // Mirror boundary faces patch by patch
249 
250 
251  labelList newToOldPatch(boundary().size(), -1);
252  labelList newPatchSizes(boundary().size(), -1);
253  labelList newPatchStarts(boundary().size(), -1);
254  label nNewPatches = 0;
255 
256  forAll(boundaryMesh(), patchi)
257  {
258  const label curPatchSize = boundaryMesh()[patchi].size();
259  const label curPatchStart = boundaryMesh()[patchi].start();
260  const boolList& curInserted = insertedBouFace[patchi];
261 
262  newPatchStarts[nNewPatches] = nNewFaces;
263 
264  // Master side
265  for (label facei = 0; facei < curPatchSize; facei++)
266  {
267  // Check if the face has already been added. If not, add it and
268  // insert the numbering details.
269  if (!curInserted[facei])
270  {
271  newFaces[nNewFaces] = oldFaces[curPatchStart + facei];
272 
273  masterFaceLookup[curPatchStart + facei] = nNewFaces;
274  nNewFaces++;
275  }
276  }
277 
278  // Mirror side
279  for (label facei = 0; facei < curPatchSize; facei++)
280  {
281  // Check if the face has already been added. If not, add it and
282  // insert the numbering details.
283  if (!curInserted[facei])
284  {
285  const face& oldFace = oldFaces[curPatchStart + facei];
286  face& nf = newFaces[nNewFaces];
287  nf.setSize(oldFace.size());
288 
289  nf[0] = mirrorPointLookup[oldFace[0]];
290 
291  for (label i = 1; i < oldFace.size(); i++)
292  {
293  nf[i] = mirrorPointLookup[oldFace[oldFace.size() - i]];
294  }
295 
296  mirrorFaceLookup[curPatchStart + facei] = nNewFaces;
297  nNewFaces++;
298  }
299  else
300  {
301  // Grab the index of the master face for the mirror side
302  mirrorFaceLookup[curPatchStart + facei] =
303  masterFaceLookup[curPatchStart + facei];
304  }
305  }
306 
307  // If patch exists, grab the name and type of the original patch
308  if (nNewFaces > newPatchStarts[nNewPatches])
309  {
310  newToOldPatch[nNewPatches] = patchi;
311 
312  newPatchSizes[nNewPatches] =
313  nNewFaces - newPatchStarts[nNewPatches];
314 
315  nNewPatches++;
316  }
317  }
318 
319  // Tidy up the lists
320  newFaces.setSize(nNewFaces);
321  Info<< " New faces: " << nNewFaces << endl;
322 
323  newToOldPatch.setSize(nNewPatches);
324  newPatchSizes.setSize(nNewPatches);
325  newPatchStarts.setSize(nNewPatches);
326 
327  Info<< "Mirroring patches. Old patches: " << boundary().size()
328  << " New patches: " << nNewPatches << endl;
329 
330  Info<< "Mirroring cells. Old cells: " << oldCells.size()
331  << " New cells: " << 2*oldCells.size() << endl;
332 
333  cellList newCells(2*oldCells.size());
334  label nNewCells = 0;
335 
336  // Construct new to old cell map
337  cellMapPtr_.reset(new labelList(newCells.size()));
338  labelList& cellMap = cellMapPtr_();
339 
340  // Grab the original cells. Take care of face renumbering.
341  forAll(oldCells, celli)
342  {
343  const cell& oc = oldCells[celli];
344 
345  cell& nc = newCells[nNewCells];
346  nc.setSize(oc.size());
347 
348  forAll(oc, i)
349  {
350  nc[i] = masterFaceLookup[oc[i]];
351  }
352 
353  cellMap[nNewCells] = celli;
354 
355  nNewCells++;
356  }
357 
358  // Mirror the cells
359  forAll(oldCells, celli)
360  {
361  const cell& oc = oldCells[celli];
362 
363  cell& nc = newCells[nNewCells];
364  nc.setSize(oc.size());
365 
366  forAll(oc, i)
367  {
368  nc[i] = mirrorFaceLookup[oc[i]];
369  }
370 
371  cellMap[nNewCells] = celli;
372 
373  nNewCells++;
374  }
375 
376  // Mirror the cell shapes
377  Info<< "Mirroring cell shapes." << endl;
378 
379  Info<< nl << "Creating new mesh" << endl;
380  mirrorMeshPtr_.reset
381  (
382  new fvMesh
383  (
384  io,
385  std::move(newPoints),
386  std::move(newFaces),
387  std::move(newCells)
388  )
389  );
390 
391  fvMesh& pMesh = mirrorMeshPtr_();
392 
393  // Add the boundary patches
394  List<polyPatch*> p(newPatchSizes.size());
395 
396  forAll(p, patchi)
397  {
398  p[patchi] = boundaryMesh()[newToOldPatch[patchi]].clone
399  (
400  pMesh.boundaryMesh(),
401  patchi,
402  newPatchSizes[patchi],
403  newPatchStarts[patchi]
404  ).ptr();
405  }
406 
407  pMesh.addPatches(p);
408 }
409 
410 
411 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
412 
414 {}
415 
416 
417 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
void setSize(const label)
Reset size of List.
Definition: List.C:281
mirrorFvMesh(const IOobject &io, const IOobject &dictIO)
Construct from IOobject.
~mirrorFvMesh()
Destructor.
label patchi
const pointField & points
const cellShapeList & cells
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar e
Elementary charge.
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
List< cell > cellList
list of cells
Definition: cellList.H:42
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
messageStream Info
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
dimensioned< scalar > mag(const dimensioned< Type > &)
List< List< bool > > boolListList
Definition: boolList.H:51
PtrList< polyPatch > polyPatchList
container classes for polyPatch
Definition: polyPatchList.H:45
List< face > faceList
Definition: faceListFwd.H:43
UList< label > labelUList
Definition: UList.H:65
static const char nl
Definition: Ostream.H:260
faceListList boundary(nPatches)
volScalarField & p