| blob | 8b83eff6bd56874a34a393f384734d3160bfaf48 |
1 /*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | foam-extend: Open Source CFD
4 \\ / O peration | Version: 3.2
5 \\ / A nd | Web: http://www.foam-extend.org
6 \\/ M anipulation | For copyright notice see file Copyright
7 -------------------------------------------------------------------------------
8 License
9 This file is part of foam-extend.
11 foam-extend is free software: you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation, either version 3 of the License, or (at your
14 option) any later version.
16 foam-extend is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 General Public License for more details.
21 You should have received a copy of the GNU General Public License
22 along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
24 Class
25 dynamicTopoFvMesh
27 Description
28 An implementation of dynamic changes to mesh-topology
30 Author
31 Sandeep Menon
32 University of Massachusetts Amherst
33 All rights reserved
35 \*---------------------------------------------------------------------------*/
37 #include "dynamicTopoFvMesh.H"
38 #include "addToRunTimeSelectionTable.H"
40 #include "eMesh.H"
41 #include "Stack.H"
42 #include "triFace.H"
43 #include "changeMap.H"
44 #include "clockTime.H"
45 #include "volFields.H"
46 #include "MeshObject.H"
47 #include "topoMapper.H"
48 #include "coupledInfo.H"
49 #include "mapPolyMesh.H"
50 #include "MapFvFields.H"
51 #include "SortableList.H"
52 #include "motionSolver.H"
53 #include "fvPatchFields.H"
54 #include "fvsPatchFields.H"
55 #include "subMeshLduAddressing.H"
56 #include "lengthScaleEstimator.H"
57 #include "conservativeMapFields.H"
59 namespace Foam
60 {
62 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
64 defineTypeNameAndDebug(dynamicTopoFvMesh,0);
66 addToRunTimeSelectionTable(dynamicFvMesh, dynamicTopoFvMesh, IOobject);
68 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
70 // Construct from IOobject
71 dynamicTopoFvMesh::dynamicTopoFvMesh(const IOobject& io)
72 :
73 dynamicFvMesh(io),
74 baseMesh_(*this),
75 nPatches_(boundaryMesh().size()),
76 topoChangeFlag_(false),
77 isSubMesh_(false),
78 dict_
79 (
80 IOobject
81 (
82 "dynamicMeshDict",
83 polyMesh::time().constant(),
84 (*this),
85 IOobject::MUST_READ,
86 IOobject::NO_WRITE,
87 false
88 )
89 ),
90 mandatory_
91 (
92 dict_.subDict("dynamicTopoFvMesh").lookup("allOptionsMandatory")
93 ),
94 twoDMesh_(polyMesh::nGeometricD() == 2 ? true : false),
95 edgeRefinement_
96 (
97 dict_.subDict("dynamicTopoFvMesh").lookup("edgeRefinement")
98 ),
99 loadMotionSolver_(true),
100 bandWidthReduction_(false),
101 coupledModification_(false),
102 lduPtr_(NULL),
103 interval_(1),
104 eMeshPtr_(NULL),
105 mapper_(NULL),
106 motionSolver_(NULL),
107 lengthEstimator_(NULL),
108 oldPoints_(polyMesh::points()),
109 points_(polyMesh::points()),
110 faces_(polyMesh::faces()),
111 owner_(polyMesh::faceOwner()),
112 neighbour_(polyMesh::faceNeighbour()),
113 cells_(primitiveMesh::cells()),
114 oldPatchSizes_(nPatches_, 0),
115 patchSizes_(nPatches_, 0),
116 oldPatchStarts_(nPatches_, -1),
117 patchStarts_(nPatches_, -1),
118 oldEdgePatchSizes_(nPatches_, 0),
119 edgePatchSizes_(nPatches_, 0),
120 oldEdgePatchStarts_(nPatches_, -1),
121 edgePatchStarts_(nPatches_, -1),
122 oldPatchNMeshPoints_(nPatches_, -1),
123 patchNMeshPoints_(nPatches_, -1),
124 nOldPoints_(primitiveMesh::nPoints()),
125 nPoints_(primitiveMesh::nPoints()),
126 nOldEdges_(0),
127 nEdges_(0),
128 nOldFaces_(primitiveMesh::nFaces()),
129 nFaces_(primitiveMesh::nFaces()),
130 nOldCells_(primitiveMesh::nCells()),
131 nCells_(primitiveMesh::nCells()),
132 nOldInternalFaces_(primitiveMesh::nInternalFaces()),
133 nInternalFaces_(primitiveMesh::nInternalFaces()),
134 nOldInternalEdges_(0),
135 nInternalEdges_(0),
136 maxModifications_(-1),
137 statistics_(TOTAL_OP_TYPES, 0),
138 sliverThreshold_(0.1),
139 slicePairs_(0),
140 maxTetsPerEdge_(-1),
141 swapDeviation_(0.0),
142 allowTableResize_(false)
143 {
144 // Check the size of owner/neighbour
145 if (owner_.size() != neighbour_.size())
146 {
147 // Size up to number of faces
148 neighbour_.setSize(nFaces_, -1);
149 }
151 const polyBoundaryMesh& boundary = polyMesh::boundaryMesh();
153 // Initialize the multiThreading environment
154 initializeThreadingEnvironment();
156 // Read optional parameters.
157 readOptionalParameters();
159 // Initialize patch-size information
160 for (label i = 0; i < nPatches_; i++)
161 {
162 patchNMeshPoints_[i] = boundary[i].meshPoints().size();
163 oldPatchSizes_[i] = patchSizes_[i] = boundary[i].size();
164 oldPatchStarts_[i] = patchStarts_[i] = boundary[i].start();
165 oldPatchNMeshPoints_[i] = patchNMeshPoints_[i];
166 }
168 // Open the tetMetric dynamic-link library (for 3D only)
169 loadMetric();
171 // Initialize edge-related connectivity structures
172 initEdges();
174 // Load the mesh-motion solver
175 loadMotionSolver();
177 // Load the field-mapper
178 loadFieldMapper();
180 // Load the length-scale estimator,
181 // and read refinement options
182 loadLengthScaleEstimator();
183 }
186 //- Construct from components. Used for subMeshes only.
187 dynamicTopoFvMesh::dynamicTopoFvMesh
188 (
189 const dynamicTopoFvMesh& mesh,
190 const IOobject& io,
191 const Xfer<pointField>& points,
192 const Xfer<pointField>& oldPoints,
193 const Xfer<edgeList>& edges,
194 const Xfer<faceList>& faces,
195 const Xfer<labelListList>& faceEdges,
196 const Xfer<labelList>& owner,
197 const Xfer<labelList>& neighbour,
198 const labelList& faceStarts,
199 const labelList& faceSizes,
200 const labelList& edgeStarts,
201 const labelList& edgeSizes,
202 const wordList& patchNames,
203 const dictionary& patchDict
204 )
205 :
206 dynamicFvMesh
207 (
208 io,
209 oldPoints,
210 faces,
211 owner,
212 neighbour,
213 false
214 ),
215 baseMesh_(mesh),
216 nPatches_(faceStarts.size()),
217 topoChangeFlag_(false),
218 isSubMesh_(true),
219 dict_(mesh.dict_),
220 mandatory_(mesh.mandatory_),
221 twoDMesh_(mesh.twoDMesh_),
222 edgeRefinement_(mesh.edgeRefinement_),
223 loadMotionSolver_(mesh.loadMotionSolver_),
224 bandWidthReduction_(mesh.bandWidthReduction_),
225 coupledModification_(false),
226 lduPtr_(NULL),
227 interval_(1),
228 eMeshPtr_(NULL),
229 mapper_(NULL),
230 motionSolver_(NULL),
231 lengthEstimator_(NULL),
232 oldPoints_(polyMesh::points()),
233 points_(points()),
234 faces_(polyMesh::faces()),
235 cells_(polyMesh::cells()),
236 edges_(edges),
237 faceEdges_(faceEdges),
238 oldPatchSizes_(faceSizes),
239 patchSizes_(faceSizes),
240 oldPatchStarts_(faceStarts),
241 patchStarts_(faceStarts),
242 oldEdgePatchSizes_(edgeSizes),
243 edgePatchSizes_(edgeSizes),
244 oldEdgePatchStarts_(edgeStarts),
245 edgePatchStarts_(edgeStarts),
246 oldPatchNMeshPoints_(faceStarts.size(), -1),
247 patchNMeshPoints_(faceStarts.size(), -1),
248 nOldPoints_(points_.size()),
249 nPoints_(points_.size()),
250 nOldEdges_(edges_.size()),
251 nEdges_(edges_.size()),
252 nOldFaces_(faces_.size()),
253 nFaces_(faces_.size()),
254 nOldCells_(cells_.size()),
255 nCells_(cells_.size()),
256 nOldInternalFaces_(faceStarts[0]),
257 nInternalFaces_(faceStarts[0]),
258 nOldInternalEdges_(edgeStarts[0]),
259 nInternalEdges_(edgeStarts[0]),
260 maxModifications_(mesh.maxModifications_),
261 statistics_(TOTAL_OP_TYPES, 0),
262 sliverThreshold_(mesh.sliverThreshold_),
263 slicePairs_(0),
264 maxTetsPerEdge_(mesh.maxTetsPerEdge_),
265 swapDeviation_(mesh.swapDeviation_),
266 allowTableResize_(mesh.allowTableResize_),
267 tetMetric_(mesh.tetMetric_)
268 {
269 // Initialize owner and neighbour
270 owner_.setSize(faces_.size(), -1);
271 neighbour_.setSize(faces_.size(), -1);
273 // Set owner and neighbour from polyMesh
274 const labelList& own = polyMesh::faceOwner();
275 const labelList& nei = polyMesh::faceNeighbour();
277 forAll(own, faceI)
278 {
279 owner_[faceI] = own[faceI];
280 }
282 forAll(nei, faceI)
283 {
284 neighbour_[faceI] = nei[faceI];
285 }
287 // Initialize the multiThreading environment.
288 // Force to single-threaded.
289 initializeThreadingEnvironment(1);
291 const polyBoundaryMesh& boundary = polyMesh::boundaryMesh();
293 // Add a boundary patches for polyMesh.
294 List<polyPatch*> patches(faceStarts.size());
296 forAll(patches, patchI)
297 {
298 // Set the pointer
299 patches[patchI] =
300 (
301 polyPatch::New
302 (
303 patchNames[patchI],
304 patchDict.subDict
305 (
306 patchNames[patchI]
307 ),
308 patchI,
309 boundary
310 ).ptr()
311 );
312 }
314 // Add patches, but don't calculate geometry, etc
315 fvMesh::addFvPatches(patches, false);
317 // Set sizes for the reverse maps
318 reversePointMap_.setSize(nPoints_, -7);
319 reverseEdgeMap_.setSize(nEdges_, -7);
320 reverseFaceMap_.setSize(nFaces_, -7);
321 reverseCellMap_.setSize(nCells_, -7);
323 // Now build edgeFaces and pointEdges information.
324 edgeFaces_ = invertManyToMany<labelList, labelList>(nEdges_, faceEdges_);
326 if (is3D())
327 {
328 pointEdges_ = invertManyToMany<edge, labelList>(nPoints_, edges_);
329 }
330 }
333 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
335 dynamicTopoFvMesh::~dynamicTopoFvMesh()
336 {
337 deleteDemandDrivenData(lduPtr_);
338 }
341 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
343 // Insert the specified cell to the mesh.
344 label dynamicTopoFvMesh::insertCell
345 (
346 const cell& newCell,
347 const scalar lengthScale,
348 const label zoneID
349 )
350 {
351 label newCellIndex = cells_.size();
353 if (debug > 2)
354 {
355 Pout<< "Inserting cell: "
356 << newCellIndex << ": "
357 << newCell << endl;
358 }
360 cells_.append(newCell);
362 if (edgeRefinement_)
363 {
364 lengthScale_.append(lengthScale);
365 }
367 // Add to the zone if necessary
368 if (zoneID >= 0)
369 {
370 addedCellZones_.insert(newCellIndex, zoneID);
371 }
373 nCells_++;
375 return newCellIndex;
376 }
379 // Remove the specified cell from the mesh
380 void dynamicTopoFvMesh::removeCell
381 (
382 const label cIndex
383 )
384 {
385 if (debug > 2)
386 {
387 Pout<< "Removing cell: "
388 << cIndex << ": "
389 << cells_[cIndex]
390 << endl;
391 }
393 cells_[cIndex].clear();
395 if (edgeRefinement_)
396 {
397 lengthScale_[cIndex] = -1.0;
398 }
400 // Update the number of cells, and the reverse cell map
401 nCells_--;
403 if (cIndex < nOldCells_)
404 {
405 reverseCellMap_[cIndex] = -1;
406 }
407 else
408 {
409 // Store this information for the reOrdering stage
410 deletedCells_.insert(cIndex);
411 }
413 // Check if this cell was added to a zone
414 Map<label>::iterator it = addedCellZones_.find(cIndex);
416 if (it != addedCellZones_.end())
417 {
418 addedCellZones_.erase(it);
419 }
421 // Check if the cell was added in the current morph, and delete
422 forAll(cellsFromPoints_, indexI)
423 {
424 if (cellsFromPoints_[indexI].index() == cIndex)
425 {
426 // Remove entry from the list
427 meshOps::removeIndex(indexI, cellsFromPoints_);
428 break;
429 }
430 }
432 forAll(cellsFromEdges_, indexI)
433 {
434 if (cellsFromEdges_[indexI].index() == cIndex)
435 {
436 // Remove entry from the list
437 meshOps::removeIndex(indexI, cellsFromEdges_);
438 break;
439 }
440 }
442 forAll(cellsFromFaces_, indexI)
443 {
444 if (cellsFromFaces_[indexI].index() == cIndex)
445 {
446 // Remove entry from the list
447 meshOps::removeIndex(indexI, cellsFromFaces_);
448 break;
449 }
450 }
452 forAll(cellsFromCells_, indexI)
453 {
454 if (cellsFromCells_[indexI].index() == cIndex)
455 {
456 // Remove entry from the list
457 meshOps::removeIndex(indexI, cellsFromCells_);
458 break;
459 }
460 }
461 }
464 // Utility method for face-insertion
465 label dynamicTopoFvMesh::insertFace
466 (
467 const label patch,
468 const face& newFace,
469 const label newOwner,
470 const label newNeighbour,
471 const labelList& newFaceEdges,
472 const label zoneID
473 )
474 {
475 // Append the specified face to each face-related list.
476 // Reordering is performed after all pending changes
477 // (flips, bisections, contractions, etc) have been made to the mesh
478 label newFaceIndex = faces_.size();
480 faces_.append(newFace);
481 owner_.append(newOwner);
482 neighbour_.append(newNeighbour);
483 faceEdges_.append(newFaceEdges);
485 if (debug > 2)
486 {
487 Pout<< "Inserting face: "
488 << newFaceIndex << ": "
489 << newFace
490 << " Owner: " << newOwner
491 << " Neighbour: " << newNeighbour
492 << " faceEdges: " << newFaceEdges;
494 const polyBoundaryMesh& boundary = boundaryMesh();
496 Pout<< " Patch: ";
498 if (patch == -1)
499 {
500 Pout<< "Internal" << endl;
501 }
502 else
503 if (patch < boundary.size())
504 {
505 Pout<< boundary[patch].name() << endl;
506 }
507 else
508 {
509 Pout<< " New patch: " << patch << endl;
510 }
511 }
513 // Keep track of added boundary faces in a separate hash-table
514 // This information will be required at the reordering stage
515 addedFacePatches_.insert(newFaceIndex, patch);
517 if (newNeighbour == -1)
518 {
519 // Modify patch information for this boundary face
520 patchSizes_[patch]++;
522 for (label i = (patch + 1); i < nPatches_; i++)
523 {
524 patchStarts_[i]++;
525 }
526 }
527 else
528 {
529 // Increment the number of internal faces,
530 // and subsequent patch-starts
531 nInternalFaces_++;
533 for (label i = 0; i < nPatches_; i++)
534 {
535 patchStarts_[i]++;
536 }
537 }
539 // Add to the zone if explicitly specified
540 if (zoneID >= 0)
541 {
542 addedFaceZones_.insert(newFaceIndex, zoneID);
543 }
544 else
545 {
546 // No zone was specified. Check if this
547 // face is added to a coupled patch associated
548 // with a faceZone.
549 forAll(patchCoupling_, patchI)
550 {
551 if (patchCoupling_(patchI))
552 {
553 if (patchI == patch)
554 {
555 if (patchCoupling_[patchI].masterFaceZone() > -1)
556 {
557 addedFaceZones_.insert
558 (
559 newFaceIndex,
560 patchCoupling_[patchI].masterFaceZone()
561 );
562 }
564 break;
565 }
567 if (patchCoupling_[patchI].map().slaveIndex() == patch)
568 {
569 if (patchCoupling_[patchI].slaveFaceZone() > -1)
570 {
571 addedFaceZones_.insert
572 (
573 newFaceIndex,
574 patchCoupling_[patchI].slaveFaceZone()
575 );
576 }
578 break;
579 }
580 }
581 }
582 }
584 // Increment the total face count
585 nFaces_++;
587 return newFaceIndex;
588 }
591 // Remove the specified face from the mesh
592 void dynamicTopoFvMesh::removeFace
593 (
594 const label fIndex
595 )
596 {
597 // Identify the patch for this face
598 label patch = whichPatch(fIndex);
600 if (debug > 2)
601 {
602 Pout<< "Removing face: "
603 << fIndex << ": "
604 << faces_[fIndex];
606 const polyBoundaryMesh& boundary = boundaryMesh();
608 Pout<< " Patch: ";
610 if (patch == -1)
611 {
612 Pout<< "Internal" << endl;
613 }
614 else
615 if (patch < boundary.size())
616 {
617 Pout<< boundary[patch].name() << endl;
618 }
619 else
620 {
621 Pout<< " New patch: " << patch << endl;
622 }
623 }
625 if (patch >= 0)
626 {
627 // Modify patch information for this boundary face
628 patchSizes_[patch]--;
630 for (label i = (patch + 1); i < nPatches_; i++)
631 {
632 patchStarts_[i]--;
633 }
634 }
635 else
636 {
637 // Decrement the internal face count, and subsequent patch-starts
638 nInternalFaces_--;
640 forAll(patchStarts_, patchI)
641 {
642 patchStarts_[patchI]--;
643 }
644 }
646 // Clear entities.
647 faces_[fIndex].clear();
648 owner_[fIndex] = -1;
649 neighbour_[fIndex] = -1;
650 faceEdges_[fIndex].clear();
652 // Entity won't be removed from the stack for efficiency
653 // It will be discarded on access instead.
655 // Update coupled face maps, if necessary.
656 forAll(patchCoupling_, patchI)
657 {
658 if (patchCoupling_(patchI))
659 {
660 const coupleMap& cMap = patchCoupling_[patchI].map();
662 cMap.removeMaster(coupleMap::FACE, fIndex);
663 cMap.removeSlave(coupleMap::FACE, fIndex);
664 }
665 }
667 // Update the reverse face-map, but only if this is a face that existed
668 // at time [n]. Added faces which are deleted during the topology change
669 // needn't be updated.
670 if (fIndex < nOldFaces_)
671 {
672 reverseFaceMap_[fIndex] = -1;
673 }
674 else
675 {
676 // Store this information for the reOrdering stage
677 deletedFaces_.insert(fIndex);
678 }
680 // Check and remove from the list of added face patches
681 Map<label>::iterator fpit = addedFacePatches_.find(fIndex);
683 if (fpit != addedFacePatches_.end())
684 {
685 addedFacePatches_.erase(fpit);
686 }
688 // Check if this face was added to a zone
689 Map<label>::iterator fzit = addedFaceZones_.find(fIndex);
691 if (fzit != addedFaceZones_.end())
692 {
693 addedFaceZones_.erase(fzit);
694 }
696 // Check if the face was added in the current morph, and delete
697 forAll(facesFromPoints_, indexI)
698 {
699 if (facesFromPoints_[indexI].index() == fIndex)
700 {
701 // Remove entry from the list
702 meshOps::removeIndex
703 (
704 indexI,
705 facesFromPoints_
706 );
708 break;
709 }
710 }
712 forAll(facesFromEdges_, indexI)
713 {
714 if (facesFromEdges_[indexI].index() == fIndex)
715 {
716 // Remove entry from the list
717 meshOps::removeIndex
718 (
719 indexI,
720 facesFromEdges_
721 );
723 break;
724 }
725 }
727 forAll(facesFromFaces_, indexI)
728 {
729 if (facesFromFaces_[indexI].index() == fIndex)
730 {
731 // Remove entry from the list
732 meshOps::removeIndex
733 (
734 indexI,
735 facesFromFaces_
736 );
738 break;
739 }
740 }
742 // Remove from the flipFaces list, if necessary
743 labelHashSet::iterator ffit = flipFaces_.find(fIndex);
745 if (ffit != flipFaces_.end())
746 {
747 flipFaces_.erase(ffit);
748 }
750 // Decrement the total face-count
751 nFaces_--;
752 }
755 // Insert the specified edge to the mesh
756 label dynamicTopoFvMesh::insertEdge
757 (
758 const label patch,
759 const edge& newEdge,
760 const labelList& edgeFaces
761 )
762 {
763 label newEdgeIndex = edges_.size();
765 edges_.append(newEdge);
766 edgeFaces_.append(edgeFaces);
768 if (debug > 2)
769 {
770 Pout<< "Inserting edge: "
771 << newEdgeIndex << ": "
772 << newEdge;
774 const polyBoundaryMesh& boundary = boundaryMesh();
776 Pout<< " Patch: ";
778 if (patch == -1)
779 {
780 Pout<< "Internal" << endl;
781 }
782 else
783 if (patch < boundary.size())
784 {
785 Pout<< boundary[patch].name() << endl;
786 }
787 else
788 {
789 Pout<< " New patch: " << patch << endl;
790 }
791 }
793 // Keep track of added edges in a separate hash-table
794 // This information will be required at the reordering stage
795 addedEdgePatches_.insert(newEdgeIndex,patch);
797 if (patch >= 0)
798 {
799 // Modify patch information for this boundary edge
800 edgePatchSizes_[patch]++;
802 for (label i = (patch + 1); i < nPatches_; i++)
803 {
804 edgePatchStarts_[i]++;
805 }
806 }
807 else
808 {
809 // Increment the number of internal edges, and subsequent patch-starts
810 nInternalEdges_++;
812 for (label i = 0; i < nPatches_; i++)
813 {
814 edgePatchStarts_[i]++;
815 }
816 }
818 // Size-up the pointEdges list
819 if (is3D())
820 {
821 meshOps::sizeUpList(newEdgeIndex, pointEdges_[newEdge[0]]);
822 meshOps::sizeUpList(newEdgeIndex, pointEdges_[newEdge[1]]);
823 }
825 // Increment the total edge count
826 nEdges_++;
828 return newEdgeIndex;
829 }
832 // Remove the specified edge from the mesh
833 void dynamicTopoFvMesh::removeEdge
834 (
835 const label eIndex
836 )
837 {
838 if (is3D())
839 {
840 const edge& rEdge = edges_[eIndex];
842 // Size-down the pointEdges list
843 if (pointEdges_[rEdge[0]].size())
844 {
845 meshOps::sizeDownList(eIndex, pointEdges_[rEdge[0]]);
846 }
848 if (pointEdges_[rEdge[1]].size())
849 {
850 meshOps::sizeDownList(eIndex, pointEdges_[rEdge[1]]);
851 }
853 // Entity won't be removed from the stack for efficiency
854 // It will be discarded on access instead.
856 // Update coupled edge maps, if necessary.
857 forAll(patchCoupling_, patchI)
858 {
859 if (patchCoupling_(patchI))
860 {
861 const coupleMap& cMap = patchCoupling_[patchI].map();
863 cMap.removeMaster(coupleMap::EDGE, eIndex);
864 cMap.removeSlave(coupleMap::EDGE, eIndex);
865 }
866 }
867 }
869 // Identify the patch for this edge
870 label patch = whichEdgePatch(eIndex);
872 if (debug > 2)
873 {
874 Pout<< "Removing edge: "
875 << eIndex << ": "
876 << edges_[eIndex];
878 const polyBoundaryMesh& boundary = boundaryMesh();
880 Pout<< " Patch: ";
882 if (patch == -1)
883 {
884 Pout<< "Internal" << endl;
885 }
886 else
887 if (patch < boundary.size())
888 {
889 Pout<< boundary[patch].name() << endl;
890 }
891 else
892 {
893 Pout<< " New patch: " << patch << endl;
894 }
895 }
897 // Invalidate edge and edgeFaces
898 edges_[eIndex] = edge(-1, -1);
899 edgeFaces_[eIndex].clear();
901 if (patch >= 0)
902 {
903 // Modify patch information for this boundary edge
904 edgePatchSizes_[patch]--;
906 for (label i = (patch + 1); i < nPatches_; i++)
907 {
908 edgePatchStarts_[i]--;
909 }
910 }
911 else
912 {
913 // Decrement the internal edge count, and subsequent patch-starts
914 nInternalEdges_--;
916 forAll(edgePatchStarts_, patchI)
917 {
918 edgePatchStarts_[patchI]--;
919 }
920 }
922 // Update reverse edge-map, but only if this is an edge that existed
923 // at time [n]. Added edges which are deleted during the topology change
924 // needn't be updated.
925 if (eIndex < nOldEdges_)
926 {
927 reverseEdgeMap_[eIndex] = -1;
928 }
929 else
930 {
931 // Store this information for the reOrdering stage
932 deletedEdges_.insert(eIndex);
933 }
935 // Check and remove from the list of added edge patches
936 Map<label>::iterator it = addedEdgePatches_.find(eIndex);
938 if (it != addedEdgePatches_.end())
939 {
940 addedEdgePatches_.erase(it);
941 }
943 // Decrement the total edge-count
944 nEdges_--;
945 }
948 // Insert the specified point to the mesh
949 label dynamicTopoFvMesh::insertPoint
950 (
951 const point& newPoint,
952 const point& oldPoint,
953 const labelList& mapPoints,
954 const label zoneID
955 )
956 {
957 // Add a new point to the end of the list
958 label newPointIndex = points_.size();
960 points_.append(newPoint);
961 oldPoints_.append(oldPoint);
963 if (debug > 2)
964 {
965 Pout<< "Inserting new point: "
966 << newPointIndex << ": "
967 << newPoint
968 << " and old point: "
969 << oldPoint
970 << " Mapped from: "
971 << mapPoints << endl;
972 }
974 // Make a pointsFromPoints entry
975 meshOps::sizeUpList
976 (
977 objectMap(newPointIndex, mapPoints),
978 pointsFromPoints_
979 );
981 // Add an empty entry to pointEdges as well.
982 // This entry can be sized-up appropriately at a later stage.
983 if (is3D())
984 {
985 pointEdges_.append(labelList(0));
986 }
988 // Add to the zone if necessary
989 if (zoneID >= 0)
990 {
991 addedPointZones_.insert(newPointIndex, zoneID);
992 }
994 nPoints_++;
996 return newPointIndex;
997 }
1000 // Remove the specified point from the mesh
1001 void dynamicTopoFvMesh::removePoint
1002 (
1003 const label pIndex
1004 )
1005 {
1006 if (debug > 2)
1007 {
1008 Pout<< "Removing point: " << pIndex
1009 << " :: new: " << points_[pIndex]
1010 << " :: old: " << oldPoints_[pIndex]
1011 << endl;
1012 }
1014 // Remove the point
1015 // (or just make sure that it's never used anywhere else)
1016 // points_[pIndex] = point();
1017 // oldPoints_[pIndex] = point();
1019 // Remove pointEdges as well
1020 if (is3D())
1021 {
1022 pointEdges_[pIndex].clear();
1023 }
1025 // Update the reverse point map
1026 if (pIndex < nOldPoints_)
1027 {
1028 reversePointMap_[pIndex] = -1;
1029 }
1030 else
1031 {
1032 deletedPoints_.insert(pIndex);
1033 }
1035 // Check if this point was added to a zone
1036 Map<label>::iterator pzit = addedPointZones_.find(pIndex);
1038 if (pzit != addedPointZones_.end())
1039 {
1040 addedPointZones_.erase(pzit);
1041 }
1043 // Update coupled point maps, if necessary.
1044 forAll(patchCoupling_, patchI)
1045 {
1046 if (patchCoupling_(patchI))
1047 {
1048 // Obtain references
1049 const coupleMap & cMap = patchCoupling_[patchI].map();
1051 Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
1052 Map<label>& rPointMap = cMap.reverseEntityMap(coupleMap::POINT);
1054 Map<label>::iterator pmit = pointMap.find(pIndex);
1056 if (pmit != pointMap.end())
1057 {
1058 // Erase the reverse map first
1059 rPointMap.erase(pmit());
1061 // Update pointMap
1062 pointMap.erase(pmit);
1063 }
1064 }
1065 }
1067 // Check if the point was added in the current morph, and delete
1068 forAll(pointsFromPoints_, indexI)
1069 {
1070 if (pointsFromPoints_[indexI].index() == pIndex)
1071 {
1072 // Remove entry from the list
1073 meshOps::removeIndex(indexI, pointsFromPoints_);
1074 break;
1075 }
1076 }
1078 // Decrement the total point-count
1079 nPoints_--;
1080 }
1083 // Utility method to build vertexHull for an edge [3D].
1084 // Assumes that edgeFaces information is consistent.
1085 void dynamicTopoFvMesh::buildVertexHull
1086 (
1087 const label eIndex,
1088 labelList& vertexHull,
1089 const label checkIndex
1090 ) const
1091 {
1092 bool found = false;
1093 label faceIndex = -1, cellIndex = -1;
1094 label otherPoint = -1, nextPoint = -1;
1095 label failMode = 0;
1097 // Obtain references
1098 const edge& edgeToCheck = edges_[eIndex];
1099 const labelList& eFaces = edgeFaces_[eIndex];
1101 // Re-size the list first
1102 vertexHull.clear();
1103 vertexHull.setSize(eFaces.size(), -1);
1105 if (whichEdgePatch(eIndex) == -1)
1106 {
1107 // Internal edge.
1108 // Pick the first face and start with that
1109 faceIndex = eFaces[0];
1110 }
1111 else
1112 {
1113 // Need to find a properly oriented start-face
1114 forAll(eFaces, faceI)
1115 {
1116 if (whichPatch(eFaces[faceI]) > -1)
1117 {
1118 meshOps::findIsolatedPoint
1119 (
1120 faces_[eFaces[faceI]],
1121 edgeToCheck,
1122 otherPoint,
1123 nextPoint
1124 );
1126 if (nextPoint == edgeToCheck[checkIndex])
1127 {
1128 faceIndex = eFaces[faceI];
1129 break;
1130 }
1131 }
1132 }
1133 }
1135 if (faceIndex == -1)
1136 {
1137 // Define failure-mode
1138 failMode = 1;
1139 }
1140 else
1141 {
1142 // Shuffle vertices to appear in CCW order
1143 forAll(vertexHull, indexI)
1144 {
1145 meshOps::findIsolatedPoint
1146 (
1147 faces_[faceIndex],
1148 edgeToCheck,
1149 otherPoint,
1150 nextPoint
1151 );
1153 // Add the isolated point
1154 vertexHull[indexI] = otherPoint;
1156 // Figure out how this edge is oriented.
1157 if (nextPoint == edgeToCheck[checkIndex])
1158 {
1159 // Counter-clockwise. Pick the owner.
1160 cellIndex = owner_[faceIndex];
1161 }
1162 else
1163 if (whichPatch(faceIndex) == -1)
1164 {
1165 // Clockwise. Pick the neighbour.
1166 cellIndex = neighbour_[faceIndex];
1167 }
1168 else
1169 if (indexI != (vertexHull.size() - 1))
1170 {
1171 // This could be a pinched manifold edge
1172 // (situation with more than two boundary faces)
1173 // Pick another (properly oriented) boundary face.
1174 faceIndex = -1;
1176 forAll(eFaces, faceI)
1177 {
1178 if (whichPatch(eFaces[faceI]) > -1)
1179 {
1180 meshOps::findIsolatedPoint
1181 (
1182 faces_[eFaces[faceI]],
1183 edgeToCheck,
1184 otherPoint,
1185 nextPoint
1186 );
1188 if
1189 (
1190 (nextPoint == edgeToCheck[checkIndex]) &&
1191 (findIndex(vertexHull, otherPoint) == -1)
1192 )
1193 {
1194 faceIndex = eFaces[faceI];
1195 break;
1196 }
1197 }
1198 }
1200 if (faceIndex == -1)
1201 {
1202 failMode = 2;
1203 break;
1204 }
1205 else
1206 {
1207 continue;
1208 }
1209 }
1210 else
1211 {
1212 // Looks like we've hit the last boundary face. Break out.
1213 break;
1214 }
1216 const cell& cellToCheck = cells_[cellIndex];
1218 found = false;
1220 // Assuming tet-cells,
1221 // Loop through edgeFaces and get the next face
1222 forAll(eFaces, faceI)
1223 {
1224 if
1225 (
1226 eFaces[faceI] != faceIndex
1227 && eFaces[faceI] == cellToCheck[0]
1228 )
1229 {
1230 faceIndex = cellToCheck[0];
1231 found = true; break;
1232 }
1234 if
1235 (
1236 eFaces[faceI] != faceIndex
1237 && eFaces[faceI] == cellToCheck[1]
1238 )
1239 {
1240 faceIndex = cellToCheck[1];
1241 found = true; break;
1242 }
1244 if
1245 (
1246 eFaces[faceI] != faceIndex
1247 && eFaces[faceI] == cellToCheck[2]
1248 )
1249 {
1250 faceIndex = cellToCheck[2];
1251 found = true; break;
1252 }
1254 if
1255 (
1256 eFaces[faceI] != faceIndex
1257 && eFaces[faceI] == cellToCheck[3]
1258 )
1259 {
1260 faceIndex = cellToCheck[3];
1261 found = true; break;
1262 }
1263 }
1265 if (!found)
1266 {
1267 failMode = 3;
1268 break;
1269 }
1270 }
1271 }
1273 if (debug > 2 && !failMode)
1274 {
1275 // Check for invalid indices
1276 if (findIndex(vertexHull, -1) > -1)
1277 {
1278 failMode = 4;
1279 }
1281 // Check for duplicate labels
1282 if (!failMode)
1283 {
1284 labelHashSet uniquePoints;
1286 forAll(vertexHull, pointI)
1287 {
1288 bool inserted = uniquePoints.insert(vertexHull[pointI]);
1290 if (!inserted)
1291 {
1292 Pout<< " vertexHull for edge: "
1293 << eIndex << "::" << edgeToCheck
1294 << " contains identical vertex labels: "
1295 << vertexHull << endl;
1297 failMode = 5;
1298 }
1299 }
1300 }
1301 }
1303 if (failMode)
1304 {
1305 // Prepare edgeCells
1306 dynamicLabelList eCells(10);
1308 forAll(eFaces, faceI)
1309 {
1310 label own = owner_[eFaces[faceI]];
1311 label nei = neighbour_[eFaces[faceI]];
1313 if (findIndex(eCells, own) == -1)
1314 {
1315 eCells.append(own);
1316 }
1318 if (nei > -1)
1319 {
1320 if (findIndex(eCells, nei) == -1)
1321 {
1322 eCells.append(nei);
1323 }
1324 }
1325 }
1327 // Write out for post-processing
1328 label eI = eIndex, epI = whichEdgePatch(eIndex);
1329 word epName(epI < 0 ? "Internal" : boundaryMesh()[epI].name());
1331 writeVTK("vRingEdge_" + Foam::name(eI), eIndex, 1, false, true);
1332 writeVTK("vRingEdgeFaces_" + Foam::name(eI), eFaces, 2, false, true);
1333 writeVTK("vRingEdgeCells_" + Foam::name(eI), eCells, 3, false, true);
1335 Pout<< "edgeFaces: " << endl;
1337 forAll(eFaces, faceI)
1338 {
1339 label fpI = whichPatch(eFaces[faceI]);
1340 word fpName(fpI < 0 ? "Internal" : boundaryMesh()[fpI].name());
1342 Pout<< " Face: " << eFaces[faceI]
1343 << ":: " << faces_[eFaces[faceI]]
1344 << " Owner: " << owner_[eFaces[faceI]]
1345 << " Neighbour: " << neighbour_[eFaces[faceI]]
1346 << " Patch: " << fpName
1347 << endl;
1348 }
1350 FatalErrorIn
1351 (
1352 "\n"
1353 "void dynamicTopoFvMesh::buildVertexHull\n"
1354 "(\n"
1355 " const label eIndex,\n"
1356 " labelList& vertexHull,\n"
1357 " const label checkIndex\n"
1358 ")\n"
1359 )
1360 << " Failed to determine a vertex ring. " << nl
1361 << " Failure mode: " << failMode << nl
1362 << " Edge: " << eIndex << ":: " << edgeToCheck << nl
1363 << " edgeFaces: " << eFaces << nl
1364 << " Patch: " << epName << nl
1365 << " Current vertexHull: " << vertexHull
1366 << abort(FatalError);
1367 }
1368 }
1371 void dynamicTopoFvMesh::handleMeshSlicing()
1372 {
1373 if (slicePairs_.empty())
1374 {
1375 return;
1376 }
1378 if (lengthEstimator().holdOff())
1379 {
1380 // Clear out data.
1381 slicePairs_.clear();
1382 lengthEstimator().clearBoxes();
1384 // Hold-off mesh slicing for a few time-steps.
1385 lengthEstimator().decrementHoldOff();
1387 return;
1388 }
1390 Info<< "Slicing Mesh...";
1392 // Loop through candidates and weed-out invalid points
1393 forAll(slicePairs_, pairI)
1394 {
1395 const labelPair& pairToCheck = slicePairs_[pairI];
1397 bool available = true;
1399 if (is2D())
1400 {
1401 forAll(pairToCheck, indexI)
1402 {
1403 if (deletedFaces_.found(pairToCheck[indexI]))
1404 {
1405 available = false;
1406 break;
1407 }
1409 if (pairToCheck[indexI] < nOldFaces_)
1410 {
1411 if (reverseFaceMap_[pairToCheck[indexI]] == -1)
1412 {
1413 available = false;
1414 break;
1415 }
1416 }
1417 }
1418 }
1419 else
1420 {
1421 forAll(pairToCheck, indexI)
1422 {
1423 if (deletedPoints_.found(pairToCheck[indexI]))
1424 {
1425 available = false;
1426 break;
1427 }
1429 if (pairToCheck[indexI] < nOldPoints_)
1430 {
1431 if (reversePointMap_[pairToCheck[indexI]] == -1)
1432 {
1433 available = false;
1434 break;
1435 }
1436 }
1437 }
1438 }
1440 if (available)
1441 {
1442 // Slice the mesh at this point.
1443 sliceMesh(pairToCheck);
1444 }
1445 }
1447 Info<< "Done." << endl;
1449 // Clear out data.
1450 slicePairs_.clear();
1451 lengthEstimator().clearBoxes();
1453 // Set the sliceHoldOff value
1454 lengthEstimator().setHoldOff(50);
1455 }
1458 // Handle layer addition / removal events
1459 void dynamicTopoFvMesh::handleLayerAdditionRemoval()
1460 {
1461 const dictionary& meshSubDict = dict_.subDict("dynamicTopoFvMesh");
1463 if (meshSubDict.found("layering") || mandatory_)
1464 {
1465 // Bail out if no entries are present
1466 if (meshSubDict.subDict("layering").empty())
1467 {
1468 return;
1469 }
1470 }
1471 else
1472 {
1473 // No dictionary present.
1474 return;
1475 }
1477 // Fetch layering patches
1478 const polyBoundaryMesh& boundary = boundaryMesh();
1479 const dictionary& layeringDict = meshSubDict.subDict("layering");
1481 wordList layeringPatches = layeringDict.toc();
1483 forAll(layeringPatches, wordI)
1484 {
1485 word pName = layeringPatches[wordI];
1486 label patchID = boundary.findPatchID(pName);
1488 // If this patch has no faces, move on
1489 if (boundary[patchID].empty())
1490 {
1491 continue;
1492 }
1494 // Use first face to determine layer thickness
1495 scalar magSf = mag(boundary[patchID].faceAreas()[0]);
1496 scalar Vc = cellVolumes()[boundary[patchID].faceCells()[0]];
1498 // Define thickness
1499 scalar thickness = (Vc / (magSf + VSMALL));
1501 // Fetch min / max thickness
1502 scalar minThickness =
1503 (
1504 readScalar(layeringDict.subDict(pName).lookup("minThickness"))
1505 );
1507 scalar maxThickness =
1508 (
1509 readScalar(layeringDict.subDict(pName).lookup("maxThickness"))
1510 );
1512 if (debug)
1513 {
1514 Pout<< " Patch: " << pName << nl
1515 << " Face area: " << magSf << nl
1516 << " Cell volume: " << Vc << nl
1517 << " Layer thickness: " << thickness << nl
1518 << " Min thickness: " << minThickness << nl
1519 << " Max thickness: " << maxThickness << nl
1520 << endl;
1521 }
1523 if (thickness > maxThickness)
1524 {
1525 // Add cell layer above patch
1526 addCellLayer(patchID);
1527 }
1528 else
1529 if (thickness < minThickness)
1530 {
1531 // Remove cell layer above patch
1532 removeCellLayer(patchID);
1533 }
1534 }
1535 }
1538 // Test an edge / face for proximity with other faces on proximity patches
1539 // and return the scalar distance to an oppositely-oriented face.
1540 scalar dynamicTopoFvMesh::testProximity
1541 (
1542 const label index,
1543 label& proximityFace
1544 ) const
1545 {
1546 scalar proxDistance = GREAT, testStep = 0.0;
1547 vector gCentre = vector::zero, gNormal = vector::zero;
1549 if (is2D())
1550 {
1551 // Obtain the face-normal.
1552 gNormal = faces_[index].normal(points_);
1554 // Obtain the face centre.
1555 gCentre = faces_[index].centre(points_);
1557 // Fetch the edge
1558 const edge& edgeToCheck = edges_[getTriBoundaryEdge(index)];
1560 // Calculate a test step-size
1561 testStep =
1562 (
1563 linePointRef
1564 (
1565 points_[edgeToCheck.start()],
1566 points_[edgeToCheck.end()]
1567 ).mag()
1568 );
1569 }
1570 else
1571 {
1572 const edge& edgeToCheck = edges_[index];
1573 const labelList& eFaces = edgeFaces_[index];
1575 linePointRef lpr
1576 (
1577 points_[edgeToCheck.start()],
1578 points_[edgeToCheck.end()]
1579 );
1581 // Obtain the edge centre.
1582 gCentre = lpr.centre();
1584 // Obtain the edge-normal
1585 forAll(eFaces, faceI)
1586 {
1587 if (neighbour_[eFaces[faceI]] == -1)
1588 {
1589 // Obtain the normal.
1590 gNormal += faces_[eFaces[faceI]].normal(points_);
1591 }
1592 }
1594 // Calculate a test step-size
1595 testStep = lpr.mag();
1596 }
1598 // Normalize
1599 gNormal /= (mag(gNormal) + VSMALL);
1601 // Test for proximity, and mark slice-pairs
1602 // is the distance is below threshold.
1603 if
1604 (
1605 lengthEstimator().testProximity
1606 (
1607 gCentre,
1608 gNormal,
1609 testStep,
1610 proximityFace,
1611 proxDistance
1612 )
1613 )
1614 {
1615 labelPair proxPoints(-1, -1);
1616 bool foundPoint = false;
1618 if (is2D())
1619 {
1620 proxPoints.first() = index;
1621 proxPoints.second() = proximityFace;
1623 if
1624 (
1625 (faces_[index].size() == 4) &&
1626 (polyMesh::faces()[proximityFace].size() == 4)
1627 )
1628 {
1629 foundPoint = true;
1630 }
1631 }
1632 else
1633 {
1634 const edge& thisEdge = edges_[index];
1635 const face& proxFace = polyMesh::faces()[proximityFace];
1637 // Check if any points on this face are still around.
1638 // If yes, mark one of them as the end point
1639 // for Dijkstra's algorithm. The start point will be a point
1640 // on this edge.
1641 proxPoints.first() = thisEdge[0];
1643 forAll(proxFace, pointI)
1644 {
1645 if (reversePointMap_[proxFace[pointI]] != -1)
1646 {
1647 proxPoints.second() = proxFace[pointI];
1648 foundPoint = true;
1649 break;
1650 }
1651 }
1652 }
1654 if (foundPoint)
1655 {
1656 // Lock the edge mutex
1657 entityMutex_[1].lock();
1659 label newSize = slicePairs_.size() + 1;
1661 // Const-cast slicePairs for resize
1662 List<labelPair>& sP = const_cast<List<labelPair>&>(slicePairs_);
1664 // Add this entry as a candidate for mesh slicing.
1665 sP.setSize(newSize, proxPoints);
1667 // Unlock the edge mutex
1668 entityMutex_[1].unlock();
1669 }
1670 }
1672 return proxDistance;
1673 }
1676 // Calculate the edge length-scale for the mesh
1677 void dynamicTopoFvMesh::calculateLengthScale(bool dump)
1678 {
1679 if (!edgeRefinement_)
1680 {
1681 return;
1682 }
1684 Switch dumpLengthScale(false);
1686 const dictionary& meshDict = dict_.subDict("dynamicTopoFvMesh");
1688 if (meshDict.found("dumpLengthScale") || mandatory_)
1689 {
1690 dumpLengthScale = readBool(meshDict.lookup("dumpLengthScale"));
1691 }
1693 autoPtr<volScalarField> lsfPtr(NULL);
1695 if (dumpLengthScale && time().outputTime() && dump)
1696 {
1697 lsfPtr.set
1698 (
1699 new volScalarField
1700 (
1701 IOobject
1702 (
1703 "lengthScale",
1704 time().timeName(),
1705 *this,
1706 IOobject::NO_READ,
1707 IOobject::NO_WRITE,
1708 false
1709 ),
1710 *this,
1711 dimensionedScalar("scalar", dimLength, 0)
1712 )
1713 );
1714 }
1716 // Bail-out if a dumping was not requested in dictionary.
1717 if (dump && !dumpLengthScale)
1718 {
1719 return;
1720 }
1722 // Size the field and calculate length-scale
1723 lengthScale_.setSize(nCells_, 0.0);
1725 lengthEstimator().calculateLengthScale(lengthScale_);
1727 // Check if length-scale is to be dumped to disk.
1728 if (dumpLengthScale && time().outputTime() && dump)
1729 {
1730 // Obtain length-scale values from the mesh
1731 lsfPtr->internalField() = lengthScale_;
1733 lsfPtr->write();
1734 }
1735 }
1738 // Read optional dictionary parameters
1739 void dynamicTopoFvMesh::readOptionalParameters(bool reRead)
1740 {
1741 // Read from disk
1742 dict_.readIfModified();
1744 const dictionary& meshSubDict = dict_.subDict("dynamicTopoFvMesh");
1746 // Enable/disable run-time debug level
1747 if (meshSubDict.found("debug") || mandatory_)
1748 {
1749 // Set an attachment point for debugging
1750 if
1751 (
1752 !reRead &&
1753 meshSubDict.found("parDebugAttach") &&
1754 readBool(meshSubDict.lookup("parDebugAttach"))
1755 )
1756 {
1757 label wait;
1759 if (Pstream::master())
1760 {
1761 Info<< nl << "Waiting for debugger attachment..." << endl;
1763 cin >> wait;
1765 if (Pstream::parRun())
1766 {
1767 for (label i = 1; i < Pstream::nProcs(); i++)
1768 {
1769 meshOps::pWrite(i, wait);
1770 }
1771 }
1772 }
1773 else
1774 {
1775 meshOps::pRead(Pstream::masterNo(), wait);
1776 }
1777 }
1779 // Look for a processor-specific debug flag
1780 if
1781 (
1782 Pstream::parRun() &&
1783 meshSubDict.found("parDebugProcs")
1784 )
1785 {
1786 if (!reRead)
1787 {
1788 labelList procs(meshSubDict.lookup("parDebugProcs"));
1789 label pdebug = readLabel(meshSubDict.lookup("debug"));
1791 if (findIndex(procs, Pstream::myProcNo()) > -1)
1792 {
1793 debug = pdebug;
1794 }
1795 else
1796 if (pdebug)
1797 {
1798 // Minimal debug information,
1799 // so that synchronizations are successful
1800 debug = 1;
1801 }
1802 }
1803 }
1804 else
1805 {
1806 debug = readLabel(meshSubDict.lookup("debug"));
1807 }
1808 }
1809 else
1810 {
1811 debug = 0;
1812 }
1814 // Set debug option for underlying classes as well.
1815 if (debug > 2)
1816 {
1817 fvMesh::debug = true;
1818 polyMesh::debug = true;
1819 }
1821 // Re-read edge-refinement options, if necessary
1822 if (edgeRefinement_ && reRead)
1823 {
1824 lengthEstimator().readRefinementOptions
1825 (
1826 meshSubDict.subDict("refinementOptions"),
1827 true,
1828 mandatory_
1829 );
1830 }
1832 // Read re-mesh interval
1833 if (meshSubDict.found("interval") || mandatory_)
1834 {
1835 interval_ = readLabel(meshSubDict.lookup("interval"));
1836 }
1837 else
1838 {
1839 interval_ = 1;
1840 }
1842 // Check if an external mesh-motion solver is used
1843 if (meshSubDict.found("loadMotionSolver") || mandatory_)
1844 {
1845 loadMotionSolver_.readIfPresent("loadMotionSolver", meshSubDict);
1846 }
1848 // Update bandwidth reduction switch
1849 if (meshSubDict.found("bandwidthReduction") || mandatory_)
1850 {
1851 bandWidthReduction_.readIfPresent("bandwidthReduction", meshSubDict);
1852 }
1854 // Update threshold for sliver cells
1855 if (meshSubDict.found("sliverThreshold") || mandatory_)
1856 {
1857 sliverThreshold_ = readScalar(meshSubDict.lookup("sliverThreshold"));
1859 if (sliverThreshold_ > 1.0 || sliverThreshold_ < 0.0)
1860 {
1861 FatalErrorIn("void dynamicTopoFvMesh::readOptionalParameters()")
1862 << " Sliver threshold out of range [0..1]"
1863 << abort(FatalError);
1864 }
1865 }
1867 // Update limit for max number of bisections / collapses
1868 if (meshSubDict.found("maxModifications") || mandatory_)
1869 {
1870 maxModifications_ = readLabel(meshSubDict.lookup("maxModifications"));
1871 }
1873 // Update limit for swap on curved surfaces
1874 if (meshSubDict.found("swapDeviation") || mandatory_)
1875 {
1876 swapDeviation_ = readScalar(meshSubDict.lookup("swapDeviation"));
1878 if (swapDeviation_ > 1.0 || swapDeviation_ < 0.0)
1879 {
1880 FatalErrorIn("void dynamicTopoFvMesh::readOptionalParameters()")
1881 << " Swap deviation out of range [0..1]"
1882 << abort(FatalError);
1883 }
1884 }
1886 // For tetrahedral meshes...
1887 if (is3D())
1888 {
1889 // Check if swapping is to be avoided on any patches
1890 if (meshSubDict.found("noSwapPatches") || mandatory_)
1891 {
1892 wordList noSwapPatches =
1893 (
1894 meshSubDict.subDict("noSwapPatches").toc()
1895 );
1897 noSwapPatchIDs_.setSize(noSwapPatches.size());
1899 label indexI = 0;
1901 forAll(noSwapPatches, wordI)
1902 {
1903 const word& patchName = noSwapPatches[wordI];
1905 noSwapPatchIDs_[indexI++] =
1906 (
1907 boundaryMesh().findPatchID(patchName)
1908 );
1909 }
1910 }
1912 // Check if a limit has been imposed on maxTetsPerEdge
1913 if (meshSubDict.found("maxTetsPerEdge") || mandatory_)
1914 {
1915 maxTetsPerEdge_ = readLabel(meshSubDict.lookup("maxTetsPerEdge"));
1916 }
1917 else
1918 {
1919 maxTetsPerEdge_ = 7;
1920 }
1922 // Check if programming tables can be resized at runtime
1923 if (meshSubDict.found("allowTableResize") || mandatory_)
1924 {
1925 allowTableResize_ =
1926 (
1927 readBool(meshSubDict.lookup("allowTableResize"))
1928 );
1929 }
1930 else
1931 {
1932 allowTableResize_ = false;
1933 }
1934 }
1936 // Check for load-balancing in parallel
1937 if (reRead && (meshSubDict.found("loadBalancing") || mandatory_))
1938 {
1939 // Fetch non-const reference
1940 dictionary& balanceDict =
1941 (
1942 dict_.subDict("dynamicTopoFvMesh").subDict("loadBalancing")
1943 );
1945 // Execute balancing
1946 executeLoadBalancing(balanceDict);
1947 }
1948 }
1951 // Initialize edge related connectivity lists
1952 void dynamicTopoFvMesh::initEdges()
1953 {
1954 // Initialize eMesh, and copy to local lists
1955 eMeshPtr_.set(new eMesh(*this));
1957 // Obtain information
1958 nEdges_ = eMeshPtr_->nEdges();
1959 nInternalEdges_ = eMeshPtr_->nInternalEdges();
1960 edgePatchSizes_ = eMeshPtr_->boundary().patchSizes();
1961 edgePatchStarts_ = eMeshPtr_->boundary().patchStarts();
1963 // Set old edge information
1964 nOldEdges_ = nEdges_;
1965 nOldInternalEdges_ = nInternalEdges_;
1966 oldEdgePatchSizes_ = edgePatchSizes_;
1967 oldEdgePatchStarts_ = edgePatchStarts_;
1969 // Set local lists with edge connectivity information
1970 edges_ = eMeshPtr_->edges();
1971 edgeFaces_ = eMeshPtr_->edgeFaces();
1972 faceEdges_ = eMeshPtr_->faceEdges();
1974 if (is3D())
1975 {
1976 // Invert edges to obtain pointEdges
1977 pointEdges_ = invertManyToMany<edge, labelList>(nPoints_, edges_);
1978 }
1980 // Clear out unwanted eMesh connectivity
1981 eMeshPtr_->clearOut();
1983 // Read coupled patch information from dictionary.
1984 readCoupledPatches();
1985 }
1988 // Load the mesh-quality metric from the library
1989 void dynamicTopoFvMesh::loadMetric()
1990 {
1991 if (is2D())
1992 {
1993 return;
1994 }
1996 // Specify the dictionary we would be looking in...
1997 const dictionary& meshDict = dict_.subDict("dynamicTopoFvMesh");
1999 // Select an appropriate metric
2000 tetMetric_ = tetMetric::New(meshDict, meshDict.lookup("tetMetric"));
2001 }
2004 // Load the mesh-motion solver
2005 void dynamicTopoFvMesh::loadMotionSolver()
2006 {
2007 if (motionSolver_.valid())
2008 {
2009 FatalErrorIn
2010 (
2011 "void dynamicTopoFvMesh::loadMotionSolver() "
2012 ) << nl << " Motion solver already loaded. "
2013 << abort(FatalError);
2014 }
2015 else
2016 if (dict_.found("solver") && loadMotionSolver_)
2017 {
2018 motionSolver_ = motionSolver::New(*this);
2019 }
2020 }
2023 // Load the field mapper
2024 void dynamicTopoFvMesh::loadFieldMapper()
2025 {
2026 if (mapper_.valid())
2027 {
2028 FatalErrorIn
2029 (
2030 "void dynamicTopoFvMesh::loadFieldMapper() "
2031 ) << nl << " Field mapper already loaded. "
2032 << abort(FatalError);
2033 }
2034 else
2035 {
2036 mapper_.set(new topoMapper(*this, dict_));
2037 }
2038 }
2041 // Load the length scale estimator
2042 void dynamicTopoFvMesh::loadLengthScaleEstimator()
2043 {
2044 if (edgeRefinement_)
2045 {
2046 if (lengthEstimator_.valid())
2047 {
2048 FatalErrorIn
2049 (
2050 "void dynamicTopoFvMesh::loadLengthScaleEstimator() "
2051 ) << nl << " Length scale estimator already loaded. "
2052 << abort(FatalError);
2053 }
2054 else
2055 {
2056 lengthEstimator_.set
2057 (
2058 new lengthScaleEstimator(*this)
2059 );
2060 }
2062 // Read options
2063 lengthEstimator().readRefinementOptions
2064 (
2065 dict_.subDict("dynamicTopoFvMesh").subDict("refinementOptions"),
2066 false,
2067 mandatory_
2068 );
2070 // Set coupled patch options, if available
2071 if (dict_.found("coupledPatches") || mandatory_)
2072 {
2073 lengthEstimator().setCoupledPatches
2074 (
2075 dict_.subDict("coupledPatches")
2076 );
2077 }
2078 }
2079 }
2082 // Initialize the threading environment.
2083 // - Provides an override option to avoid reading from the dictionary.
2084 void dynamicTopoFvMesh::initializeThreadingEnvironment
2085 (
2086 const label specThreads
2087 )
2088 {
2089 // Initialize an IOobject for the IOmultiThreader object
2090 IOobject io
2091 (
2092 "threader",
2093 this->time().timeName(),
2094 (*this),
2095 IOobject::NO_READ,
2096 IOobject::NO_WRITE,
2097 true
2098 );
2100 if (specThreads > 0)
2101 {
2102 threader_.set(new IOmultiThreader(io, specThreads));
2103 }
2104 else
2105 {
2106 if (dict_.subDict("dynamicTopoFvMesh").found("threads") || mandatory_)
2107 {
2108 threader_.set
2109 (
2110 new IOmultiThreader
2111 (
2112 io,
2113 readLabel
2114 (
2115 dict_.subDict("dynamicTopoFvMesh").lookup("threads")
2116 )
2117 )
2118 );
2119 }
2120 else
2121 {
2122 threader_.set(new IOmultiThreader(io, 1));
2123 }
2124 }
2126 // Get the number of threads and allocate threadHandlers
2127 label nThreads = threader_->getNumThreads();
2129 if (nThreads == 1)
2130 {
2131 handlerPtr_.setSize(1);
2133 handlerPtr_.set
2134 (
2135 0,
2136 new meshHandler(*this, threader())
2137 );
2139 handlerPtr_[0].setMaster();
2141 // Size the stacks
2142 entityStack_.setSize(1);
2143 }
2144 else
2145 {
2146 // Index '0' is master, rest are slaves
2147 handlerPtr_.setSize(nThreads + 1);
2149 // Size the stacks
2150 entityStack_.setSize(nThreads + 1);
2152 forAll(handlerPtr_, threadI)
2153 {
2154 handlerPtr_.set
2155 (
2156 threadI,
2157 new meshHandler(*this, threader())
2158 );
2160 if (threadI == 0)
2161 {
2162 handlerPtr_[0].setMaster();
2163 }
2164 else
2165 {
2166 handlerPtr_[threadI].setID(threader_->getID(threadI - 1));
2167 handlerPtr_[threadI].setSlave();
2168 }
2169 }
2170 }
2171 }
2174 // 2D Edge-swapping engine
2175 void dynamicTopoFvMesh::swap2DEdges(void *argument)
2176 {
2177 // Recast the argument
2178 meshHandler *thread = static_cast<meshHandler*>(argument);
2180 if (thread->slave())
2181 {
2182 thread->sendSignal(meshHandler::START);
2183 }
2185 dynamicTopoFvMesh& mesh = thread->reference();
2187 // Figure out which thread this is...
2188 label tIndex = mesh.self();
2190 // Set the timer
2191 clockTime sTimer;
2193 bool reported = false;
2194 label stackSize = mesh.stack(tIndex).size();
2195 scalar interval = mesh.reportInterval(), oIndex = 0.0, nIndex = 0.0;
2197 oIndex = ::floor(sTimer.elapsedTime() / interval);
2199 // Pick items off the stack
2200 while (!mesh.stack(tIndex).empty())
2201 {
2202 // Report progress
2203 if (thread->master())
2204 {
2205 // Update the index, if its changed
2206 nIndex = ::floor(sTimer.elapsedTime() / interval);
2208 if ((nIndex - oIndex) > VSMALL)
2209 {
2210 oIndex = nIndex;
2212 scalar percent =
2213 (
2214 100.0 -
2215 (
2216 (100.0 * mesh.stack(tIndex).size())
2217 / (stackSize + VSMALL)
2218 )
2219 );
2221 Info<< " Swap Progress: " << percent << "% :"
2222 << " Total: " << mesh.status(TOTAL_SWAPS)
2223 << " \r"
2224 << flush;
2226 reported = true;
2227 }
2228 }
2230 // Retrieve the index for this face
2231 label fIndex = mesh.stack(tIndex).pop();
2233 // Perform a Delaunay test and check if a flip is necesary.
2234 bool failed = mesh.testDelaunay(fIndex);
2236 if (failed)
2237 {
2238 if (thread->master())
2239 {
2240 // Swap this face.
2241 mesh.swapQuadFace(fIndex);
2242 }
2243 else
2244 {
2245 // Push this on to the master stack
2246 mesh.stack(0).push(fIndex);
2247 }
2248 }
2249 }
2251 if (thread->slave())
2252 {
2253 thread->sendSignal(meshHandler::STOP);
2254 }
2256 if (reported)
2257 {
2258 Info<< " Swap Progress: 100% :"
2259 << " Total: " << mesh.status(TOTAL_SWAPS)
2260 << " \r"
2261 << endl;
2262 }
2263 }
2266 // 3D Edge-swapping engine
2267 void dynamicTopoFvMesh::swap3DEdges
2268 (
2269 void *argument
2270 )
2271 {
2272 // Recast the argument
2273 meshHandler *thread = static_cast<meshHandler*>(argument);
2275 if (thread->slave())
2276 {
2277 thread->sendSignal(meshHandler::START);
2278 }
2280 dynamicTopoFvMesh& mesh = thread->reference();
2282 // Figure out which thread this is...
2283 label tIndex = mesh.self();
2285 // Dynamic programming variables
2286 labelList m;
2287 PtrList<scalarListList> Q;
2288 PtrList<labelListList> K, triangulations;
2290 // Hull vertices information
2291 labelList hullV;
2293 // Allocate dynamic programming tables
2294 mesh.initTables(m, Q, K, triangulations);
2296 // Set the timer
2297 clockTime sTimer;
2299 bool reported = false;
2300 label stackSize = mesh.stack(tIndex).size();
2301 scalar interval = mesh.reportInterval(), oIndex = 0.0, nIndex = 0.0;
2303 oIndex = ::floor(sTimer.elapsedTime() / interval);
2305 // Pick edges off the stack
2306 while (!mesh.stack(tIndex).empty())
2307 {
2308 // Report progress
2309 if (thread->master())
2310 {
2311 // Update the index, if its changed
2312 nIndex = ::floor(sTimer.elapsedTime() / interval);
2314 if ((nIndex - oIndex) > VSMALL)
2315 {
2316 oIndex = nIndex;
2318 scalar percent =
2319 (
2320 100.0 -
2321 (
2322 (100.0 * mesh.stack(tIndex).size())
2323 / (stackSize + VSMALL)
2324 )
2325 );
2327 Info<< " Swap Progress: " << percent << "% :"
2328 << " Surface: " << mesh.status(SURFACE_SWAPS)
2329 << ", Total: " << mesh.status(TOTAL_SWAPS)
2330 << " \r"
2331 << flush;
2333 reported = true;
2334 }
2335 }
2337 // Retrieve an edge from the stack
2338 label eIndex = mesh.stack(tIndex).pop();
2340 // Compute the minimum quality of cells around this edge
2341 scalar minQuality = mesh.computeMinQuality(eIndex, hullV);
2343 // Check if this edge is on a bounding curve
2344 // (Override purity check for processor edges)
2345 if (mesh.checkBoundingCurve(eIndex, true))
2346 {
2347 continue;
2348 }
2350 // Fill the dynamic programming tables
2351 if (mesh.fillTables(eIndex, minQuality, m, hullV, Q, K, triangulations))
2352 {
2353 // Check if edge-swapping is required.
2354 if (mesh.checkQuality(eIndex, m, Q, minQuality))
2355 {
2356 if (thread->master())
2357 {
2358 // Remove this edge according to the swap sequence
2359 mesh.removeEdgeFlips
2360 (
2361 eIndex,
2362 minQuality,
2363 hullV,
2364 Q,
2365 K,
2366 triangulations
2367 );
2368 }
2369 else
2370 {
2371 // Push this on to the master stack
2372 mesh.stack(0).push(eIndex);
2373 }
2374 }
2375 }
2376 }
2378 if (thread->slave())
2379 {
2380 thread->sendSignal(meshHandler::STOP);
2381 }
2383 if (reported)
2384 {
2385 Info<< " Swap Progress: 100% :"
2386 << " Surface: " << mesh.status(SURFACE_SWAPS)
2387 << ", Total: " << mesh.status(TOTAL_SWAPS)
2388 << " \r"
2389 << endl;
2390 }
2391 }
2394 // Edge refinement engine
2395 void dynamicTopoFvMesh::edgeRefinementEngine
2396 (
2397 void *argument
2398 )
2399 {
2400 // Loop through all edges and bisect/collapse by the criterion:
2401 // Bisect when edge-length > ratioMax_*lengthScale
2402 // Collapse when edge-length < ratioMin_*lengthScale
2404 // Recast the argument
2405 meshHandler *thread = static_cast<meshHandler*>(argument);
2407 if (thread->slave())
2408 {
2409 thread->sendSignal(meshHandler::START);
2410 }
2412 dynamicTopoFvMesh& mesh = thread->reference();
2414 // Figure out which thread this is...
2415 label tIndex = mesh.self();
2417 // Set the timer
2418 clockTime sTimer;
2420 bool reported = false;
2421 label stackSize = mesh.stack(tIndex).size();
2422 scalar interval = mesh.reportInterval(), oIndex = 0.0, nIndex = 0.0;
2424 oIndex = ::floor(sTimer.elapsedTime() / interval);
2426 while (!mesh.stack(tIndex).empty())
2427 {
2428 // Update the index, if its changed
2429 // Report progress
2430 if (thread->master())
2431 {
2432 nIndex = ::floor(sTimer.elapsedTime() / interval);
2434 if ((nIndex - oIndex) > VSMALL)
2435 {
2436 oIndex = nIndex;
2438 scalar percent =
2439 (
2440 100.0 -
2441 (
2442 (100.0 * mesh.stack(tIndex).size())
2443 / (stackSize + VSMALL)
2444 )
2445 );
2447 Info<< " Refinement Progress: " << percent << "% :"
2448 << " Bisections: " << mesh.status(TOTAL_BISECTIONS)
2449 << ", Collapses: " << mesh.status(TOTAL_COLLAPSES)
2450 << ", Total: " << mesh.status(TOTAL_MODIFICATIONS)
2451 << " \r"
2452 << flush;
2454 reported = true;
2455 }
2456 }
2458 // Retrieve an entity from the stack
2459 label eIndex = mesh.stack(tIndex).pop();
2461 if (mesh.checkBisection(eIndex))
2462 {
2463 if (thread->master())
2464 {
2465 // Bisect this edge
2466 mesh.bisectEdge(eIndex);
2467 }
2468 else
2469 {
2470 // Push this on to the master stack
2471 mesh.stack(0).push(eIndex);
2472 }
2473 }
2474 else
2475 if (mesh.checkCollapse(eIndex))
2476 {
2477 if (thread->master())
2478 {
2479 // Collapse this edge
2480 mesh.collapseEdge(eIndex);
2481 }
2482 else
2483 {
2484 // Push this on to the master stack
2485 mesh.stack(0).push(eIndex);
2486 }
2487 }
2488 }
2490 if (thread->slave())
2491 {
2492 thread->sendSignal(meshHandler::STOP);
2493 }
2495 if (reported)
2496 {
2497 Info<< " Refinement Progress: 100% :"
2498 << " Bisections: " << mesh.status(TOTAL_BISECTIONS)
2499 << ", Collapses: " << mesh.status(TOTAL_COLLAPSES)
2500 << ", Total: " << mesh.status(TOTAL_MODIFICATIONS)
2501 << " \r"
2502 << endl;
2503 }
2504 }
2507 // Remove 2D sliver cells from the mesh
2508 void dynamicTopoFvMesh::remove2DSlivers()
2509 {
2510 // If coupled patches exist, set the flag
2511 if (patchCoupling_.size() || procIndices_.size())
2512 {
2513 setCoupledModification();
2514 }
2516 // Sort by sliver-quality.
2517 labelList cIndices(thresholdSlivers_.toc());
2518 SortableList<scalar> values(cIndices.size());
2520 // Fill-in values to sort by...
2521 forAll(cIndices, indexI)
2522 {
2523 values[indexI] = thresholdSlivers_[cIndices[indexI]];
2524 }
2526 // Explicitly sort by quality value.
2527 values.sort();
2529 const labelList& indices = values.indices();
2531 if (debug)
2532 {
2533 if (thresholdSlivers_.size())
2534 {
2535 Pout<< "Sliver list: " << endl;
2537 forAll(indices, indexI)
2538 {
2539 label cIndex = cIndices[indices[indexI]];
2541 Pout<< " Cell: " << cIndex
2542 << " Quality: " << thresholdSlivers_[cIndex]
2543 << endl;
2544 }
2546 if (debug > 1)
2547 {
2548 writeVTK("sliverCells", cIndices, 3);
2549 }
2550 }
2551 }
2553 forAll(indices, indexI)
2554 {
2555 // Fetch the cell
2556 label cIndex = cIndices[indices[indexI]];
2558 // Ensure that this cell actually exists.
2559 if (cells_[cIndex].empty())
2560 {
2561 continue;
2562 }
2564 const cell& cellToCheck = cells_[cIndex];
2566 // Find an appropriate quad-face
2567 label fIndex = -1;
2569 forAll(cellToCheck, faceI)
2570 {
2571 if (faces_[cellToCheck[faceI]].size() == 4)
2572 {
2573 fIndex = cellToCheck[faceI];
2574 break;
2575 }
2576 }
2578 // Find the prism faces
2579 FixedList<label,2> cBdyIndex(-1), cIntIndex(-1);
2580 FixedList<face,2> cBdyFace, cIntFace;
2582 meshOps::findPrismFaces
2583 (
2584 fIndex,
2585 cIndex,
2586 faces_,
2587 cells_,
2588 neighbour_,
2589 cBdyFace,
2590 cBdyIndex,
2591 cIntFace,
2592 cIntIndex
2593 );
2595 if (debug > 1)
2596 {
2597 InfoIn("void dynamicTopoFvMesh::remove2DSlivers()")
2598 << nl
2599 << " Considering cell: " << cIndex
2600 << ":: " << cells_[cIndex]
2601 << " for sliver removal."
2602 << " Using face: " << fIndex
2603 << ":: " << faces_[fIndex]
2604 << endl;
2605 }
2607 label triFace = cBdyIndex[0], triEdge = -1;
2609 // Find the common-edge between quad-tri faces
2610 meshOps::findCommonEdge
2611 (
2612 fIndex,
2613 triFace,
2614 faceEdges_,
2615 triEdge
2616 );
2618 // Find the isolated point.
2619 label ptIndex = -1, nextPtIndex = -1;
2621 const edge& edgeToCheck = edges_[triEdge];
2623 meshOps::findIsolatedPoint
2624 (
2625 faces_[triFace],
2626 edgeToCheck,
2627 ptIndex,
2628 nextPtIndex
2629 );
2631 // Determine the interior faces connected to each edge-point.
2632 label firstFace = -1, secondFace = -1;
2634 if (cIntFace[0].which(edgeToCheck[0]) > -1)
2635 {
2636 firstFace = cIntIndex[0];
2637 secondFace = cIntIndex[1];
2638 }
2639 else
2640 {
2641 firstFace = cIntIndex[1];
2642 secondFace = cIntIndex[0];
2643 }
2645 point ec =
2646 (
2647 linePointRef
2648 (
2649 points_[edgeToCheck[0]],
2650 points_[edgeToCheck[1]]
2651 ).centre()
2652 );
2654 FixedList<vector, 2> p(vector::zero), q(vector::zero);
2655 FixedList<scalar, 2> proj(0.0);
2657 // Find the projection on the edge.
2658 forAll(edgeToCheck, pointI)
2659 {
2660 p[pointI] = (points_[ptIndex] - points_[edgeToCheck[pointI]]);
2661 q[pointI] = (ec - points_[edgeToCheck[pointI]]);
2663 q[pointI] /= (mag(q[pointI]) + VSMALL);
2665 proj[pointI] = (p[pointI] & q[pointI]);
2666 }
2668 // Take action based on the magnitude of the projection.
2669 if (mag(proj[0]) < VSMALL)
2670 {
2671 if (collapseQuadFace(firstFace).type() > 0)
2672 {
2673 status(TOTAL_SLIVERS)++;
2674 }
2676 continue;
2677 }
2679 if (mag(proj[1]) < VSMALL)
2680 {
2681 if (collapseQuadFace(secondFace).type() > 0)
2682 {
2683 status(TOTAL_SLIVERS)++;
2684 }
2686 continue;
2687 }
2689 if (proj[0] > 0.0 && proj[1] < 0.0)
2690 {
2691 changeMap map = bisectQuadFace(firstFace, changeMap());
2693 // Loop through added faces, and collapse
2694 // the appropriate one
2695 const List<objectMap>& aF = map.addedFaceList();
2697 forAll(aF, faceI)
2698 {
2699 if
2700 (
2701 (owner_[aF[faceI].index()] == cIndex) &&
2702 (aF[faceI].index() != firstFace)
2703 )
2704 {
2705 if (collapseQuadFace(aF[faceI].index()).type() > 0)
2706 {
2707 status(TOTAL_SLIVERS)++;
2708 }
2710 break;
2711 }
2712 }
2714 continue;
2715 }
2717 if (proj[0] < 0.0 && proj[1] > 0.0)
2718 {
2719 changeMap map = bisectQuadFace(secondFace, changeMap());
2721 // Loop through added faces, and collapse
2722 // the appropriate one
2723 const List<objectMap>& aF = map.addedFaceList();
2725 forAll(aF, faceI)
2726 {
2727 if
2728 (
2729 (owner_[aF[faceI].index()] == cIndex) &&
2730 (aF[faceI].index() != secondFace)
2731 )
2732 {
2733 if (collapseQuadFace(aF[faceI].index()).type() > 0)
2734 {
2735 status(TOTAL_SLIVERS)++;
2736 }
2738 break;
2739 }
2740 }
2742 continue;
2743 }
2745 if (proj[0] > 0.0 && proj[1] > 0.0)
2746 {
2747 changeMap map = bisectQuadFace(fIndex, changeMap());
2749 // Loop through added faces, and collapse
2750 // the appropriate one
2751 const List<objectMap>& aF = map.addedFaceList();
2753 forAll(aF, faceI)
2754 {
2755 if
2756 (
2757 (owner_[aF[faceI].index()] == cIndex) &&
2758 (aF[faceI].index() != fIndex)
2759 )
2760 {
2761 if (collapseQuadFace(aF[faceI].index()).type() > 0)
2762 {
2763 status(TOTAL_SLIVERS)++;
2764 }
2766 break;
2767 }
2768 }
2770 continue;
2771 }
2772 }
2774 // Clear out the list
2775 thresholdSlivers_.clear();
2777 // If coupled patches exist, reset the flag
2778 if (patchCoupling_.size() || procIndices_.size())
2779 {
2780 unsetCoupledModification();
2781 }
2782 }
2785 // Identify the sliver type in 3D
2786 const changeMap dynamicTopoFvMesh::identifySliverType
2787 (
2788 const label cIndex
2789 ) const
2790 {
2791 changeMap map;
2793 // Ensure that this cell actually exists.
2794 if (cells_[cIndex].empty())
2795 {
2796 return map;
2797 }
2799 label fourthPoint = -1;
2800 scalar minDistance = GREAT;
2801 face tFace(3), testFace(3), faceToCheck(3);
2802 FixedList<edge, 2> edgeToCheck(edge(-1,-1));
2804 const cell& cellToCheck = cells_[cIndex];
2806 // Find the point-face pair with minimum perpendicular distance
2807 forAll(cellToCheck, faceI)
2808 {
2809 label isolatedPoint = -1;
2810 label nextFaceI = cellToCheck.fcIndex(faceI);
2812 // Pick two faces from this cell.
2813 const face& currFace = faces_[cellToCheck[faceI]];
2814 const face& nextFace = faces_[cellToCheck[nextFaceI]];
2816 // Get the fourth point
2817 forAll(nextFace, pointI)
2818 {
2819 if
2820 (
2821 nextFace[pointI] != currFace[0]
2822 && nextFace[pointI] != currFace[1]
2823 && nextFace[pointI] != currFace[2]
2824 )
2825 {
2826 isolatedPoint = nextFace[pointI];
2828 // Configure a triangular face with correct orientation.
2829 if (owner_[cellToCheck[faceI]] == cIndex)
2830 {
2831 testFace[0] = currFace[2];
2832 testFace[1] = currFace[1];
2833 testFace[2] = currFace[0];
2834 }
2835 else
2836 {
2837 testFace[0] = currFace[0];
2838 testFace[1] = currFace[1];
2839 testFace[2] = currFace[2];
2840 }
2842 break;
2843 }
2844 }
2846 // Obtain the unit normal.
2847 vector testNormal = testFace.normal(points_);
2849 testNormal /= (mag(testNormal) + VSMALL);
2851 // Project the isolated point onto the face.
2852 vector p = points_[isolatedPoint] - points_[testFace[0]];
2853 vector q = p - ((p & testNormal)*testNormal);
2855 // Compute the distance
2856 scalar distance = mag(p - q);
2858 // Is it the least so far?
2859 if (distance < minDistance)
2860 {
2861 // Use this point-face pair.
2862 fourthPoint = isolatedPoint;
2863 tFace = testFace;
2864 minDistance = distance;
2865 }
2866 }
2868 // Obtain the face-normal.
2869 vector refArea = tFace.normal(points_);
2871 // Normalize it.
2872 vector n = refArea/mag(refArea);
2874 // Define edge-vectors.
2875 vector r1 = points_[tFace[1]] - points_[tFace[0]];
2876 vector r2 = points_[tFace[2]] - points_[tFace[1]];
2877 vector r3 = points_[tFace[0]] - points_[tFace[2]];
2879 // Project the fourth point onto the face.
2880 vector r4 = points_[fourthPoint] - points_[tFace[0]];
2882 r4 = r4 - ((r4 & n)*n);
2884 // Define the two other vectors.
2885 vector r5 = r4 - r1;
2886 vector r6 = r5 - r2;
2888 // Calculate three signed triangle areas, using tFace[0] as the origin.
2889 scalar t1 = n & (0.5 * (r1 ^ r4));
2890 scalar t2 = n & (0.5 * (r2 ^ r5));
2891 scalar t3 = n & (0.5 * (r3 ^ r6));
2893 // Determine sliver types based on are magnitudes.
2894 if (t1 > 0 && t2 > 0 && t3 > 0)
2895 {
2896 // Region R0: Cap cell.
2897 map.type() = 2;
2898 map.add("apexPoint", fourthPoint);
2900 faceToCheck[0] = tFace[0];
2901 faceToCheck[1] = tFace[1];
2902 faceToCheck[2] = tFace[2];
2903 }
2905 if (t1 < 0 && t2 > 0 && t3 > 0)
2906 {
2907 // Region R1: Sliver cell.
2908 map.type() = 1;
2910 edgeToCheck[0][0] = tFace[2];
2911 edgeToCheck[0][1] = fourthPoint;
2912 edgeToCheck[1][0] = tFace[0];
2913 edgeToCheck[1][1] = tFace[1];
2914 }
2916 if (t1 > 0 && t2 < 0 && t3 > 0)
2917 {
2918 // Region R2: Sliver cell.
2919 map.type() = 1;
2921 edgeToCheck[0][0] = tFace[0];
2922 edgeToCheck[0][1] = fourthPoint;
2923 edgeToCheck[1][0] = tFace[1];
2924 edgeToCheck[1][1] = tFace[2];
2925 }
2927 if (t1 > 0 && t2 > 0 && t3 < 0)
2928 {
2929 // Region R3: Sliver cell.
2930 map.type() = 1;
2932 edgeToCheck[0][0] = tFace[1];
2933 edgeToCheck[0][1] = fourthPoint;
2934 edgeToCheck[1][0] = tFace[2];
2935 edgeToCheck[1][1] = tFace[0];
2936 }
2938 if (t1 < 0 && t2 > 0 && t3 < 0)
2939 {
2940 // Region R4: Cap cell.
2941 map.type() = 2;
2942 map.add("apexPoint", tFace[0]);
2944 faceToCheck[0] = tFace[1];
2945 faceToCheck[1] = tFace[2];
2946 faceToCheck[2] = fourthPoint;
2947 }
2949 if (t1 < 0 && t2 < 0 && t3 > 0)
2950 {
2951 // Region R5: Cap cell.
2952 map.type() = 2;
2953 map.add("apexPoint", tFace[1]);
2955 faceToCheck[0] = tFace[2];
2956 faceToCheck[1] = tFace[0];
2957 faceToCheck[2] = fourthPoint;
2958 }
2960 if (t1 > 0 && t2 < 0 && t3 < 0)
2961 {
2962 // Region R6: Cap cell.
2963 map.type() = 2;
2964 map.add("apexPoint", tFace[2]);
2966 faceToCheck[0] = tFace[0];
2967 faceToCheck[1] = tFace[1];
2968 faceToCheck[2] = fourthPoint;
2969 }
2971 // See if an over-ride to wedge/spade is necessary.
2972 // Obtain a reference area magnitude.
2973 scalar refMag = 0.1*(refArea & n);
2975 if (mag(t1) < refMag)
2976 {
2977 if (mag(t3) < refMag)
2978 {
2979 // Wedge case: Too close to point [0]
2980 map.type() = 4;
2982 edgeToCheck[0][0] = tFace[0];
2983 edgeToCheck[0][1] = fourthPoint;
2984 }
2985 else
2986 if (mag(t2) < refMag)
2987 {
2988 // Wedge case: Too close to point [1]
2989 map.type() = 4;
2991 edgeToCheck[0][0] = tFace[1];
2992 edgeToCheck[0][1] = fourthPoint;
2993 }
2994 else
2995 if ((mag(t2) > refMag) && (mag(t3) > refMag))
2996 {
2997 // Spade case: Too close to edge vector r1
2998 map.type() = 3;
2999 map.add("apexPoint", fourthPoint);
3001 edgeToCheck[0][0] = tFace[0];
3002 edgeToCheck[0][1] = tFace[1];
3003 }
3004 }
3006 if (mag(t2) < refMag)
3007 {
3008 if (mag(t3) < refMag)
3009 {
3010 // Wedge case: Too close to point [2]
3011 map.type() = 4;
3013 edgeToCheck[0][0] = tFace[2];
3014 edgeToCheck[0][1] = fourthPoint;
3015 }
3016 else
3017 if ((mag(t1) > refMag) && (mag(t3) > refMag))
3018 {
3019 // Spade case: Too close to edge vector r2
3020 map.type() = 3;
3021 map.add("apexPoint", fourthPoint);
3023 edgeToCheck[0][0] = tFace[1];
3024 edgeToCheck[0][1] = tFace[2];
3025 }
3026 }
3028 if (mag(t3) < refMag)
3029 {
3030 if ((mag(t1) > refMag) && (mag(t2) > refMag))
3031 {
3032 // Spade case: Too close to edge vector r3
3033 map.type() = 3;
3034 map.add("apexPoint", fourthPoint);
3036 edgeToCheck[0][0] = tFace[2];
3037 edgeToCheck[0][1] = tFace[0];
3038 }
3039 }
3041 // Determine appropriate information for sliver exudation.
3042 switch (map.type())
3043 {
3044 case 1:
3045 {
3046 FixedList<bool, 2> foundEdge(false);
3048 // Search the cell-faces for first and second edges.
3049 forAll(cellToCheck, faceI)
3050 {
3051 const labelList& fEdges = faceEdges_[cellToCheck[faceI]];
3053 forAll(fEdges, edgeI)
3054 {
3055 const edge& thisEdge = edges_[fEdges[edgeI]];
3057 if (thisEdge == edgeToCheck[0])
3058 {
3059 map.add("firstEdge", fEdges[edgeI]);
3061 foundEdge[0] = true;
3062 }
3064 if (thisEdge == edgeToCheck[1])
3065 {
3066 map.add("secondEdge", fEdges[edgeI]);
3068 foundEdge[1] = true;
3069 }
3070 }
3072 if (foundEdge[0] && foundEdge[1])
3073 {
3074 break;
3075 }
3076 }
3078 break;
3079 }
3081 case 2:
3082 {
3083 // Search the cell-faces for opposing faces.
3084 forAll(cellToCheck, faceI)
3085 {
3086 const face& thisFace = faces_[cellToCheck[faceI]];
3088 if (triFace::compare(triFace(thisFace), triFace(faceToCheck)))
3089 {
3090 map.add("opposingFace", cellToCheck[faceI]);
3092 break;
3093 }
3094 }
3096 break;
3097 }
3099 case 3:
3100 case 4:
3101 {
3102 bool foundEdge = false;
3104 // Search the cell-faces for first edge.
3105 forAll(cellToCheck, faceI)
3106 {
3107 const labelList& fEdges = faceEdges_[cellToCheck[faceI]];
3109 forAll(fEdges, edgeI)
3110 {
3111 const edge& thisEdge = edges_[fEdges[edgeI]];
3113 if (thisEdge == edgeToCheck[0])
3114 {
3115 map.add("firstEdge", fEdges[edgeI]);
3117 foundEdge = true;
3118 }
3119 }
3121 if (foundEdge)
3122 {
3123 break;
3124 }
3125 }
3127 break;
3128 }
3130 default:
3131 {
3132 WarningIn
3133 (
3134 "void dynamicTopoFvMesh::identifySliverType"
3135 "(const label cIndex) const"
3136 )
3137 << nl << "Could not identify sliver type for cell: "
3138 << cIndex
3139 << endl;
3140 }
3141 }
3143 if (debug > 2)
3144 {
3145 Pout<< "Cell: " << cIndex
3146 << " Identified sliver type as: "
3147 << map.type() << endl;
3148 }
3150 // Return the result.
3151 return map;
3152 }
3155 // Remove sliver cells
3156 void dynamicTopoFvMesh::removeSlivers()
3157 {
3158 if (!edgeRefinement_)
3159 {
3160 return;
3161 }
3163 // Check if a removeSlivers entry was found in the dictionary
3164 if (dict_.subDict("dynamicTopoFvMesh").found("removeSlivers"))
3165 {
3166 Switch rs =
3167 (
3168 dict_.subDict("dynamicTopoFvMesh").lookup("removeSlivers")
3169 );
3171 if (!rs)
3172 {
3173 return;
3174 }
3175 }
3177 // Invoke the 2D sliver removal routine
3178 if (is2D())
3179 {
3180 remove2DSlivers();
3181 return;
3182 }
3184 // If coupled patches exist, set the flag
3185 if (patchCoupling_.size() || procIndices_.size())
3186 {
3187 setCoupledModification();
3188 }
3190 // Sort by sliver-quality.
3191 labelList cIndices(thresholdSlivers_.toc());
3192 SortableList<scalar> values(cIndices.size());
3194 // Fill-in values to sort by...
3195 forAll(cIndices, indexI)
3196 {
3197 values[indexI] = thresholdSlivers_[cIndices[indexI]];
3198 }
3200 // Explicitly sort by quality value.
3201 values.sort();
3203 const labelList& indices = values.indices();
3205 if (debug && thresholdSlivers_.size())
3206 {
3207 Pout<< "Sliver list: " << endl;
3209 forAll(indices, indexI)
3210 {
3211 label cIndex = cIndices[indices[indexI]];
3213 Pout<< " Cell: " << cIndex
3214 << " Quality: " << thresholdSlivers_[cIndex]
3215 << endl;
3216 }
3218 if (debug > 1)
3219 {
3220 writeVTK("sliverCells", cIndices, 3);
3221 }
3222 }
3224 forAll(indices, indexI)
3225 {
3226 // Fetch the cell index
3227 label cIndex = cIndices[indices[indexI]];
3229 // First check if this sliver cell is handled elsewhere.
3230 if (procIndices_.size())
3231 {
3232 bool foundInSubMesh = false;
3234 forAll(procIndices_, procI)
3235 {
3236 label proc = procIndices_[procI];
3238 if (priority(proc, lessOp<label>(), Pstream::myProcNo()))
3239 {
3240 Map<label>& rCellMap =
3241 (
3242 sendMeshes_[proc].map().reverseEntityMap
3243 (
3244 coupleMap::CELL
3245 )
3246 );
3248 if (rCellMap.found(cIndex))
3249 {
3250 // This cell was sent to another sub-domain.
3251 foundInSubMesh = true;
3252 break;
3253 }
3254 }
3255 }
3257 if (foundInSubMesh)
3258 {
3259 continue;
3260 }
3261 }
3263 // Identify the sliver type.
3264 changeMap map = identifySliverType(cIndex);
3266 if (debug)
3267 {
3268 InfoIn("void dynamicTopoFvMesh::removeSlivers()")
3269 << nl << "Removing Cell: " << cIndex
3270 << " of sliver type: " << map.type()
3271 << " with quality: " << thresholdSlivers_[cIndex]
3272 << endl;
3273 }
3275 bool success = false;
3277 // Take action based on the type of sliver.
3278 switch (map.type())
3279 {
3280 case -1:
3281 {
3282 // Sliver cell was removed by a prior operation.
3283 // Nothing needs to be done.
3284 break;
3285 }
3287 case 1:
3288 {
3289 // Sliver cell.
3290 // Determine which edges need to be bisected.
3291 label firstEdge = readLabel(map.lookup("firstEdge"));
3292 label secondEdge = readLabel(map.lookup("secondEdge"));
3294 // Force bisection on both edges.
3295 changeMap firstMap = bisectEdge(firstEdge, false, true);
3296 changeMap secondMap = bisectEdge(secondEdge, false, true);
3298 // Collapse the intermediate edge.
3299 // Since we don't know which edge it is, search
3300 // through recently added edges and compare.
3301 edge edgeToCheck
3302 (
3303 firstMap.addedPointList()[0].index(),
3304 secondMap.addedPointList()[0].index()
3305 );
3307 bool foundCollapseEdge = false;
3309 const List<objectMap>& firstMapEdges =
3310 (
3311 firstMap.addedEdgeList()
3312 );
3314 const List<objectMap>& secondMapEdges =
3315 (
3316 secondMap.addedEdgeList()
3317 );
3319 // Loop through the first list.
3320 forAll(firstMapEdges, edgeI)
3321 {
3322 const edge& thisEdge =
3323 (
3324 edges_[firstMapEdges[edgeI].index()]
3325 );
3327 if (thisEdge == edgeToCheck)
3328 {
3329 // Collapse this edge.
3330 if
3331 (
3332 collapseEdge
3333 (
3334 firstMapEdges[edgeI].index(),
3335 -1,
3336 false,
3337 true
3338 ).type() > 0
3339 )
3340 {
3341 success = true;
3342 }
3344 foundCollapseEdge = true;
3345 break;
3346 }
3347 }
3349 // Loop through the second list.
3350 if (!foundCollapseEdge)
3351 {
3352 forAll(secondMapEdges, edgeI)
3353 {
3354 const edge& thisEdge =
3355 (
3356 edges_[secondMapEdges[edgeI].index()]
3357 );
3359 if (thisEdge == edgeToCheck)
3360 {
3361 // Collapse this edge.
3362 if
3363 (
3364 collapseEdge
3365 (
3366 secondMapEdges[edgeI].index(),
3367 -1,
3368 false,
3369 true
3370 ).type() > 0
3371 )
3372 {
3373 success = true;
3374 }
3376 break;
3377 }
3378 }
3379 }
3381 break;
3382 }
3384 case 2:
3385 {
3386 // Cap cell.
3387 //label opposingFace = readLabel(map.lookup("opposingFace"));
3389 // Force trisection of the opposing face.
3390 changeMap faceMap;
3391 //(
3392 // trisectFace(opposingFace, false, true)
3393 //);
3395 // Collapse the intermediate edge.
3396 // Since we don't know which edge it is, search
3397 // through recently added edges and compare.
3398 edge edgeToCheck
3399 (
3400 readLabel(map.lookup("apexPoint")),
3401 faceMap.addedPointList()[0].index()
3402 );
3404 const List<objectMap>& faceMapEdges =
3405 (
3406 faceMap.addedEdgeList()
3407 );
3409 forAll(faceMapEdges, edgeI)
3410 {
3411 const edge& thisEdge = edges_[faceMapEdges[edgeI].index()];
3413 if (thisEdge == edgeToCheck)
3414 {
3415 // Collapse this edge.
3416 if
3417 (
3418 collapseEdge
3419 (
3420 faceMapEdges[edgeI].index(),
3421 -1,
3422 false,
3423 true
3424 ).type() > 0
3425 )
3426 {
3427 success = true;
3428 }
3430 break;
3431 }
3432 }
3434 break;
3435 }
3437 case 3:
3438 {
3439 // Spade cell.
3441 // Force bisection on the first edge.
3442 changeMap firstMap =
3443 (
3444 bisectEdge
3445 (
3446 readLabel(map.lookup("firstEdge")),
3447 false,
3448 true
3449 )
3450 );
3452 // Collapse the intermediate edge.
3453 // Since we don't know which edge it is, search
3454 // through recently added edges and compare.
3455 edge edgeToCheck
3456 (
3457 readLabel(map.lookup("apexPoint")),
3458 firstMap.addedPointList()[0].index()
3459 );
3461 const List<objectMap>& firstMapEdges =
3462 (
3463 firstMap.addedEdgeList()
3464 );
3466 // Loop through the first list.
3467 forAll(firstMapEdges, edgeI)
3468 {
3469 const edge& thisEdge = edges_[firstMapEdges[edgeI].index()];
3471 if (thisEdge == edgeToCheck)
3472 {
3473 // Collapse this edge.
3474 if
3475 (
3476 collapseEdge
3477 (
3478 firstMapEdges[edgeI].index(),
3479 -1,
3480 false,
3481 true
3482 ).type() > 0
3483 )
3484 {
3485 success = true;
3486 }
3488 break;
3489 }
3490 }
3492 break;
3493 }
3495 case 4:
3496 {
3497 // Wedge cell.
3499 // Collapse the first edge.
3500 if
3501 (
3502 collapseEdge
3503 (
3504 readLabel(map.lookup("firstEdge")),
3505 -1,
3506 false,
3507 true
3508 ).type() > 0
3509 )
3510 {
3511 success = true;
3512 }
3514 break;
3515 }
3517 default:
3518 {
3519 WarningIn("void dynamicTopoFvMesh::removeSlivers()")
3520 << nl << "Could not identify sliver type for cell: "
3521 << cIndex
3522 << endl;
3523 }
3524 }
3526 // Increment the count for successful sliver removal
3527 if (success)
3528 {
3529 status(TOTAL_SLIVERS)++;
3530 }
3531 }
3533 // Clear out the list
3534 thresholdSlivers_.clear();
3536 // If coupled patches exist, reset the flag
3537 if (patchCoupling_.size() || procIndices_.size())
3538 {
3539 unsetCoupledModification();
3540 }
3541 }
3544 // Given an input quad-face, determine checkEdges from mesh
3545 // - Does not refer to member data directly,
3546 // since this is also used by subMeshes
3547 void dynamicTopoFvMesh::getCheckEdges
3548 (
3549 const label fIndex,
3550 const dynamicTopoFvMesh& mesh,
3551 changeMap& map,
3552 FixedList<edge,4>& checkEdge,
3553 FixedList<label,4>& checkEdgeIndex
3554 )
3555 {
3556 checkEdgeIndex[0] = mesh.getTriBoundaryEdge(fIndex);
3557 checkEdge[0] = mesh.edges_[checkEdgeIndex[0]];
3559 const labelList& fEdges = mesh.faceEdges_[fIndex];
3561 forAll(fEdges, edgeI)
3562 {
3563 if (checkEdgeIndex[0] != fEdges[edgeI])
3564 {
3565 const edge& thisEdge = mesh.edges_[fEdges[edgeI]];
3567 if
3568 (
3569 checkEdge[0].start() == thisEdge[0] ||
3570 checkEdge[0].start() == thisEdge[1]
3571 )
3572 {
3573 checkEdgeIndex[1] = fEdges[edgeI];
3574 checkEdge[1] = thisEdge;
3576 // Update the map
3577 map.add("firstEdge", checkEdgeIndex[1]);
3578 }
3579 else
3580 if
3581 (
3582 checkEdge[0].end() == thisEdge[0] ||
3583 checkEdge[0].end() == thisEdge[1]
3584 )
3585 {
3586 checkEdgeIndex[2] = fEdges[edgeI];
3587 checkEdge[2] = thisEdge;
3589 // Update the map
3590 map.add("secondEdge", checkEdgeIndex[2]);
3591 }
3592 else
3593 {
3594 checkEdgeIndex[3] = fEdges[edgeI];
3595 checkEdge[3] = thisEdge;
3596 }
3597 }
3598 }
3599 }
3602 // Return length-scale at an face-location in the mesh [2D]
3603 scalar dynamicTopoFvMesh::faceLengthScale
3604 (
3605 const label fIndex
3606 ) const
3607 {
3608 // Reset the scale first
3609 scalar scale = 0.0;
3611 label facePatch = whichPatch(fIndex);
3613 // Determine whether the face is internal
3614 if (facePatch < 0)
3615 {
3616 scale =
3617 (
3618 0.5 *
3619 (
3620 lengthScale_[owner_[fIndex]]
3621 + lengthScale_[neighbour_[fIndex]]
3622 )
3623 );
3624 }
3625 else
3626 {
3627 // Fetch the fixed-length scale
3628 scale = lengthEstimator().fixedLengthScale(fIndex, facePatch);
3630 // If this is a floating face, pick the owner length-scale
3631 if (lengthEstimator().isFreePatch(facePatch))
3632 {
3633 scale = lengthScale_[owner_[fIndex]];
3634 }
3636 // If proximity-based refinement is requested,
3637 // test the proximity to the nearest non-neighbour.
3638 if (lengthEstimator().isProximityPatch(facePatch))
3639 {
3640 label proximityFace = -1;
3642 // Perform a proximity-check.
3643 scalar distance = testProximity(fIndex, proximityFace);
3645 if (debug > 3 && self() == 0)
3646 {
3647 if
3648 (
3649 (proximityFace > -1) &&
3650 ((distance / 5.0) < scale)
3651 )
3652 {
3653 Pout<< " Closest opposing face detected for face: " << nl
3654 << '\t' << fIndex
3655 << " :: " << faces_[fIndex]
3656 << " was face:\n"
3657 << '\t' << proximityFace
3658 << " :: " << polyMesh::faces()[proximityFace] << nl
3659 << " with distance: " << distance
3660 << endl;
3661 }
3662 }
3664 scale =
3665 (
3666 Foam::min
3667 (
3668 scale,
3669 ((distance / 3.0) - SMALL)/lengthEstimator().ratioMax()
3670 )
3671 );
3672 }
3674 // If this face lies on a processor patch,
3675 // fetch lengthScale info from patchSubMeshes
3676 if (processorCoupledEntity(fIndex))
3677 {
3678 scale = processorLengthScale(fIndex);
3679 }
3681 // Limit scales if necessary
3682 lengthEstimator().limitScale(scale);
3683 }
3685 return scale;
3686 }
3689 // Compute length-scale at an edge-location in the mesh [3D]
3690 scalar dynamicTopoFvMesh::edgeLengthScale
3691 (
3692 const label eIndex
3693 ) const
3694 {
3695 // Reset the scale first
3696 scalar scale = 0.0;
3698 const labelList& eFaces = edgeFaces_[eIndex];
3700 label edgePatch = whichEdgePatch(eIndex);
3702 // Determine whether the edge is internal
3703 if (edgePatch < 0)
3704 {
3705 forAll(eFaces, faceI)
3706 {
3707 scale += lengthScale_[owner_[eFaces[faceI]]];
3708 scale += lengthScale_[neighbour_[eFaces[faceI]]];
3709 }
3711 scale /= (2.0*eFaces.size());
3712 }
3713 else
3714 {
3715 // Search for boundary faces, and average their scale
3716 forAll(eFaces, faceI)
3717 {
3718 label facePatch = whichPatch(eFaces[faceI]);
3720 // Skip internal faces
3721 if (facePatch == -1)
3722 {
3723 continue;
3724 }
3726 // If this is a floating face, pick the owner length-scale
3727 if (lengthEstimator().isFreePatch(facePatch))
3728 {
3729 scale += lengthScale_[owner_[eFaces[faceI]]];
3730 }
3731 else
3732 {
3733 // Fetch fixed length-scale
3734 scale +=
3735 (
3736 lengthEstimator().fixedLengthScale
3737 (
3738 eFaces[faceI],
3739 facePatch
3740 )
3741 );
3742 }
3743 }
3745 scale *= 0.5;
3747 // If proximity-based refinement is requested,
3748 // test the proximity to the nearest non-neighbour.
3749 if (lengthEstimator().isProximityPatch(edgePatch))
3750 {
3751 label proximityFace = -1;
3753 // Perform a proximity-check.
3754 scalar distance = testProximity(eIndex, proximityFace);
3756 if (debug > 3 && self() == 0)
3757 {
3758 if
3759 (
3760 (proximityFace > -1) &&
3761 ((distance / 5.0) < scale)
3762 )
3763 {
3764 Pout<< " Closest opposing face detected for edge: " << nl
3765 << '\t' << eIndex
3766 << " :: " << edges_[eIndex]
3767 << " was face:\n"
3768 << '\t' << proximityFace
3769 << " :: " << polyMesh::faces()[proximityFace] << nl
3770 << " with distance: " << distance
3771 << endl;
3772 }
3773 }
3775 scale =
3776 (
3777 Foam::min
3778 (
3779 scale,
3780 ((distance / 3.0) - SMALL)/lengthEstimator().ratioMax()
3781 )
3782 );
3783 }
3785 // If curvature-based refinement is requested,
3786 // test the variation in face-normal directions.
3787 if (lengthEstimator().isCurvaturePatch(edgePatch))
3788 {
3789 // Obtain face-normals for both faces.
3790 label count = 0;
3791 FixedList<vector, 2> fNorm;
3793 forAll(eFaces, faceI)
3794 {
3795 if (neighbour_[eFaces[faceI]] == -1)
3796 {
3797 // Obtain the normal.
3798 fNorm[count] = faces_[eFaces[faceI]].normal(points_);
3800 // Normalize it.
3801 fNorm[count] /= mag(fNorm[count]);
3803 count++;
3804 }
3805 }
3807 scalar deviation = (fNorm[0] & fNorm[1]);
3808 scalar refDeviation = lengthEstimator().curvatureDeviation();
3810 if (mag(deviation) < refDeviation)
3811 {
3812 // Fetch the edge
3813 const edge& edgeToCheck = edges_[eIndex];
3815 // Get the edge-length.
3816 scalar length =
3817 (
3818 linePointRef
3819 (
3820 points_[edgeToCheck.start()],
3821 points_[edgeToCheck.end()]
3822 ).mag()
3823 );
3825 if (debug > 3 && self() == 0)
3826 {
3827 Pout<< "Deviation: " << deviation << nl
3828 << "curvatureDeviation: " << refDeviation
3829 << ", Edge: " << eIndex << ", Length: " << length
3830 << ", Scale: " << scale << nl
3831 << " Half-length: " << (0.5*length) << nl
3832 << " MinRatio: "
3833 << (lengthEstimator().ratioMin()*scale)
3834 << endl;
3835 }
3837 scale =
3838 (
3839 Foam::min
3840 (
3841 scale,
3842 ((length - SMALL)/lengthEstimator().ratioMax())
3843 )
3844 );
3845 }
3846 }
3848 // If this edge lies on a processor patch,
3849 // fetch lengthScale info from patchSubMeshes
3850 if (processorCoupledEntity(eIndex))
3851 {
3852 scale = processorLengthScale(eIndex);
3853 }
3855 // Limit scales if necessary
3856 lengthEstimator().limitScale(scale);
3857 }
3859 return scale;
3860 }
3863 // MultiThreaded topology modifier
3864 void dynamicTopoFvMesh::threadedTopoModifier()
3865 {
3866 // Set sizes for the reverse maps
3867 reversePointMap_.setSize(nPoints_, -7);
3868 reverseEdgeMap_.setSize(nEdges_, -7);
3869 reverseFaceMap_.setSize(nFaces_, -7);
3870 reverseCellMap_.setSize(nCells_, -7);
3872 // Remove sliver cells first.
3873 removeSlivers();
3875 // Coupled entities to avoid during normal modification
3876 labelHashSet entities;
3878 // Handle coupled patches.
3879 handleCoupledPatches(entities);
3881 // Handle layer addition / removal
3882 handleLayerAdditionRemoval();
3884 // Set the thread scheduling sequence
3885 labelList topoSequence(threader_->getNumThreads());
3887 // Linear sequence from 1 to nThreads
3888 forAll(topoSequence, indexI)
3889 {
3890 topoSequence[indexI] = indexI + 1;
3891 }
3893 if (edgeRefinement_)
3894 {
3895 // Initialize stacks
3896 initStacks(entities);
3898 // Execute threads
3899 if (threader_->multiThreaded())
3900 {
3901 executeThreads(topoSequence, handlerPtr_, &edgeRefinementEngine);
3902 }
3904 // Set the master thread to implement modifications
3905 edgeRefinementEngine(&(handlerPtr_[0]));
3907 // Handle mesh slicing events, if necessary
3908 handleMeshSlicing();
3910 if (debug)
3911 {
3912 Info<< nl << "Edge Bisection/Collapse complete." << endl;
3913 }
3914 }
3916 // Re-Initialize stacks
3917 initStacks(entities);
3919 // Execute threads
3920 if (threader_->multiThreaded())
3921 {
3922 if (is2D())
3923 {
3924 executeThreads(topoSequence, handlerPtr_, &swap2DEdges);
3925 }
3926 else
3927 {
3928 executeThreads(topoSequence, handlerPtr_, &swap3DEdges);
3929 }
3930 }
3932 // Set the master thread to implement modifications
3933 if (is2D())
3934 {
3935 swap2DEdges(&(handlerPtr_[0]));
3936 }
3937 else
3938 {
3939 swap3DEdges(&(handlerPtr_[0]));
3940 }
3942 if (debug)
3943 {
3944 Info<< nl << "Edge Swapping complete." << endl;
3945 }
3947 // Synchronize coupled patches
3948 syncCoupledPatches(entities);
3949 }
3952 // Reset the mesh and generate mapping information
3953 // - Return true if topology changes were made.
3954 // - Return false otherwise (motion only)
3955 bool dynamicTopoFvMesh::resetMesh()
3956 {
3957 // Reduce across processors.
3958 reduce(topoChangeFlag_, orOp<bool>());
3960 const dictionary& meshSubDict = dict_.subDict("dynamicTopoFvMesh");
3962 if (topoChangeFlag_)
3963 {
3964 // Write out statistics
3965 if (Pstream::parRun())
3966 {
3967 if (debug)
3968 {
3969 Pout<< " Bisections :: Total: " << status(TOTAL_BISECTIONS)
3970 << ", Surface: " << status(SURFACE_BISECTIONS) << nl
3971 << " Collapses :: Total: " << status(TOTAL_COLLAPSES)
3972 << ", Surface: " << status(SURFACE_COLLAPSES) << nl
3973 << " Swaps :: Total: " << status(TOTAL_SWAPS)
3974 << ", Surface: " << status(SURFACE_SWAPS) << endl;
3975 }
3977 reduce(statistics_, sumOp<labelList>());
3978 }
3980 Info<< " Bisections :: Total: " << status(TOTAL_BISECTIONS)
3981 << ", Surface: " << status(SURFACE_BISECTIONS) << nl
3982 << " Collapses :: Total: " << status(TOTAL_COLLAPSES)
3983 << ", Surface: " << status(SURFACE_COLLAPSES) << nl
3984 << " Swaps :: Total: " << status(TOTAL_SWAPS)
3985 << ", Surface: " << status(SURFACE_SWAPS) << endl;
3987 if (status(TOTAL_SLIVERS))
3988 {
3989 Pout<< " Slivers :: " << status(TOTAL_SLIVERS) << endl;
3990 }
3992 // Fetch reference to mapper
3993 const topoMapper& fieldMapper = mapper_();
3995 // Set information for the mapping stage
3996 // - Must be done prior to field-transfers and mesh reset
3997 fieldMapper.storeMeshInformation();
3999 // Obtain the number of patches before
4000 // any possible boundary reset
4001 label nOldPatches = boundaryMesh().size();
4003 // If the number of patches have changed
4004 // at run-time, reset boundaries first
4005 if (nPatches_ != nOldPatches)
4006 {
4007 resetBoundaries();
4008 }
4010 // Set up field-transfers before dealing with mapping
4011 wordList fieldTypes;
4012 List<wordList> fieldNames;
4013 List<List<char> > sendBuffer, recvBuffer;
4015 // Subset fields and transfer
4016 initFieldTransfers
4017 (
4018 fieldTypes,
4019 fieldNames,
4020 sendBuffer,
4021 recvBuffer
4022 );
4024 // Set sizes for mapping
4025 faceWeights_.setSize(facesFromFaces_.size(), scalarField(0));
4026 faceCentres_.setSize(facesFromFaces_.size(), vectorField(0));
4027 cellWeights_.setSize(cellsFromCells_.size(), scalarField(0));
4028 cellCentres_.setSize(cellsFromCells_.size(), vectorField(0));
4030 // Determine if mapping is to be skipped
4031 // Optionally skip mapping for remeshing-only / pre-processing
4032 bool skipMapping = false;
4034 if (meshSubDict.found("skipMapping") || mandatory_)
4035 {
4036 skipMapping = readBool(meshSubDict.lookup("skipMapping"));
4037 }
4039 // Fetch the tolerance for mapping
4040 scalar mapTol = 1e-10;
4042 if (meshSubDict.found("mappingTol") || mandatory_)
4043 {
4044 mapTol = readScalar(meshSubDict.lookup("mappingTol"));
4045 }
4047 // Check if outputs are enabled on failure
4048 bool mappingOutput = false;
4050 if (meshSubDict.found("mappingOutput") || mandatory_)
4051 {
4052 mappingOutput = readBool(meshSubDict.lookup("mappingOutput"));
4053 }
4055 clockTime mappingTimer;
4057 // Compute mapping weights for modified entities
4058 threadedMapping(mapTol, skipMapping, mappingOutput);
4060 // Print out stats
4061 Info<< " Mapping time: "
4062 << mappingTimer.elapsedTime() << " s"
4063 << endl;
4065 // Synchronize field transfers prior to the reOrdering stage
4066 syncFieldTransfers
4067 (
4068 fieldTypes,
4069 fieldNames,
4070 recvBuffer
4071 );
4073 // Obtain references to zones, if any
4074 pointZoneMesh& pointZones = polyMesh::pointZones();
4075 faceZoneMesh& faceZones = polyMesh::faceZones();
4076 cellZoneMesh& cellZones = polyMesh::cellZones();
4078 // Allocate temporary lists for mesh-reset
4079 pointField points(nPoints_);
4080 pointField preMotionPoints(nPoints_);
4081 edgeList edges(nEdges_);
4082 faceList faces(nFaces_);
4083 labelList owner(nFaces_);
4084 labelList neighbour(nInternalFaces_);
4085 labelListList faceEdges(nFaces_);
4086 labelListList edgeFaces(nEdges_);
4087 labelList oldPatchStarts(oldPatchStarts_);
4088 labelList oldPatchNMeshPoints(oldPatchNMeshPoints_);
4089 labelListList pointZoneMap(pointZones.size());
4090 labelListList faceZonePointMap(faceZones.size());
4091 labelListList faceZoneFaceMap(faceZones.size());
4092 labelListList cellZoneMap(cellZones.size());
4094 // Obtain faceZone point maps before reordering
4095 List<Map<label> > oldFaceZonePointMaps(faceZones.size());
4097 forAll(faceZones, fzI)
4098 {
4099 oldFaceZonePointMaps[fzI] = faceZones[fzI]().meshPointMap();
4100 }
4102 clockTime reOrderingTimer;
4104 // Reorder the mesh and obtain current topological information
4105 reOrderMesh
4106 (
4107 points,
4108 preMotionPoints,
4109 edges,
4110 faces,
4111 owner,
4112 neighbour,
4113 faceEdges,
4114 edgeFaces,
4115 pointZoneMap,
4116 faceZoneFaceMap,
4117 cellZoneMap
4118 );
4120 // Print out stats
4121 Info<< " Reordering time: "
4122 << reOrderingTimer.elapsedTime() << " s"
4123 << endl;
4125 // Obtain the patch-point maps before resetting the mesh
4126 List<Map<label> > oldPatchPointMaps(nOldPatches);
4128 forAll(oldPatchPointMaps, patchI)
4129 {
4130 oldPatchPointMaps[patchI] = boundaryMesh()[patchI].meshPointMap();
4131 }
4133 // Set weighting information.
4134 // This takes over the weight data.
4135 fieldMapper.setFaceWeights
4136 (
4137 xferMove(faceWeights_),
4138 xferMove(faceCentres_)
4139 );
4141 fieldMapper.setCellWeights
4142 (
4143 xferMove(cellWeights_),
4144 xferMove(cellCentres_)
4145 );
4147 // Reset the mesh, and specify a non-valid
4148 // boundary to avoid globalData construction
4149 polyMesh::resetPrimitives
4150 (
4151 xferCopy(preMotionPoints),
4152 xferMove(faces),
4153 xferMove(owner),
4154 xferMove(neighbour),
4155 patchSizes_,
4156 patchStarts_,
4157 false
4158 );
4160 // Check the dictionary to determine whether
4161 // edge-connectivity is to be stored on disk.
4162 // This usually benefits restart-time for large
4163 // cases, at the expense of disk-space.
4164 bool storePrimitives = false;
4166 if (meshSubDict.found("storeEdgePrimitives") || mandatory_)
4167 {
4168 storePrimitives =
4169 (
4170 readBool(meshSubDict.lookup("storeEdgePrimitives"))
4171 );
4172 }
4174 // Reset the edge mesh
4175 eMeshPtr_->resetPrimitives
4176 (
4177 edges,
4178 faceEdges,
4179 edgeFaces,
4180 edgePatchSizes_,
4181 edgePatchStarts_,
4182 true,
4183 (time().outputTime() && storePrimitives)
4184 );
4186 // Generate mapping for points on boundary patches
4187 labelListList patchPointMap(nPatches_);
4189 for (label i = 0; i < nPatches_; i++)
4190 {
4191 // Obtain new patch mesh points after reset.
4192 const labelList& meshPointLabels = boundaryMesh()[i].meshPoints();
4194 patchNMeshPoints_[i] = meshPointLabels.size();
4196 patchPointMap[i].setSize(meshPointLabels.size(), -1);
4198 // Skip for newly introduced patches
4199 if (i >= nOldPatches)
4200 {
4201 continue;
4202 }
4204 forAll(meshPointLabels, pointI)
4205 {
4206 label oldIndex = pointMap_[meshPointLabels[pointI]];
4208 // Check if the point existed before...
4209 if (oldIndex > -1)
4210 {
4211 // Look for the old position on this patch.
4212 Map<label>::const_iterator oIter =
4213 (
4214 oldPatchPointMaps[i].find(oldIndex)
4215 );
4217 // Add an entry if the point was found
4218 if (oIter != oldPatchPointMaps[i].end())
4219 {
4220 patchPointMap[i][pointI] = oIter();
4221 }
4222 }
4223 }
4224 }
4226 // Generate mapping for faceZone points
4227 forAll(faceZones, fzI)
4228 {
4229 // Obtain new face zone mesh points after reset.
4230 const labelList& meshPointLabels = faceZones[fzI]().meshPoints();
4232 faceZonePointMap[fzI].setSize(meshPointLabels.size(), -1);
4234 forAll(meshPointLabels, pointI)
4235 {
4236 label oldIndex = pointMap_[meshPointLabels[pointI]];
4238 // Check if the point existed before...
4239 if (oldIndex > -1)
4240 {
4241 // Look for the old position on this patch.
4242 Map<label>::const_iterator oIter =
4243 (
4244 oldFaceZonePointMaps[fzI].find(oldIndex)
4245 );
4247 // Add an entry if the point was found
4248 if (oIter != oldFaceZonePointMaps[fzI].end())
4249 {
4250 faceZonePointMap[fzI][pointI] = oIter();
4251 }
4252 }
4253 }
4254 }
4256 // Generate new mesh mapping information
4257 mapPolyMesh mpm
4258 (
4259 (*this),
4260 nOldPoints_,
4261 nOldFaces_,
4262 nOldCells_,
4263 pointMap_,
4264 pointsFromPoints_,
4265 faceMap_,
4266 facesFromPoints_,
4267 facesFromEdges_,
4268 facesFromFaces_,
4269 cellMap_,
4270 cellsFromPoints_,
4271 cellsFromEdges_,
4272 cellsFromFaces_,
4273 cellsFromCells_,
4274 reversePointMap_,
4275 reverseFaceMap_,
4276 reverseCellMap_,
4277 flipFaces_,
4278 patchPointMap,
4279 pointZoneMap,
4280 faceZonePointMap,
4281 faceZoneFaceMap,
4282 cellZoneMap,
4283 preMotionPoints,
4284 oldPatchStarts,
4285 oldPatchNMeshPoints,
4286 true
4287 );
4289 // Update the underlying mesh, and map all related fields
4290 updateMesh(mpm);
4292 if (mpm.hasMotionPoints())
4293 {
4294 // Perform a dummy movePoints to force V0 creation
4295 movePoints(mpm.preMotionPoints());
4297 // Reset old-volumes
4298 resetMotion();
4299 setV0();
4300 }
4302 // Correct volume fluxes on the old mesh
4303 fieldMapper.correctFluxes();
4305 // Clear mapper after use
4306 fieldMapper.clear();
4308 if (mpm.hasMotionPoints())
4309 {
4310 // Now move mesh to new points and
4311 // compute correct mesh-fluxes.
4312 movePoints(points);
4313 }
4315 // Update the mesh-motion solver
4316 if (motionSolver_.valid())
4317 {
4318 motionSolver_->updateMesh(mpm);
4319 }
4321 // Clear unwanted member data
4322 addedFacePatches_.clear();
4323 addedEdgePatches_.clear();
4324 addedPointZones_.clear();
4325 addedFaceZones_.clear();
4326 addedCellZones_.clear();
4327 faceParents_.clear();
4328 cellParents_.clear();
4330 // Clear the deleted entity map
4331 deletedPoints_.clear();
4332 deletedEdges_.clear();
4333 deletedFaces_.clear();
4334 deletedCells_.clear();
4336 // Clear flipFaces
4337 flipFaces_.clear();
4339 // Clear reverse maps
4340 reversePointMap_.clear();
4341 reverseEdgeMap_.clear();
4342 reverseFaceMap_.clear();
4343 reverseCellMap_.clear();
4345 // Update "old" information
4346 nOldPoints_ = nPoints_;
4347 nOldEdges_ = nEdges_;
4348 nOldFaces_ = nFaces_;
4349 nOldCells_ = nCells_;
4350 nOldInternalFaces_ = nInternalFaces_;
4351 nOldInternalEdges_ = nInternalEdges_;
4353 for (label i = 0; i < nPatches_; i++)
4354 {
4355 oldPatchSizes_[i] = patchSizes_[i];
4356 oldPatchStarts_[i] = patchStarts_[i];
4357 oldEdgePatchSizes_[i] = edgePatchSizes_[i];
4358 oldEdgePatchStarts_[i] = edgePatchStarts_[i];
4359 oldPatchNMeshPoints_[i] = patchNMeshPoints_[i];
4360 }
4362 // Clear parallel structures
4363 if (Pstream::parRun())
4364 {
4365 procIndices_.clear();
4366 procPriority_.clear();
4367 sendMeshes_.clear();
4368 recvMeshes_.clear();
4369 }
4371 bool checkCplBoundaries = false;
4373 if (meshSubDict.found("checkCoupledBoundaries") || mandatory_)
4374 {
4375 checkCplBoundaries =
4376 (
4377 readBool(meshSubDict.lookup("checkCoupledBoundaries"))
4378 );
4379 }
4381 if (checkCplBoundaries)
4382 {
4383 bool failed = checkCoupledBoundaries();
4385 if (failed)
4386 {
4387 FatalErrorIn("bool dynamicTopoFvMesh::resetMesh()")
4388 << " Coupled boundary check failed on processor: "
4389 << Pstream::myProcNo()
4390 << abort(FatalError);
4391 }
4392 }
4394 // Reset statistics
4395 statistics_ = 0;
4396 }
4397 else
4398 {
4399 // No topology changes were made.
4400 // Only execute mesh-motion.
4401 if (motionSolver_.valid())
4402 {
4403 movePoints(points_);
4404 }
4405 }
4407 // Obtain mesh stats after topo-changes
4408 meshQuality(true);
4410 // Dump length-scale to disk, if requested.
4411 calculateLengthScale(true);
4413 // Optionally check if decomposition is to be written out
4414 if
4415 (
4416 Pstream::parRun() &&
4417 time().outputTime() &&
4418 meshSubDict.found("writeDecomposition") &&
4419 readBool(meshSubDict.lookup("writeDecomposition"))
4420 )
4421 {
4422 volScalarField
4423 (
4424 IOobject
4425 (
4426 "decomposition",
4427 time().timeName(),
4428 *this,
4429 IOobject::NO_READ,
4430 IOobject::NO_WRITE,
4431 false
4432 ),
4433 *this,
4434 dimensionedScalar("rank", dimless, Pstream::myProcNo())
4435 ).write();
4436 }
4438 // Reset and return flag
4439 if (topoChangeFlag_)
4440 {
4441 topoChangeFlag_ = false;
4442 return true;
4443 }
4445 // No changes were made.
4446 return false;
4447 }
4450 // Return custom lduAddressing
4451 const lduAddressing& dynamicTopoFvMesh::lduAddr() const
4452 {
4453 // For subMeshes, avoid globalData construction
4454 // due to lduAddressing::patchSchedule, using a
4455 // customized approach.
4456 if (isSubMesh_)
4457 {
4458 if (!lduPtr_)
4459 {
4460 lduPtr_ = new subMeshLduAddressing(*this);
4461 }
4463 return *lduPtr_;
4464 }
4466 return fvMesh::lduAddr();
4467 }
4470 // Update mesh corresponding to the given map
4471 void dynamicTopoFvMesh::updateMesh(const mapPolyMesh& mpm)
4472 {
4473 if (coupledModification_)
4474 {
4475 // This bit gets called only during the load-balancing
4476 // stage, since the fvMesh::updateMesh is a bit different
4477 fvMesh::updateMesh(mpm);
4478 return;
4479 }
4481 // Delete oldPoints in polyMesh
4482 polyMesh::resetMotion();
4484 // Update fvMesh and polyMesh, and map all fields
4485 fvMesh::updateMesh(mpm);
4486 }
4489 // Map all fields in time using a customized mapper
4490 void dynamicTopoFvMesh::mapFields(const mapPolyMesh& mpm) const
4491 {
4492 if (coupledModification_)
4493 {
4494 // This bit gets called only during the load-balancing
4495 // stage, since the fvMesh::mapFields is a bit different
4496 fvMesh::mapFields(mpm);
4497 return;
4498 }
4500 if (debug)
4501 {
4502 Info<< "void dynamicTopoFvMesh::mapFields(const mapPolyMesh&) const: "
4503 << "Mapping fv fields."
4504 << endl;
4505 }
4507 const topoMapper& fieldMapper = mapper_();
4509 // Set the mapPolyMesh object in the mapper
4510 fieldMapper.setMapper(mpm);
4512 // Deregister gradient fields and centres,
4513 // but retain for mapping
4514 fieldMapper.deregisterMeshInformation();
4516 // Conservatively map scalar/vector volFields
4517 conservativeMapVolFields<scalar>(fieldMapper);
4518 conservativeMapVolFields<vector>(fieldMapper);
4520 // Map all the volFields in the objectRegistry
4521 MapGeometricFields<sphericalTensor,fvPatchField,topoMapper,volMesh>
4522 (fieldMapper);
4523 MapGeometricFields<symmTensor,fvPatchField,topoMapper,volMesh>
4524 (fieldMapper);
4525 MapGeometricFields<tensor,fvPatchField,topoMapper,volMesh>
4526 (fieldMapper);
4528 // Conservatively map scalar/vector surfaceFields
4529 conservativeMapSurfaceFields<scalar>(fieldMapper);
4530 conservativeMapSurfaceFields<vector>(fieldMapper);
4532 // Map all the surfaceFields in the objectRegistry
4533 MapGeometricFields<sphericalTensor,fvsPatchField,topoMapper,surfaceMesh>
4534 (fieldMapper);
4535 MapGeometricFields<symmTensor,fvsPatchField,topoMapper,surfaceMesh>
4536 (fieldMapper);
4537 MapGeometricFields<tensor,fvsPatchField,topoMapper,surfaceMesh>
4538 (fieldMapper);
4539 }
4542 // Update the mesh for motion / topology changes.
4543 // - Return true if topology changes have occurred
4544 bool dynamicTopoFvMesh::update()
4545 {
4546 // Re-read options, in case they have been modified at run-time
4547 readOptionalParameters(true);
4549 // Set old point positions
4550 oldPoints_ = polyMesh::points();
4552 // Invoke mesh-motion solver and store new points
4553 if (motionSolver_.valid())
4554 {
4555 points_ = motionSolver_->newPoints()();
4556 }
4557 else
4558 {
4559 // Set point positions from mesh
4560 points_ = polyMesh::points();
4561 }
4563 // Obtain mesh stats before topo-changes
4564 bool noSlivers = meshQuality(true);
4566 // Return if the interval is invalid,
4567 // not at re-mesh interval, or slivers are absent.
4568 // Handy while using only mesh-motion.
4569 if (interval_ < 0 || ((time().timeIndex() % interval_ != 0) && noSlivers))
4570 {
4571 return resetMesh();
4572 }
4574 // Calculate the edge length-scale for the mesh
4575 calculateLengthScale();
4577 // Track mesh topology modification time
4578 clockTime topoTimer;
4580 // Invoke the threaded topoModifier
4581 threadedTopoModifier();
4583 Info<< " Topo modifier time: "
4584 << topoTimer.elapsedTime() << " s"
4585 << endl;
4587 // Apply all topology changes (if any) and reset mesh.
4588 return resetMesh();
4589 }
4592 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
4594 void dynamicTopoFvMesh::operator=(const dynamicTopoFvMesh& rhs)
4595 {
4596 // Check for assignment to self
4597 if (this == &rhs)
4598 {
4599 FatalErrorIn
4600 (
4601 "dynamicTopoFvMesh::operator=(const dynamicTopoFvMesh&)"
4602 )
4603 << "Attempted assignment to self"
4604 << abort(FatalError);
4605 }
4606 }
4609 } // End namespace Foam
4611 // ************************************************************************* //