FreeFOAM The Cross-Platform CFD Toolkit
modifyMesh.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 Application
25  modifyMesh
26 
27 Description
28  Manipulates mesh elements.
29 
30  Actions are:
31  (boundary)points:
32  - move
33 
34  (boundary)edges:
35  - split and move introduced point
36 
37  (boundary)faces:
38  - split(triangulate) and move introduced point
39 
40  edges:
41  - collapse
42 
43  cells:
44  - split into polygonal base pyramids around newly introduced mid
45  point
46 
47  Is a bit of a loose collection of mesh change drivers.
48 
49 Usage
50 
51  - modifyMesh [OPTIONS]
52 
53  @param -overwrite \n
54  Overwrite existing data.
55 
56  @param -case <dir>\n
57  Case directory.
58 
59  @param -parallel \n
60  Run in parallel.
61 
62  @param -help \n
63  Display help message.
64 
65  @param -doc \n
66  Display Doxygen API documentation page for this application.
67 
68  @param -srcDoc \n
69  Display Doxygen source documentation page for this application.
70 
71 \*---------------------------------------------------------------------------*/
72 
73 #include <OpenFOAM/argList.H>
74 #include <OpenFOAM/Time.H>
75 #include <OpenFOAM/polyMesh.H>
77 #include <OpenFOAM/mapPolyMesh.H>
79 #include "cellSplitter.H"
81 #include <meshTools/meshTools.H>
82 #include <OpenFOAM/Pair.H>
83 
84 using namespace Foam;
85 
86 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
87 
88 // Locate point on patch. Returns (mesh) point label.
89 label findPoint(const primitivePatch& pp, const point& nearPoint)
90 {
91  const pointField& points = pp.points();
92  const labelList& meshPoints = pp.meshPoints();
93 
94  // Find nearest and next nearest
95  scalar minDistSqr = GREAT;
96  label minI = -1;
97 
98  scalar almostMinDistSqr = GREAT;
99  label almostMinI = -1;
100 
101  forAll(meshPoints, i)
102  {
103  label pointI = meshPoints[i];
104 
105  scalar distSqr = magSqr(nearPoint - points[pointI]);
106 
107  if (distSqr < minDistSqr)
108  {
109  almostMinDistSqr = minDistSqr;
110  almostMinI = minI;
111 
112  minDistSqr = distSqr;
113  minI = pointI;
114  }
115  else if (distSqr < almostMinDistSqr)
116  {
117  almostMinDistSqr = distSqr;
118  almostMinI = pointI;
119  }
120  }
121 
122 
123  // Decide if nearPoint unique enough.
124  Info<< "Found to point " << nearPoint << nl
125  << " nearest point : " << minI
126  << " distance " << Foam::sqrt(minDistSqr)
127  << " at " << points[minI] << nl
128  << " next nearest point : " << almostMinI
129  << " distance " << Foam::sqrt(almostMinDistSqr)
130  << " at " << points[almostMinI] << endl;
131 
132  if (almostMinDistSqr < 4*minDistSqr)
133  {
134  Info<< "Next nearest too close to nearest. Aborting" << endl;
135 
136  return -1;
137  }
138  else
139  {
140  return minI;
141  }
142 }
143 
144 
145 // Locate edge on patch. Return mesh edge label.
146 label findEdge
147 (
148  const primitiveMesh& mesh,
149  const primitivePatch& pp,
150  const point& nearPoint
151 )
152 {
153  const pointField& localPoints = pp.localPoints();
154  const pointField& points = pp.points();
155  const labelList& meshPoints = pp.meshPoints();
156  const edgeList& edges = pp.edges();
157 
158  // Find nearest and next nearest
159  scalar minDist = GREAT;
160  label minI = -1;
161 
162  scalar almostMinDist = GREAT;
163  label almostMinI = -1;
164 
165  forAll(edges, edgeI)
166  {
167  const edge& e = edges[edgeI];
168 
169  pointHit pHit(e.line(localPoints).nearestDist(nearPoint));
170 
171  if (pHit.hit())
172  {
173  if (pHit.distance() < minDist)
174  {
175  almostMinDist = minDist;
176  almostMinI = minI;
177 
178  minDist = pHit.distance();
179  minI = meshTools::findEdge
180  (
181  mesh,
182  meshPoints[e[0]],
183  meshPoints[e[1]]
184  );
185  }
186  else if (pHit.distance() < almostMinDist)
187  {
188  almostMinDist = pHit.distance();
189  almostMinI = meshTools::findEdge
190  (
191  mesh,
192  meshPoints[e[0]],
193  meshPoints[e[1]]
194  );
195  }
196  }
197  }
198 
199  if (minI == -1)
200  {
201  Info<< "Did not find edge close to point " << nearPoint << endl;
202 
203  return -1;
204  }
205 
206 
207  // Decide if nearPoint unique enough.
208  Info<< "Found to point " << nearPoint << nl
209  << " nearest edge : " << minI
210  << " distance " << minDist << " endpoints "
211  << mesh.edges()[minI].line(points) << nl
212  << " next nearest edge : " << almostMinI
213  << " distance " << almostMinDist << " endpoints "
214  << mesh.edges()[almostMinI].line(points) << nl
215  << endl;
216 
217  if (almostMinDist < 2*minDist)
218  {
219  Info<< "Next nearest too close to nearest. Aborting" << endl;
220 
221  return -1;
222  }
223  else
224  {
225  return minI;
226  }
227 }
228 
229 
230 // Find face on patch. Return mesh face label.
231 label findFace
232 (
233  const primitiveMesh& mesh,
234  const primitivePatch& pp,
235  const point& nearPoint
236 )
237 {
238  const pointField& points = pp.points();
239 
240  // Find nearest and next nearest
241  scalar minDist = GREAT;
242  label minI = -1;
243 
244  scalar almostMinDist = GREAT;
245  label almostMinI = -1;
246 
247  forAll(pp, patchFaceI)
248  {
249  pointHit pHit(pp[patchFaceI].nearestPoint(nearPoint, points));
250 
251  if (pHit.hit())
252  {
253  if (pHit.distance() < minDist)
254  {
255  almostMinDist = minDist;
256  almostMinI = minI;
257 
258  minDist = pHit.distance();
259  minI = patchFaceI + mesh.nInternalFaces();
260  }
261  else if (pHit.distance() < almostMinDist)
262  {
263  almostMinDist = pHit.distance();
264  almostMinI = patchFaceI + mesh.nInternalFaces();
265  }
266  }
267  }
268 
269  if (minI == -1)
270  {
271  Info<< "Did not find face close to point " << nearPoint << endl;
272 
273  return -1;
274  }
275 
276 
277  // Decide if nearPoint unique enough.
278  Info<< "Found to point " << nearPoint << nl
279  << " nearest face : " << minI
280  << " distance " << minDist
281  << " to face centre " << mesh.faceCentres()[minI] << nl
282  << " next nearest face : " << almostMinI
283  << " distance " << almostMinDist
284  << " to face centre " << mesh.faceCentres()[almostMinI] << nl
285  << endl;
286 
287  if (almostMinDist < 2*minDist)
288  {
289  Info<< "Next nearest too close to nearest. Aborting" << endl;
290 
291  return -1;
292  }
293  else
294  {
295  return minI;
296  }
297 }
298 
299 
300 // Find cell with cell centre close to given point.
301 label findCell(const primitiveMesh& mesh, const point& nearPoint)
302 {
303  label cellI = mesh.findCell(nearPoint);
304 
305  if (cellI != -1)
306  {
307  scalar distToCcSqr = magSqr(nearPoint - mesh.cellCentres()[cellI]);
308 
309  const labelList& cPoints = mesh.cellPoints()[cellI];
310 
311  label minI = -1;
312  scalar minDistSqr = GREAT;
313 
314  forAll(cPoints, i)
315  {
316  label pointI = cPoints[i];
317 
318  scalar distSqr = magSqr(nearPoint - mesh.points()[pointI]);
319 
320  if (distSqr < minDistSqr)
321  {
322  minDistSqr = distSqr;
323  minI = pointI;
324  }
325  }
326 
327  // Decide if nearPoint unique enough.
328  Info<< "Found to point " << nearPoint << nl
329  << " nearest cell : " << cellI
330  << " distance " << Foam::sqrt(distToCcSqr)
331  << " to cell centre " << mesh.cellCentres()[cellI] << nl
332  << " nearest mesh point : " << minI
333  << " distance " << Foam::sqrt(minDistSqr)
334  << " to " << mesh.points()[minI] << nl
335  << endl;
336 
337  if (minDistSqr < 4*distToCcSqr)
338  {
339  Info<< "Mesh point too close to nearest cell centre. Aborting"
340  << endl;
341 
342  cellI = -1;
343  }
344  }
345 
346  return cellI;
347 }
348 
349 
350 
351 // Main program:
352 
353 int main(int argc, char *argv[])
354 {
355  argList::validOptions.insert("overwrite", "");
356 
357 # include <OpenFOAM/setRootCase.H>
358 # include <OpenFOAM/createTime.H>
359  runTime.functionObjects().off();
360 # include <OpenFOAM/createPolyMesh.H>
361  const word oldInstance = mesh.pointsInstance();
362 
363  bool overwrite = args.optionFound("overwrite");
364 
365  Info<< "Reading modifyMeshDict\n" << endl;
366 
367  IOdictionary dict
368  (
369  IOobject
370  (
371  "modifyMeshDict",
372  runTime.system(),
373  mesh,
376  )
377  );
378 
379  // Read all from the dictionary.
380  List<Pair<point> > pointsToMove(dict.lookup("pointsToMove"));
381  List<Pair<point> > edgesToSplit(dict.lookup("edgesToSplit"));
382  List<Pair<point> > facesToTriangulate
383  (
384  dict.lookup("facesToTriangulate")
385  );
386 
387  bool cutBoundary =
388  (
389  pointsToMove.size()
390  || edgesToSplit.size()
391  || facesToTriangulate.size()
392  );
393 
394  List<Pair<point> > edgesToCollapse(dict.lookup("edgesToCollapse"));
395 
396  bool collapseEdge = edgesToCollapse.size();
397 
398  List<Pair<point> > cellsToPyramidise(dict.lookup("cellsToSplit"));
399 
400  bool cellsToSplit = cellsToPyramidise.size();
401 
402  //List<Tuple<pointField,point> >
403  // cellsToCreate(dict.lookup("cellsToCreate"));
404 
405  Info<< "Read from " << dict.name() << nl
406  << " Boundary cutting module:" << nl
407  << " points to move :" << pointsToMove.size() << nl
408  << " edges to split :" << edgesToSplit.size() << nl
409  << " faces to triangulate:" << facesToTriangulate.size() << nl
410  << " Cell splitting module:" << nl
411  << " cells to split :" << cellsToPyramidise.size() << nl
412  << " Edge collapsing module:" << nl
413  << " edges to collapse :" << edgesToCollapse.size() << nl
414  //<< " cells to create :" << cellsToCreate.size() << nl
415  << endl;
416 
417  if
418  (
419  (cutBoundary && collapseEdge)
420  || (cutBoundary && cellsToSplit)
421  || (collapseEdge && cellsToSplit)
422  )
423  {
425  << "Used more than one mesh modifying module "
426  << "(boundary cutting, cell splitting, edge collapsing)" << nl
427  << "Please do them in separate passes." << exit(FatalError);
428  }
429 
430 
431 
432  // Get calculating engine for all of outside
433  const SubList<face> outsideFaces
434  (
435  mesh.faces(),
436  mesh.nFaces() - mesh.nInternalFaces(),
437  mesh.nInternalFaces()
438  );
439 
440  primitivePatch allBoundary(outsideFaces, mesh.points());
441 
442 
443  // Look up mesh labels and convert to input for boundaryCutter.
444 
445  bool validInputs = true;
446 
447 
448  Info<< nl << "Looking up points to move ..." << nl << endl;
449  Map<point> pointToPos(pointsToMove.size());
450  forAll(pointsToMove, i)
451  {
452  const Pair<point>& pts = pointsToMove[i];
453 
454  label pointI = findPoint(allBoundary, pts.first());
455 
456  if (pointI == -1 || !pointToPos.insert(pointI, pts.second()))
457  {
458  Info<< "Could not insert mesh point " << pointI
459  << " for input point " << pts.first() << nl
460  << "Perhaps the point is already marked for moving?" << endl;
461  validInputs = false;
462  }
463  }
464 
465 
466  Info<< nl << "Looking up edges to split ..." << nl << endl;
467  Map<List<point> > edgeToCuts(edgesToSplit.size());
468  forAll(edgesToSplit, i)
469  {
470  const Pair<point>& pts = edgesToSplit[i];
471 
472  label edgeI = findEdge(mesh, allBoundary, pts.first());
473 
474  if
475  (
476  edgeI == -1
477  || !edgeToCuts.insert(edgeI, List<point>(1, pts.second()))
478  )
479  {
480  Info<< "Could not insert mesh edge " << edgeI
481  << " for input point " << pts.first() << nl
482  << "Perhaps the edge is already marked for cutting?" << endl;
483 
484  validInputs = false;
485  }
486  }
487 
488 
489  Info<< nl << "Looking up faces to triangulate ..." << nl << endl;
490  Map<point> faceToDecompose(facesToTriangulate.size());
491  forAll(facesToTriangulate, i)
492  {
493  const Pair<point>& pts = facesToTriangulate[i];
494 
495  label faceI = findFace(mesh, allBoundary, pts.first());
496 
497  if (faceI == -1 || !faceToDecompose.insert(faceI, pts.second()))
498  {
499  Info<< "Could not insert mesh face " << faceI
500  << " for input point " << pts.first() << nl
501  << "Perhaps the face is already marked for splitting?" << endl;
502 
503  validInputs = false;
504  }
505  }
506 
507 
508 
509  Info<< nl << "Looking up cells to convert to pyramids around"
510  << " cell centre ..." << nl << endl;
511  Map<point> cellToPyrCentre(cellsToPyramidise.size());
512  forAll(cellsToPyramidise, i)
513  {
514  const Pair<point>& pts = cellsToPyramidise[i];
515 
516  label cellI = findCell(mesh, pts.first());
517 
518  if (cellI == -1 || !cellToPyrCentre.insert(cellI, pts.second()))
519  {
520  Info<< "Could not insert mesh cell " << cellI
521  << " for input point " << pts.first() << nl
522  << "Perhaps the cell is already marked for splitting?" << endl;
523 
524  validInputs = false;
525  }
526  }
527 
528 
529  Info<< nl << "Looking up edges to collapse ..." << nl << endl;
530  Map<point> edgeToPos(edgesToCollapse.size());
531  forAll(edgesToCollapse, i)
532  {
533  const Pair<point>& pts = edgesToCollapse[i];
534 
535  label edgeI = findEdge(mesh, allBoundary, pts.first());
536 
537  if (edgeI == -1 || !edgeToPos.insert(edgeI, pts.second()))
538  {
539  Info<< "Could not insert mesh edge " << edgeI
540  << " for input point " << pts.first() << nl
541  << "Perhaps the edge is already marked for collaping?" << endl;
542 
543  validInputs = false;
544  }
545  }
546 
547 
548 
549  if (!validInputs)
550  {
551  Info<< nl << "There was a problem in one of the inputs in the"
552  << " dictionary. Not modifying mesh." << endl;
553  }
554  else if (cellToPyrCentre.size())
555  {
556  Info<< nl << "All input cells located. Modifying mesh." << endl;
557 
558  // Mesh change engine
559  cellSplitter cutter(mesh);
560 
561  // Topo change container
562  polyTopoChange meshMod(mesh);
563 
564  // Insert commands into meshMod
565  cutter.setRefinement(cellToPyrCentre, meshMod);
566 
567  // Do changes
568  autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
569 
570  if (morphMap().hasMotionPoints())
571  {
572  mesh.movePoints(morphMap().preMotionPoints());
573  }
574 
575  cutter.updateMesh(morphMap());
576 
577  if (!overwrite)
578  {
579  runTime++;
580  }
581  else
582  {
583  mesh.setInstance(oldInstance);
584  }
585 
586  // Write resulting mesh
587  Info << "Writing modified mesh to time " << runTime.timeName() << endl;
588  mesh.write();
589  }
590  else if (edgeToPos.size())
591  {
592  Info<< nl << "All input edges located. Modifying mesh." << endl;
593 
594  // Mesh change engine
595  edgeCollapser cutter(mesh);
596 
597  pointField newPoints(mesh.points());
598 
599  // Get new positions and construct collapse network
600  forAllConstIter(Map<point>, edgeToPos, iter)
601  {
602  label edgeI = iter.key();
603  const edge& e = mesh.edges()[edgeI];
604 
605  cutter.collapseEdge(edgeI, e[0]);
606  newPoints[e[0]] = iter();
607  }
608 
609  // Move master point to destination.
610  mesh.movePoints(newPoints);
611 
612  // Topo change container
613  polyTopoChange meshMod(mesh);
614 
615  // Insert
616  cutter.setRefinement(meshMod);
617 
618  // Do changes
619  autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
620 
621  if (morphMap().hasMotionPoints())
622  {
623  mesh.movePoints(morphMap().preMotionPoints());
624  }
625 
626  // Not implemented yet:
627  //cutter.updateMesh(morphMap());
628 
629 
630  if (!overwrite)
631  {
632  runTime++;
633  }
634  else
635  {
636  mesh.setInstance(oldInstance);
637  }
638 
639  // Write resulting mesh
640  Info << "Writing modified mesh to time " << runTime.timeName() << endl;
641  mesh.write();
642  }
643  else
644  {
645  Info<< nl << "All input points located. Modifying mesh." << endl;
646 
647  // Mesh change engine
648  boundaryCutter cutter(mesh);
649 
650  // Topo change container
651  polyTopoChange meshMod(mesh);
652 
653  // Insert commands into meshMod
654  cutter.setRefinement
655  (
656  pointToPos,
657  edgeToCuts,
658  Map<labelPair>(0), // Faces to split diagonally
659  faceToDecompose, // Faces to triangulate
660  meshMod
661  );
662 
663  // Do changes
664  autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
665 
666  if (morphMap().hasMotionPoints())
667  {
668  mesh.movePoints(morphMap().preMotionPoints());
669  }
670 
671  cutter.updateMesh(morphMap());
672 
673  if (!overwrite)
674  {
675  runTime++;
676  }
677  else
678  {
679  mesh.setInstance(oldInstance);
680  }
681 
682  // Write resulting mesh
683  Info << "Writing modified mesh to time " << runTime.timeName() << endl;
684  mesh.write();
685  }
686 
687 
688  Info << nl << "End" << endl;
689 
690  return 0;
691 }
692 
693 
694 // ************************ vim: set sw=4 sts=4 et: ************************ //