diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverLayers.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverLayers.C
deleted file mode 100644
index 3ebf41390cbb3c52986e2b41bdb4e8647b06ee3e..0000000000000000000000000000000000000000
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverLayers.C
+++ /dev/null
@@ -1,2961 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software; you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by the
- Free Software Foundation; either version 2 of the License, or (at your
- option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM; if not, write to the Free Software Foundation,
- Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
- All to do with adding cell layers
-
-\*----------------------------------------------------------------------------*/
-
-#include "autoHexMeshDriver.H"
-#include "fvMesh.H"
-#include "Time.H"
-#include "combineFaces.H"
-#include "removePoints.H"
-#include "pointFields.H"
-#include "motionSmoother.H"
-#include "mathematicalConstants.H"
-#include "pointSet.H"
-#include "faceSet.H"
-#include "cellSet.H"
-#include "polyTopoChange.H"
-#include "mapPolyMesh.H"
-#include "addPatchCellLayer.H"
-#include "mapDistributePolyMesh.H"
-#include "OFstream.H"
-#include "layerParameters.H"
-
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
-Foam::label Foam::autoHexMeshDriver::mergePatchFacesUndo
-(
- const scalar minCos,
- const scalar concaveCos,
- const dictionary& motionDict
-)
-{
- // Patch face merging engine
- combineFaces faceCombiner(mesh_, true);
-
- // Pick up all candidate cells on boundary
- labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
-
- {
- labelList patchIDs(meshRefinement::addedPatches(globalToPatch_));
-
- const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
- forAll(patchIDs, i)
- {
- label patchI = patchIDs[i];
-
- const polyPatch& patch = patches[patchI];
-
- if (!patch.coupled())
- {
- forAll(patch, i)
- {
- boundaryCells.insert(mesh_.faceOwner()[patch.start()+i]);
- }
- }
- }
- }
-
- // Get all sets of faces that can be merged
- labelListList allFaceSets
- (
- faceCombiner.getMergeSets
- (
- minCos,
- concaveCos,
- boundaryCells
- )
- );
-
- label nFaceSets = returnReduce(allFaceSets.size(), sumOp<label>());
-
- Info<< "Merging " << nFaceSets << " sets of faces." << nl << endl;
-
- if (nFaceSets > 0)
- {
- if (debug_)
- {
- faceSet allSets(mesh_, "allFaceSets", allFaceSets.size());
- forAll(allFaceSets, setI)
- {
- forAll(allFaceSets[setI], i)
- {
- allSets.insert(allFaceSets[setI][i]);
- }
- }
- Pout<< "Writing all faces to be merged to set "
- << allSets.objectPath() << endl;
- allSets.write();
- }
-
-
- // Topology changes container
- polyTopoChange meshMod(mesh_);
-
- // Merge all faces of a set into the first face of the set.
- faceCombiner.setRefinement(allFaceSets, meshMod);
-
- // Experimental: store data for all the points that have been deleted
- meshRefinerPtr_().storeData
- (
- faceCombiner.savedPointLabels(), // points to store
- labelList(0), // faces to store
- labelList(0) // cells to store
- );
-
- // Change the mesh (no inflation)
- autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
-
- // Update fields
- mesh_.updateMesh(map);
-
- // Move mesh (since morphing does not do this)
- if (map().hasMotionPoints())
- {
- mesh_.movePoints(map().preMotionPoints());
- }
-
- faceCombiner.updateMesh(map);
-
- meshRefinerPtr_().updateMesh(map, labelList(0));
-
-
-
- for (label iteration = 0; iteration < 100; iteration++)
- {
- Info<< nl
- << "Undo iteration " << iteration << nl
- << "----------------" << endl;
-
-
- // Check mesh for errors
- // ~~~~~~~~~~~~~~~~~~~~~
-
- faceSet errorFaces
- (
- mesh_,
- "errorFaces",
- mesh_.nFaces()-mesh_.nInternalFaces()
- );
- bool hasErrors = motionSmoother::checkMesh
- (
- false, // report
- mesh_,
- motionDict,
- errorFaces
- );
- //if (checkEdgeConnectivity)
- //{
- // Info<< "Checking edge-face connectivity (duplicate faces"
- // << " or non-consecutive shared vertices)" << endl;
- //
- // label nOldSize = errorFaces.size();
- //
- // hasErrors =
- // mesh_.checkFaceFaces
- // (
- // false,
- // &errorFaces
- // )
- // || hasErrors;
- //
- // Info<< "Detected additional "
- // << returnReduce(errorFaces.size()-nOldSize, sumOp<label>())
- // << " faces with illegal face-face connectivity"
- // << endl;
- //}
-
- if (!hasErrors)
- {
- break;
- }
-
-
- if (debug_)
- {
- Pout<< "Writing all faces in error to faceSet "
- << errorFaces.objectPath() << nl << endl;
- errorFaces.write();
- }
-
-
- // Check any master cells for using any of the error faces
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- DynamicList<label> mastersToRestore(allFaceSets.size());
-
- forAll(allFaceSets, setI)
- {
- label masterFaceI = faceCombiner.masterFace()[setI];
-
- if (masterFaceI != -1)
- {
- label masterCellII = mesh_.faceOwner()[masterFaceI];
-
- const cell& cFaces = mesh_.cells()[masterCellII];
-
- forAll(cFaces, i)
- {
- if (errorFaces.found(cFaces[i]))
- {
- mastersToRestore.append(masterFaceI);
- break;
- }
- }
- }
- }
- mastersToRestore.shrink();
-
- label nRestore = returnReduce
- (
- mastersToRestore.size(),
- sumOp<label>()
- );
-
- Info<< "Masters that need to be restored:"
- << nRestore << endl;
-
- if (debug_)
- {
- faceSet restoreSet(mesh_, "mastersToRestore", mastersToRestore);
- Pout<< "Writing all " << mastersToRestore.size()
- << " masterfaces to be restored to set "
- << restoreSet.objectPath() << endl;
- restoreSet.write();
- }
-
-
- if (nRestore == 0)
- {
- break;
- }
-
-
- // Undo
- // ~~~~
-
- // Topology changes container
- polyTopoChange meshMod(mesh_);
-
- // Merge all faces of a set into the first face of the set.
- // Experimental:mark all points/faces/cells that have been restored.
- Map<label> restoredPoints(0);
- Map<label> restoredFaces(0);
- Map<label> restoredCells(0);
-
- faceCombiner.setUnrefinement
- (
- mastersToRestore,
- meshMod,
- restoredPoints,
- restoredFaces,
- restoredCells
- );
-
- // Change the mesh (no inflation)
- autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
-
- // Update fields
- mesh_.updateMesh(map);
-
- // Move mesh (since morphing does not do this)
- if (map().hasMotionPoints())
- {
- mesh_.movePoints(map().preMotionPoints());
- }
-
- faceCombiner.updateMesh(map);
-
- // Renumber restore maps
- inplaceMapKey(map().reversePointMap(), restoredPoints);
- inplaceMapKey(map().reverseFaceMap(), restoredFaces);
- inplaceMapKey(map().reverseCellMap(), restoredCells);
-
- // Experimental:restore all points/face/cells in maps
- meshRefinerPtr_().updateMesh
- (
- map,
- labelList(0), // changedFaces
- restoredPoints,
- restoredFaces,
- restoredCells
- );
-
- Info<< endl;
- }
-
- if (debug_)
- {
- Pout<< "Writing merged-faces mesh to time "
- << mesh_.time().timeName() << nl << endl;
- mesh_.write();
- }
- }
- else
- {
- Info<< "No faces merged ..." << endl;
- }
-
- return nFaceSets;
-}
-
-
-// Remove points. pointCanBeDeleted is parallel synchronised.
-Foam::autoPtr<Foam::mapPolyMesh> Foam::autoHexMeshDriver::doRemovePoints
-(
- removePoints& pointRemover,
- const boolList& pointCanBeDeleted
-)
-{
- // Topology changes container
- polyTopoChange meshMod(mesh_);
-
- pointRemover.setRefinement(pointCanBeDeleted, meshMod);
-
- // Change the mesh (no inflation)
- autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
-
- // Update fields
- mesh_.updateMesh(map);
-
- // Move mesh (since morphing does not do this)
- if (map().hasMotionPoints())
- {
- mesh_.movePoints(map().preMotionPoints());
- }
-
- pointRemover.updateMesh(map);
- meshRefinerPtr_().updateMesh(map, labelList(0));
-
- return map;
-}
-
-
-// Restore faces (which contain removed points)
-Foam::autoPtr<Foam::mapPolyMesh> Foam::autoHexMeshDriver::doRestorePoints
-(
- removePoints& pointRemover,
- const labelList& facesToRestore
-)
-{
- // Topology changes container
- polyTopoChange meshMod(mesh_);
-
- // Determine sets of points and faces to restore
- labelList localFaces, localPoints;
- pointRemover.getUnrefimentSet
- (
- facesToRestore,
- localFaces,
- localPoints
- );
-
- // Undo the changes on the faces that are in error.
- pointRemover.setUnrefinement
- (
- localFaces,
- localPoints,
- meshMod
- );
-
- // Change the mesh (no inflation)
- autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
-
- // Update fields
- mesh_.updateMesh(map);
-
- // Move mesh (since morphing does not do this)
- if (map().hasMotionPoints())
- {
- mesh_.movePoints(map().preMotionPoints());
- }
-
- pointRemover.updateMesh(map);
- meshRefinerPtr_().updateMesh(map, labelList(0));
-
- return map;
-}
-
-
-// Collect all faces that are both in candidateFaces and in the set.
-// If coupled face also collects the coupled face.
-Foam::labelList Foam::autoHexMeshDriver::collectFaces
-(
- const labelList& candidateFaces,
- const labelHashSet& set
-) const
-{
- // Has face been selected?
- boolList selected(mesh_.nFaces(), false);
-
- forAll(candidateFaces, i)
- {
- label faceI = candidateFaces[i];
-
- if (set.found(faceI))
- {
- selected[faceI] = true;
- }
- }
- syncTools::syncFaceList
- (
- mesh_,
- selected,
- orEqOp<bool>(), // combine operator
- false // separation
- );
-
- labelList selectedFaces(findIndices(selected, true));
-
- return selectedFaces;
-}
-
-
-// Pick up faces of cells of faces in set.
-Foam::labelList Foam::autoHexMeshDriver::growFaceCellFace
-(
- const labelHashSet& set
-) const
-{
- boolList selected(mesh_.nFaces(), false);
-
- forAllConstIter(faceSet, set, iter)
- {
- label faceI = iter.key();
-
- label own = mesh_.faceOwner()[faceI];
-
- const cell& ownFaces = mesh_.cells()[own];
- forAll(ownFaces, ownFaceI)
- {
- selected[ownFaces[ownFaceI]] = true;
- }
-
- if (mesh_.isInternalFace(faceI))
- {
- label nbr = mesh_.faceNeighbour()[faceI];
-
- const cell& nbrFaces = mesh_.cells()[nbr];
- forAll(nbrFaces, nbrFaceI)
- {
- selected[nbrFaces[nbrFaceI]] = true;
- }
- }
- }
- syncTools::syncFaceList
- (
- mesh_,
- selected,
- orEqOp<bool>(), // combine operator
- false // separation
- );
- return findIndices(selected, true);
-}
-
-
-// Remove points not used by any face or points used by only two faces where
-// the edges are in line
-Foam::label Foam::autoHexMeshDriver::mergeEdgesUndo
-(
- const scalar minCos,
- const dictionary& motionDict
-)
-{
- Info<< nl
- << "Merging all points on surface that" << nl
- << "- are used by only two boundary faces and" << nl
- << "- make an angle with a cosine of more than " << minCos
- << "." << nl << endl;
-
- // Point removal analysis engine with undo
- removePoints pointRemover(mesh_, true);
-
- // Count usage of points
- boolList pointCanBeDeleted;
- label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
-
- if (nRemove > 0)
- {
- Info<< "Removing " << nRemove
- << " straight edge points ..." << nl << endl;
-
- // Remove points
- // ~~~~~~~~~~~~~
-
- doRemovePoints(pointRemover, pointCanBeDeleted);
-
-
- for (label iteration = 0; iteration < 100; iteration++)
- {
- Info<< nl
- << "Undo iteration " << iteration << nl
- << "----------------" << endl;
-
-
- // Check mesh for errors
- // ~~~~~~~~~~~~~~~~~~~~~
-
- faceSet errorFaces
- (
- mesh_,
- "errorFaces",
- mesh_.nFaces()-mesh_.nInternalFaces()
- );
- bool hasErrors = motionSmoother::checkMesh
- (
- false, // report
- mesh_,
- motionDict,
- errorFaces
- );
- //if (checkEdgeConnectivity)
- //{
- // Info<< "Checking edge-face connectivity (duplicate faces"
- // << " or non-consecutive shared vertices)" << endl;
- //
- // label nOldSize = errorFaces.size();
- //
- // hasErrors =
- // mesh_.checkFaceFaces
- // (
- // false,
- // &errorFaces
- // )
- // || hasErrors;
- //
- // Info<< "Detected additional "
- // << returnReduce(errorFaces.size()-nOldSize, sumOp<label>())
- // << " faces with illegal face-face connectivity"
- // << endl;
- //}
-
- if (!hasErrors)
- {
- break;
- }
-
- if (debug_)
- {
- Pout<< "**Writing all faces in error to faceSet "
- << errorFaces.objectPath() << nl << endl;
- errorFaces.write();
- }
-
- labelList masterErrorFaces
- (
- collectFaces
- (
- pointRemover.savedFaceLabels(),
- errorFaces
- )
- );
-
- label n = returnReduce(masterErrorFaces.size(), sumOp<label>());
-
- Info<< "Detected " << n
- << " error faces on boundaries that have been merged."
- << " These will be restored to their original faces." << nl
- << endl;
-
- if (n == 0)
- {
- if (hasErrors)
- {
- Info<< "Detected "
- << returnReduce(errorFaces.size(), sumOp<label>())
- << " error faces in mesh."
- << " Restoring neighbours of faces in error." << nl
- << endl;
-
- labelList expandedErrorFaces
- (
- growFaceCellFace
- (
- errorFaces
- )
- );
-
- doRestorePoints(pointRemover, expandedErrorFaces);
- }
-
- break;
- }
-
- doRestorePoints(pointRemover, masterErrorFaces);
- }
-
- if (debug_)
- {
- Pout<< "Writing merged-edges mesh to time "
- << mesh_.time().timeName() << nl << endl;
- mesh_.write();
- }
- }
- else
- {
- Info<< "No straight edges simplified and no points removed ..." << endl;
- }
-
- return nRemove;
-}
-
-
-// For debugging: Dump displacement to .obj files
-void Foam::autoHexMeshDriver::dumpDisplacement
-(
- const fileName& prefix,
- const indirectPrimitivePatch& pp,
- const vectorField& patchDisp,
- const List<extrudeMode>& extrudeStatus
-)
-{
- OFstream dispStr(prefix + "_disp.obj");
- Info<< "Writing all displacements to " << dispStr.name() << nl << endl;
-
- label vertI = 0;
-
- forAll(patchDisp, patchPointI)
- {
- const point& pt = pp.localPoints()[patchPointI];
-
- meshTools::writeOBJ(dispStr, pt); vertI++;
- meshTools::writeOBJ(dispStr, pt + patchDisp[patchPointI]); vertI++;
-
- dispStr << "l " << vertI-1 << ' ' << vertI << nl;
- }
-
-
- OFstream illStr(prefix + "_illegal.obj");
- Info<< "Writing invalid displacements to " << illStr.name() << nl << endl;
-
- vertI = 0;
-
- forAll(patchDisp, patchPointI)
- {
- if (extrudeStatus[patchPointI] != EXTRUDE)
- {
- const point& pt = pp.localPoints()[patchPointI];
-
- meshTools::writeOBJ(illStr, pt); vertI++;
- meshTools::writeOBJ(illStr, pt + patchDisp[patchPointI]); vertI++;
-
- illStr << "l " << vertI-1 << ' ' << vertI << nl;
- }
- }
-}
-
-
-// Check that primitivePatch is not multiply connected. Collect non-manifold
-// points in pointSet.
-void Foam::autoHexMeshDriver::checkManifold
-(
- const indirectPrimitivePatch& fp,
- pointSet& nonManifoldPoints
-)
-{
- // Check for non-manifold points (surface pinched at point)
- fp.checkPointManifold(false, &nonManifoldPoints);
-
- // Check for edge-faces (surface pinched at edge)
- const labelListList& edgeFaces = fp.edgeFaces();
-
- forAll(edgeFaces, edgeI)
- {
- const labelList& eFaces = edgeFaces[edgeI];
-
- if (eFaces.size() > 2)
- {
- const edge& e = fp.edges()[edgeI];
-
- nonManifoldPoints.insert(fp.meshPoints()[e[0]]);
- nonManifoldPoints.insert(fp.meshPoints()[e[1]]);
- }
- }
-}
-
-
-void Foam::autoHexMeshDriver::checkMeshManifold() const
-{
- Info<< nl << "Checking mesh manifoldness ..." << endl;
-
- // Get all outside faces
- labelList outsideFaces(mesh_.nFaces() - mesh_.nInternalFaces());
-
- for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
- {
- outsideFaces[faceI - mesh_.nInternalFaces()] = faceI;
- }
-
- pointSet nonManifoldPoints
- (
- mesh_,
- "nonManifoldPoints",
- mesh_.nPoints() / 100
- );
-
- // Build primitivePatch out of faces and check it for problems.
- checkManifold
- (
- indirectPrimitivePatch
- (
- IndirectList<face>(mesh_.faces(), outsideFaces),
- mesh_.points()
- ),
- nonManifoldPoints
- );
-
- label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
-
- if (nNonManif > 0)
- {
- Info<< "Outside of mesh is multiply connected across edges or"
- << " points." << nl
- << "This is not a fatal error but might cause some unexpected"
- << " behaviour." << nl
- << "Writing " << nNonManif
- << " points where this happens to pointSet "
- << nonManifoldPoints.name() << endl;
-
- nonManifoldPoints.write();
- }
- Info<< endl;
-}
-
-
-
-// Unset extrusion on point. Returns true if anything unset.
-bool Foam::autoHexMeshDriver::unmarkExtrusion
-(
- const label patchPointI,
- pointField& patchDisp,
- labelList& patchNLayers,
- List<extrudeMode>& extrudeStatus
-)
-{
- if (extrudeStatus[patchPointI] == EXTRUDE)
- {
- extrudeStatus[patchPointI] = NOEXTRUDE;
- patchNLayers[patchPointI] = 0;
- patchDisp[patchPointI] = vector::zero;
- return true;
- }
- else if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
- {
- extrudeStatus[patchPointI] = NOEXTRUDE;
- patchNLayers[patchPointI] = 0;
- patchDisp[patchPointI] = vector::zero;
- return true;
- }
- else
- {
- return false;
- }
-}
-
-
-// Unset extrusion on face. Returns true if anything unset.
-bool Foam::autoHexMeshDriver::unmarkExtrusion
-(
- const face& localFace,
- pointField& patchDisp,
- labelList& patchNLayers,
- List<extrudeMode>& extrudeStatus
-)
-{
- bool unextruded = false;
-
- forAll(localFace, fp)
- {
- if
- (
- unmarkExtrusion
- (
- localFace[fp],
- patchDisp,
- patchNLayers,
- extrudeStatus
- )
- )
- {
- unextruded = true;
- }
- }
- return unextruded;
-}
-
-
-// No extrusion at non-manifold points.
-void Foam::autoHexMeshDriver::handleNonManifolds
-(
- const indirectPrimitivePatch& pp,
- const labelList& meshEdges,
- pointField& patchDisp,
- labelList& patchNLayers,
- List<extrudeMode>& extrudeStatus
-) const
-{
- Info<< nl << "Handling non-manifold points ..." << endl;
-
- // Detect non-manifold points
- Info<< nl << "Checking patch manifoldness ..." << endl;
-
- pointSet nonManifoldPoints(mesh_, "nonManifoldPoints", pp.nPoints());
-
- // 1. Local check
- checkManifold(pp, nonManifoldPoints);
-
- label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
-
- Info<< "Outside of local patch is multiply connected across edges or"
- << " points at " << nNonManif << " points." << endl;
-
-
- const labelList& meshPoints = pp.meshPoints();
-
-
- // 2. Parallel check
- // For now disable extrusion at any shared edges.
- const labelHashSet sharedEdgeSet(mesh_.globalData().sharedEdgeLabels());
-
- forAll(pp.edges(), edgeI)
- {
- if (sharedEdgeSet.found(meshEdges[edgeI]))
- {
- const edge& e = mesh_.edges()[meshEdges[edgeI]];
-
- Pout<< "Disabling extrusion at edge "
- << mesh_.points()[e[0]] << mesh_.points()[e[1]]
- << " since it is non-manifold across coupled patches."
- << endl;
-
- nonManifoldPoints.insert(e[0]);
- nonManifoldPoints.insert(e[1]);
- }
- }
-
- // 3b. extrusion can produce multiple faces between 2 cells
- // across processor boundary
- // This occurs when a coupled face shares more than 1 edge with a
- // non-coupled boundary face.
- // This is now correctly handled by addPatchCellLayer in that it
- // extrudes a single face from the stringed up edges.
-
-
- nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
-
- if (nNonManif > 0)
- {
- Info<< "Outside of patches is multiply connected across edges or"
- << " points at " << nNonManif << " points." << nl
- << "Writing " << nNonManif
- << " points where this happens to pointSet "
- << nonManifoldPoints.name() << nl
- << "and setting layer thickness to zero on these points."
- << endl;
-
- nonManifoldPoints.write();
-
- forAll(meshPoints, patchPointI)
- {
- if (nonManifoldPoints.found(meshPoints[patchPointI]))
- {
- unmarkExtrusion
- (
- patchPointI,
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
- }
- }
- }
-
- Info<< "Set displacement to zero for all " << nNonManif
- << " non-manifold points" << endl;
-}
-
-
-// Parallel feature edge detection. Assumes non-manifold edges already handled.
-void Foam::autoHexMeshDriver::handleFeatureAngle
-(
- const indirectPrimitivePatch& pp,
- const labelList& meshEdges,
- const scalar minCos,
- pointField& patchDisp,
- labelList& patchNLayers,
- List<extrudeMode>& extrudeStatus
-) const
-{
- Info<< nl << "Handling feature edges ..." << endl;
-
- if (minCos < 1-SMALL)
- {
- // Normal component of normals of connected faces.
- vectorField edgeNormal(mesh_.nEdges(), wallPoint::greatPoint);
-
- const labelListList& edgeFaces = pp.edgeFaces();
-
- forAll(edgeFaces, edgeI)
- {
- const labelList& eFaces = pp.edgeFaces()[edgeI];
-
- label meshEdgeI = meshEdges[edgeI];
-
- forAll(eFaces, i)
- {
- nomalsCombine()
- (
- edgeNormal[meshEdgeI],
- pp.faceNormals()[eFaces[i]]
- );
- }
- }
-
- syncTools::syncEdgeList
- (
- mesh_,
- edgeNormal,
- nomalsCombine(),
- wallPoint::greatPoint, // null value
- false // no separation
- );
-
- label vertI = 0;
- autoPtr<OFstream> str;
- if (debug_)
- {
- str.reset(new OFstream(mesh_.time().path()/"featureEdges.obj"));
- Info<< "Writing feature edges to " << str().name() << endl;
- }
-
- label nFeats = 0;
-
- // Now on coupled edges the edgeNormal will have been truncated and
- // only be still be the old value where two faces have the same normal
- forAll(edgeFaces, edgeI)
- {
- const labelList& eFaces = pp.edgeFaces()[edgeI];
-
- label meshEdgeI = meshEdges[edgeI];
-
- const vector& n = edgeNormal[meshEdgeI];
-
- if (n != wallPoint::greatPoint)
- {
- scalar cos = n & pp.faceNormals()[eFaces[0]];
-
- if (cos < minCos)
- {
- const edge& e = pp.edges()[edgeI];
-
- unmarkExtrusion
- (
- e[0],
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
- unmarkExtrusion
- (
- e[1],
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
-
- nFeats++;
-
- if (str.valid())
- {
- meshTools::writeOBJ(str(), pp.localPoints()[e[0]]);
- vertI++;
- meshTools::writeOBJ(str(), pp.localPoints()[e[1]]);
- vertI++;
- str()<< "l " << vertI-1 << ' ' << vertI << nl;
- }
- }
- }
- }
-
- Info<< "Set displacement to zero for points on "
- << returnReduce(nFeats, sumOp<label>())
- << " feature edges" << endl;
- }
-}
-
-
-// No extrusion on cells with warped faces. Calculates the thickness of the
-// layer and compares it to the space the warped face takes up. Disables
-// extrusion if layer thickness is more than faceRatio of the thickness of
-// the face.
-void Foam::autoHexMeshDriver::handleWarpedFaces
-(
- const indirectPrimitivePatch& pp,
- const scalar faceRatio,
- const scalar edge0Len,
- const labelList& cellLevel,
- pointField& patchDisp,
- labelList& patchNLayers,
- List<extrudeMode>& extrudeStatus
-) const
-{
- Info<< nl << "Handling cells with warped patch faces ..." << nl;
-
- const pointField& points = mesh_.points();
-
- label nWarpedFaces = 0;
-
- forAll(pp, i)
- {
- const face& f = pp[i];
-
- if (f.size() > 3)
- {
- label faceI = pp.addressing()[i];
-
- label ownLevel = cellLevel[mesh_.faceOwner()[faceI]];
- scalar edgeLen = edge0Len/(1<<ownLevel);
-
- // Normal distance to face centre plane
- const point& fc = mesh_.faceCentres()[faceI];
- const vector& fn = pp.faceNormals()[i];
-
- scalarField vProj(f.size());
-
- forAll(f, fp)
- {
- vector n = points[f[fp]] - fc;
- vProj[fp] = (n & fn);
- }
-
- // Get normal 'span' of face
- scalar minVal = min(vProj);
- scalar maxVal = max(vProj);
-
- if ((maxVal - minVal) > faceRatio * edgeLen)
- {
- if
- (
- unmarkExtrusion
- (
- pp.localFaces()[i],
- patchDisp,
- patchNLayers,
- extrudeStatus
- )
- )
- {
- nWarpedFaces++;
- }
- }
- }
- }
-
- Info<< "Set displacement to zero on "
- << returnReduce(nWarpedFaces, sumOp<label>())
- << " warped faces since layer would be > " << faceRatio
- << " of the size of the bounding box." << endl;
-}
-
-
-
-//// No extrusion on cells with multiple patch faces. There ususally is a reason
-//// why combinePatchFaces hasn't succeeded.
-//void Foam::autoHexMeshDriver::handleMultiplePatchFaces
-//(
-// const indirectPrimitivePatch& pp,
-// pointField& patchDisp,
-// labelList& patchNLayers,
-// List<extrudeMode>& extrudeStatus
-//) const
-//{
-// Info<< nl << "Handling cells with multiple patch faces ..." << nl;
-//
-// const labelListList& pointFaces = pp.pointFaces();
-//
-// // Cells that should not get an extrusion layer
-// cellSet multiPatchCells(mesh_, "multiPatchCells", pp.size());
-//
-// // Detect points that use multiple faces on same cell.
-// forAll(pointFaces, patchPointI)
-// {
-// const labelList& pFaces = pointFaces[patchPointI];
-//
-// labelHashSet pointCells(pFaces.size());
-//
-// forAll(pFaces, i)
-// {
-// label cellI = mesh_.faceOwner()[pp.addressing()[pFaces[i]]];
-//
-// if (!pointCells.insert(cellI))
-// {
-// // Second or more occurrence of cell so cell has two or more
-// // pp faces connected to this point.
-// multiPatchCells.insert(cellI);
-// }
-// }
-// }
-//
-// label nMultiPatchCells = returnReduce
-// (
-// multiPatchCells.size(),
-// sumOp<label>()
-// );
-//
-// Info<< "Detected " << nMultiPatchCells
-// << " cells with multiple (connected) patch faces." << endl;
-//
-// label nChanged = 0;
-//
-// if (nMultiPatchCells > 0)
-// {
-// Info<< "Writing " << nMultiPatchCells
-// << " cells with multiple (connected) patch faces to cellSet "
-// << multiPatchCells.objectPath() << endl;
-// multiPatchCells.write();
-//
-//
-// // Go through all points and remove extrusion on any cell in
-// // multiPatchCells
-// // (has to be done in separate loop since having one point on
-// // multipatches has to reset extrusion on all points of cell)
-//
-// forAll(pointFaces, patchPointI)
-// {
-// if (extrudeStatus[patchPointI] != NOEXTRUDE)
-// {
-// const labelList& pFaces = pointFaces[patchPointI];
-//
-// forAll(pFaces, i)
-// {
-// label cellI =
-// mesh_.faceOwner()[pp.addressing()[pFaces[i]]];
-//
-// if (multiPatchCells.found(cellI))
-// {
-// if
-// (
-// unmarkExtrusion
-// (
-// patchPointI,
-// patchDisp,
-// patchNLayers,
-// extrudeStatus
-// )
-// )
-// {
-// nChanged++;
-// }
-// }
-// }
-// }
-// }
-//
-// reduce(nChanged, sumOp<label>());
-// }
-//
-// Info<< "Prevented extrusion on " << nChanged
-// << " points due to multiple patch faces." << nl << endl;
-//}
-
-
-// No extrusion on faces with differing number of layers for points
-void Foam::autoHexMeshDriver::setNumLayers
-(
- const labelList& patchToNLayers,
- const labelList& patchIDs,
- const indirectPrimitivePatch& pp,
- pointField& patchDisp,
- labelList& patchNLayers,
- List<extrudeMode>& extrudeStatus
-) const
-{
- Info<< nl << "Handling points with inconsistent layer specification ..."
- << endl;
-
- // Get for every point (really only nessecary on patch external points)
- // the max and min of any patch faces using it.
- labelList maxLayers(patchNLayers.size(), labelMin);
- labelList minLayers(patchNLayers.size(), labelMax);
-
- forAll(patchIDs, i)
- {
- label patchI = patchIDs[i];
-
- const labelList& meshPoints = mesh_.boundaryMesh()[patchI].meshPoints();
-
- label wantedLayers = patchToNLayers[patchI];
-
- forAll(meshPoints, patchPointI)
- {
- label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
-
- maxLayers[ppPointI] = max(wantedLayers, maxLayers[ppPointI]);
- minLayers[ppPointI] = min(wantedLayers, minLayers[ppPointI]);
- }
- }
-
- syncTools::syncPointList
- (
- mesh_,
- pp.meshPoints(),
- maxLayers,
- maxEqOp<label>(),
- labelMin, // null value
- false // no separation
- );
- syncTools::syncPointList
- (
- mesh_,
- pp.meshPoints(),
- minLayers,
- minEqOp<label>(),
- labelMax, // null value
- false // no separation
- );
-
- // Unmark any point with different min and max
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- //label nConflicts = 0;
-
- forAll(maxLayers, i)
- {
- if (maxLayers[i] == labelMin || minLayers[i] == labelMax)
- {
- FatalErrorIn("setNumLayers(..)")
- << "Patchpoint:" << i << " coord:" << pp.localPoints()[i]
- << " maxLayers:" << maxLayers
- << " minLayers:" << minLayers
- << abort(FatalError);
- }
- else if (maxLayers[i] == minLayers[i])
- {
- // Ok setting.
- patchNLayers[i] = maxLayers[i];
- }
- else
- {
- // Inconsistent num layers between patch faces using point
- //if
- //(
- // unmarkExtrusion
- // (
- // i,
- // patchDisp,
- // patchNLayers,
- // extrudeStatus
- // )
- //)
- //{
- // nConflicts++;
- //}
- patchNLayers[i] = maxLayers[i];
- }
- }
-
- //reduce(nConflicts, sumOp<label>());
- //
- //Info<< "Set displacement to zero for " << nConflicts
- // << " points due to points being on multiple regions"
- // << " with inconsistent nLayers specification." << endl;
-}
-
-
-// Grow no-extrusion layer
-void Foam::autoHexMeshDriver::growNoExtrusion
-(
- const indirectPrimitivePatch& pp,
- pointField& patchDisp,
- labelList& patchNLayers,
- List<extrudeMode>& extrudeStatus
-)
-{
- Info<< nl << "Growing non-extrusion points by one layer ..." << endl;
-
- List<extrudeMode> grownExtrudeStatus(extrudeStatus);
-
- const faceList& localFaces = pp.localFaces();
-
- label nGrown = 0;
-
- forAll(localFaces, faceI)
- {
- const face& f = localFaces[faceI];
-
- bool hasSqueeze = false;
- forAll(f, fp)
- {
- if (extrudeStatus[f[fp]] == NOEXTRUDE)
- {
- hasSqueeze = true;
- break;
- }
- }
-
- if (hasSqueeze)
- {
- // Squeeze all points of face
- forAll(f, fp)
- {
- if
- (
- extrudeStatus[f[fp]] == NOEXTRUDE
- && grownExtrudeStatus[f[fp]] != NOEXTRUDE
- )
- {
- grownExtrudeStatus[f[fp]] = NOEXTRUDE;
- nGrown++;
- }
- }
- }
- }
-
- extrudeStatus.transfer(grownExtrudeStatus);
-
- forAll(extrudeStatus, patchPointI)
- {
- if (extrudeStatus[patchPointI] == NOEXTRUDE)
- {
- patchDisp[patchPointI] = vector::zero;
- patchNLayers[patchPointI] = 0;
- }
- }
-
- reduce(nGrown, sumOp<label>());
-
- Info<< "Set displacement to zero for an additional " << nGrown
- << " points." << endl;
-}
-
-
-void Foam::autoHexMeshDriver::calculateLayerThickness
-(
- const indirectPrimitivePatch& pp,
- const labelList& patchIDs,
- const scalarField& patchExpansionRatio,
- const scalarField& patchFinalLayerRatio,
- const scalarField& patchRelMinThickness,
- const labelList& cellLevel,
- const labelList& patchNLayers,
- const scalar edge0Len,
-
- scalarField& thickness,
- scalarField& minThickness,
- scalarField& expansionRatio
-) const
-{
- const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
- if (min(patchRelMinThickness) < 0 || max(patchRelMinThickness) > 2)
- {
- FatalErrorIn("calculateLayerThickness(..)")
- << "Thickness should be factor of local undistorted cell size."
- << " Valid values are [0..2]." << nl
- << " minThickness:" << patchRelMinThickness
- << exit(FatalError);
- }
-
-
- // Per point the max cell level of connected cells
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- labelList maxPointLevel(pp.nPoints(), labelMin);
-
- {
- forAll(pp, i)
- {
- label ownLevel = cellLevel[mesh_.faceOwner()[pp.addressing()[i]]];
-
- const face& f = pp.localFaces()[i];
-
- forAll(f, fp)
- {
- maxPointLevel[f[fp]] = max(maxPointLevel[f[fp]], ownLevel);
- }
- }
-
- syncTools::syncPointList
- (
- mesh_,
- pp.meshPoints(),
- maxPointLevel,
- maxEqOp<label>(),
- labelMin, // null value
- false // no separation
- );
- }
-
-
- // Rework patch-wise layer parameters into minimum per point
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- expansionRatio.setSize(pp.nPoints());
- expansionRatio = GREAT;
- scalarField finalLayerRatio(pp.nPoints(), GREAT);
- scalarField relMinThickness(pp.nPoints(), GREAT);
-
- {
- forAll(patchIDs, i)
- {
- label patchI = patchIDs[i];
-
- const labelList& meshPoints = patches[patchI].meshPoints();
-
- forAll(meshPoints, patchPointI)
- {
- label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
-
- expansionRatio[ppPointI] = min
- (
- expansionRatio[ppPointI],
- patchExpansionRatio[patchI]
- );
- finalLayerRatio[ppPointI] = min
- (
- finalLayerRatio[ppPointI],
- patchFinalLayerRatio[patchI]
- );
- relMinThickness[ppPointI] = min
- (
- relMinThickness[ppPointI],
- patchRelMinThickness[patchI]
- );
- }
- }
-
- syncTools::syncPointList
- (
- mesh_,
- pp.meshPoints(),
- expansionRatio,
- minEqOp<scalar>(),
- GREAT, // null value
- false // no separation
- );
- syncTools::syncPointList
- (
- mesh_,
- pp.meshPoints(),
- finalLayerRatio,
- minEqOp<scalar>(),
- GREAT, // null value
- false // no separation
- );
- syncTools::syncPointList
- (
- mesh_,
- pp.meshPoints(),
- relMinThickness,
- minEqOp<scalar>(),
- GREAT, // null value
- false // no separation
- );
- }
-
-
-
- // Per mesh point the expansion parameters
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- thickness.setSize(pp.nPoints());
- minThickness.setSize(pp.nPoints());
-
- forAll(maxPointLevel, pointI)
- {
- // Find undistorted edge size for this level.
- scalar edgeLen = edge0Len/(1<<maxPointLevel[pointI]);
-
- // Calculate layer thickness based on expansion ratio
- // and final layer height
- if (expansionRatio[pointI] == 1)
- {
- thickness[pointI] =
- finalLayerRatio[pointI]
- * patchNLayers[pointI]
- * edgeLen;
- minThickness[pointI] = relMinThickness[pointI]*edgeLen;
- }
- else
- {
- scalar invExpansion = 1.0 / expansionRatio[pointI];
- label nLay = patchNLayers[pointI];
- thickness[pointI] =
- finalLayerRatio[pointI]
- * edgeLen
- * (1.0 - pow(invExpansion, nLay))
- / (1.0 - invExpansion);
- minThickness[pointI] = relMinThickness[pointI]*edgeLen;
- }
- }
-
- Info<< "calculateLayerThickness : min:" << gMin(thickness)
- << " max:" << gMax(thickness) << endl;
-}
-
-
-// Synchronize displacement among coupled patches.
-void Foam::autoHexMeshDriver::syncPatchDisplacement
-(
- const motionSmoother& meshMover,
- const scalarField& minThickness,
- pointField& patchDisp,
- labelList& patchNLayers,
- List<extrudeMode>& extrudeStatus
-) const
-{
- const labelList& meshPoints = meshMover.patch().meshPoints();
-
- label nChangedTotal = 0;
-
- while (true)
- {
- label nChanged = 0;
-
- // Sync displacement (by taking min)
- syncTools::syncPointList
- (
- mesh_,
- meshPoints,
- patchDisp,
- minEqOp<vector>(),
- wallPoint::greatPoint, // null value
- false // no separation
- );
-
- // Unmark if displacement too small
- forAll(patchDisp, i)
- {
- if (mag(patchDisp[i]) < minThickness[i])
- {
- if
- (
- unmarkExtrusion
- (
- i,
- patchDisp,
- patchNLayers,
- extrudeStatus
- )
- )
- {
- nChanged++;
- }
- }
- }
-
- labelList syncPatchNLayers(patchNLayers);
-
- syncTools::syncPointList
- (
- mesh_,
- meshPoints,
- syncPatchNLayers,
- minEqOp<label>(),
- labelMax, // null value
- false // no separation
- );
-
- // Reset if differs
- forAll(syncPatchNLayers, i)
- {
- if (syncPatchNLayers[i] != patchNLayers[i])
- {
- if
- (
- unmarkExtrusion
- (
- i,
- patchDisp,
- patchNLayers,
- extrudeStatus
- )
- )
- {
- nChanged++;
- }
- }
- }
-
- syncTools::syncPointList
- (
- mesh_,
- meshPoints,
- syncPatchNLayers,
- maxEqOp<label>(),
- labelMin, // null value
- false // no separation
- );
-
- // Reset if differs
- forAll(syncPatchNLayers, i)
- {
- if (syncPatchNLayers[i] != patchNLayers[i])
- {
- if
- (
- unmarkExtrusion
- (
- i,
- patchDisp,
- patchNLayers,
- extrudeStatus
- )
- )
- {
- nChanged++;
- }
- }
- }
- nChangedTotal += nChanged;
-
- if (!returnReduce(nChanged, sumOp<label>()))
- {
- break;
- }
- }
-
- Info<< "Prevented extrusion on "
- << returnReduce(nChangedTotal, sumOp<label>())
- << " coupled patch points during syncPatchDisplacement." << endl;
-}
-
-
-// Calculate displacement vector for all patch points. Uses pointNormal.
-// Checks that displaced patch point would be visible from all centres
-// of the faces using it.
-// extrudeStatus is both input and output and gives the status of each
-// patch point.
-void Foam::autoHexMeshDriver::getPatchDisplacement
-(
- const motionSmoother& meshMover,
- const scalarField& thickness,
- const scalarField& minThickness,
- pointField& patchDisp,
- labelList& patchNLayers,
- List<extrudeMode>& extrudeStatus
-) const
-{
- Info<< nl << "Determining displacement for added points"
- << " according to pointNormal ..." << endl;
-
- const indirectPrimitivePatch& pp = meshMover.patch();
- const vectorField& faceNormals = pp.faceNormals();
- const labelListList& pointFaces = pp.pointFaces();
- const pointField& localPoints = pp.localPoints();
- const labelList& meshPoints = pp.meshPoints();
-
- // Determine pointNormal
- // ~~~~~~~~~~~~~~~~~~~~~
-
- pointField pointNormals(pp.nPoints(), vector::zero);
- {
- labelList nPointFaces(pp.nPoints(), 0);
-
- forAll(faceNormals, faceI)
- {
- const face& f = pp.localFaces()[faceI];
-
- forAll(f, fp)
- {
- pointNormals[f[fp]] += faceNormals[faceI];
- nPointFaces[f[fp]] ++;
- }
- }
-
- syncTools::syncPointList
- (
- mesh_,
- meshPoints,
- pointNormals,
- plusEqOp<vector>(),
- vector::zero, // null value
- false // no separation
- );
-
- syncTools::syncPointList
- (
- mesh_,
- meshPoints,
- nPointFaces,
- plusEqOp<label>(),
- 0, // null value
- false // no separation
- );
-
- forAll(pointNormals, i)
- {
- pointNormals[i] /= nPointFaces[i];
- }
- }
-
-
- // Determine local length scale on patch
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- // Start off from same thickness everywhere (except where no extrusion)
- patchDisp = thickness*pointNormals;
-
- // Check if no extrude possible.
- forAll(pointNormals, patchPointI)
- {
- label meshPointI = pp.meshPoints()[patchPointI];
-
- if (extrudeStatus[patchPointI] == NOEXTRUDE)
- {
- // Do not use unmarkExtrusion; forcibly set to zero extrusion.
- patchNLayers[patchPointI] = 0;
- patchDisp[patchPointI] = vector::zero;
- }
- else
- {
- // Get normal
- const vector& n = pointNormals[patchPointI];
-
- if (!meshTools::visNormal(n, faceNormals, pointFaces[patchPointI]))
- {
- Pout<< "No valid normal for point " << meshPointI
- << ' ' << pp.points()[meshPointI]
- << "; setting displacement to " << patchDisp[patchPointI]
- << endl;
-
- extrudeStatus[patchPointI] = EXTRUDEREMOVE;
- }
- }
- }
-
- // At illegal points make displacement average of new neighbour positions
- forAll(extrudeStatus, patchPointI)
- {
- if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
- {
- point avg(vector::zero);
- label nPoints = 0;
-
- const labelList& pEdges = pp.pointEdges()[patchPointI];
-
- forAll(pEdges, i)
- {
- label edgeI = pEdges[i];
-
- label otherPointI = pp.edges()[edgeI].otherVertex(patchPointI);
-
- if (extrudeStatus[otherPointI] != NOEXTRUDE)
- {
- avg += localPoints[otherPointI] + patchDisp[otherPointI];
- nPoints++;
- }
- }
-
- if (nPoints > 0)
- {
- Pout<< "Displacement at illegal point "
- << localPoints[patchPointI]
- << " set to " << (avg / nPoints - localPoints[patchPointI])
- << endl;
-
- patchDisp[patchPointI] =
- avg / nPoints
- - localPoints[patchPointI];
- }
- }
- }
-
- // Make sure displacement is equal on both sides of coupled patches.
- syncPatchDisplacement
- (
- meshMover,
- minThickness,
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
-
- Info<< endl;
-}
-
-
-// Truncates displacement
-// - for all patchFaces in the faceset displacement gets set to zero
-// - all displacement < minThickness gets set to zero
-Foam::label Foam::autoHexMeshDriver::truncateDisplacement
-(
- const motionSmoother& meshMover,
- const scalarField& minThickness,
- const faceSet& illegalPatchFaces,
- pointField& patchDisp,
- labelList& patchNLayers,
- List<extrudeMode>& extrudeStatus
-) const
-{
- const polyMesh& mesh = meshMover.mesh();
- const indirectPrimitivePatch& pp = meshMover.patch();
-
- label nChanged = 0;
-
- const Map<label>& meshPointMap = pp.meshPointMap();
-
- forAllConstIter(faceSet, illegalPatchFaces, iter)
- {
- label faceI = iter.key();
-
- if (mesh.isInternalFace(faceI))
- {
- FatalErrorIn("truncateDisplacement(..)")
- << "Faceset " << illegalPatchFaces.name()
- << " contains internal face " << faceI << nl
- << "It should only contain patch faces" << abort(FatalError);
- }
-
- const face& f = mesh.faces()[faceI];
-
-
- forAll(f, fp)
- {
- if (meshPointMap.found(f[fp]))
- {
- label patchPointI = meshPointMap[f[fp]];
-
- if (extrudeStatus[patchPointI] != NOEXTRUDE)
- {
- unmarkExtrusion
- (
- patchPointI,
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
- nChanged++;
- }
- }
- }
- }
-
- forAll(patchDisp, patchPointI)
- {
- if (mag(patchDisp[patchPointI]) < minThickness[patchPointI])
- {
- if
- (
- unmarkExtrusion
- (
- patchPointI,
- patchDisp,
- patchNLayers,
- extrudeStatus
- )
- )
- {
- nChanged++;
- }
- }
- else if (extrudeStatus[patchPointI] == NOEXTRUDE)
- {
- // Make sure displacement is 0. Should already be so but ...
- patchDisp[patchPointI] = vector::zero;
- patchNLayers[patchPointI] = 0;
- }
- }
-
-
- const faceList& localFaces = pp.localFaces();
-
- while (true)
- {
- syncPatchDisplacement
- (
- meshMover,
- minThickness,
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
-
- // Make sure that a face doesn't have two non-consecutive areas
- // not extruded (e.g. quad where vertex 0 and 2 are not extruded
- // but 1 and 3 are) since this gives topological errors.
-
- label nPinched = 0;
-
- forAll(localFaces, i)
- {
- const face& localF = localFaces[i];
-
- // Count number of transitions from unsnapped to snapped.
- label nTrans = 0;
-
- extrudeMode prevMode = extrudeStatus[localF.prevLabel(0)];
-
- forAll(localF, fp)
- {
- extrudeMode fpMode = extrudeStatus[localF[fp]];
-
- if (prevMode == NOEXTRUDE && fpMode != NOEXTRUDE)
- {
- nTrans++;
- }
- prevMode = fpMode;
- }
-
- if (nTrans > 1)
- {
- // Multiple pinches. Reset whole face as unextruded.
- if
- (
- unmarkExtrusion
- (
- localF,
- patchDisp,
- patchNLayers,
- extrudeStatus
- )
- )
- {
- nPinched++;
- nChanged++;
- }
- }
- }
-
- reduce(nPinched, sumOp<label>());
-
- Info<< "truncateDisplacement : Unextruded " << nPinched
- << " faces due to non-consecutive vertices being extruded." << endl;
-
-
- // Make sure that a face has consistent number of layers for all
- // its vertices.
-
- label nDiffering = 0;
-
- //forAll(localFaces, i)
- //{
- // const face& localF = localFaces[i];
- //
- // label numLayers = -1;
- //
- // forAll(localF, fp)
- // {
- // if (patchNLayers[localF[fp]] > 0)
- // {
- // if (numLayers == -1)
- // {
- // numLayers = patchNLayers[localF[fp]];
- // }
- // else if (numLayers != patchNLayers[localF[fp]])
- // {
- // // Differing number of layers
- // if
- // (
- // unmarkExtrusion
- // (
- // localF,
- // patchDisp,
- // patchNLayers,
- // extrudeStatus
- // )
- // )
- // {
- // nDiffering++;
- // nChanged++;
- // }
- // break;
- // }
- // }
- // }
- //}
- //
- //reduce(nDiffering, sumOp<label>());
- //
- //Info<< "truncateDisplacement : Unextruded " << nDiffering
- // << " faces due to having differing number of layers." << endl;
-
- if (nPinched+nDiffering == 0)
- {
- break;
- }
- }
-
- return nChanged;
-}
-
-
-// Setup layer information (at points and faces) to modify mesh topology in
-// regions where layer mesh terminates.
-void Foam::autoHexMeshDriver::setupLayerInfoTruncation
-(
- const motionSmoother& meshMover,
- const labelList& patchNLayers,
- const List<extrudeMode>& extrudeStatus,
- const label nBufferCellsNoExtrude,
- labelList& nPatchPointLayers,
- labelList& nPatchFaceLayers
-) const
-{
- Info<< nl << "Setting up information for layer truncation ..." << endl;
-
- const indirectPrimitivePatch& pp = meshMover.patch();
- const polyMesh& mesh = meshMover.mesh();
-
- if (nBufferCellsNoExtrude < 0)
- {
- Info<< nl << "Performing no layer truncation."
- << " nBufferCellsNoExtrude set to less than 0 ..." << endl;
-
- // Face layers if any point get extruded
- forAll(pp.localFaces(), patchFaceI)
- {
- const face& f = pp.localFaces()[patchFaceI];
-
- forAll(f, fp)
- {
- if (patchNLayers[f[fp]] > 0)
- {
- nPatchFaceLayers[patchFaceI] = patchNLayers[f[fp]];
- break;
- }
- }
- }
- nPatchPointLayers = patchNLayers;
- }
- else
- {
- // Determine max point layers per face.
- labelList maxLevel(pp.size(), 0);
-
- forAll(pp.localFaces(), patchFaceI)
- {
- const face& f = pp.localFaces()[patchFaceI];
-
- // find patch faces where layer terminates (i.e contains extrude
- // and noextrude points).
-
- bool noExtrude = false;
- label mLevel = 0;
-
- forAll(f, fp)
- {
- if (extrudeStatus[f[fp]] == NOEXTRUDE)
- {
- noExtrude = true;
- }
- mLevel = max(mLevel, patchNLayers[f[fp]]);
- }
-
- if (mLevel > 0)
- {
- // So one of the points is extruded. Check if all are extruded
- // or is a mix.
-
- if (noExtrude)
- {
- nPatchFaceLayers[patchFaceI] = 1;
- maxLevel[patchFaceI] = mLevel;
- }
- else
- {
- maxLevel[patchFaceI] = mLevel;
- }
- }
- }
-
- // We have the seed faces (faces with nPatchFaceLayers != maxLevel)
- // Now do a meshwave across the patch where we pick up neighbours
- // of seed faces.
- // Note: quite inefficient. Could probably be coded better.
-
- const labelListList& pointFaces = pp.pointFaces();
-
- label nLevels = gMax(patchNLayers);
-
- // flag neighbouring patch faces with number of layers to grow
- for (label ilevel = 1; ilevel < nLevels; ilevel++)
- {
- label nBuffer;
-
- if (ilevel == 1)
- {
- nBuffer = nBufferCellsNoExtrude - 1;
- }
- else
- {
- nBuffer = nBufferCellsNoExtrude;
- }
-
- for (label ibuffer = 0; ibuffer < nBuffer + 1; ibuffer++)
- {
- labelList tempCounter(nPatchFaceLayers);
-
- boolList foundNeighbour(pp.nPoints(), false);
-
- forAll(pp.meshPoints(), patchPointI)
- {
- forAll(pointFaces[patchPointI], pointFaceI)
- {
- label faceI = pointFaces[patchPointI][pointFaceI];
-
- if
- (
- nPatchFaceLayers[faceI] != -1
- && maxLevel[faceI] > 0
- )
- {
- foundNeighbour[patchPointI] = true;
- break;
- }
- }
- }
-
- syncTools::syncPointList
- (
- mesh,
- pp.meshPoints(),
- foundNeighbour,
- orEqOp<bool>(),
- false, // null value
- false // no separation
- );
-
- forAll(pp.meshPoints(), patchPointI)
- {
- if (foundNeighbour[patchPointI])
- {
- forAll(pointFaces[patchPointI], pointFaceI)
- {
- label faceI = pointFaces[patchPointI][pointFaceI];
- if
- (
- nPatchFaceLayers[faceI] == -1
- && maxLevel[faceI] > 0
- && ilevel < maxLevel[faceI]
- )
- {
- tempCounter[faceI] = ilevel;
- }
- }
- }
- }
- nPatchFaceLayers = tempCounter;
- }
- }
-
- forAll(pp.localFaces(), patchFaceI)
- {
- if (nPatchFaceLayers[patchFaceI] == -1)
- {
- nPatchFaceLayers[patchFaceI] = maxLevel[patchFaceI];
- }
- }
-
- forAll(pp.meshPoints(), patchPointI)
- {
- if (extrudeStatus[patchPointI] != NOEXTRUDE)
- {
- forAll(pointFaces[patchPointI], pointFaceI)
- {
- label face = pointFaces[patchPointI][pointFaceI];
- nPatchPointLayers[patchPointI] = max
- (
- nPatchPointLayers[patchPointI],
- nPatchFaceLayers[face]
- );
- }
- }
- else
- {
- nPatchPointLayers[patchPointI] = 0;
- }
- }
- syncTools::syncPointList
- (
- mesh,
- pp.meshPoints(),
- nPatchPointLayers,
- maxEqOp<label>(),
- 0, // null value
- false // no separation
- );
- }
-}
-
-
-// Does any of the cells use a face from faces?
-bool Foam::autoHexMeshDriver::cellsUseFace
-(
- const polyMesh& mesh,
- const labelList& cellLabels,
- const labelHashSet& faces
-)
-{
- forAll(cellLabels, i)
- {
- const cell& cFaces = mesh.cells()[cellLabels[i]];
-
- forAll(cFaces, cFaceI)
- {
- if (faces.found(cFaces[cFaceI]))
- {
- return true;
- }
- }
- }
- return false;
-}
-
-
-// Checks the newly added cells and locally unmarks points so they
-// will not get extruded next time round. Returns global number of unmarked
-// points (0 if all was fine)
-Foam::label Foam::autoHexMeshDriver::checkAndUnmark
-(
- const addPatchCellLayer& addLayer,
- const dictionary& motionDict,
- const indirectPrimitivePatch& pp,
- const fvMesh& newMesh,
-
- pointField& patchDisp,
- labelList& patchNLayers,
- List<extrudeMode>& extrudeStatus
-)
-{
- // Check the resulting mesh for errors
- Info<< nl << "Checking mesh with layer ..." << endl;
- faceSet wrongFaces(newMesh, "wrongFaces", newMesh.nFaces()/1000);
- motionSmoother::checkMesh(false, newMesh, motionDict, wrongFaces);
- Info<< "Detected " << returnReduce(wrongFaces.size(), sumOp<label>())
- << " illegal faces"
- << " (concave, zero area or negative cell pyramid volume)"
- << endl;
-
- // Undo local extrusion if
- // - any of the added cells in error
-
- label nChanged = 0;
-
- // Get all cells in the layer.
- labelListList addedCells
- (
- addPatchCellLayer::addedCells
- (
- newMesh,
- addLayer.layerFaces()
- )
- );
-
- // Check if any of the faces in error uses any face of an added cell
- forAll(addedCells, oldPatchFaceI)
- {
- // Get the cells (in newMesh labels) per old patch face (in mesh
- // labels)
- const labelList& fCells = addedCells[oldPatchFaceI];
-
- if (cellsUseFace(newMesh, fCells, wrongFaces))
- {
- // Unmark points on old mesh
- if
- (
- unmarkExtrusion
- (
- pp.localFaces()[oldPatchFaceI],
- patchDisp,
- patchNLayers,
- extrudeStatus
- )
- )
- {
- nChanged++;
- }
- }
- }
-
- return returnReduce(nChanged, sumOp<label>());
-}
-
-
-//- Count global number of extruded faces
-Foam::label Foam::autoHexMeshDriver::countExtrusion
-(
- const indirectPrimitivePatch& pp,
- const List<extrudeMode>& extrudeStatus
-)
-{
- // Count number of extruded patch faces
- label nExtruded = 0;
- {
- const faceList& localFaces = pp.localFaces();
-
- forAll(localFaces, i)
- {
- const face& localFace = localFaces[i];
-
- forAll(localFace, fp)
- {
- if (extrudeStatus[localFace[fp]] != NOEXTRUDE)
- {
- nExtruded++;
- break;
- }
- }
- }
- }
-
- return returnReduce(nExtruded, sumOp<label>());
-}
-
-
-// Collect layer faces and layer cells into bools for ease of handling
-void Foam::autoHexMeshDriver::getLayerCellsFaces
-(
- const polyMesh& mesh,
- const addPatchCellLayer& addLayer,
- boolList& flaggedCells,
- boolList& flaggedFaces
-)
-{
- flaggedCells.setSize(mesh.nCells());
- flaggedCells = false;
- flaggedFaces.setSize(mesh.nFaces());
- flaggedFaces = false;
-
- // Mark all faces in the layer
- const labelListList& layerFaces = addLayer.layerFaces();
-
- // Mark all cells in the layer.
- labelListList addedCells(addPatchCellLayer::addedCells(mesh, layerFaces));
-
- forAll(addedCells, oldPatchFaceI)
- {
- const labelList& added = addedCells[oldPatchFaceI];
-
- forAll(added, i)
- {
- flaggedCells[added[i]] = true;
- }
- }
-
- forAll(layerFaces, oldPatchFaceI)
- {
- const labelList& layer = layerFaces[oldPatchFaceI];
-
- if (layer.size() > 0)
- {
- for (label i = 1; i < layer.size()-1; i++)
- {
- flaggedFaces[layer[i]] = true;
- }
- }
- }
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-void Foam::autoHexMeshDriver::mergePatchFacesUndo
-(
- const layerParameters& layerParams,
- const dictionary& motionDict
-)
-{
- scalar minCos = Foam::cos
- (
- layerParams.featureAngle()
- * mathematicalConstant::pi/180.0
- );
-
- scalar concaveCos = Foam::cos
- (
- layerParams.concaveAngle()
- * mathematicalConstant::pi/180.0
- );
-
- Info<< nl
- << "Merging all faces of a cell" << nl
- << "---------------------------" << nl
- << " - which are on the same patch" << nl
- << " - which make an angle < " << layerParams.featureAngle()
- << " degrees"
- << nl
- << " (cos:" << minCos << ')' << nl
- << " - as long as the resulting face doesn't become concave"
- << " by more than "
- << layerParams.concaveAngle()
- << " degrees (0=straight, 180=fully concave)" << nl
- << endl;
-
- label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict);
-
- nChanged += mergeEdgesUndo(minCos, motionDict);
-}
-
-
-void Foam::autoHexMeshDriver::addLayers
-(
- const layerParameters& layerParams,
- const dictionary& motionDict,
- const label nAllowableErrors,
- motionSmoother& meshMover
-)
-{
- const indirectPrimitivePatch& pp = meshMover.patch();
- const labelList& meshPoints = pp.meshPoints();
-
- // Precalculate mesh edge labels for patch edges
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- labelList meshEdges(pp.nEdges());
- forAll(pp.edges(), edgeI)
- {
- const edge& ppEdge = pp.edges()[edgeI];
- label v0 = meshPoints[ppEdge[0]];
- label v1 = meshPoints[ppEdge[1]];
- meshEdges[edgeI] = meshTools::findEdge
- (
- mesh_.edges(),
- mesh_.pointEdges()[v0],
- v0,
- v1
- );
- }
-
- // Displacement for all pp.localPoints.
- vectorField patchDisp(pp.nPoints(), vector::one);
-
- // Number of layers for all pp.localPoints. Note: only valid if
- // extrudeStatus = EXTRUDE.
- labelList patchNLayers(pp.nPoints(), 0);
-
- // Whether to add edge for all pp.localPoints.
- List<extrudeMode> extrudeStatus(pp.nPoints(), EXTRUDE);
-
-
- // Get number of layer per point from number of layers per patch
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- setNumLayers
- (
- layerParams.numLayers(), // per patch the num layers
- meshMover.adaptPatchIDs(), // patches that are being moved
- pp, // indirectpatch for all faces moving
-
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
-
- // Disable extrusion on non-manifold points
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- handleNonManifolds
- (
- pp,
- meshEdges,
-
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
-
- // Disable extrusion on feature angles
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- handleFeatureAngle
- (
- pp,
- meshEdges,
- layerParams.featureAngle()*mathematicalConstant::pi/180.0,
-
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
-
- // Disable extrusion on warped faces
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- // Undistorted edge length
- const scalar edge0Len = meshRefinerPtr_().meshCutter().level0EdgeLength();
- const labelList& cellLevel = meshRefinerPtr_().meshCutter().cellLevel();
-
- handleWarpedFaces
- (
- pp,
- layerParams.maxFaceThicknessRatio(),
- edge0Len,
- cellLevel,
-
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
-
- //// Disable extrusion on cells with multiple patch faces
- //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- //
- //handleMultiplePatchFaces
- //(
- // pp,
- //
- // patchDisp,
- // patchNLayers,
- // extrudeStatus
- //);
-
-
- // Grow out region of non-extrusion
- for (label i = 0; i < layerParams.nGrow(); i++)
- {
- growNoExtrusion
- (
- pp,
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
- }
-
- // Determine (wanted) point-wise layer thickness and expansion ratio
- scalarField thickness(pp.nPoints());
- scalarField minThickness(pp.nPoints());
- scalarField expansionRatio(pp.nPoints());
- calculateLayerThickness
- (
- pp,
- meshMover.adaptPatchIDs(),
- layerParams.expansionRatio(),
- layerParams.finalLayerRatio(),
- layerParams.minThickness(),
- cellLevel,
- patchNLayers,
- edge0Len,
-
- thickness,
- minThickness,
- expansionRatio
- );
-
- // Calculate wall to medial axis distance for smoothing displacement
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- pointScalarField pointMedialDist
- (
- IOobject
- (
- "pointMedialDist",
- mesh_.time().timeName(),
- mesh_,
- IOobject::NO_READ,
- IOobject::NO_WRITE,
- false
- ),
- meshMover.pMesh(),
- dimensionedScalar("pointMedialDist", dimless, 0.0)
- );
-
- pointVectorField dispVec
- (
- IOobject
- (
- "dispVec",
- mesh_.time().timeName(),
- mesh_,
- IOobject::NO_READ,
- IOobject::NO_WRITE,
- false
- ),
- meshMover.pMesh(),
- dimensionedVector("dispVec", dimless, vector::zero)
- );
-
- pointScalarField medialRatio
- (
- IOobject
- (
- "medialRatio",
- mesh_.time().timeName(),
- mesh_,
- IOobject::NO_READ,
- IOobject::NO_WRITE,
- false
- ),
- meshMover.pMesh(),
- dimensionedScalar("medialRatio", dimless, 0.0)
- );
-
- // Setup information for medial axis smoothing. Calculates medial axis
- // and a smoothed displacement direction.
- // - pointMedialDist : distance to medial axis
- // - dispVec : normalised direction of nearest displacement
- // - medialRatio : ratio of medial distance to wall distance.
- // (1 at wall, 0 at medial axis)
- medialAxisSmoothingInfo
- (
- meshMover,
- layerParams.nSmoothNormals(),
- layerParams.nSmoothSurfaceNormals(),
- layerParams.minMedianAxisAngleCos(),
-
- dispVec,
- medialRatio,
- pointMedialDist
- );
-
-
-
- // Saved old points
- pointField oldPoints(mesh_.points());
-
- // Last set of topology changes. (changing mesh clears out polyTopoChange)
- polyTopoChange savedMeshMod(mesh_.boundaryMesh().size());
-
- boolList flaggedCells;
- boolList flaggedFaces;
-
- while (true)
- {
- // Make sure displacement is equal on both sides of coupled patches.
- syncPatchDisplacement
- (
- meshMover,
- minThickness,
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
-
- // Displacement acc. to pointnormals
- getPatchDisplacement
- (
- meshMover,
- thickness,
- minThickness,
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
-
- // Shrink mesh by displacement value first.
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- {
- pointField oldPatchPos(pp.localPoints());
-
- //// Laplacian displacement shrinking.
- //shrinkMeshDistance
- //(
- // debug,
- // meshMover,
- // -patchDisp, // Shrink in opposite direction of addedPoints
- // layerParams.nSmoothDisp(),
- // layerParams.nSnap()
- //);
-
- // Medial axis based shrinking
- shrinkMeshMedialDistance
- (
- meshMover,
-
- layerParams.nSmoothThickness(),
- layerParams.maxThicknessToMedialRatio(),
- nAllowableErrors,
- layerParams.nSnap(),
- layerParams.layerTerminationCos(),
-
- thickness,
- minThickness,
-
- dispVec,
- medialRatio,
- pointMedialDist,
-
- extrudeStatus,
- patchDisp,
- patchNLayers
- );
-
- // Update patchDisp (since not all might have been honoured)
- patchDisp = oldPatchPos - pp.localPoints();
- }
-
- // Truncate displacements that are too small (this will do internal
- // ones, coupled ones have already been truncated by syncPatch)
- faceSet dummySet(mesh_, "wrongPatchFaces", 0);
- truncateDisplacement
- (
- meshMover,
- minThickness,
- dummySet,
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
-
-
- // Dump to .obj file for debugging.
- if (debug_)
- {
- dumpDisplacement
- (
- mesh_.time().path()/"layer",
- pp,
- patchDisp,
- extrudeStatus
- );
-
- const_cast<Time&>(mesh_.time())++;
- Info<< "Writing shrunk mesh to " << mesh_.time().timeName() << endl;
- mesh_.write();
- }
-
-
- // Mesh topo change engine
- polyTopoChange meshMod(mesh_);
-
- // Grow layer of cells on to patch. Handles zero sized displacement.
- addPatchCellLayer addLayer(mesh_);
-
- // Determine per point/per face number of layers to extrude. Also
- // handles the slow termination of layers when going switching layers
-
- labelList nPatchPointLayers(pp.nPoints(),-1);
- labelList nPatchFaceLayers(pp.localFaces().size(),-1);
- setupLayerInfoTruncation
- (
- meshMover,
- patchNLayers,
- extrudeStatus,
- layerParams.nBufferCellsNoExtrude(),
- nPatchPointLayers,
- nPatchFaceLayers
- );
-
- // Calculate displacement for first layer for addPatchLayer
- vectorField firstDisp(patchNLayers.size(), vector::zero);
-
- forAll(patchNLayers, i)
- {
- if (patchNLayers[i] > 0)
- {
- if (expansionRatio[i] == 1.0)
- {
- firstDisp[i] = patchDisp[i]/nPatchPointLayers[i];
- }
- else
- {
- label nLay = nPatchPointLayers[i];
- scalar h =
- pow(expansionRatio[i], nLay - 1)
- * (mag(patchDisp[i])*(1.0 - expansionRatio[i]))
- / (1.0 - pow(expansionRatio[i], nLay));
- firstDisp[i] = h/mag(patchDisp[i])*patchDisp[i];
- }
- }
- }
-
- scalarField invExpansionRatio = 1.0 / expansionRatio;
-
- // Add topo regardless of whether extrudeStatus is extruderemove.
- // Not add layer if patchDisp is zero.
- addLayer.setRefinement
- (
- invExpansionRatio,
- pp,
- nPatchFaceLayers, // layers per face
- nPatchPointLayers, // layers per point
- firstDisp, // thickness of first layer
- meshMod
- );
-
- if (debug_)
- {
- const_cast<Time&>(mesh_.time())++;
- }
-
- // Store mesh changes for if mesh is correct.
- savedMeshMod = meshMod;
-
-
- // With the stored topo changes we create a new mesh so we can
- // undo if neccesary.
-
- autoPtr<fvMesh> newMeshPtr;
- autoPtr<mapPolyMesh> map = meshMod.makeMesh
- (
- newMeshPtr,
- IOobject
- (
- //mesh.name()+"_layer",
- mesh_.name(),
- static_cast<polyMesh&>(mesh_).instance(),
- mesh_.db(),
- static_cast<polyMesh&>(mesh_).readOpt(),
- static_cast<polyMesh&>(mesh_).writeOpt()
- ), // io params from original mesh but new name
- mesh_, // original mesh
- true // parallel sync
- );
- fvMesh& newMesh = newMeshPtr();
-
- //?neccesary? Update fields
- newMesh.updateMesh(map);
-
- // Update numbering on addLayer:
- // - cell/point labels to be newMesh.
- // - patchFaces to remain in oldMesh order.
- addLayer.updateMesh
- (
- map,
- identity(pp.size()),
- identity(pp.nPoints())
- );
-
- // Collect layer faces and cells for outside loop.
- getLayerCellsFaces
- (
- newMesh,
- addLayer,
- flaggedCells,
- flaggedFaces
- );
-
-
- if (debug_)
- {
- Info<< "Writing layer mesh to " << mesh_.time().timeName() << endl;
- newMesh.write();
- cellSet addedCellSet
- (
- newMesh,
- "addedCells",
- findIndices(flaggedCells, true)
- );
- Info<< "Writing "
- << returnReduce(addedCellSet.size(), sumOp<label>())
- << " added cells to cellSet "
- << addedCellSet.objectPath()
- << endl;
- addedCellSet.write();
-
- faceSet layerFacesSet
- (
- newMesh,
- "layerFaces",
- findIndices(flaggedCells, true)
- );
- Info<< "Writing "
- << returnReduce(layerFacesSet.size(), sumOp<label>())
- << " faces inside added layer to faceSet "
- << layerFacesSet.objectPath()
- << endl;
- layerFacesSet.write();
- }
-
- // Unset the extrusion at the pp.
- label nTotChanged = checkAndUnmark
- (
- addLayer,
- motionDict,
- pp,
- newMesh,
-
- patchDisp,
- patchNLayers,
- extrudeStatus
- );
-
- Info<< "Extruding " << countExtrusion(pp, extrudeStatus)
- << " out of " << returnReduce(pp.size(), sumOp<label>())
- << " faces. Removed extrusion at " << nTotChanged << " faces."
- << endl;
-
- if (nTotChanged == 0)
- {
- break;
- }
-
- // Reset mesh points and start again
- meshMover.movePoints(oldPoints);
- meshMover.correct();
- }
-
-
- // At this point we have a (shrunk) mesh and a set of topology changes
- // which will make a valid mesh with layer. Apply these changes to the
- // current mesh.
-
- // Apply the stored topo changes to the current mesh.
- autoPtr<mapPolyMesh> map = savedMeshMod.changeMesh(mesh_, false);
-
- // Update fields
- mesh_.updateMesh(map);
-
- // Move mesh (since morphing does not do this)
- if (map().hasMotionPoints())
- {
- mesh_.movePoints(map().preMotionPoints());
- }
-
- meshRefinerPtr_().updateMesh(map, labelList(0));
-
-
- // Do final balancing
- // ~~~~~~~~~~~~~~~~~~
-
- if (Pstream::parRun())
- {
- // Balance. No restriction on face zones and baffles.
- autoPtr<mapDistributePolyMesh> map = balance(false, false);
-
- // Re-distribute flag of layer faces and cells
- map().distributeCellData(flaggedCells);
- map().distributeFaceData(flaggedFaces);
- }
-
-
- // Write mesh
- // ~~~~~~~~~~
-
- //writeMesh("Layer mesh");
- cellSet addedCellSet(mesh_, "addedCells", findIndices(flaggedCells, true));
- Info<< "Writing "
- << returnReduce(addedCellSet.size(), sumOp<label>())
- << " added cells to cellSet "
- << addedCellSet.objectPath()
- << endl;
- addedCellSet.write();
-
- faceSet layerFacesSet(mesh_, "layerFaces", findIndices(flaggedFaces, true));
- Info<< "Writing "
- << returnReduce(layerFacesSet.size(), sumOp<label>())
- << " faces inside added layer to faceSet "
- << layerFacesSet.objectPath()
- << endl;
- layerFacesSet.write();
-}
-
-
-// ************************************************************************* //
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverShrink.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverShrink.C
deleted file mode 100644
index 0c66073e84c422135bcdb5aa82a2b6314641d61b..0000000000000000000000000000000000000000
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverShrink.C
+++ /dev/null
@@ -1,1151 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software; you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by the
- Free Software Foundation; either version 2 of the License, or (at your
- option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM; if not, write to the Free Software Foundation,
- Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
- Shrinking mesh (part of adding cell layers)
-
-\*----------------------------------------------------------------------------*/
-
-#include "autoHexMeshDriver.H"
-#include "fvMesh.H"
-#include "Time.H"
-#include "pointFields.H"
-#include "motionSmoother.H"
-#include "pointData.H"
-#include "PointEdgeWave.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
-// Calculate inverse sum of edge weights (currently always 1.0)
-void Foam::autoHexMeshDriver::sumWeights
-(
- const PackedList<1>& isMasterEdge,
- const labelList& meshEdges,
- const labelList& meshPoints,
- const edgeList& edges,
- scalarField& invSumWeight
-) const
-{
- invSumWeight = 0;
-
- forAll(edges, edgeI)
- {
- if (isMasterEdge.get(meshEdges[edgeI]) == 1)
- {
- const edge& e = edges[edgeI];
- //scalar eWeight = edgeWeights[edgeI];
- scalar eWeight = 1.0;
-
- invSumWeight[e[0]] += eWeight;
- invSumWeight[e[1]] += eWeight;
- }
- }
-
- syncTools::syncPointList
- (
- mesh_,
- meshPoints,
- invSumWeight,
- plusEqOp<scalar>(),
- scalar(0.0), // null value
- false // no separation
- );
-
- forAll(invSumWeight, pointI)
- {
- scalar w = invSumWeight[pointI];
-
- if (w > 0.0)
- {
- invSumWeight[pointI] = 1.0/w;
- }
- }
-}
-
-
-// Smooth field on moving patch
-void Foam::autoHexMeshDriver::smoothField
-(
- const motionSmoother& meshMover,
- const PackedList<1>& isMasterEdge,
- const labelList& meshEdges,
- const scalarField& fieldMin,
- const label& nSmoothDisp,
- scalarField& field
-) const
-{
- const indirectPrimitivePatch& pp = meshMover.patch();
- const edgeList& edges = pp.edges();
- const labelList& meshPoints = pp.meshPoints();
-
- scalarField invSumWeight(pp.nPoints());
- sumWeights
- (
- isMasterEdge,
- meshEdges,
- meshPoints,
- edges,
- invSumWeight
- );
-
- // Get smoothly varying patch field.
- Info<< "shrinkMeshDistance : Smoothing field ..." << endl;
-
- for (label iter = 0; iter < nSmoothDisp; iter++)
- {
- scalarField average(pp.nPoints());
- averageNeighbours
- (
- isMasterEdge,
- meshEdges,
- meshPoints,
- pp.edges(),
- invSumWeight,
- field,
- average
- );
-
- // Transfer to field
- forAll(field, pointI)
- {
- //full smoothing neighbours + point value
- average[pointI] = 0.5*(field[pointI]+average[pointI]);
-
- // perform monotonic smoothing
- if
- (
- average[pointI] < field[pointI]
- && average[pointI] >= fieldMin[pointI]
- )
- {
- field[pointI] = average[pointI];
- }
- }
-
- // Do residual calculation every so often.
- if ((iter % 10) == 0)
- {
- Info<< " Iteration " << iter << " residual "
- << gSum(mag(field-average))
- /returnReduce(average.size(), sumOp<label>())
- << endl;
- }
- }
-}
-
-
-// Smooth normals on moving patch.
-void Foam::autoHexMeshDriver::smoothPatchNormals
-(
- const motionSmoother& meshMover,
- const PackedList<1>& isMasterEdge,
- const labelList& meshEdges,
- const label nSmoothDisp,
- pointField& normals
-) const
-{
- const indirectPrimitivePatch& pp = meshMover.patch();
- const edgeList& edges = pp.edges();
- const labelList& meshPoints = pp.meshPoints();
-
- // Calculate inverse sum of weights
-
- scalarField invSumWeight(pp.nPoints());
- sumWeights
- (
- isMasterEdge,
- meshEdges,
- meshPoints,
- edges,
- invSumWeight
- );
-
- // Get smoothly varying internal normals field.
- Info<< "shrinkMeshDistance : Smoothing normals ..." << endl;
-
- for (label iter = 0; iter < nSmoothDisp; iter++)
- {
- vectorField average(pp.nPoints());
- averageNeighbours
- (
- isMasterEdge,
- meshEdges,
- meshPoints,
- pp.edges(),
- invSumWeight,
- normals,
- average
- );
-
- // Do residual calculation every so often.
- if ((iter % 10) == 0)
- {
- Info<< " Iteration " << iter << " residual "
- << gSum(mag(normals-average))
- /returnReduce(average.size(), sumOp<label>())
- << endl;
- }
-
- // Transfer to normals vector field
- forAll(average, pointI)
- {
- // full smoothing neighbours + point value
- average[pointI] = 0.5*(normals[pointI]+average[pointI]);
- normals[pointI] = average[pointI];
- normals[pointI] /= mag(normals[pointI]) + VSMALL;
- }
- }
-}
-
-
-// Smooth normals in interior.
-void Foam::autoHexMeshDriver::smoothNormals
-(
- const label nSmoothDisp,
- const PackedList<1>& isMasterEdge,
- const labelList& fixedPoints,
- pointVectorField& normals
-) const
-{
- // Get smoothly varying internal normals field.
- Info<< "shrinkMeshDistance : Smoothing normals ..." << endl;
-
- const edgeList& edges = mesh_.edges();
-
- // Points that do not change.
- PackedList<1> isFixedPoint(mesh_.nPoints(), 0);
-
- // Internal points that are fixed
- forAll(fixedPoints, i)
- {
- label meshPointI = fixedPoints[i];
- isFixedPoint.set(meshPointI, 1);
- }
-
- // Correspondence between local edges/points and mesh edges/points
- const labelList meshEdges(identity(mesh_.nEdges()));
- const labelList meshPoints(identity(mesh_.nPoints()));
-
- // Calculate inverse sum of weights
-
- scalarField invSumWeight(mesh_.nPoints(), 0);
- sumWeights
- (
- isMasterEdge,
- meshEdges,
- meshPoints,
- edges,
- invSumWeight
- );
-
- Info<< "shrinkMeshDistance : Smoothing normals in interior ..." << endl;
-
- for (label iter = 0; iter < nSmoothDisp; iter++)
- {
- vectorField average(mesh_.nPoints());
- averageNeighbours
- (
- isMasterEdge,
- meshEdges,
- meshPoints,
- edges,
- invSumWeight,
- normals,
- average
- );
-
- // Do residual calculation every so often.
- if ((iter % 10) == 0)
- {
- Info<< " Iteration " << iter << " residual "
- << gSum(mag(normals-average))
- /returnReduce(average.size(), sumOp<label>())
- << endl;
- }
-
-
- // Transfer to normals vector field
- forAll(average, pointI)
- {
- if (isFixedPoint.get(pointI) == 0)
- {
- //full smoothing neighbours + point value
- average[pointI] = 0.5*(normals[pointI]+average[pointI]);
- normals[pointI] = average[pointI];
- normals[pointI] /= mag(normals[pointI]) + VSMALL;
- }
- }
- }
-}
-
-
-// Tries and find a medial axis point. Done by comparing vectors to nearest
-// wall point for both vertices of edge.
-bool Foam::autoHexMeshDriver::isMaxEdge
-(
- const List<pointData>& pointWallDist,
- const label edgeI,
- const scalar minCos
-) const
-{
- const pointField& points = mesh_.points();
-
- // Do not mark edges with one side on moving wall.
-
- const edge& e = mesh_.edges()[edgeI];
-
- vector v0(points[e[0]] - pointWallDist[e[0]].origin());
- scalar magV0(mag(v0));
-
- if (magV0 < SMALL)
- {
- return false;
- }
-
- vector v1(points[e[1]] - pointWallDist[e[1]].origin());
- scalar magV1(mag(v1));
-
- if (magV1 < SMALL)
- {
- return false;
- }
-
- v0 /= magV0;
- v1 /= magV1;
-
- // Test angle.
- if ((v0 & v1) < minCos)
- {
- return true;
- }
- else
- {
- return false;
- }
-}
-
-
-// Stop layer growth where mesh wraps around edge with a
-// large feature angle
-void Foam::autoHexMeshDriver::handleFeatureAngleLayerTerminations
-(
- const indirectPrimitivePatch& pp,
- const scalar minCos,
- List<extrudeMode>& extrudeStatus,
- pointField& patchDisp,
- labelList& patchNLayers,
- label& nPointCounter
-) const
-{
- // Mark faces that have all points extruded
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- boolList extrudedFaces(pp.size(), true);
-
- forAll(pp.localFaces(), faceI)
- {
- const face& f = pp.localFaces()[faceI];
-
- forAll(f, fp)
- {
- if (extrudeStatus[f[fp]] == NOEXTRUDE)
- {
- extrudedFaces[faceI] = false;
- break;
- }
- }
- }
-
-
- // Detect situation where two featureedge-neighbouring faces are partly or
- // not extruded and the edge itself is extruded. In this case unmark the
- // edge for extrusion.
-
- forAll(pp.edgeFaces(), edgeI)
- {
- const labelList& eFaces = pp.edgeFaces()[edgeI];
-
- if (eFaces.size() == 2)
- {
- const edge& e = pp.edges()[edgeI];
- label v0 = e[0];
- label v1 = e[1];
-
- if
- (
- extrudeStatus[v0] != NOEXTRUDE
- || extrudeStatus[v1] != NOEXTRUDE
- )
- {
- if (!extrudedFaces[eFaces[0]] || !extrudedFaces[eFaces[1]])
- {
- const vector& n0 = pp.faceNormals()[eFaces[0]];
- const vector& n1 = pp.faceNormals()[eFaces[1]];
-
- if ((n0 & n1) < minCos)
- {
- if
- (
- unmarkExtrusion
- (
- v0,
- patchDisp,
- patchNLayers,
- extrudeStatus
- )
- )
- {
- nPointCounter++;
- }
- if
- (
- unmarkExtrusion
- (
- v1,
- patchDisp,
- patchNLayers,
- extrudeStatus
- )
- )
- {
- nPointCounter++;
- }
- }
- }
- }
- }
- }
-}
-
-
-// Find isolated islands (points, edges and faces and layer terminations)
-// in the layer mesh and stop any layer growth at these points.
-void Foam::autoHexMeshDriver::findIsolatedRegions
-(
- const indirectPrimitivePatch& pp,
- const PackedList<1>& isMasterEdge,
- const labelList& meshEdges,
- const scalar minCosLayerTermination,
- scalarField& field,
- List<extrudeMode>& extrudeStatus,
- pointField& patchDisp,
- labelList& patchNLayers
-) const
-{
- Info<< "shrinkMeshDistance : Removing isolated regions ..." << endl;
-
- // Keep count of number of points unextruded
- label nPointCounter = 0;
-
- while(true)
- {
- // Stop layer growth where mesh wraps around edge with a
- // large feature angle
- handleFeatureAngleLayerTerminations
- (
- pp,
- minCosLayerTermination,
-
- extrudeStatus,
- patchDisp,
- patchNLayers,
- nPointCounter
- );
-
-
- // Do not extrude from point where all neighbouring
- // faces are not grown
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- boolList extrudedFaces(pp.size(), true);
- forAll(pp.localFaces(), faceI)
- {
- const face& f = pp.localFaces()[faceI];
- forAll(f, fp)
- {
- if (extrudeStatus[f[fp]] == NOEXTRUDE)
- {
- extrudedFaces[faceI] = false;
- break;
- }
- }
- }
-
- const labelListList& pointFaces = pp.pointFaces();
-
- boolList keptPoints(pp.nPoints(), false);
- forAll(keptPoints, patchPointI)
- {
- const labelList& pFaces = pointFaces[patchPointI];
-
- forAll(pFaces, i)
- {
- label faceI = pFaces[i];
- if (extrudedFaces[faceI])
- {
- keptPoints[patchPointI] = true;
- break;
- }
- }
- }
-
- syncTools::syncPointList
- (
- mesh_,
- pp.meshPoints(),
- keptPoints,
- orEqOp<bool>(),
- false, // null value
- false // no separation
- );
-
- label nChanged = 0;
-
- forAll(keptPoints, patchPointI)
- {
- if (!keptPoints[patchPointI])
- {
- if
- (
- unmarkExtrusion
- (
- patchPointI,
- patchDisp,
- patchNLayers,
- extrudeStatus
- )
- )
- {
- nPointCounter++;
- nChanged++;
- }
- }
- }
-
-
- if (returnReduce(nChanged, sumOp<label>()) == 0)
- {
- break;
- }
- }
-
- const edgeList& edges = pp.edges();
-
-
- // Count number of mesh edges using a point
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- labelList isolatedPoint(pp.nPoints(),0);
-
- forAll(edges, edgeI)
- {
- if (isMasterEdge.get(meshEdges[edgeI]) == 1)
- {
- const edge& e = edges[edgeI];
-
- label v0 = e[0];
- label v1 = e[1];
-
- if (extrudeStatus[v1] != NOEXTRUDE)
- {
- isolatedPoint[v0] += 1;
- }
- if (extrudeStatus[v0] != NOEXTRUDE)
- {
- isolatedPoint[v1] += 1;
- }
- }
- }
-
- syncTools::syncPointList
- (
- mesh_,
- pp.meshPoints(),
- isolatedPoint,
- plusEqOp<label>(),
- 0, // null value
- false // no separation
- );
-
- // stop layer growth on isolated faces
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- forAll(pp, faceI)
- {
- const face& f = pp.localFaces()[faceI];
- bool failed = false;
- forAll(f, fp)
- {
- if (isolatedPoint[f[fp]] > 2)
- {
- failed = true;
- break;
- }
- }
- bool allPointsExtruded = true;
- if (!failed)
- {
- forAll(f, fp)
- {
- if (extrudeStatus[f[fp]] == NOEXTRUDE)
- {
- allPointsExtruded = false;
- break;
- }
- }
-
- if (allPointsExtruded)
- {
- forAll(f, fp)
- {
- if
- (
- unmarkExtrusion
- (
- f[fp],
- patchDisp,
- patchNLayers,
- extrudeStatus
- )
- )
- {
- field[f[fp]] = 0.0;
- nPointCounter++;
- }
- }
- }
- }
- }
-
- reduce(nPointCounter, sumOp<label>());
- Info<< "Number isolated points extrusion stopped : "<< nPointCounter<< endl;
-}
-
-
-// Calculates medial axis fields:
-// dispVec : normal of nearest wall. Where this normal changes direction
-// defines the medial axis
-// medialDist : distance to medial axis
-// medialRatio : ratio of medial distance to wall distance.
-// (1 at wall, 0 at medial axis)
-void Foam::autoHexMeshDriver::medialAxisSmoothingInfo
-(
- const motionSmoother& meshMover,
- const label nSmoothNormals,
- const label nSmoothSurfaceNormals,
- const scalar minMedianAxisAngleCos,
-
- pointVectorField& dispVec,
- pointScalarField& medialRatio,
- pointScalarField& medialDist
-) const
-{
-
- Info<< "medialAxisSmoothingInfo :"
- << " Calculate distance to Medial Axis ..." << endl;
-
- const pointField& points = mesh_.points();
- const pointMesh& pMesh = meshMover.pMesh();
-
- const indirectPrimitivePatch& pp = meshMover.patch();
- const vectorField& faceNormals = pp.faceNormals();
- const labelList& meshPoints = pp.meshPoints();
-
- // Predetermine mesh edges
- // ~~~~~~~~~~~~~~~~~~~~~~~
-
- // Precalulate master edge (only relevant for shared edges)
- PackedList<1> isMasterEdge(syncTools::getMasterEdges(mesh_));
- // Precalculate meshEdge per pp edge
- labelList meshEdges(pp.nEdges());
-
- forAll(meshEdges, patchEdgeI)
- {
- const edge& e = pp.edges()[patchEdgeI];
-
- label v0 = pp.meshPoints()[e[0]];
- label v1 = pp.meshPoints()[e[1]];
- meshEdges[patchEdgeI] = meshTools::findEdge
- (
- mesh_.edges(),
- mesh_.pointEdges()[v0],
- v0,
- v1
- );
- }
-
-
- // Determine pointNormal
- // ~~~~~~~~~~~~~~~~~~~~~
-
- pointField pointNormals(pp.nPoints(), vector::zero);
- {
- labelList nPointFaces(pp.nPoints(), 0);
-
- forAll(faceNormals, faceI)
- {
- const face& f = pp.localFaces()[faceI];
-
- forAll(f, fp)
- {
- pointNormals[f[fp]] += faceNormals[faceI];
- nPointFaces[f[fp]] ++;
- }
- }
-
- syncTools::syncPointList
- (
- mesh_,
- meshPoints,
- pointNormals,
- plusEqOp<vector>(),
- vector::zero, // null value
- false // no separation
- );
-
- syncTools::syncPointList
- (
- mesh_,
- meshPoints,
- nPointFaces,
- plusEqOp<label>(),
- 0, // null value
- false // no separation
- );
-
- forAll(pointNormals, i)
- {
- pointNormals[i] /= nPointFaces[i];
- }
- }
-
- // Smooth patch normal vectors
- smoothPatchNormals
- (
- meshMover,
- isMasterEdge,
- meshEdges,
- nSmoothSurfaceNormals,
- pointNormals
- );
-
-
- // Calculate distance to pp points
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- // Distance to wall
- List<pointData> pointWallDist(mesh_.nPoints());
-
-
- // 1. Calculate distance to points where displacement is specified.
- {
- // Seed data.
- List<pointData> wallInfo(meshPoints.size());
-
- forAll(meshPoints, patchPointI)
- {
- label pointI = meshPoints[patchPointI];
- wallInfo[patchPointI] = pointData
- (
- points[pointI],
- 0.0,
- pointI, // passive scalar
- pointNormals[patchPointI] // surface normals
- );
- }
-
- // Do all calculations
- List<pointData> edgeWallDist(mesh_.nEdges());
- PointEdgeWave<pointData> wallDistCalc
- (
- pMesh,
- meshPoints,
- wallInfo,
- pointWallDist,
- edgeWallDist,
- mesh_.nPoints() // max iterations
- );
- }
-
- // 2. Find points with max distance and transport information back to
- // wall.
- {
- List<pointData> pointMedialDist(mesh_.nPoints());
- List<pointData> edgeMedialDist(mesh_.nEdges());
-
- // Seed point data.
- DynamicList<pointData> maxInfo(meshPoints.size());
- DynamicList<label> maxPoints(meshPoints.size());
-
- // 1. Medial axis points
-
- const edgeList& edges = mesh_.edges();
-
- forAll(edges, edgeI)
- {
- if (isMaxEdge(pointWallDist, edgeI, minMedianAxisAngleCos))
- {
- // Both end points of edge have very different nearest wall
- // point. Mark both points as medial axis points.
- const edge& e = edges[edgeI];
-
- forAll(e, ep)
- {
- label pointI = e[ep];
-
- if (!pointMedialDist[pointI].valid())
- {
- maxPoints.append(pointI);
- maxInfo.append
- (
- pointData
- (
- points[pointI],
- 0.0,
- pointI, // passive data
- vector::zero // passive data
- )
- );
- pointMedialDist[pointI] = maxInfo[maxInfo.size()-1];
- }
- }
- }
- }
-
-
- // 2. Seed non-adapt patches
- const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
- labelHashSet adaptPatches(meshMover.adaptPatchIDs());
-
- forAll(patches, patchI)
- {
- const polyPatch& pp = patches[patchI];
-
- if
- (
- !pp.coupled()
- && !isA<emptyPolyPatch>(pp)
- && !adaptPatches.found(patchI)
- )
- {
- Info<< "Inserting points on patch " << pp.name() << endl;
-
- const labelList& meshPoints = pp.meshPoints();
-
- forAll(meshPoints, i)
- {
- label pointI = meshPoints[i];
-
- if (!pointMedialDist[pointI].valid())
- {
- maxPoints.append(pointI);
- maxInfo.append
- (
- pointData
- (
- points[pointI],
- 0.0,
- pointI, // passive data
- vector::zero // passive data
- )
- );
- pointMedialDist[pointI] = maxInfo[maxInfo.size()-1];
- }
- }
- }
- }
-
- maxInfo.shrink();
- maxPoints.shrink();
-
- // Do all calculations
- PointEdgeWave<pointData> medialDistCalc
- (
- pMesh,
- maxPoints,
- maxInfo,
-
- pointMedialDist,
- edgeMedialDist,
- mesh_.nPoints() // max iterations
- );
-
- // Extract medial axis distance as pointScalarField
- forAll(pointMedialDist, pointI)
- {
- medialDist[pointI] = Foam::sqrt(pointMedialDist[pointI].distSqr());
- }
- }
-
- // Extract transported surface normals as pointVectorField
- forAll(dispVec, i)
- {
- dispVec[i] = pointWallDist[i].v();
- }
-
- // Smooth normal vectors. Do not change normals on pp.meshPoints
- smoothNormals(nSmoothNormals, isMasterEdge, meshPoints, dispVec);
-
- // Calculate ratio point medial distance to point wall distance
- forAll(medialRatio, pointI)
- {
- scalar wDist2 = pointWallDist[pointI].distSqr();
- scalar mDist = medialDist[pointI];
-
- if (wDist2 < sqr(SMALL) && mDist < SMALL)
- {
- medialRatio[pointI] = 0.0;
- }
- else
- {
- medialRatio[pointI] = mDist / (Foam::sqrt(wDist2) + mDist);
- }
- }
-
- if (debug_)
- {
- Info<< "medialAxisSmoothingInfo :"
- << " Writing:" << nl
- << " " << dispVec.name()
- << " : normalised direction of nearest displacement" << nl
- << " " << medialDist.name()
- << " : distance to medial axis" << nl
- << " " << medialRatio.name()
- << " : ratio of medial distance to wall distance" << nl
- << endl;
- dispVec.write();
- medialDist.write();
- medialRatio.write();
- }
-}
-
-
-void Foam::autoHexMeshDriver::shrinkMeshMedialDistance
-(
- motionSmoother& meshMover,
- const label nSmoothThickness,
- const scalar maxThicknessToMedialRatio,
- const label nAllowableErrors,
- const label nSnap,
- const scalar minCosLayerTermination,
-
- const scalarField& layerThickness,
- const scalarField& minThickness,
-
- const pointVectorField& dispVec,
- const pointScalarField& medialRatio,
- const pointScalarField& medialDist,
-
- List<extrudeMode>& extrudeStatus,
- pointField& patchDisp,
- labelList& patchNLayers
-) const
-{
- Info<< "shrinkMeshMedialDistance : Smoothing using Medial Axis ..." << endl;
-
- const polyMesh& mesh = meshMover.mesh();
- const pointMesh& pMesh = meshMover.pMesh();
-
- const indirectPrimitivePatch& pp = meshMover.patch();
- const labelList& meshPoints = pp.meshPoints();
-
- // Precalulate master edge (only relevant for shared edges)
- PackedList<1> isMasterEdge(syncTools::getMasterEdges(mesh_));
- // Precalculate meshEdge per pp edge
- labelList meshEdges(pp.nEdges());
-
- forAll(meshEdges, patchEdgeI)
- {
- const edge& e = pp.edges()[patchEdgeI];
-
- label v0 = pp.meshPoints()[e[0]];
- label v1 = pp.meshPoints()[e[1]];
- meshEdges[patchEdgeI] = meshTools::findEdge
- (
- mesh_.edges(),
- mesh_.pointEdges()[v0],
- v0,
- v1
- );
- }
-
-
- scalarField thickness(layerThickness.size());
-
- thickness = mag(patchDisp);
-
- forAll(thickness, patchPointI)
- {
- if (extrudeStatus[patchPointI] == NOEXTRUDE)
- {
- thickness[patchPointI] = 0.0;
- }
- }
-
- label numThicknessRatioExclude = 0;
-
- // reduce thickness where thickness/medial axis distance large
- forAll(meshPoints, patchPointI)
- {
- if (extrudeStatus[patchPointI] != NOEXTRUDE)
- {
- label pointI = meshPoints[patchPointI];
-
- scalar mDist = medialDist[pointI];
-
- scalar thicknessRatio = thickness[patchPointI]/(mDist+VSMALL);
-
- if (thicknessRatio > maxThicknessToMedialRatio)
- {
- // Truncate thickness.
- thickness[patchPointI] =
- 0.5*(minThickness[patchPointI]+thickness[patchPointI]);
-
- patchDisp[patchPointI] =
- thickness[patchPointI]
- * patchDisp[patchPointI]
- / (mag(patchDisp[patchPointI]) + VSMALL);
- numThicknessRatioExclude++;
- }
- }
- }
-
- reduce(numThicknessRatioExclude, sumOp<label>());
- Info<< "shrinkMeshMedialDistance : Reduce layer thickness at "
- << numThicknessRatioExclude
- << " nodes where thickness to medial axis distance is large " << endl;
-
- // find points where layer growth isolated to a lone point, edge or face
-
- findIsolatedRegions
- (
- pp,
- isMasterEdge,
- meshEdges,
- minCosLayerTermination,
- thickness,
- extrudeStatus,
- patchDisp,
- patchNLayers
- );
-
- // smooth layer thickness on moving patch
- smoothField
- (
- meshMover,
- isMasterEdge,
- meshEdges,
- minThickness,
- nSmoothThickness,
- thickness
- );
-
- List<pointData> pointWallDist(mesh.nPoints());
-
- const pointField& points = mesh.points();
- // 1. Calculate distance to points where displacement is specified.
- // This wave is used to transport layer thickness
- {
- // Distance to wall and medial axis on edges.
- List<pointData> edgeWallDist(mesh.nEdges());
- labelList wallPoints(meshPoints.size());
-
- // Seed data.
- List<pointData> wallInfo(meshPoints.size());
-
- forAll(meshPoints, patchPointI)
- {
- label pointI = meshPoints[patchPointI];
- wallPoints[patchPointI] = pointI;
- wallInfo[patchPointI] = pointData
- (
- points[pointI],
- 0.0,
- thickness[patchPointI], // transport layer thickness
- vector::zero // passive vector
- );
- }
-
- // Do all calculations
- PointEdgeWave<pointData> wallDistCalc
- (
- pMesh,
- wallPoints,
- wallInfo,
- pointWallDist,
- edgeWallDist,
- mesh.nPoints() // max iterations
- );
- }
-
- // Calculate scaled displacement vector
- pointVectorField& displacement = meshMover.displacement();
-
- forAll(displacement, pointI)
- {
- // 1) displacement on nearest wall point, scaled by medialRatio
- // (wall distance / medial distance)
- // 2) pointWallDist[pointI].s() is layer thickness transported
- // from closest wall point.
- // 3) shrink in opposite direction of addedPoints
- displacement[pointI] =
- -medialRatio[pointI]
- *pointWallDist[pointI].s()
- *dispVec[pointI];
- }
-
- // Current faces to check. Gets modified in meshMover.scaleMesh
- labelList checkFaces(identity(mesh.nFaces()));
-
- Info<< "shrinkMeshMedialDistance : Moving mesh ..." << endl;
- scalar oldErrorReduction = -1;
-
- for (label iter = 0; iter < 2*nSnap ; iter++)
- {
- Info<< "Iteration " << iter << endl;
- if (iter == nSnap)
- {
- Info<< "Displacement scaling for error reduction set to 0." << endl;
- oldErrorReduction = meshMover.setErrorReduction(0.0);
- }
-
- if (meshMover.scaleMesh(checkFaces, true, nAllowableErrors))
- {
- Info<< "shrinkMeshMedialDistance : Successfully moved mesh" << endl;
- break;
- }
- }
-
- if (oldErrorReduction >= 0)
- {
- meshMover.setErrorReduction(oldErrorReduction);
- }
-
- Info<< "shrinkMeshMedialDistance : Finished moving mesh ..." << endl;
-}
-
-
-// ************************************************************************* //
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverSnap.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverSnap.C
deleted file mode 100644
index 8e074cd05e2506da9f76c2fa877cbda00ce1246e..0000000000000000000000000000000000000000
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverSnap.C
+++ /dev/null
@@ -1,1245 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software; you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by the
- Free Software Foundation; either version 2 of the License, or (at your
- option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM; if not, write to the Free Software Foundation,
- Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
- All to do with snapping to the surface
-
-\*----------------------------------------------------------------------------*/
-
-#include "autoHexMeshDriver.H"
-#include "syncTools.H"
-#include "fvMesh.H"
-#include "Time.H"
-#include "OFstream.H"
-#include "mapPolyMesh.H"
-#include "motionSmoother.H"
-#include "pointEdgePoint.H"
-#include "PointEdgeWave.H"
-#include "mergePoints.H"
-#include "snapParameters.H"
-
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
-void Foam::autoHexMeshDriver::getZonedSurfaces
-(
- labelList& zonedSurfaces,
- labelList& unzonedSurfaces
-) const
-{ // Surfaces with zone information
- const wordList& faceZoneNames = surfaces().faceZoneNames();
-
- zonedSurfaces.setSize(faceZoneNames.size());
- label zonedI = 0;
- unzonedSurfaces.setSize(faceZoneNames.size());
- label unzonedI = 0;
-
- forAll(faceZoneNames, surfI)
- {
- if (faceZoneNames[surfI].size() > 0)
- {
- zonedSurfaces[zonedI++] = surfI;
- }
- else
- {
- unzonedSurfaces[unzonedI++] = surfI;
- }
- }
- zonedSurfaces.setSize(zonedI);
- unzonedSurfaces.setSize(unzonedI);
-}
-
-
-// Get faces to repatch. Returns map from face to patch.
-Foam::Map<Foam::label> Foam::autoHexMeshDriver::getZoneBafflePatches
-(
- const bool allowBoundary
-) const
-{
- Map<label> bafflePatch(mesh_.nFaces()/1000);
-
- const wordList& faceZoneNames = surfaces().faceZoneNames();
- const faceZoneMesh& fZones = mesh_.faceZones();
-
- forAll(faceZoneNames, surfI)
- {
- if (faceZoneNames[surfI].size() > 0)
- {
- // Get zone
- label zoneI = fZones.findZoneID(faceZoneNames[surfI]);
-
- const faceZone& fZone = fZones[zoneI];
-
- //// Get patch allocated for zone
- //label patchI = surfaceToCyclicPatch_[surfI];
- // Get patch of (first region) of surface
- label patchI = globalToPatch_[surfaces().globalRegion(surfI, 0)];
-
- Info<< "For surface "
- << surfaces().names()[surfI]
- << " found faceZone " << fZone.name()
- << " and patch " << mesh_.boundaryMesh()[patchI].name()
- << endl;
-
-
- forAll(fZone, i)
- {
- label faceI = fZone[i];
-
- if (allowBoundary || mesh_.isInternalFace(faceI))
- {
- if (!bafflePatch.insert(faceI, patchI))
- {
- label oldPatchI = bafflePatch[faceI];
-
- if (oldPatchI != patchI)
- {
- FatalErrorIn("getZoneBafflePatches(const bool)")
- << "Face " << faceI
- << " fc:" << mesh_.faceCentres()[faceI]
- << " is in faceZone "
- << mesh_.boundaryMesh()[oldPatchI].name()
- << " and in faceZone "
- << mesh_.boundaryMesh()[patchI].name()
- << abort(FatalError);
- }
- }
- }
- }
- }
- }
- return bafflePatch;
-}
-
-
-// Calculate geometrically collocated points, Requires PackedList to be
-// sizes and initalised!
-Foam::label Foam::autoHexMeshDriver::getCollocatedPoints
-(
- const scalar tol,
- const pointField& points,
- PackedList<1>& isCollocatedPoint
-)
-{
- labelList pointMap;
- pointField newPoints;
- bool hasMerged = mergePoints
- (
- points, // points
- tol, // mergeTol
- false, // verbose
- pointMap,
- newPoints
- );
-
- if (!returnReduce(hasMerged, orOp<bool>()))
- {
- return 0;
- }
-
- // Determine which newPoints are referenced more than once
- label nCollocated = 0;
-
- // Per old point the newPoint. Or -1 (not set yet) or -2 (already seen
- // twice)
- labelList firstOldPoint(newPoints.size(), -1);
- forAll(pointMap, oldPointI)
- {
- label newPointI = pointMap[oldPointI];
-
- if (firstOldPoint[newPointI] == -1)
- {
- // First use of oldPointI. Store.
- firstOldPoint[newPointI] = oldPointI;
- }
- else if (firstOldPoint[newPointI] == -2)
- {
- // Third or more reference of oldPointI -> non-manifold
- isCollocatedPoint.set(oldPointI, 1u);
- nCollocated++;
- }
- else
- {
- // Second reference of oldPointI -> non-manifold
- isCollocatedPoint.set(firstOldPoint[newPointI], 1u);
- nCollocated++;
-
- isCollocatedPoint.set(oldPointI, 1u);
- nCollocated++;
-
- // Mark with special value to save checking next time round
- firstOldPoint[newPointI] = -2;
- }
- }
- return returnReduce(nCollocated, sumOp<label>());
-}
-
-
-// Calculate displacement as average of patch points.
-Foam::pointField Foam::autoHexMeshDriver::smoothPatchDisplacement
-(
- const motionSmoother& meshMover
-) const
-{
- const indirectPrimitivePatch& pp = meshMover.patch();
-
- // Calculate geometrically non-manifold points on the patch to be moved.
- PackedList<1> nonManifoldPoint(pp.nPoints());
- label nNonManifoldPoints = getCollocatedPoints
- (
- SMALL,
- pp.localPoints(),
- nonManifoldPoint
- );
- Info<< "Found " << nNonManifoldPoints << " non-mainfold point(s)."
- << endl;
-
-
- // Average points
- // ~~~~~~~~~~~~~~
-
- // We determine three points:
- // - average of (centres of) connected patch faces
- // - average of (centres of) connected internal mesh faces
- // - as fallback: centre of any connected cell
- // so we can do something moderately sensible for non/manifold points.
-
- // Note: the averages are calculated properly parallel. This is
- // necessary to get the points shared by processors correct.
-
-
- const labelListList& pointFaces = pp.pointFaces();
- const labelList& meshPoints = pp.meshPoints();
- const pointField& points = pp.points();
- const polyMesh& mesh = meshMover.mesh();
-
-
- // Get average position of boundary face centres
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- vectorField avgBoundary(pointFaces.size(), vector::zero);
- labelList nBoundary(pointFaces.size(), 0);
-
- forAll(pointFaces, patchPointI)
- {
- const labelList& pFaces = pointFaces[patchPointI];
-
- forAll(pFaces, pfI)
- {
- avgBoundary[patchPointI] += pp[pFaces[pfI]].centre(points);
- }
- nBoundary[patchPointI] = pFaces.size();
- }
-
- syncTools::syncPointList
- (
- mesh,
- pp.meshPoints(),
- avgBoundary,
- plusEqOp<point>(), // combine op
- vector::zero, // null value
- false // no separation
- );
- syncTools::syncPointList
- (
- mesh,
- pp.meshPoints(),
- nBoundary,
- plusEqOp<label>(), // combine op
- 0, // null value
- false // no separation
- );
-
- forAll(avgBoundary, i)
- {
- avgBoundary[i] /= nBoundary[i];
- }
-
-
- // Get average position of internal face centres
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- vectorField avgInternal;
- labelList nInternal;
- {
- vectorField globalSum(mesh.nPoints(), vector::zero);
- labelList globalNum(mesh.nPoints(), 0);
-
- // Note: no use of pointFaces
- const faceList& faces = mesh.faces();
-
- for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
- {
- const face& f = faces[faceI];
- const point& fc = mesh.faceCentres()[faceI];
-
- forAll(f, fp)
- {
- globalSum[f[fp]] += fc;
- globalNum[f[fp]]++;
- }
- }
-
- // Count coupled faces as internal ones (but only once)
- const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
- forAll(patches, patchI)
- {
- if (Pstream::parRun() && isA<processorPolyPatch>(patches[patchI]))
- {
- const processorPolyPatch& pp =
- refCast<const processorPolyPatch>(patches[patchI]);
-
- if (pp.myProcNo() < pp.neighbProcNo())
- {
- const vectorField::subField faceCentres = pp.faceCentres();
-
- forAll(pp, i)
- {
- const face& f = pp[i];
- const point& fc = faceCentres[i];
-
- forAll(f, fp)
- {
- globalSum[f[fp]] += fc;
- globalNum[f[fp]]++;
- }
- }
- }
- }
- else if (isA<cyclicPolyPatch>(patches[patchI]))
- {
- const cyclicPolyPatch& pp =
- refCast<const cyclicPolyPatch>(patches[patchI]);
-
- const vectorField::subField faceCentres = pp.faceCentres();
-
- for (label i = 0; i < pp.size()/2; i++)
- {
- const face& f = pp[i];
- const point& fc = faceCentres[i];
-
- forAll(f, fp)
- {
- globalSum[f[fp]] += fc;
- globalNum[f[fp]]++;
- }
- }
- }
- }
-
- syncTools::syncPointList
- (
- mesh,
- globalSum,
- plusEqOp<vector>(), // combine op
- vector::zero, // null value
- false // no separation
- );
- syncTools::syncPointList
- (
- mesh,
- globalNum,
- plusEqOp<label>(), // combine op
- 0, // null value
- false // no separation
- );
-
- avgInternal.setSize(meshPoints.size());
- nInternal.setSize(meshPoints.size());
-
- forAll(avgInternal, patchPointI)
- {
- label meshPointI = meshPoints[patchPointI];
-
- nInternal[patchPointI] = globalNum[meshPointI];
-
- if (nInternal[patchPointI] == 0)
- {
- avgInternal[patchPointI] = globalSum[meshPointI];
- }
- else
- {
- avgInternal[patchPointI] =
- globalSum[meshPointI]
- / nInternal[patchPointI];
- }
- }
- }
-
-
- // Precalculate any cell using mesh point (replacement of pointCells()[])
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- labelList anyCell(mesh.nPoints(), -1);
- forAll(mesh.faceNeighbour(), faceI)
- {
- label own = mesh.faceOwner()[faceI];
- const face& f = mesh.faces()[faceI];
-
- forAll(f, fp)
- {
- anyCell[f[fp]] = own;
- }
- }
- for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
- {
- label own = mesh.faceOwner()[faceI];
-
- const face& f = mesh.faces()[faceI];
-
- forAll(f, fp)
- {
- anyCell[f[fp]] = own;
- }
- }
-
-
- // Displacement to calculate.
- pointField patchDisp(meshPoints.size(), vector::zero);
-
- forAll(pointFaces, i)
- {
- label meshPointI = meshPoints[i];
- const point& currentPos = pp.points()[meshPointI];
-
- // Now we have the two average points: avgBoundary and avgInternal
- // and how many boundary/internal faces connect to the point
- // (nBoundary, nInternal)
- // Do some blending between the two.
- // Note: the following section has some reasoning behind it but the
- // blending factors can be experimented with.
-
- point newPos;
-
- if (nonManifoldPoint.get(i) == 0u)
- {
- // Points that are manifold. Weight the internal and boundary
- // by their number of faces and blend with
- scalar internalBlend = 0.1;
- scalar blend = 0.1;
-
- point avgPos =
- (
- internalBlend*nInternal[i]*avgInternal[i]
- +(1-internalBlend)*nBoundary[i]*avgBoundary[i]
- )
- / (internalBlend*nInternal[i]+(1-internalBlend)*nBoundary[i]);
-
- newPos = (1-blend)*avgPos + blend*currentPos;
- }
- else if (nInternal[i] == 0)
- {
- // Non-manifold without internal faces. Use any connected cell
- // as internal point instead. Use precalculated any cell to avoid
- // e.g. pointCells()[meshPointI][0]
-
- const point& cc = mesh.cellCentres()[anyCell[meshPointI]];
-
- scalar cellCBlend = 0.8;
- scalar blend = 0.1;
-
- point avgPos = (1-cellCBlend)*avgBoundary[i] + cellCBlend*cc;
-
- newPos = (1-blend)*avgPos + blend*currentPos;
- }
- else
- {
- // Non-manifold point with internal faces connected to them
- scalar internalBlend = 0.9;
- scalar blend = 0.1;
-
- point avgPos =
- internalBlend*avgInternal[i]
- + (1-internalBlend)*avgBoundary[i];
-
- newPos = (1-blend)*avgPos + blend*currentPos;
- }
-
- patchDisp[i] = newPos - currentPos;
- }
-
- return patchDisp;
-}
-
-
-Foam::tmp<Foam::scalarField> Foam::autoHexMeshDriver::edgePatchDist
-(
- const pointMesh& pMesh,
- const indirectPrimitivePatch& pp
-)
-{
- const polyMesh& mesh = pMesh();
-
- // Set initial changed points to all the patch points
- List<pointEdgePoint> wallInfo(pp.nPoints());
-
- forAll(pp.localPoints(), ppI)
- {
- wallInfo[ppI] = pointEdgePoint(pp.localPoints()[ppI], 0.0);
- }
-
- // Current info on points
- List<pointEdgePoint> allPointInfo(mesh.nPoints());
-
- // Current info on edges
- List<pointEdgePoint> allEdgeInfo(mesh.nEdges());
-
- PointEdgeWave<pointEdgePoint> wallCalc
- (
- pMesh,
- pp.meshPoints(),
- wallInfo,
-
- allPointInfo,
- allEdgeInfo,
- mesh.globalData().nTotalPoints() // max iterations
- );
-
- // Copy edge values into scalarField
- tmp<scalarField> tedgeDist(new scalarField(mesh.nEdges()));
- scalarField& edgeDist = tedgeDist();
-
- forAll(allEdgeInfo, edgeI)
- {
- edgeDist[edgeI] = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
- }
-
-
- //{
- // // For debugging: dump to file
- // pointScalarField pointDist
- // (
- // IOobject
- // (
- // "pointDist",
- // mesh.DB().timeName(),
- // mesh.DB(),
- // IOobject::NO_READ,
- // IOobject::AUTO_WRITE
- // ),
- // pMesh,
- // dimensionedScalar("pointDist", dimless, 0.0)
- // );
- //
- // forAll(allEdgeInfo, edgeI)
- // {
- // scalar d = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
- //
- // const edge& e = mesh.edges()[edgeI];
- //
- // pointDist[e[0]] += d;
- // pointDist[e[1]] += d;
- // }
- // forAll(pointDist, pointI)
- // {
- // pointDist[pointI] /= mesh.pointEdges()[pointI].size();
- // }
- // Info<< "Writing patch distance to " << pointDist.name()
- // << " at time " << mesh.DB().timeName() << endl;
- //
- // pointDist.write();
- //}
-
- return tedgeDist;
-}
-
-
-void Foam::autoHexMeshDriver::dumpMove
-(
- const fileName& fName,
- const pointField& meshPts,
- const pointField& surfPts
-)
-{
- // Dump direction of growth into file
- Pout<< nl << "Dumping move direction to " << fName << nl
- << "View this Lightwave-OBJ file with e.g. javaview" << nl
- << endl;
-
- OFstream nearestStream(fName);
-
- label vertI = 0;
-
- forAll(meshPts, ptI)
- {
- meshTools::writeOBJ(nearestStream, meshPts[ptI]);
- vertI++;
-
- meshTools::writeOBJ(nearestStream, surfPts[ptI]);
- vertI++;
-
- nearestStream<< "l " << vertI-1 << ' ' << vertI << nl;
- }
-}
-
-
-// Check whether all displacement vectors point outwards of patch. Return true
-// if so.
-bool Foam::autoHexMeshDriver::outwardsDisplacement
-(
- const indirectPrimitivePatch& pp,
- const vectorField& patchDisp
-)
-{
- const vectorField& faceNormals = pp.faceNormals();
- const labelListList& pointFaces = pp.pointFaces();
-
- forAll(pointFaces, pointI)
- {
- const labelList& pFaces = pointFaces[pointI];
-
- vector disp(patchDisp[pointI]);
-
- scalar magDisp = mag(disp);
-
- if (magDisp > SMALL)
- {
- disp /= magDisp;
-
- bool outwards = meshTools::visNormal(disp, faceNormals, pFaces);
-
- if (!outwards)
- {
- Warning<< "Displacement " << patchDisp[pointI]
- << " at mesh point " << pp.meshPoints()[pointI]
- << " coord " << pp.points()[pp.meshPoints()[pointI]]
- << " points through the surrounding patch faces" << endl;
- return false;
- }
- }
- else
- {
- //? Displacement small but in wrong direction. Would probably be ok.
- }
- }
- return true;
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-Foam::autoPtr<Foam::mapPolyMesh> Foam::autoHexMeshDriver::createZoneBaffles
-(
- List<labelPair>& baffles
-)
-{
- labelList zonedSurfaces;
- labelList unzonedSurfaces;
- getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
-
- autoPtr<mapPolyMesh> map;
-
- // No need to sync; all processors will have all same zonedSurfaces.
- if (zonedSurfaces.size() > 0)
- {
- // Split internal faces on interface surfaces
- Info<< "Converting zoned faces into baffles ..." << endl;
-
- // Get faces (internal only) to be baffled. Map from face to patch
- // label.
- Map<label> faceToPatch(getZoneBafflePatches(false));
-
- label nZoneFaces = returnReduce(faceToPatch.size(), sumOp<label>());
- if (nZoneFaces > 0)
- {
- // Convert into labelLists
- labelList ownPatch(mesh_.nFaces(), -1);
- forAllConstIter(Map<label>, faceToPatch, iter)
- {
- ownPatch[iter.key()] = iter();
- }
-
- // Create baffles. both sides same patch.
- map = meshRefinerPtr_().createBaffles(ownPatch, ownPatch);
-
- // Get pairs of faces created.
- // Just loop over faceMap and store baffle if we encounter a slave
- // face.
-
- baffles.setSize(faceToPatch.size());
- label baffleI = 0;
-
- const labelList& faceMap = map().faceMap();
- const labelList& reverseFaceMap = map().reverseFaceMap();
-
- forAll(faceMap, faceI)
- {
- label oldFaceI = faceMap[faceI];
-
- // Does face originate from face-to-patch
- Map<label>::const_iterator iter = faceToPatch.find(oldFaceI);
-
- if (iter != faceToPatch.end())
- {
- label masterFaceI = reverseFaceMap[oldFaceI];
- if (faceI != masterFaceI)
- {
- baffles[baffleI++] = labelPair(masterFaceI, faceI);
- }
- }
- }
-
- if (baffleI != faceToPatch.size())
- {
- FatalErrorIn("autoHexMeshDriver::createZoneBaffles(..)")
- << "Had " << faceToPatch.size() << " patches to create "
- << " but encountered " << baffleI
- << " slave faces originating from patcheable faces."
- << abort(FatalError);
- }
-
- if (debug_)
- {
- const_cast<Time&>(mesh_.time())++;
- Pout<< "Writing baffled mesh to time "
- << mesh_.time().timeName() << endl;
- mesh_.write();
- }
- }
- Info<< "Created " << nZoneFaces << " baffles in = "
- << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
- }
- return map;
-}
-
-
-Foam::autoPtr<Foam::mapPolyMesh> Foam::autoHexMeshDriver::mergeZoneBaffles
-(
- const List<labelPair>& baffles
-)
-{
- labelList zonedSurfaces;
- labelList unzonedSurfaces;
- getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
-
- autoPtr<mapPolyMesh> map;
-
- // No need to sync; all processors will have all same zonedSurfaces.
- label nBaffles = returnReduce(baffles.size(), sumOp<label>());
- if (zonedSurfaces.size() > 0 && nBaffles > 0)
- {
- // Merge any baffles
- Info<< "Converting " << nBaffles << " baffles back into zoned faces ..."
- << endl;
-
- map = meshRefinerPtr_().mergeBaffles(baffles);
-
- Info<< "Converted baffles in = "
- << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
- }
-
- return map;
-}
-
-
-Foam::scalarField Foam::autoHexMeshDriver::calcSnapDistance
-(
- const snapParameters& snapParams,
- const indirectPrimitivePatch& pp
-) const
-{
- const edgeList& edges = pp.edges();
- const labelListList& pointEdges = pp.pointEdges();
- const pointField& localPoints = pp.localPoints();
-
- scalarField maxEdgeLen(localPoints.size(), -GREAT);
-
- forAll(pointEdges, pointI)
- {
- const labelList& pEdges = pointEdges[pointI];
-
- forAll(pEdges, pEdgeI)
- {
- const edge& e = edges[pEdges[pEdgeI]];
-
- scalar len = e.mag(localPoints);
-
- maxEdgeLen[pointI] = max(maxEdgeLen[pointI], len);
- }
- }
-
- syncTools::syncPointList
- (
- mesh_,
- pp.meshPoints(),
- maxEdgeLen,
- maxEqOp<scalar>(), // combine op
- -GREAT, // null value
- false // no separation
- );
-
- return snapParams.snapTol()*maxEdgeLen;
-}
-
-
-// Invert globalToPatch_ to get the patches related to surfaces.
-Foam::labelList Foam::autoHexMeshDriver::getSurfacePatches() const
-{
- // Set of patches originating from surface
- labelHashSet surfacePatchSet(globalToPatch_.size());
-
- forAll(globalToPatch_, i)
- {
- if (globalToPatch_[i] != -1)
- {
- surfacePatchSet.insert(globalToPatch_[i]);
- }
- }
-
- DynamicList<label> surfacePatches(surfacePatchSet.size());
-
- for (label patchI = 0; patchI < mesh_.boundaryMesh().size(); patchI++)
- {
- if (surfacePatchSet.found(patchI))
- {
- surfacePatches.append(patchI);
- }
- }
- return surfacePatches.shrink();
-}
-
-
-void Foam::autoHexMeshDriver::preSmoothPatch
-(
- const snapParameters& snapParams,
- const label nInitErrors,
- const List<labelPair>& baffles,
- motionSmoother& meshMover
-) const
-{
- labelList checkFaces;
-
- Info<< "Smoothing patch points ..." << endl;
- for
- (
- label smoothIter = 0;
- smoothIter < snapParams.nSmoothPatch();
- smoothIter++
- )
- {
- Info<< "Smoothing iteration " << smoothIter << endl;
- checkFaces.setSize(mesh_.nFaces());
- forAll(checkFaces, faceI)
- {
- checkFaces[faceI] = faceI;
- }
-
- pointField patchDisp(smoothPatchDisplacement(meshMover));
-
- // The current mesh is the starting mesh to smooth from.
- meshMover.setDisplacement(patchDisp);
- meshMover.correct();
-
- scalar oldErrorReduction = -1;
-
- for (label snapIter = 0; snapIter < 2*snapParams.nSnap(); snapIter++)
- {
- Info<< nl << "Scaling iteration " << snapIter << endl;
-
- if (snapIter == snapParams.nSnap())
- {
- Info<< "Displacement scaling for error reduction set to 0."
- << endl;
- oldErrorReduction = meshMover.setErrorReduction(0.0);
- }
-
- // Try to adapt mesh to obtain displacement by smoothly
- // decreasing displacement at error locations.
- if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
- {
- Info<< "Successfully moved mesh" << endl;
- break;
- }
- }
-
- if (oldErrorReduction >= 0)
- {
- meshMover.setErrorReduction(oldErrorReduction);
- }
- Info<< endl;
- }
-
-
- // The current mesh is the starting mesh to smooth from.
- meshMover.correct();
-
- if (debug_)
- {
- const_cast<Time&>(mesh_.time())++;
- Pout<< "Writing patch smoothed mesh to time " << mesh_.time().timeName()
- << endl;
- mesh_.write();
- }
-
- Info<< "Patch points smoothed in = "
- << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
-}
-
-
-// Get (pp-local) indices of points that are both on zone and on patched surface
-Foam::labelList Foam::autoHexMeshDriver::getZoneSurfacePoints
-(
- const indirectPrimitivePatch& pp,
- const word& zoneName
-) const
-{
- label zoneI = mesh_.faceZones().findZoneID(zoneName);
-
- if (zoneI == -1)
- {
- FatalErrorIn
- (
- "autoHexMeshDriver::getZoneSurfacePoints"
- "(const indirectPrimitivePatch&, const word&)"
- ) << "Cannot find zone " << zoneName
- << exit(FatalError);
- }
-
- const faceZone& fZone = mesh_.faceZones()[zoneI];
-
-
- // Could use PrimitivePatch & localFaces to extract points but might just
- // as well do it ourselves.
-
- boolList pointOnZone(pp.nPoints(), false);
-
- forAll(fZone, i)
- {
- const face& f = mesh_.faces()[fZone[i]];
-
- forAll(f, fp)
- {
- label meshPointI = f[fp];
-
- Map<label>::const_iterator iter =
- pp.meshPointMap().find(meshPointI);
-
- if (iter != pp.meshPointMap().end())
- {
- label pointI = iter();
- pointOnZone[pointI] = true;
- }
- }
- }
-
- return findIndices(pointOnZone, true);
-}
-
-
-Foam::vectorField Foam::autoHexMeshDriver::calcNearestSurface
-(
- const scalarField& snapDist,
- motionSmoother& meshMover
-) const
-{
- Info<< "Calculating patchDisplacement as distance to nearest surface"
- << " point ..." << endl;
-
- const indirectPrimitivePatch& pp = meshMover.patch();
- const pointField& localPoints = pp.localPoints();
-
- // Divide surfaces into zoned and unzoned
- labelList zonedSurfaces;
- labelList unzonedSurfaces;
- getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
-
- // Displacement per patch point
- vectorField patchDisp(localPoints.size(), vector::zero);
-
-
- // 1. All points to non-interface surfaces
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- {
- List<pointIndexHit> hitInfo;
- labelList hitSurface;
- surfaces().findNearest
- (
- unzonedSurfaces,
- localPoints,
- sqr(4*snapDist), // sqr of attract distance
- hitSurface,
- hitInfo
- );
-
- forAll(hitInfo, pointI)
- {
- if (hitInfo[pointI].hit())
- {
- patchDisp[pointI] =
- hitInfo[pointI].hitPoint()
- - localPoints[pointI];
- }
- //else
- //{
- // WarningIn("autoHexMeshDriver::calcNearestSurface(..)")
- // << "For point:" << pointI
- // << " coordinate:" << localPoints[pointI]
- // << " did not find any surface within:"
- // << 4*snapDist[pointI]
- // << " meter." << endl;
- //}
- }
- }
-
-
-
- // 2. All points on zones to their respective surface
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- // Surfaces with zone information
- const wordList& faceZoneNames = surfaces().faceZoneNames();
-
- forAll(zonedSurfaces, i)
- {
- label zoneSurfI = zonedSurfaces[i];
-
- const labelList surfacesToTest(1, zoneSurfI);
-
- // Get indices of points both on faceZone and on pp.
- labelList zonePointIndices
- (
- getZoneSurfacePoints
- (
- pp,
- faceZoneNames[zoneSurfI]
- )
- );
-
- pointField zonePoints(zonePointIndices.size());
- forAll(zonePointIndices, i)
- {
- zonePoints[i] = localPoints[zonePointIndices[i]];
- }
-
- // Find nearest for points both on faceZone and pp.
- List<pointIndexHit> hitInfo;
- labelList hitSurface;
- surfaces().findNearest
- (
- labelList(1, zoneSurfI),
- zonePoints,
- sqr(4*snapDist),
- hitSurface,
- hitInfo
- );
-
- forAll(hitInfo, pointI)
- {
- if (hitInfo[pointI].hit())
- {
- patchDisp[pointI] =
- hitInfo[pointI].hitPoint()
- - localPoints[pointI];
- }
- else
- {
- WarningIn("autoHexMeshDriver::calcNearestSurface(..)")
- << "For point:" << pointI
- << " coordinate:" << localPoints[pointI]
- << " did not find any surface within:"
- << 4*snapDist[pointI]
- << " meter." << endl;
- }
- }
- }
-
-
- {
- scalarField magDisp(mag(patchDisp));
-
- Info<< "Wanted displacement : average:"
- << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
- << " min:" << gMin(magDisp)
- << " max:" << gMax(magDisp) << endl;
- }
-
- Info<< "Calculated surface displacement in = "
- << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
-
-
- // Limit amount of movement.
- forAll(patchDisp, patchPointI)
- {
- scalar magDisp = mag(patchDisp[patchPointI]);
-
- if (magDisp > snapDist[patchPointI])
- {
- patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp;
-
- Pout<< "Limiting displacement for " << patchPointI
- << " from " << magDisp << " to " << snapDist[patchPointI]
- << endl;
- }
- }
-
- // Points on zones in one domain but only present as point on other
- // will not do condition 2 on all. Sync explicitly.
- syncTools::syncPointList
- (
- mesh_,
- pp.meshPoints(),
- patchDisp,
- minMagEqOp(), // combine op
- vector(GREAT, GREAT, GREAT), // null value
- false // no separation
- );
-
-
- // Check for displacement being outwards.
- outwardsDisplacement(pp, patchDisp);
-
- // Set initial distribution of displacement field (on patches) from
- // patchDisp and make displacement consistent with b.c. on displacement
- // pointVectorField.
- meshMover.setDisplacement(patchDisp);
-
- if (debug_)
- {
- dumpMove
- (
- mesh_.time().path()/"patchDisplacement.obj",
- pp.localPoints(),
- pp.localPoints() + patchDisp
- );
- }
-
- return patchDisp;
-}
-
-
-void Foam::autoHexMeshDriver::smoothDisplacement
-(
- const snapParameters& snapParams,
- motionSmoother& meshMover
-) const
-{
- const pointMesh& pMesh = meshMover.pMesh();
- const indirectPrimitivePatch& pp = meshMover.patch();
-
- Info<< "Smoothing displacement ..." << endl;
-
- // Set edge diffusivity as inverse of distance to patch
- scalarField edgeGamma(1.0/(edgePatchDist(pMesh, pp) + SMALL));
- //scalarField edgeGamma(mesh_.nEdges(), 1.0);
- //scalarField edgeGamma(wallGamma(mesh, pp, 10, 1));
-
- // Get displacement field
- pointVectorField& disp = meshMover.displacement();
-
- for (label iter = 0; iter < snapParams.nSmoothDispl(); iter++)
- {
- if ((iter % 10) == 0)
- {
- Info<< "Iteration " << iter << endl;
- }
- pointVectorField oldDisp(disp);
-
- meshMover.smooth(oldDisp, edgeGamma, false, disp);
- }
- Info<< "Displacement smoothed in = "
- << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
-
- if (debug_)
- {
- const_cast<Time&>(mesh_.time())++;
- Pout<< "Writing smoothed mesh to time " << mesh_.time().timeName()
- << endl;
- mesh_.write();
-
- Pout<< "Writing displacement field ..." << endl;
- disp.write();
- tmp<pointScalarField> magDisp(mag(disp));
- magDisp().write();
-
- Pout<< "Writing actual patch displacement ..." << endl;
- vectorField actualPatchDisp
- (
- IndirectList<point>(disp, pp.meshPoints())()
- );
- dumpMove
- (
- mesh_.time().path()/"actualPatchDisplacement.obj",
- pp.localPoints(),
- pp.localPoints() + actualPatchDisp
- );
- }
-}
-
-
-void Foam::autoHexMeshDriver::scaleMesh
-(
- const snapParameters& snapParams,
- const label nInitErrors,
- const List<labelPair>& baffles,
- motionSmoother& meshMover
-)
-{
- // Relax displacement until correct mesh
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- labelList checkFaces(identity(mesh_.nFaces()));
-
- scalar oldErrorReduction = -1;
-
- Info<< "Moving mesh ..." << endl;
- for (label iter = 0; iter < 2*snapParams.nSnap(); iter++)
- {
- Info<< nl << "Iteration " << iter << endl;
-
- if (iter == snapParams.nSnap())
- {
- Info<< "Displacement scaling for error reduction set to 0." << endl;
- oldErrorReduction = meshMover.setErrorReduction(0.0);
- }
-
- if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
- {
- Info<< "Successfully moved mesh" << endl;
-
- break;
- }
- if (debug_)
- {
- const_cast<Time&>(mesh_.time())++;
- Pout<< "Writing scaled mesh to time " << mesh_.time().timeName()
- << endl;
- mesh_.write();
-
- Pout<< "Writing displacement field ..." << endl;
- meshMover.displacement().write();
- tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
- magDisp().write();
- }
- }
-
- if (oldErrorReduction >= 0)
- {
- meshMover.setErrorReduction(oldErrorReduction);
- }
- Info<< "Moved mesh in = "
- << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
-}
-
-
-// ************************************************************************* //
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverTemplates.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverTemplates.C
deleted file mode 100644
index 31807f0dee77b8e43cd5d9e70684391adcf566b3..0000000000000000000000000000000000000000
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverTemplates.C
+++ /dev/null
@@ -1,75 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software; you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by the
- Free Software Foundation; either version 2 of the License, or (at your
- option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM; if not, write to the Free Software Foundation,
- Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "autoHexMeshDriver.H"
-#include "syncTools.H"
-
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
-template<class Type>
-void Foam::autoHexMeshDriver::averageNeighbours
-(
- const PackedList<1>& isMasterEdge,
- const labelList& meshEdges,
- const labelList& meshPoints,
- const edgeList& edges,
- const scalarField& invSumWeight,
- const Field<Type>& data,
- Field<Type>& average
-) const
-{
- average = pTraits<Type>::zero;
-
- forAll(edges, edgeI)
- {
- if (isMasterEdge.get(meshEdges[edgeI]) == 1)
- {
- const edge& e = edges[edgeI];
- //scalar eWeight = edgeWeights[edgeI];
- scalar eWeight = 1.0;
- label v0 = e[0];
- label v1 = e[1];
-
- average[v0] += eWeight*data[v1];
- average[v1] += eWeight*data[v0];
- }
- }
-
- syncTools::syncPointList
- (
- mesh_,
- meshPoints,
- average,
- plusEqOp<Type>(),
- pTraits<Type>::zero, // null value
- false // no separation
- );
-
- average *= invSumWeight;
-}
-
-
-// ************************************************************************* //
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
index d764c3fa6f2e996b42023eaff80023488ea79194..2db82c8eed8954d395b596580688cf3402c36362 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
@@ -361,7 +361,7 @@ public:
// Member Functions
- //- Helper funntion: count cells per processor in wanted distribution
+ //- Helper function: count cells per processor in wanted distribution
static labelList countCells(const labelList&);
//- Send cells to neighbours according to distribution