Commit f55da541 authored by franjo_j@hotmail.com's avatar franjo_j@hotmail.com
Browse files

2D cartesian templete working


git-svn-id: https://pl5.projectlocker.com/igui/meshGeneration/svn@42 fdcce57e-7e00-11e2-b579-49867b4cea03
parent 88a13ccc
......@@ -9,6 +9,7 @@ meshSurfaceCheckInvertedVertices = utilities/surfaceTools/meshSurfaceCheckInvert
meshSurfaceCheckEdgeTypes = utilities/surfaceTools/meshSurfaceCheckEdgeTypes
meshSurfaceCutter = utilities/surfaceTools/meshSurfaceCutter
meshSurfaceMapper = utilities/surfaceTools/meshSurfaceMapper
meshSurfaceMapper2D = utilities/surfaceTools/meshSurfaceMapper2D
edgeExtraction = utilities/surfaceTools/edgeExtraction
edgeExtractor = $(edgeExtraction)/edgeExtractor
meshSurfaceEdgeExtractor = utilities/surfaceTools/meshSurfaceEdgeExtractor
......@@ -82,9 +83,11 @@ triSurfaceDetectFeatureEdges = utilities/triSurfaceTools/triSurfaceDetectFeature
triSurfaceClassifyEdges = utilities/triSurfaceTools/triSurfaceClassifyEdges
triSurfacePatchManipulator = utilities/triSurfaceTools/triSurfacePatchManipulator
triSurfaceRemoveFacets = utilities/triSurfaceTools/triSurfaceRemoveFacets
triSurfaceExtrude2DEdges = utilities/triSurfaceTools/triSurfaceExtrude2DEdges
polyMeshGen = utilities/meshes/polyMeshGen
writePatch = utilities/meshes/polyMeshGen/writePatch
polyMeshGen2DEngine = utilities/meshes/polyMeshGen2DEngine
polyMeshGenModifier = utilities/meshes/polyMeshGenModifier
polyMeshGenAddressing = utilities/meshes/polyMeshGenAddressing
polyMeshGenChecks = utilities/meshes/polyMeshGenChecks
......@@ -148,6 +151,8 @@ $(polyMeshGen)/polyMeshGenPoints.C
$(polyMeshGen)/polyMeshGenFaces.C
$(polyMeshGen)/polyMeshGenCells.C
$(polyMeshGen2DEngine)/polyMeshGen2DEngine.C
$(writePatch)/writePatchBase.C
$(writePatch)/writePatch.C
$(writePatch)/writeProcessorPatch.C
......@@ -299,6 +304,10 @@ $(meshSurfaceMapper)/meshSurfaceMapperMapVertices.C
$(meshSurfaceMapper)/meshSurfaceMapperCornersAndEdges.C
$(meshSurfaceMapper)/meshSurfaceMapperPremapVertices.C
$(meshSurfaceMapper2D)/meshSurfaceMapper2D.C
$(meshSurfaceMapper2D)/meshSurfaceMapper2DMapVertices.C
$(meshSurfaceMapper2D)/meshSurfaceMapper2DPremapVertices.C
$(edgeExtractor)/edgeExtractor.C
$(meshSurfaceEdgeExtractor)/meshSurfaceEdgeExtractor.C
......@@ -439,6 +448,8 @@ $(triSurfacePatchManipulator)/triSurfacePatchManipulatorFunctions.C
$(triSurfaceRemoveFacets)/triSurfaceRemoveFacets.C
$(triSurfaceRemoveFacets)/triSurfaceRemoveFacetsFunctions.C
$(triSurfaceExtrude2DEdges)/triSurfaceExtrude2DEdges.C
$(cartesianMeshExtractor)/cartesianMeshExtractor.C
$(cartesianMeshExtractor)/cartesianMeshExtractorPointsAndAddressing.C
$(cartesianMeshExtractor)/cartesianMeshExtractorPolyMesh.C
......
......@@ -52,12 +52,12 @@ namespace Foam
void cartesianMeshExtractor::createPolyMesh()
{
Info << "Creating polyMesh from octree" << endl;
const meshOctree& octree = octreeCheck_.octree();
//- give labels to cubes which will be used as mesh cells
const List<direction>& cType = octreeCheck_.boxType();
labelList& leafCellLabel = *leafCellLabelPtr_;
label nCells(0);
forAll(cType, leafI)
......@@ -67,7 +67,7 @@ void cartesianMeshExtractor::createPolyMesh()
(octree.returnLeaf(leafI).procNo() != Pstream::myProcNo())
)
continue;
if( cType[leafI] & meshOctreeAddressing::MESHCELL )
{
leafCellLabel[leafI] = nCells++;
......@@ -78,32 +78,32 @@ void cartesianMeshExtractor::createPolyMesh()
polyMeshGenModifier meshModifier(mesh_);
faceListPMG& faces = meshModifier.facesAccess();
cellListPMG& cells = meshModifier.cellsAccess();
//- start creating octree mesh
cells.setSize(nCells);
List<direction> nFacesInCell(nCells, direction(0));
label nFaces(0);
const VRWGraph& octreeFaces = octreeCheck_.octreeFaces();
const labelListPMG& owner = octreeCheck_.octreeFaceOwner();
const labelListPMG& neighbour = octreeCheck_.octreeFaceNeighbour();
//- map storing box label and a direction for each processor face
//- The map stores data in the same order on both sides of processor
//- boundaries. This is a consequence of Morton ordering of
//- leaf boxes in the octree.
std::map<label, labelListPMG> procFaces;
forAll(octreeFaces, faceI)
{
const label own = owner[faceI];
const label nei = neighbour[faceI];
const label ownLabel = leafCellLabel[own];
label neiLabel(-1);
if( nei != -1 )
neiLabel = leafCellLabel[nei];
if( (ownLabel != -1) && (neiLabel != -1) )
{
++nFaces;
......@@ -114,11 +114,11 @@ void cartesianMeshExtractor::createPolyMesh()
{
++nFaces;
++nFacesInCell[ownLabel];
if( (nei != -1) && (cType[nei] & meshOctreeAddressing::MESHCELL) )
{
const label procNo = octree.returnLeaf(nei).procNo();
procFaces[procNo].append(faceI);
}
}
......@@ -126,41 +126,41 @@ void cartesianMeshExtractor::createPolyMesh()
{
++nFaces;
++nFacesInCell[neiLabel];
if( (own != -1) && (cType[own] & meshOctreeAddressing::MESHCELL) )
{
const label procNo = octree.returnLeaf(own).procNo();
procFaces[procNo].append(faceI);
}
}
}
//- case is a serial run
faces.setSize(nFaces);
forAll(cells, cI)
cells[cI].setSize(nFacesInCell[cI]);
nFacesInCell = 0;
//- calculate faces in processor patches
if( Pstream::parRun() )
{
PtrList<writeProcessorPatch>& procBoundaries =
meshModifier.procBoundariesAccess();
//- set the number of procBoundaries
procBoundaries.setSize(procFaces.size());
std::ostringstream ss;
ss << Pstream::myProcNo();
const word name("processor"+ss.str()+"to");
label nProcBoundaries(nFaces), patchI(0);
//- allocate memory for processor patches
std::map<label, labelListPMG>::const_iterator iter;
for(iter=procFaces.begin();iter!=procFaces.end();++iter)
{
const label procI = iter->first;
std::ostringstream ssNei;
ssNei << procI;
procBoundaries.set
......@@ -176,7 +176,7 @@ void cartesianMeshExtractor::createPolyMesh()
procI
)
);
nProcBoundaries -= iter->second.size();
++patchI;
}
......@@ -187,15 +187,15 @@ void cartesianMeshExtractor::createPolyMesh()
for(iter=procFaces.begin();iter!=procFaces.end();++iter)
{
procBoundaries[patchI].patchStart() = nProcBoundaries;
const labelListPMG& patchFaces = iter->second;
forAll(patchFaces, pfI)
{
const label fLabel = patchFaces[pfI];
const label own = owner[fLabel];
const label nei = neighbour[fLabel];
const label curCell = leafCellLabel[own];
label neiCell(-1);
if( nei != -1 )
......@@ -229,7 +229,7 @@ void cartesianMeshExtractor::createPolyMesh()
<< abort(FatalError);
}
}
if( procBoundaries[patchI].patchSize() !=
(nProcBoundaries - procBoundaries[patchI].patchStart())
)
......@@ -238,30 +238,30 @@ void cartesianMeshExtractor::createPolyMesh()
"cartesianMeshExtractor::createPolyMesh()"
) << "Invalid patch size!" << Pstream::myProcNo()
<< abort(FatalError);
++patchI;
}
}
nFaces = 0;
forAll(octreeFaces, faceI)
{
const label own = owner[faceI];
const label nei = neighbour[faceI];
const label ownLabel = leafCellLabel[own];
label neiLabel(-1);
if( nei != -1 )
neiLabel = leafCellLabel[nei];
if( (ownLabel != -1) && (neiLabel != -1) )
{
//- internal face
faces[nFaces].setSize(octreeFaces.sizeOfRow(faceI));
forAllRow(octreeFaces, faceI, pI)
faces[nFaces][pI] = octreeFaces(faceI, pI);
cells[ownLabel][nFacesInCell[ownLabel]++] = nFaces;
cells[neiLabel][nFacesInCell[neiLabel]++] = nFaces;
++nFaces;
......@@ -273,12 +273,12 @@ void cartesianMeshExtractor::createPolyMesh()
//- face at a parallel boundary
continue;
}
//- boundary face
faces[nFaces].setSize(octreeFaces.sizeOfRow(faceI));
forAllRow(octreeFaces, faceI, pI)
faces[nFaces][pI] = octreeFaces(faceI, pI);
cells[ownLabel][nFacesInCell[ownLabel]++] = nFaces;
++nFaces;
}
......@@ -289,24 +289,24 @@ void cartesianMeshExtractor::createPolyMesh()
//- face at a parallel boundary
continue;
}
//- boundary face
faces[nFaces].setSize(octreeFaces.sizeOfRow(faceI));
faces[nFaces][0] = octreeFaces(faceI, 0);
for(label pI=octreeFaces.sizeOfRow(faceI)-1;pI>0;--pI)
faces[nFaces][octreeFaces.sizeOfRow(faceI)-pI] =
octreeFaces(faceI, pI);
cells[neiLabel][nFacesInCell[neiLabel]++] = nFaces;
++nFaces;
}
}
# ifdef DEBUGMesh
label nProcBoundaries(0);
forAll(procBoundaries, patchI)
nProcBoundaries += procBoundaries[patchI].patchSize();
if( faces.size() != (nProcBoundaries + nFaces) )
{
Serr << "Number of faces " << faces.size() << endl;
......@@ -318,7 +318,7 @@ void cartesianMeshExtractor::createPolyMesh()
) << Pstream::myProcNo() << "This mesh is invalid!"
<< abort(FatalError);
}
vectorField closedness(cells.size(), vector::zero);
const labelList& owner = mesh_.owner();
const labelList& neighbour = mesh_.neighbour();
......@@ -354,7 +354,53 @@ void cartesianMeshExtractor::createPolyMesh()
# endif
meshModifier.reorderBoundaryFaces();
if( octree.isQuadtree() )
{
//- generate empty patches
//- search for faces with a dominant z coordinate and store them
//- into an empty patch
meshSurfaceEngine mse(mesh_);
const vectorField& fNormals = mse.faceNormals();
const faceList::subList& bFaces = mse.boundaryFaces();
const labelList& fOwner = mse.faceOwners();
wordList patchNames(2);
patchNames[0] = "defaultFaces";
patchNames[1] = "unusedFaces";
VRWGraph boundaryFaces;
labelListPMG newFaceOwner;
labelListPMG newFacePatch;
forAll(fNormals, bfI)
{
//- store the face and its owner
boundaryFaces.appendList(bFaces[bfI]);
newFaceOwner.append(fOwner[bfI]);
const vector& fNormal = fNormals[bfI];
if( Foam::mag(fNormal.z()) > Foam::mag(fNormal.x() + fNormal.y()) )
{
newFacePatch.append(1);
}
else
{
newFacePatch.append(0);
}
}
//- replace the boundary with faces in correct patches
meshModifier.replaceBoundary
(
patchNames,
boundaryFaces,
newFaceOwner,
newFacePatch
);
}
Info << "Finished creating polyMesh" << endl;
}
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2005-2007 Franjo Juretic
\\/ 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Description
\*---------------------------------------------------------------------------*/
#include "polyMeshGen2DEngine.H"
#include "polyMeshGenAddressing.H"
#include "demandDrivenData.H"
# ifdef USE_OMP
#include <omp.h>
# endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void polyMeshGen2DEngine::findActiveFacesAndPoints()
{
const faceListPMG& faces = mesh_.faces();
const pointFieldPMG& points = mesh_.points();
//- classify points
const scalar tZ = 0.05 * (bb_.max().z() - bb_.min().z());
activePoint_.setSize(points.size());
label counter(0);
# ifdef USE_OMP
# pragma omp parallel for schedule(dynamic, 50) reduction(+ : counter)
# endif
forAll(points, pointI)
{
if( Foam::mag(points[pointI].z() - bb_.min().z()) < tZ )
{
activePoint_[pointI] = true;
++counter;
}
else
{
activePoint_[pointI] = false;
}
}
if( 2 * counter != points.size() )
{
FatalErrorIn
(
"void polyMeshGen2DEngine::findActiveFacesAndPoints()"
) << "The number of points at smallest x coordinate is"
<< " not half of the total number of points."
<< " This is not a 2D mesh or is not aligned with the z axis"
<< exit(FatalError);
}
activePoints_.setSize(counter);
offsetPoints_.setSize(counter);
counter = 0;
forAll(activePoint_, pointI)
if( activePoint_[pointI] )
activePoints_[counter++] = pointI;
const VRWGraph& pointPoints = mesh_.addressingData().pointPoints();
# ifdef USE_OMP
# pragma omp parallel for schedule(dynamic, 50)
# endif
forAll(activePoints_, pI)
{
const label pointI = activePoints_[pI];
label nInactive(0), offsetPoint(-1);
forAllRow(pointPoints, pointI, ppI)
{
if( !activePoint_[pointPoints(pointI, ppI)] )
{
++nInactive;
offsetPoint = pointPoints(pointI, ppI);
}
}
if( nInactive == 1 )
{
offsetPoints_[pI] = offsetPoint;
}
else
{
FatalErrorIn
(
"void polyMeshGen2DEngine::findActiveFacesAndPoints()"
) << "This cannot be a 2D mesh" << exit(FatalError);
}
}
//- find active faces
activeFace_.setSize(faces.size());
counter = 0;
# ifdef USE_OMP
# pragma omp parallel for schedule(dynamic, 50) reduction(+ : counter)
# endif
forAll(faces, faceI)
{
const face& f = faces[faceI];
bool allActive(true);
forAll(f, pI)
{
if( !activePoint_[f[pI]] )
{
allActive = false;
break;
}
}
if( allActive )
{
activeFace_[faceI] = true;
++counter;
}
else
{
activeFace_[faceI] = false;
}
}
//- set active faces
activeFaces_.setSize(counter);
forAll(activeFace_, faceI)
{
if( activeFace_[faceI] )
activeFaces_[counter++] = faceI;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
polyMeshGen2DEngine::polyMeshGen2DEngine(polyMeshGen& mesh)
:
mesh_(mesh),
bb_(mesh.points()),
activeFace_(),
activeFaces_(),
activePoint_(),
activePoints_(),
offsetPoints_()
{}
polyMeshGen2DEngine::~polyMeshGen2DEngine()
{}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void polyMeshGen2DEngine::correctPoints()
{
pointFieldPMG& points = mesh_.points();
# ifdef USE_OMP
# pragma omp parallel for schedule(dynamic, 50)
# endif
forAll(activePoints_, apI)
{
point& p = points[activePoints_[apI]];
point& op = points[offsetPoints_[apI]];
op.x() = p.x();
op.y() = p.y();
p.z() = bb_.min().z();
op.z() = bb_.max().z();
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2005-2007 Franjo Juretic
\\/ 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Class
polyMeshGen2DEngine
Description
A simple engine which provides topological information of a 2D mesh
and allows for maintaining consistency
SourceFiles
polyMeshGen2DEngine.C
\*---------------------------------------------------------------------------*/
#ifndef polyMeshGen2DEngine_H
#define polyMeshGen2DEngine_H
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //