From 92a66ebc61dcb3fb3fca95909f84260a0111703c Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 2 Apr 2012 10:56:59 +0100
Subject: [PATCH] ENH: cvMeshbackgroundMesh: cleaned up building
---
.../cvMesh/cvMeshBackgroundMesh/Make/files | 3 +
.../cvMesh/cvMeshBackgroundMesh/Make/options | 31 +
.../cvMeshBackgroundMesh.C | 920 ++++++++++++++++++
3 files changed, 954 insertions(+)
create mode 100644 applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/Make/files
create mode 100644 applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/Make/options
create mode 100644 applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/cvMeshBackgroundMesh.C
diff --git a/applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/Make/files b/applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/Make/files
new file mode 100644
index 00000000000..58406b931af
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/Make/files
@@ -0,0 +1,3 @@
+cvMeshBackgroundMesh.C
+
+EXE = $(FOAM_APPBIN)/cvMeshBackgroundMesh
diff --git a/applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/Make/options
new file mode 100644
index 00000000000..0e2d0c9f0ac
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/Make/options
@@ -0,0 +1,31 @@
+/*
+EXE_DEBUG = -DFULLDEBUG -g -O0
+EXE_FROUNDING_MATH = -frounding-math
+EXE_NDEBUG = -DNDEBUG
+*/
+
+include $(GENERAL_RULES)/CGAL
+
+EXE_INC = \
+ ${EXE_FROUNDING_MATH} \
+ ${EXE_NDEBUG} \
+ ${CGAL_INC} \
+ -I../conformalVoronoiMesh/lnInclude \
+ -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
+ -I$(LIB_SRC)/edgeMesh/lnInclude \
+ -I$(LIB_SRC)/sampling/lnInclude \
+ -I$(LIB_SRC)/triSurface/lnInclude \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/dynamicMesh/lnInclude
+
+EXE_LIBS = \
+ $(CGAL_LIBS) \
+ -lconformalVoronoiMesh \
+ -ldecompositionMethods /* -L$(FOAM_LIBBIN)/dummy -lscotchDecomp */ \
+ -ledgeMesh \
+ -ltriSurface \
+ -lmeshTools \
+ -ldynamicMesh \
+ -lsampling \
+ -lfiniteVolume
diff --git a/applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/cvMeshBackgroundMesh.C b/applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/cvMeshBackgroundMesh.C
new file mode 100644
index 00000000000..c2642366a92
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/cvMeshBackgroundMesh/cvMeshBackgroundMesh.C
@@ -0,0 +1,920 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Application
+ cvMeshBackGroundMesh
+
+Description
+ Writes out background mesh as constructed by cvMesh and constructs
+ distanceSurface.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PatchTools.H"
+#include "argList.H"
+#include "Time.H"
+#include "triSurface.H"
+#include "searchableSurfaces.H"
+#include "conformationSurfaces.H"
+#include "cellSizeControlSurfaces.H"
+#include "backgroundMeshDecomposition.H"
+#include "cellShape.H"
+#include "cellModeller.H"
+#include "DynamicField.H"
+#include "isoSurfaceCell.H"
+#include "vtkSurfaceWriter.H"
+#include "syncTools.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Tolerance (as fraction of the bounding box). Needs to be fairly lax since
+// usually meshes get written with limited precision (6 digits)
+static const scalar defaultMergeTol = 1E-6;
+
+// Get merging distance when matching face centres
+scalar getMergeDistance
+(
+ const argList& args,
+ const Time& runTime,
+ const boundBox& bb
+)
+{
+ scalar mergeTol = defaultMergeTol;
+ args.optionReadIfPresent("mergeTol", mergeTol);
+
+ scalar writeTol =
+ Foam::pow(scalar(10.0), -scalar(IOstream::defaultPrecision()));
+
+ Info<< "Merge tolerance : " << mergeTol << nl
+ << "Write tolerance : " << writeTol << endl;
+
+ if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
+ {
+ FatalErrorIn("getMergeDistance")
+ << "Your current settings specify ASCII writing with "
+ << IOstream::defaultPrecision() << " digits precision." << endl
+ << "Your merging tolerance (" << mergeTol << ") is finer than this."
+ << endl
+ << "Please change your writeFormat to binary"
+ << " or increase the writePrecision" << endl
+ << "or adjust the merge tolerance (-mergeTol)."
+ << exit(FatalError);
+ }
+
+ scalar mergeDist = mergeTol * bb.mag();
+
+ Info<< "Overall meshes bounding box : " << bb << nl
+ << "Relative tolerance : " << mergeTol << nl
+ << "Absolute matching distance : " << mergeDist << nl
+ << endl;
+
+ return mergeDist;
+}
+
+
+void printMeshData(const polyMesh& mesh)
+{
+ // Collect all data on master
+
+ globalIndex globalCells(mesh.nCells());
+ labelListList patchNeiProcNo(Pstream::nProcs());
+ labelListList patchSize(Pstream::nProcs());
+ const labelList& pPatches = mesh.globalData().processorPatches();
+ patchNeiProcNo[Pstream::myProcNo()].setSize(pPatches.size());
+ patchSize[Pstream::myProcNo()].setSize(pPatches.size());
+ forAll(pPatches, i)
+ {
+ const processorPolyPatch& ppp = refCast<const processorPolyPatch>
+ (
+ mesh.boundaryMesh()[pPatches[i]]
+ );
+ patchNeiProcNo[Pstream::myProcNo()][i] = ppp.neighbProcNo();
+ patchSize[Pstream::myProcNo()][i] = ppp.size();
+ }
+ Pstream::gatherList(patchNeiProcNo);
+ Pstream::gatherList(patchSize);
+
+
+ // Print stats
+
+ globalIndex globalBoundaryFaces(mesh.nFaces()-mesh.nInternalFaces());
+
+ label maxProcCells = 0;
+ label totProcFaces = 0;
+ label maxProcPatches = 0;
+ label totProcPatches = 0;
+ label maxProcFaces = 0;
+
+ for (label procI = 0; procI < Pstream::nProcs(); procI++)
+ {
+ Info<< endl
+ << "Processor " << procI << nl
+ << " Number of cells = " << globalCells.localSize(procI)
+ << endl;
+
+ label nProcFaces = 0;
+
+ const labelList& nei = patchNeiProcNo[procI];
+
+ forAll(patchNeiProcNo[procI], i)
+ {
+ Info<< " Number of faces shared with processor "
+ << patchNeiProcNo[procI][i] << " = " << patchSize[procI][i]
+ << endl;
+
+ nProcFaces += patchSize[procI][i];
+ }
+
+ Info<< " Number of processor patches = " << nei.size() << nl
+ << " Number of processor faces = " << nProcFaces << nl
+ << " Number of boundary faces = "
+ << globalBoundaryFaces.localSize(procI) << endl;
+
+ maxProcCells = max(maxProcCells, globalCells.localSize(procI));
+ totProcFaces += nProcFaces;
+ totProcPatches += nei.size();
+ maxProcPatches = max(maxProcPatches, nei.size());
+ maxProcFaces = max(maxProcFaces, nProcFaces);
+ }
+
+ // Stats
+
+ scalar avgProcCells = scalar(globalCells.size())/Pstream::nProcs();
+ scalar avgProcPatches = scalar(totProcPatches)/Pstream::nProcs();
+ scalar avgProcFaces = scalar(totProcFaces)/Pstream::nProcs();
+
+ // In case of all faces on one processor. Just to avoid division by 0.
+ if (totProcPatches == 0)
+ {
+ avgProcPatches = 1;
+ }
+ if (totProcFaces == 0)
+ {
+ avgProcFaces = 1;
+ }
+
+ Info<< nl
+ << "Number of processor faces = " << totProcFaces/2 << nl
+ << "Max number of cells = " << maxProcCells
+ << " (" << 100.0*(maxProcCells-avgProcCells)/avgProcCells
+ << "% above average " << avgProcCells << ")" << nl
+ << "Max number of processor patches = " << maxProcPatches
+ << " (" << 100.0*(maxProcPatches-avgProcPatches)/avgProcPatches
+ << "% above average " << avgProcPatches << ")" << nl
+ << "Max number of faces between processors = " << maxProcFaces
+ << " (" << 100.0*(maxProcFaces-avgProcFaces)/avgProcFaces
+ << "% above average " << avgProcFaces << ")" << nl
+ << endl;
+}
+
+
+// Return cellID
+label cellLabel
+(
+ const Vector<label>& nCells,
+ const label i,
+ const label j,
+ const label k
+)
+{
+ return i*nCells[1]*nCells[2]+j*nCells[2]+k;
+}
+
+label vtxLabel
+(
+ const Vector<label>& nCells,
+ const label i,
+ const label j,
+ const label k
+)
+{
+ Vector<label> nPoints
+ (
+ nCells[0]+1,
+ nCells[1]+1,
+ nCells[2]+1
+ );
+ return i*nPoints[1]*nPoints[2]+j*nPoints[2]+k;
+}
+
+
+//// Generate hex mesh around surface.
+//void setRefinement
+//(
+// const vector& cellSize,
+// const Vector<label>& nCells
+// polyTopoChange& meshMod
+//)
+//{
+// const treeBoundBox& bb = geometryToConformTo.globalBounds();
+//
+// // Determine the number of cells in each direction.
+// vector nScalarCells = bb.span()/cellSizeControl.defaultCellSize();
+// Vector<label> nCells
+// (
+// label(nScalarCells[0]),
+// label(nScalarCells[1]),
+// label(nScalarCells[2])
+// );
+// vector cellSize(bb.span()/nCells);
+//
+// const label nPatches = 1;
+// polyTopoChange meshMod(nPatches);
+//
+// // Generate points
+// for (label i = 0; i <= nCells[0]; i++)
+// {
+// for (label j = 0; j <= nCells[1]; j++)
+// {
+// for (label k = 0; k <= nCells[2]; k++)
+// {
+// meshMod.addPoint
+// (
+// bb.min()+i*cellSize[0]+j*cellSize[1]+k*cellSize[2],
+// -1, // masterPointID,
+// -1, // zoneID,
+// true
+// );
+// }
+// }
+// }
+//
+// // Generate cells such that cell
+// for (label i = 0; i < nCells[0]; i++)
+// {
+// for (label j = 0; j < nCells[1]; j++)
+// {
+// for (label k = 0; k < nCells[2]; k++)
+// {
+// meshMod.addCell
+// (
+// -1, // masterPointID,
+// -1, // masterEdgeID,
+// -1, // masterFaceID,
+// -1, // masterCellID,
+// -1 //zoneID
+// );
+// }
+// }
+// }
+//
+// // Generate internal faces
+// face f(4);
+//
+// // Internal faces between i and i+1
+// for (label i = 0; i < nCells[0]-1; i++)
+// {
+// for (label j = 0; j < nCells[1]-1; j++)
+// {
+// for (label k = 0; k < nCells[2]-1; k++)
+// {
+// face[0] = vertexID(nCells, i, j, k);
+// face[1] = vertexID(nCells, i, j+1, k);
+// face[2] = vertexID(nCells, i, j+1, k+1);
+// face[3] = vertexID(nCells, i, j, k+1);
+//
+// meshMod.addFace
+// (
+// f,
+// cellID(nCells, i,j, k), //own
+// cellID(nCells, i+1,j, k+1), //nei
+// -1, //masterPointID,
+// -1, //masterEdgeID,
+// -1, //masterFaceID,
+// false, //flipFaceFlux,
+// -1, //patchID,
+// -1, //zoneID,
+// false //zoneFlip
+// );
+// }
+// }
+// }
+// // Internal faces between j and j+1
+// for (label i = 0; i < nCells[0]-1; i++)
+// {
+// for (label j = 0; j < nCells[1]-1; j++)
+// {
+// for (label k = 0; k < nCells[2]-1; k++)
+// {
+// face[0] = vertexID(nCells, i, j, k);
+// face[1] = vertexID(nCells, i, j, k+1);
+// face[2] = vertexID(nCells, i+1, j, k+1);
+// face[3] = vertexID(nCells, i+1, j, k);
+//
+// meshMod.addFace
+// (
+// f,
+// cellID(nCells, i,j, k), //own
+// cellID(nCells, i,j+1, k+1), //nei
+// -1, //masterPointID,
+// -1, //masterEdgeID,
+// -1, //masterFaceID,
+// false, //flipFaceFlux,
+// -1, //patchID,
+// -1, //zoneID,
+// false //zoneFlip
+// );
+// }
+// }
+// }
+// // Internal faces between k and k+1
+// for (label i = 0; i < nCells[0]-1; i++)
+// {
+// for (label j = 0; j < nCells[1]-1; j++)
+// {
+// for (label k = 0; k < nCells[2]-1; k++)
+// {
+// face[0] = vertexID(nCells, i, j, k);
+// face[1] = vertexID(nCells, i+1, j, k);
+// face[2] = vertexID(nCells, i+1, j+1, k);
+// face[3] = vertexID(nCells, i, j+1, k);
+//
+// meshMod.addFace
+// (
+// f,
+// cellID(nCells, i,j, k), //own
+// cellID(nCells, i,j, k+1), //nei
+// -1, //masterPointID,
+// -1, //masterEdgeID,
+// -1, //masterFaceID,
+// false, //flipFaceFlux,
+// -1, //patchID,
+// -1, //zoneID,
+// false //zoneFlip
+// );
+// }
+// }
+// }
+//}
+
+
+autoPtr<polyMesh> generateHexMesh
+(
+ const IOobject& io,
+ const point& origin,
+ const vector& cellSize,
+ const Vector<label>& nCells
+)
+{
+ pointField points;
+ if (nCells[0]+nCells[1]+nCells[2] > 0)
+ {
+ points.setSize((nCells[0]+1)*(nCells[1]+1)*(nCells[2]+1));
+
+ // Generate points
+ for (label i = 0; i <= nCells[0]; i++)
+ {
+ for (label j = 0; j <= nCells[1]; j++)
+ {
+ for (label k = 0; k <= nCells[2]; k++)
+ {
+ point pt = origin;
+ pt.x() += i*cellSize[0];
+ pt.y() += j*cellSize[1];
+ pt.z() += k*cellSize[2];
+ points[vtxLabel(nCells, i, j, k)] = pt;
+ }
+ }
+ }
+ }
+
+
+ const cellModel& hex = *(cellModeller::lookup("hex"));
+ cellShapeList cellShapes(nCells[0]*nCells[1]*nCells[2]);
+
+ labelList hexPoints(8);
+ label cellI = 0;
+ for (label i = 0; i < nCells[0]; i++)
+ {
+ for (label j = 0; j < nCells[1]; j++)
+ {
+ for (label k = 0; k < nCells[2]; k++)
+ {
+ hexPoints[0] = vtxLabel(nCells, i, j, k);
+ hexPoints[1] = vtxLabel(nCells, i+1, j, k);
+ hexPoints[2] = vtxLabel(nCells, i+1, j+1, k);
+ hexPoints[3] = vtxLabel(nCells, i, j+1, k);
+ hexPoints[4] = vtxLabel(nCells, i, j, k+1);
+ hexPoints[5] = vtxLabel(nCells, i+1, j, k+1);
+ hexPoints[6] = vtxLabel(nCells, i+1, j+1, k+1);
+ hexPoints[7] = vtxLabel(nCells, i, j+1, k+1);
+ cellShapes[cellI++] = cellShape(hex, hexPoints);
+ }
+ }
+ }
+
+ faceListList boundary(0);
+ wordList patchNames(0);
+ wordList patchTypes(0);
+ word defaultFacesName = "defaultFaces";
+ word defaultFacesType = polyPatch::typeName;
+ wordList patchPhysicalTypes(0);
+
+ return autoPtr<polyMesh>
+ (
+ new polyMesh
+ (
+ io,
+ xferMoveTo<pointField>(points),
+ cellShapes,
+ boundary,
+ patchNames,
+ patchTypes,
+ defaultFacesName,
+ defaultFacesType,
+ patchPhysicalTypes
+ )
+ );
+}
+
+
+// Determine for every point a signed distance to the nearest surface
+// (outside is positive)
+tmp<scalarField> signedDistance
+(
+ const scalarField& distSqr,
+ const pointField& points,
+ const searchableSurfaces& geometry,
+ const labelList& surfaces
+)
+{
+ tmp<scalarField> tfld(new scalarField(points.size(), Foam::sqr(GREAT)));
+ scalarField& fld = tfld();
+
+ // Find nearest
+ List<pointIndexHit> nearest;
+ labelList nearestSurfaces;
+ searchableSurfacesQueries::findNearest
+ (
+ geometry,
+ surfaces,
+ points,
+ scalarField(points.size(), Foam::sqr(GREAT)),//distSqr
+ nearestSurfaces,
+ nearest
+ );
+
+ // Determine sign of nearest. Sort by surface to do this.
+ DynamicField<point> surfPoints(points.size());
+ DynamicList<label> surfIndices(points.size());
+
+ forAll(surfaces, surfI)
+ {
+ // Extract points on this surface
+ surfPoints.clear();
+ surfIndices.clear();
+ forAll(nearestSurfaces, i)
+ {
+ if (nearestSurfaces[i] == surfI)
+ {
+ surfPoints.append(points[i]);
+ surfIndices.append(i);
+ }
+ }
+
+ // Calculate sideness of these surface points
+ label geomI = surfaces[surfI];
+ List<searchableSurface::volumeType> volType;
+ geometry[geomI].getVolumeType(surfPoints, volType);
+
+ // Push back to original
+ forAll(volType, i)
+ {
+ label pointI = surfIndices[i];
+ scalar dist = mag(points[pointI] - nearest[pointI].hitPoint());
+
+ searchableSurface::volumeType vT = volType[i];
+
+ if (vT == searchableSurface::OUTSIDE)
+ {
+ fld[pointI] = dist;
+ }
+ else if (vT == searchableSurface::INSIDE)
+ {
+ fld[i] = -dist;
+ }
+ else
+ {
+ FatalErrorIn("signedDistance()")
+ << "getVolumeType failure, neither INSIDE or OUTSIDE"
+ << exit(FatalError);
+ }
+ }
+ }
+
+ return tfld;
+}
+
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+ argList::addNote
+ (
+ "Generate cvMesh-consistent representation of surfaces"
+ );
+ argList::addBoolOption
+ (
+ "writeMesh",
+ "write the resulting mesh and distance fields"
+ );
+ argList::addOption
+ (
+ "mergeTol",
+ "scalar",
+ "specify the merge distance relative to the bounding box size "
+ "(default 1E-6)"
+ );
+
+ #include "setRootCase.H"
+ #include "createTime.H"
+ runTime.functionObjects().off();
+
+ const bool writeMesh = args.optionFound("writeMesh");
+
+ if (writeMesh)
+ {
+ Info<< "Writing resulting mesh and cellDistance, pointDistance fields."
+ << nl << endl;
+ }
+
+
+ IOdictionary cvMeshDict
+ (
+ IOobject
+ (
+ "cvMeshDict",
+ runTime.system(),
+ runTime,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE
+ )
+ );
+
+ // Define/load all geometry
+ searchableSurfaces allGeometry
+ (
+ IOobject
+ (
+ "cvSearchableSurfaces",
+ runTime.constant(),
+ "triSurface",
+ runTime,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ ),
+ cvMeshDict.subDict("geometry")
+ );
+
+ Random rndGen(64293*Pstream::myProcNo());
+
+ conformationSurfaces geometryToConformTo
+ (
+ runTime,
+ rndGen,
+ allGeometry,
+ cvMeshDict.subDict("surfaceConformation")
+ );
+
+ cellSizeControlSurfaces cellSizeControl
+ (
+ allGeometry,
+ cvMeshDict.subDict("motionControl")
+ );
+
+
+ // Generate starting block mesh
+ vector cellSize;
+ {
+ const treeBoundBox& bb = geometryToConformTo.globalBounds();
+
+ // Determine the number of cells in each direction.
+ const vector span = bb.span();
+ vector nScalarCells = span/cellSizeControl.defaultCellSize();
+
+ // Calculate initial cell size to be a little bit smaller than the
+ // defaultCellSize to avoid initial refinement triggering.
+ Vector<label> nCells = Vector<label>
+ (
+ label(nScalarCells.x())+2,
+ label(nScalarCells.y())+2,
+ label(nScalarCells.z())+2
+ );
+ cellSize = vector
+ (
+ span[0]/nCells[0],
+ span[1]/nCells[1],
+ span[2]/nCells[2]
+ );
+
+ Info<< "Generating initial hex mesh with" << nl
+ << " bounding box : " << bb << nl
+ << " nCells : " << nCells << nl
+ << " cellSize : " << cellSize << nl
+ << endl;
+
+ autoPtr<polyMesh> meshPtr
+ (
+ generateHexMesh
+ (
+ IOobject
+ (
+ polyMesh::defaultRegion,
+ runTime.constant(),
+ runTime
+ ),
+ bb.min(),
+ cellSize,
+ (
+ Pstream::master()
+ ? nCells
+ : Vector<label>(0, 0, 0)
+ )
+ )
+ );
+ Info<< "Writing initial hex mesh to " << meshPtr().instance() << nl
+ << endl;
+ meshPtr().write();
+ }
+
+ // Distribute the initial mesh
+ if (Pstream::parRun())
+ {
+# include "createMesh.H"
+ Info<< "Loaded mesh:" << endl;
+ printMeshData(mesh);
+
+ // Allocate a decomposer
+ IOdictionary decompositionDict
+ (
+ IOobject
+ (
+ "decomposeParDict",
+ runTime.system(),
+ mesh,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE
+ )
+ );
+
+ autoPtr<decompositionMethod> decomposer
+ (
+ decompositionMethod::New
+ (
+ decompositionDict
+ )
+ );
+
+ labelList decomp = decomposer().decompose(mesh, mesh.cellCentres());
+
+ // Global matching tolerance
+ const scalar tolDim = getMergeDistance
+ (
+ args,
+ runTime,
+ mesh.bounds()
+ );
+
+ // Mesh distribution engine
+ fvMeshDistribute distributor(mesh, tolDim);
+
+ Info<< "Wanted distribution:"
+ << distributor.countCells(decomp) << nl << endl;
+
+ // Do actual sending/receiving of mesh
+ autoPtr<mapDistributePolyMesh> map = distributor.distribute(decomp);
+
+ // Print some statistics
+ //Info<< "After distribution:" << endl;
+ //printMeshData(mesh);
+
+ mesh.setInstance(runTime.constant());
+ Info<< "Writing redistributed mesh" << nl << endl;
+ mesh.write();
+ }
+
+
+ Info<< "Refining backgroud mesh according to cell size specification" << nl
+ << endl;
+
+ backgroundMeshDecomposition backgroundMesh
+ (
+ 1.0, //spanScale,ratio of poly cell size v.s. hex cell size
+ 0.0, //minCellSizeLimit
+ 0, //minLevels
+ 4, //volRes, check multiple points per cell
+ 20.0, //maxCellWeightCoeff
+ runTime,
+ geometryToConformTo,
+ cellSizeControl,
+ rndGen,
+ cvMeshDict
+ );
+
+ if (writeMesh)
+ {
+ runTime++;
+ Info<< "Writing mesh to " << runTime.timeName() << endl;
+ backgroundMesh.mesh().write();
+ }
+
+ const scalar tolDim = getMergeDistance
+ (
+ args,
+ runTime,
+ backgroundMesh.mesh().bounds()
+ );
+
+
+ faceList isoFaces;
+ pointField isoPoints;
+
+ {
+ // Apply a distanceSurface to it.
+ const fvMesh& fvm = backgroundMesh.mesh();
+
+ volScalarField cellDistance
+ (
+ IOobject
+ (
+ "cellDistance",
+ fvm.time().timeName(),
+ fvm.time(),
+ IOobject::NO_READ,
+ IOobject::NO_WRITE,
+ false
+ ),
+ fvm,
+ dimensionedScalar("zero", dimLength, 0)
+ );
+
+ const searchableSurfaces& geometry = geometryToConformTo.geometry();
+ const labelList& surfaces = geometryToConformTo.surfaces();
+
+
+ // Get maximum search size per cell
+ scalarField distSqr(cellDistance.size());
+
+ const labelList& cellLevel = backgroundMesh.cellLevel();
+ forAll(cellLevel, cellI)
+ {
+ // The largest edge of the cell will always be less than the
+ // span of the bounding box of the cell.
+ distSqr[cellI] = magSqr(cellSize)/pow(2, cellLevel[cellI]);
+ }
+
+ {
+ // Internal field
+ cellDistance.internalField() = signedDistance
+ (
+ distSqr,
+ fvm.C(),
+ geometry,
+ surfaces
+ );
+ // Patch fields
+ forAll(fvm.C().boundaryField(), patchI)
+ {
+ const pointField& cc = fvm.C().boundaryField()[patchI];
+ fvPatchScalarField& fld = cellDistance.boundaryField()[patchI];
+ scalarField patchDistSqr
+ (
+ fld.patch().patchInternalField(distSqr)
+ );
+ fld = signedDistance(patchDistSqr, cc, geometry, surfaces);
+ }
+
+ // On processor patches the fvm.C() will already be the cell centre
+ // on the opposite side so no need to swap cellDistance.
+
+ if (writeMesh)
+ {
+ cellDistance.write();
+ }
+ }
+
+
+ // Distance to points
+ pointScalarField pointDistance
+ (
+ IOobject
+ (
+ "pointDistance",
+ fvm.time().timeName(),
+ fvm.time(),
+ IOobject::NO_READ,
+ IOobject::NO_WRITE,
+ false
+ ),
+ pointMesh::New(fvm),
+ dimensionedScalar("zero", dimLength, 0)
+ );
+ {
+ scalarField pointDistSqr(fvm.nPoints(), -sqr(GREAT));
+ for (label faceI = 0; faceI < fvm.nInternalFaces(); faceI++)
+ {
+ label own = fvm.faceOwner()[faceI];
+ label ownDistSqr = distSqr[own];
+
+ const face& f = fvm.faces()[faceI];
+ forAll(f, fp)
+ {
+ pointDistSqr[f[fp]] = max(pointDistSqr[f[fp]], ownDistSqr);
+ }
+ }
+ syncTools::syncPointList
+ (
+ fvm,
+ pointDistSqr,
+ maxEqOp<scalar>(),
+ -sqr(GREAT) // null value
+ );
+
+ pointDistance.internalField() = signedDistance
+ (
+ pointDistSqr,
+ fvm.points(),
+ geometry,
+ surfaces
+ );
+
+ if (writeMesh)
+ {
+ pointDistance.write();
+ }
+ }
+
+ isoSurfaceCell iso
+ (
+ fvm,
+ cellDistance,
+ pointDistance,
+ 0, //distance,
+ false //regularise
+ );
+
+ isoFaces.setSize(iso.size());
+ forAll(isoFaces, i)
+ {
+ isoFaces[i] = iso[i].triFaceFace();
+ }
+ isoPoints = iso.points();
+ }
+
+
+ pointField mergedPoints;
+ faceList mergedFaces;
+ labelList pointMergeMap;
+ PatchTools::gatherAndMerge
+ (
+ tolDim,
+ primitivePatch
+ (
+ SubList<face>(isoFaces, isoFaces.size()),
+ isoPoints
+ ),
+ mergedPoints,
+ mergedFaces,
+ pointMergeMap
+ );
+
+ vtkSurfaceWriter writer;
+ writer.write
+ (
+ runTime.path(),
+ "iso",
+ mergedPoints,
+ mergedFaces
+ );
+
+ Info<< "End\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
--
GitLab