Merge branch 'cv2d'

This commit is contained in:
laurence
2011-12-09 10:47:54 +00:00
73 changed files with 7008 additions and 512 deletions

View File

@ -4,12 +4,15 @@ set -x
wmake blockMesh wmake blockMesh
wmake all extrude wmake all extrude
wmake extrude2DMesh
extrude2DMesh/Allwmake
wmake snappyHexMesh wmake snappyHexMesh
if [ -d "$CGAL_ARCH_PATH" ] if [ -d "$CGAL_ARCH_PATH" ]
then then
cd cvMesh && ./Allwmake cvMesh/Allwmake
cv2DMesh/Allwmake
fi fi
# ----------------------------------------------------------------- end-of-file # ----------------------------------------------------------------- end-of-file

View File

@ -0,0 +1,8 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wclean libso conformalVoronoi2DMesh
wclean
# ----------------------------------------------------------------- end-of-file

View File

@ -0,0 +1,8 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wmake libso conformalVoronoi2DMesh
wmake
# ----------------------------------------------------------------- end-of-file

View File

@ -0,0 +1,82 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2011 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/>.
Typedefs
CGALTriangulation2Ddefs
Description
CGAL data structures used for 2D Delaunay meshing.
Define CGAL_INEXACT to use Exact_predicates_inexact_constructions kernel
otherwise the more robust but much less efficient
Exact_predicates_exact_constructions will be used.
Define CGAL_HIERARCHY to use hierarchical Delaunay triangulation which is
faster but uses more memory than the standard Delaunay triangulation.
\*---------------------------------------------------------------------------*/
#ifndef CGALTriangulation2Ddefs_H
#define CGALTriangulation2Ddefs_H
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "CGAL/Delaunay_triangulation_2.h"
#include "indexedVertex.H"
#include "indexedFace.H"
#ifdef CGAL_INEXACT
// Fast kernel using a double as the storage type but the triangulation
// may fail
#include "CGAL/Exact_predicates_inexact_constructions_kernel.h"
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
#else
// Very robust but expensive kernel
#include "CGAL/Exact_predicates_exact_constructions_kernel.h"
typedef CGAL::Exact_predicates_exact_constructions_kernel K;
#endif
typedef CGAL::indexedVertex<K> Vb;
typedef CGAL::indexedFace<K> Fb;
#ifdef CGAL_HIERARCHY
// Data structures for hierarchical Delaunay triangulation which is more
// efficient but also uses more storage
#include "CGAL/Triangulation_hierarchy_2.h"
typedef CGAL::Triangulation_hierarchy_vertex_base_2<Vb> Vbh;
typedef CGAL::Triangulation_data_structure_2<Vbh, Fb> Tds;
typedef CGAL::Delaunay_triangulation_2<K, Tds> DT;
typedef CGAL::Triangulation_hierarchy_2<DT> Delaunay;
#else
// Data structures for standard Delaunay triangulation
typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
typedef CGAL::Delaunay_triangulation_2<K, Tds> Delaunay;
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,471 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2011 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/>.
Class
CV2D
Description
Conformal-Voronoi 2D automatic mesher with grid or read initial points
and point position relaxation with optional "squarification".
There are a substantial number of options to this mesher read from
CV2DMesherDict file e.g.:
// Min cell size used in tolerances when inserting points for
// boundary conforming.
// Also used to as the grid spacing usind in insertGrid.
minCellSize 0.05;
// Feature angle used to inser feature points
// 0 = all features, 180 = no features
featureAngle 45;
// Maximum quadrant angle allowed at a concave corner before
// additional "mitering" lines are added
maxQuadAngle 110;
// Should the mesh be square-dominated or of unbiased hexagons
squares yes;
// Near-wall region where cells are aligned with the wall specified as a
// number of cell layers
nearWallAlignedDist 3;
// Chose if the cell orientation should relax during the iterations
// or remain fixed to the x-y directions
relaxOrientation no;
// Insert near-boundary point mirror or point-pairs
insertSurfaceNearestPointPairs yes;
// Mirror near-boundary points rather than insert point-pairs
mirrorPoints no;
// Insert point-pairs vor dual-cell vertices very near the surface
insertSurfaceNearPointPairs yes;
// Choose if to randomise the initial grid created by insertGrid.
randomiseInitialGrid yes;
// Perturbation fraction, 1 = cell-size.
randomPurturbation 0.1;
// Number of relaxation iterations.
nIterations 5;
// Relaxation factor at the start of the iteration sequence.
// 0.5 is a sensible maximum and < 0.2 converges better.
relaxationFactorStart 0.8;
// Relaxation factor at the end of the iteration sequence.
// Should be <= relaxationFactorStart
relaxationFactorEnd 0;
writeInitialTriangulation no;
writeFeatureTriangulation no;
writeNearestTriangulation no;
writeInsertedPointPairs no;
writeFinalTriangulation yes;
// Maximum number of iterations used in boundaryConform.
maxBoundaryConformingIter 5;
minEdgeLenCoeff 0.5;
maxNotchLenCoeff 0.3;
minNearPointDistCoeff 0.25;
ppDistCoeff 0.05;
SourceFiles
CGALTriangulation2Ddefs.H
indexedVertex.H
indexedFace.H
CV2DI.H
CV2D.C
CV2DIO.C
tolerances.C
controls.C
insertFeaturePoints.C
insertSurfaceNearestPointPairs.C
insertSurfaceNearPointPairs.C
insertBoundaryConformPointPairs.C
\*---------------------------------------------------------------------------*/
#ifndef CV2D_H
#define CV2D_H
#define CGAL_INEXACT
#define CGAL_HIERARCHY
#include "CGALTriangulation2Ddefs.H"
#include "Time.H"
#include "point2DFieldFwd.H"
#include "dictionary.H"
#include "Switch.H"
#include "PackedBoolList.H"
#include "EdgeMap.H"
#include "cv2DControls.H"
#include "tolerances.H"
#include "meshTools.H"
#include "triSurface.H"
#include "searchableSurfaces.H"
#include "conformationSurfaces.H"
#include "cellSizeControlSurfaces.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class CV2D Declaration
\*---------------------------------------------------------------------------*/
class CV2D
:
public Delaunay
{
private:
// Private data
//- The time registry of the application
const Time& runTime_;
mutable Random rndGen_;
//- The surface to mesh
//const querySurface& qSurf_;
//- All geometry of the meshing process, including surfaces to be
// conformed to and those to be used for refinement
searchableSurfaces allGeometry_;
conformationSurfaces qSurf_;
//- Meshing controls
cv2DControls controls_;
//- The cell size control object
cellSizeControlSurfaces cellSizeControl_;
//- z-level
scalar z_;
//- Keep track of the start of the internal points
label startOfInternalPoints_;
//- Keep track of the start of the surface point-pairs
label startOfSurfacePointPairs_;
//- Keep track of the boundary conform point-pairs
// stored after the insertion of the surface point-pairs in case
// the boundary conform function is called more than once without
// removing and insertin the surface point-pairs
label startOfBoundaryConformPointPairs_;
//- Store the feature points
std::list<Vb> featurePoints_;
//- Temporary storage for a dual-cell
static const label maxNvert = 20;
mutable point2D vertices[maxNvert+1];
mutable vector2D edges[maxNvert+1];
// Private Member Functions
//- Disallow default bitwise copy construct
CV2D(const CV2D&);
//- Disallow default bitwise assignment
void operator=(const CV2D&);
//- Insert point and return it's index
inline label insertPoint
(
const point2D& pt,
const label type
);
//- Insert point and return it's index
inline label insertPoint
(
const point2D& pt,
const label index,
const label type
);
inline label insertPoint
(
const Point& p,
const label index,
const label type
);
inline bool insertMirrorPoint
(
const point2D& nearSurfPt,
const point2D& surfPt
);
//- Insert a point-pair at a distance ppDist either side of
// surface point point surfPt in the direction n
inline void insertPointPair
(
const scalar mirrorDist,
const point2D& surfPt,
const vector2D& n
);
//- Create the initial mesh from the bounding-box
void insertBoundingBox();
//- Check if a point is within a line.
bool on2DLine(const point2D& p, const linePointRef& line);
//- Insert point groups at the feature points.
void insertFeaturePoints();
//- Re-insert point groups at the feature points.
void reinsertFeaturePoints();
//- Insert point-pairs at the given set of points using the surface
// normals corresponding to the given set of surface triangles
// and write the inserted point locations to the given file.
void insertPointPairs
(
const DynamicList<point2D>& nearSurfacePoints,
const DynamicList<point2D>& surfacePoints,
const DynamicList<label>& surfaceTris,
const DynamicList<label>& surfaceHits,
const fileName fName
);
//- Check to see if dual cell specified by given vertex iterator
// intersects the boundary and hence reqires a point-pair.
bool dualCellSurfaceIntersection
(
const Triangulation::Finite_vertices_iterator& vit
) const;
//- Insert point-pairs at the nearest points on the surface to the
// control vertex of dual-cells which intersect the boundary in order
// to provide a boundary-layer mesh.
// NB: This is not guaranteed to close the boundary
void insertSurfaceNearestPointPairs();
//- Insert point-pairs at small dual-cell edges on the surface in order
// to improve the boundary-layer mesh generated by
// insertSurfaceNearestPointPairs.
void insertSurfaceNearPointPairs();
//- Insert point-pair and correcting the Finite_vertices_iterator
// to account for the additional vertices
void insertPointPair
(
Triangulation::Finite_vertices_iterator& vit,
const point2D& p,
const label trii,
const label hitSurface
);
//- Insert point-pair at the best intersection point between the lines
// from the dual-cell real centroid and it's vertices and the surface.
bool insertPointPairAtIntersection
(
Triangulation::Finite_vertices_iterator& vit,
const point2D& defVert,
const point2D vertices[],
const scalar maxProtSize
);
//- Insert point-pairs corresponding to dual-cells which intersect
// the boundary surface
label insertBoundaryConformPointPairs(const fileName& fName);
void markNearBoundaryPoints();
//- Restore the Delaunay contraint
void fast_restore_Delaunay(Vertex_handle vh);
// Flip operations used by fast_restore_Delaunay
void external_flip(Face_handle& f, int i);
bool internal_flip(Face_handle& f, int i);
//- Write all the faces and all the triangles at a particular stage.
void write(const word& stage) const;
public:
//- Runtime type information
ClassName("CV2D");
// Constructors
//- Construct for given surface
CV2D
(
const Time& runTime,
const dictionary& controlDict
);
//- Destructor
~CV2D();
// Member Functions
// Access
inline const cv2DControls& meshControls() const;
// Conversion functions between point2D, point and Point
inline const point2D& toPoint2D(const point&) const;
inline const point2DField toPoint2D(const pointField&) const;
inline point toPoint3D(const point2D&) const;
#ifdef CGAL_INEXACT
typedef const point2D& point2DFromPoint;
typedef const Point& PointFromPoint2D;
#else
typedef point2D point2DFromPoint;
typedef Point PointFromPoint2D;
#endif
inline point2DFromPoint toPoint2D(const Point&) const;
inline PointFromPoint2D toPoint(const point2D&) const;
inline point toPoint3D(const Point&) const;
// Point insertion
//- Create the initial mesh from the given internal points.
// Points must be inside the boundary by at least nearness
// otherwise they are ignored.
void insertPoints
(
const point2DField& points,
const scalar nearness
);
//- Create the initial mesh from the internal points in the given
// file. Points outside the geometry are ignored.
void insertPoints(const fileName& pointFileName);
//- Create the initial mesh as a regular grid of points.
// Points outside the geometry are ignored.
void insertGrid();
//- Insert all surface point-pairs from
// insertSurfaceNearestPointPairs and
// findIntersectionForOutsideCentroid
void insertSurfacePointPairs();
//- Insert point-pairs where there are protrusions into
// or out of the surface
void boundaryConform();
// Point removal
//- Remove the point-pairs introduced by insertSurfacePointPairs
// and boundaryConform
void removeSurfacePointPairs();
// Point motion
inline void movePoint(const Vertex_handle& vh, const Point& P);
//- Move the internal points to the given new locations and update
// the triangulation to ensure it is Delaunay
// void moveInternalPoints(const point2DField& newPoints);
//- Calculate the displacements to create the new points
void newPoints(const scalar relaxation);
//- Extract patch names and sizes.
void extractPatches
(
wordList& patchNames,
labelList& patchSizes,
EdgeMap<label>& mapEdgesRegion
) const;
// Write
//- Write internal points to .obj file
void writePoints(const fileName& fName, bool internalOnly) const;
//- Write triangles as .obj file
void writeTriangles(const fileName& fName, bool internalOnly) const;
//- Write dual faces as .obj file
void writeFaces(const fileName& fName, bool internalOnly) const;
//- Calculates dual points (circumcentres of tets) and faces
// (point-cell walk of tets).
// Returns:
// - dualPoints (in triangle ordering)
// - dualFaces (compacted)
void calcDual
(
point2DField& dualPoints,
faceList& dualFaces,
wordList& patchNames,
labelList& patchSizes,
EdgeMap<label>& mapEdgesRegion
) const;
//- Write patch
void writePatch(const fileName& fName) const;
void write() const;
};
inline bool boundaryTriangle(const CV2D::Face_handle fc);
inline bool outsideTriangle(const CV2D::Face_handle fc);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "CV2DI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,221 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2011 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/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
inline Foam::label Foam::CV2D::insertPoint
(
const point2D& p,
const label type
)
{
uint nVert = number_of_vertices();
return insertPoint(toPoint(p), nVert, type);
}
inline Foam::label Foam::CV2D::insertPoint
(
const point2D& p,
const label index,
const label type
)
{
return insertPoint(toPoint(p), index, type);
}
inline Foam::label Foam::CV2D::insertPoint
(
const Point& p,
const label index,
const label type
)
{
uint nVert = number_of_vertices();
Vertex_handle vh = insert(p);
if (nVert == number_of_vertices())
{
WarningIn("Foam::CV2D::insertPoint")
<< "Failed to insert point " << toPoint2D(p) << endl;
}
else
{
vh->index() = index;
vh->type() = type;
}
return vh->index();
}
inline bool Foam::CV2D::insertMirrorPoint
(
const point2D& nearSurfPt,
const point2D& surfPt
)
{
point2D mirrorPoint(2*surfPt - nearSurfPt);
if (qSurf_.outside(toPoint3D(mirrorPoint)))
{
insertPoint(mirrorPoint, Vb::MIRROR_POINT);
return true;
}
else
{
return false;
}
}
inline void Foam::CV2D::insertPointPair
(
const scalar ppDist,
const point2D& surfPt,
const vector2D& n
)
{
vector2D ppDistn = ppDist*n;
label master = insertPoint
(
surfPt - ppDistn,
number_of_vertices() + 1
);
insertPoint(surfPt + ppDistn, master);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const Foam::cv2DControls& Foam::CV2D::meshControls() const
{
return controls_;
}
inline const Foam::point2D& Foam::CV2D::toPoint2D(const point& p) const
{
return reinterpret_cast<const point2D&>(p);
}
inline const Foam::point2DField Foam::CV2D::toPoint2D(const pointField& p) const
{
point2DField temp(p.size());
forAll(temp, pointI)
{
temp[pointI] = point2D(p[pointI].x(), p[pointI].y());
}
return temp;
}
inline Foam::point Foam::CV2D::toPoint3D(const point2D& p) const
{
return point(p.x(), p.y(), z_);
}
#ifdef CGAL_INEXACT
inline Foam::CV2D::point2DFromPoint Foam::CV2D::toPoint2D(const Point& P) const
{
return reinterpret_cast<point2DFromPoint>(P);
}
inline Foam::CV2D::PointFromPoint2D Foam::CV2D::toPoint(const point2D& p) const
{
return reinterpret_cast<PointFromPoint2D>(p);
}
#else
inline Foam::CV2D::point2DFromPoint Foam::CV2D::toPoint2D(const Point& P) const
{
return point2D(CGAL::to_double(P.x()), CGAL::to_double(P.y()));
}
inline Foam::CV2D::PointFromPoint2D Foam::CV2D::toPoint(const point2D& p) const
{
return Point(p.x(), p.y());
}
#endif
inline Foam::point Foam::CV2D::toPoint3D(const Point& P) const
{
return point(CGAL::to_double(P.x()), CGAL::to_double(P.y()), z_);
}
inline void Foam::CV2D::movePoint(const Vertex_handle& vh, const Point& P)
{
int i = vh->index();
int t = vh->type();
remove(vh);
Vertex_handle newVh = insert(P);
newVh->index() = i;
newVh->type() = t;
// label i = vh->index();
// move(vh, P);
// vh->index() = i;
//vh->set_point(P);
//fast_restore_Delaunay(vh);
}
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
inline bool Foam::boundaryTriangle(const CV2D::Face_handle fc)
{
return boundaryTriangle
(
*fc->vertex(0),
*fc->vertex(1),
*fc->vertex(2)
);
}
inline bool Foam::outsideTriangle(const CV2D::Face_handle fc)
{
return outsideTriangle
(
*fc->vertex(0),
*fc->vertex(1),
*fc->vertex(2)
);
}
// ************************************************************************* //

View File

@ -0,0 +1,304 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2011 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/>.
\*---------------------------------------------------------------------------*/
#include "CV2D.H"
#include "OFstream.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::CV2D::writePoints(const fileName& fName, bool internalOnly) const
{
Info<< "Writing points to " << fName << nl << endl;
OFstream str(fName);
for
(
Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
++vit
)
{
if (!internalOnly || vit->internalOrBoundaryPoint())
{
meshTools::writeOBJ(str, toPoint3D(vit->point()));
}
}
}
void Foam::CV2D::writeTriangles(const fileName& fName, bool internalOnly) const
{
Info<< "Writing triangles to " << fName << nl << endl;
OFstream str(fName);
labelList vertexMap(number_of_vertices(), -2);
label verti = 0;
for
(
Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
++vit
)
{
if (!internalOnly || !vit->farPoint())
{
vertexMap[vit->index()] = verti++;
meshTools::writeOBJ(str, toPoint3D(vit->point()));
}
}
for
(
Triangulation::Finite_faces_iterator fit = finite_faces_begin();
fit != finite_faces_end();
++fit
)
{
if
(
!internalOnly
|| (
fit->vertex(0)->internalOrBoundaryPoint()
|| fit->vertex(1)->internalOrBoundaryPoint()
|| fit->vertex(2)->internalOrBoundaryPoint()
)
)
{
str << "f";
for (label i = 0; i < 3; ++i)
{
str << " " << vertexMap[fit->vertex(i)->index()] + 1;
}
str << nl;
}
}
}
void Foam::CV2D::writeFaces(const fileName& fName, bool internalOnly) const
{
Info<< "Writing dual faces to " << fName << nl << endl;
OFstream str(fName);
label dualVerti = 0;
for
(
Triangulation::Finite_faces_iterator fit = finite_faces_begin();
fit != finite_faces_end();
++fit
)
{
if
(
!internalOnly
|| (
fit->vertex(0)->internalOrBoundaryPoint()
|| fit->vertex(1)->internalOrBoundaryPoint()
|| fit->vertex(2)->internalOrBoundaryPoint()
)
)
{
fit->faceIndex() = dualVerti++;
meshTools::writeOBJ(str, toPoint3D(circumcenter(fit)));
}
else
{
fit->faceIndex() = -1;
}
}
for
(
Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
++vit
)
{
if (!internalOnly || vit->internalOrBoundaryPoint())
{
Face_circulator fcStart = incident_faces(vit);
Face_circulator fc = fcStart;
str<< 'f';
do
{
if (!is_infinite(fc))
{
if (fc->faceIndex() < 0)
{
FatalErrorIn
(
"Foam::CV2D::writeFaces"
"(const fileName& fName, bool internalOnly)"
)<< "Dual face uses vertex defined by a triangle"
" defined by an external point"
<< exit(FatalError);
}
str<< ' ' << fc->faceIndex() + 1;
}
} while (++fc != fcStart);
str<< nl;
}
}
}
void Foam::CV2D::calcDual
(
point2DField& dualPoints,
faceList& dualFaces,
wordList& patchNames,
labelList& patchSizes,
EdgeMap<label>& mapEdgesRegion
) const
{
// Dual points stored in triangle order.
dualPoints.setSize(number_of_faces());
label dualVerti = 0;
for
(
Triangulation::Finite_faces_iterator fit = finite_faces_begin();
fit != finite_faces_end();
++fit
)
{
if
(
fit->vertex(0)->internalOrBoundaryPoint()
|| fit->vertex(1)->internalOrBoundaryPoint()
|| fit->vertex(2)->internalOrBoundaryPoint()
)
{
fit->faceIndex() = dualVerti;
dualPoints[dualVerti++] = toPoint2D(circumcenter(fit));
}
else
{
fit->faceIndex() = -1;
}
}
dualPoints.setSize(dualVerti);
extractPatches(patchNames, patchSizes, mapEdgesRegion);
forAll(patchNames, patchI)
{
Info<< "Patch " << patchNames[patchI]
<< " has size " << patchSizes[patchI] << endl;
}
// Create dual faces
// ~~~~~~~~~~~~~~~~~
dualFaces.setSize(number_of_vertices());
label dualFacei = 0;
labelList faceVerts(maxNvert);
for
(
Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
++vit
)
{
if (vit->internalOrBoundaryPoint())
{
Face_circulator fcStart = incident_faces(vit);
Face_circulator fc = fcStart;
label verti = 0;
do
{
if (!is_infinite(fc))
{
if (fc->faceIndex() < 0)
{
FatalErrorIn
(
"Foam::CV2D::calcDual"
"(point2DField& dualPoints, faceList& dualFaces)"
)<< "Dual face uses vertex defined by a triangle"
" defined by an external point"
<< exit(FatalError);
}
// Look up the index of the triangle
faceVerts[verti++] = fc->faceIndex();
}
} while (++fc != fcStart);
if (faceVerts.size() > 2)
{
dualFaces[dualFacei++] =
face(labelList::subList(faceVerts, verti));
}
else
{
Info<< "From triangle point:" << vit->index()
<< " coord:" << toPoint2D(vit->point())
<< " generated illegal dualFace:" << faceVerts
<< endl;
}
}
}
dualFaces.setSize(dualFacei);
}
void Foam::CV2D::writePatch(const fileName& fName) const
{
point2DField dual2DPoints;
faceList dualFaces;
wordList patchNames;
labelList patchSizes;
EdgeMap<label> mapEdgesRegion;
calcDual(dual2DPoints, dualFaces, patchNames, patchSizes, mapEdgesRegion);
pointField dualPoints(dual2DPoints.size());
forAll(dualPoints, ip)
{
dualPoints[ip] = toPoint3D(dual2DPoints[ip]);
}
// Dump as primitive patch to be read by extrudeMesh.
OFstream str(fName);
Info<< "Writing patch to be used with extrudeMesh to file " << fName
<< endl;
str << dualPoints << nl << dualFaces << nl;
}
// ************************************************************************* //

View File

@ -0,0 +1,13 @@
#include CGAL_FILES
CV2D.C
insertFeaturePoints.C
insertSurfaceNearestPointPairs.C
insertSurfaceNearPointPairs.C
insertBoundaryConformPointPairs.C
CV2DIO.C
shortEdgeFilter2D.C
cv2DMesh.C
EXE = $(FOAM_APPBIN)/cv2DMesh

View File

@ -0,0 +1,34 @@
EXE_DEBUG = -DFULLDEBUG -g -O0
EXE_FROUNDING_MATH = -frounding-math
EXE_NDEBUG = -DNDEBUG
include $(GENERAL_RULES)/CGAL
FFLAGS = -DCGAL_FILES='"${CGAL_ARCH_PATH}/share/files"'
EXE_INC = \
${EXE_FROUNDING_MATH} \
${EXE_NDEBUG} \
${CGAL_INC} \
-I$(FOAM_APP)/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/lnInclude \
-I$(FOAM_APP)/utilities/mesh/generation/extrude/extrudeModel/lnInclude \
-IconformalVoronoi2DMesh/lnInclude \
-I$(FOAM_APP)/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/edgeMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude
EXE_LIBS = \
$(CGAL_LIBS) \
-lextrude2DMesh \
-lextrudeModel \
-lcv2DMesh \
-lconformalVoronoiMesh \
-lmeshTools \
-lsurfMesh \
-ledgeMesh \
-ltriSurface \
-ldynamicMesh \
-ldecompositionMethods

View File

@ -0,0 +1,4 @@
cv2DControls/cv2DControls.C
LIB = $(FOAM_LIBBIN)/libcv2DMesh

View File

@ -0,0 +1,3 @@
EXE_INC =
LIB_LIBS =

View File

@ -0,0 +1,149 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2011 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/>.
\*----------------------------------------------------------------------------*/
#include "cv2DControls.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::cv2DControls::cv2DControls
(
const dictionary& controlDict,
const boundBox& bb
)
:
dict_(controlDict),
motionControl_(controlDict.subDict("motionControl")),
conformationControl_(controlDict.subDict("surfaceConformation")),
minCellSize_(readScalar(motionControl_.lookup("minCellSize"))),
minCellSize2_(Foam::sqr(minCellSize_)),
maxQuadAngle_(readScalar(conformationControl_.lookup("maxQuadAngle"))),
nearWallAlignedDist_
(
readScalar(motionControl_.lookup("nearWallAlignedDist"))*minCellSize_
),
nearWallAlignedDist2_(Foam::sqr(nearWallAlignedDist_)),
insertSurfaceNearestPointPairs_
(
conformationControl_.lookup("insertSurfaceNearestPointPairs")
),
mirrorPoints_(conformationControl_.lookup("mirrorPoints")),
insertSurfaceNearPointPairs_
(
conformationControl_.lookup("insertSurfaceNearPointPairs")
),
objOutput_(motionControl_.lookupOrDefault<Switch>("objOutput", false)),
randomiseInitialGrid_(conformationControl_.lookup("randomiseInitialGrid")),
randomPerturbation_
(
readScalar(conformationControl_.lookup("randomPerturbation"))
),
maxBoundaryConformingIter_
(
readLabel(conformationControl_.lookup("maxBoundaryConformingIter"))
),
span_
(
max(mag(bb.max().x()), mag(bb.min().x()))
+ max(mag(bb.max().y()), mag(bb.min().y()))
),
span2_(Foam::sqr(span_)),
minEdgeLen_
(
readScalar(conformationControl_.lookup("minEdgeLenCoeff"))
*minCellSize_
),
minEdgeLen2_(Foam::sqr(minEdgeLen_)),
maxNotchLen_
(
readScalar(conformationControl_.lookup("maxNotchLenCoeff"))
*minCellSize_
),
maxNotchLen2_(Foam::sqr(maxNotchLen_)),
minNearPointDist_
(
readScalar(conformationControl_.lookup("minNearPointDistCoeff"))
*minCellSize_
),
minNearPointDist2_(Foam::sqr(minNearPointDist_)),
ppDist_
(
readScalar(conformationControl_.lookup("pointPairDistanceCoeff"))
*minCellSize_
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::cv2DControls::~cv2DControls()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::cv2DControls::write(Ostream& os) const
{
os.indentLevel() = 1;
os.precision(2);
os.flags(ios_base::scientific);
os << nl << "Outputting CV2D Mesher controls:" << nl
<< token::BEGIN_BLOCK << nl
<< indent << "minCellSize2_ : " << minCellSize2_ << nl
<< indent << "span_ / span2_ : " << span_ << " / " << span2_ << nl
<< indent << "maxNotchLen2_ : " << maxNotchLen2_ << nl
<< indent << "minNearPointDist2_ : " << minNearPointDist2_ << nl
<< indent << "nearWallAlignedDist2_ : " << nearWallAlignedDist2_ << nl
<< indent << "ppDist_ : " << ppDist_ << nl
<< indent << "minEdgeLen2_ : " << minEdgeLen2_ << nl
<< token::END_BLOCK << endl;
}
// * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<(Ostream& os, const cv2DControls& s)
{
s.write(os);
return os;
}
// ************************************************************************* //

View File

@ -0,0 +1,257 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 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/>.
Class
Foam::cv2DControls
Description
Controls for the 2D CV mesh generator.
SourceFiles
cv2DControls.C
cv2DControlsI.H
\*---------------------------------------------------------------------------*/
#ifndef cv2DControls_H
#define cv2DControls_H
#include "Switch.H"
#include "dictionary.H"
#include "boundBox.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class cv2DControls Declaration
\*---------------------------------------------------------------------------*/
class cv2DControls
{
// Private data
//- Description of data_
const dictionary& dict_;
const dictionary& motionControl_;
const dictionary& conformationControl_;
// Private Member Functions
//- Disallow default bitwise copy construct
cv2DControls(const cv2DControls&);
//- Disallow default bitwise assignment
void operator=(const cv2DControls&);
public:
// Controls
//- Minimum cell size below which protusions through the surface are
// not split
scalar minCellSize_;
//- Square of minCellSize
scalar minCellSize2_;
//- Maximum quadrant angle allowed at a concave corner before
// additional "mitering" lines are added
scalar maxQuadAngle_;
//- Near-wall region where cells are aligned with the wall
scalar nearWallAlignedDist_;
//- Square of nearWallAlignedDist
scalar nearWallAlignedDist2_;
//- Insert near-boundary point mirror or point-pairs
Switch insertSurfaceNearestPointPairs_;
//- Mirror near-boundary points rather than insert point-pairs
Switch mirrorPoints_;
//- Insert point-pairs vor dual-cell vertices very near the surface
Switch insertSurfaceNearPointPairs_;
Switch objOutput_;
Switch randomiseInitialGrid_;
scalar randomPerturbation_;
label maxBoundaryConformingIter_;
// Tolerances
//- Maximum cartesian span of the geometry
scalar span_;
//- Square of span
scalar span2_;
//- Minumum edge-length of the cell size below which protusions
// through the surface are not split
scalar minEdgeLen_;
//- Square of minEdgeLen
scalar minEdgeLen2_;
//- Maximum notch size below which protusions into the surface are
// not filled
scalar maxNotchLen_;
//- Square of maxNotchLen
scalar maxNotchLen2_;
//- The minimum distance alowed between a dual-cell vertex
// and the surface before a point-pair is introduced
scalar minNearPointDist_;
//- Square of minNearPoint
scalar minNearPointDist2_;
//- Distance between boundary conforming point-pairs
scalar ppDist_;
//- Square of ppDist
scalar ppDist2_;
// Constructors
cv2DControls
(
const dictionary& controlDict,
const boundBox& bb
);
//- Destructor
~cv2DControls();
// Member Functions
// Access
//- Return the minimum cell size
inline scalar minCellSize() const;
//- Return the square of the minimum cell size
inline scalar minCellSize2() const;
//- Return the maximum quadrant angle
inline scalar maxQuadAngle() const;
//- Return number of layers to align with the nearest wall
inline scalar nearWallAlignedDist() const;
//- Return square of nearWallAlignedDist
inline scalar nearWallAlignedDist2() const;
//- Return insertSurfaceNearestPointPairs Switch
inline Switch insertSurfaceNearestPointPairs() const;
//- Return mirrorPoints Switch
inline Switch mirrorPoints() const;
//- Return insertSurfaceNearPointPairs Switch
inline Switch insertSurfaceNearPointPairs() const;
//- Return the objOutput Switch
inline Switch objOutput() const;
//- Return the randomise initial point layout Switch
inline Switch randomiseInitialGrid() const;
//- Return the random perturbation factor
inline scalar randomPerturbation() const;
//- Return the maximum number of boundary conformation iterations
inline label maxBoundaryConformingIter() const;
//- Return the span
inline scalar span() const;
//- Return the span squared
inline scalar span2() const;
//- Return the minEdgeLen
inline scalar minEdgeLen() const;
//- Return the minEdgeLen squared
inline scalar minEdgeLen2() const;
//- Return the maxNotchLen
inline scalar maxNotchLen() const;
//- Return the maxNotchLen squared
inline scalar maxNotchLen2() const;
//- Return the minNearPointDist
inline scalar minNearPointDist() const;
//- Return the minNearPointDist squared
inline scalar minNearPointDist2() const;
//- Return the ppDist
inline scalar ppDist() const;
// Write
//- Write controls to output stream.
void write(Ostream& os) const;
//- Ostream Operator
friend Ostream& operator<<
(
Ostream& os,
const cv2DControls& s
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "cv2DControlsI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,153 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
\*---------------------------------------------------------------------------*/
inline Foam::scalar Foam::cv2DControls::minCellSize() const
{
return minCellSize_;
}
inline Foam::scalar Foam::cv2DControls::minCellSize2() const
{
return minCellSize2_;
}
inline Foam::scalar Foam::cv2DControls::maxQuadAngle() const
{
return maxQuadAngle_;
}
inline Foam::scalar Foam::cv2DControls::nearWallAlignedDist() const
{
return nearWallAlignedDist_;
}
inline Foam::scalar Foam::cv2DControls::nearWallAlignedDist2() const
{
return nearWallAlignedDist2_;
}
inline Foam::Switch Foam::cv2DControls::insertSurfaceNearestPointPairs() const
{
return insertSurfaceNearestPointPairs_;
}
inline Foam::Switch Foam::cv2DControls::mirrorPoints() const
{
return mirrorPoints_;
}
inline Foam::Switch Foam::cv2DControls::insertSurfaceNearPointPairs() const
{
return insertSurfaceNearPointPairs_;
}
inline Foam::Switch Foam::cv2DControls::objOutput() const
{
return objOutput_;
}
inline Foam::Switch Foam::cv2DControls::randomiseInitialGrid() const
{
return randomiseInitialGrid_;
}
inline Foam::scalar Foam::cv2DControls::randomPerturbation() const
{
return randomPerturbation_;
}
inline Foam::label Foam::cv2DControls::maxBoundaryConformingIter() const
{
return maxBoundaryConformingIter_;
}
inline Foam::scalar Foam::cv2DControls::span() const
{
return span_;
}
inline Foam::scalar Foam::cv2DControls::span2() const
{
return span2_;
}
inline Foam::scalar Foam::cv2DControls::minEdgeLen() const
{
return minEdgeLen_;
}
inline Foam::scalar Foam::cv2DControls::minEdgeLen2() const
{
return minEdgeLen2_;
}
inline Foam::scalar Foam::cv2DControls::maxNotchLen() const
{
return maxNotchLen_;
}
inline Foam::scalar Foam::cv2DControls::maxNotchLen2() const
{
return maxNotchLen2_;
}
inline Foam::scalar Foam::cv2DControls::minNearPointDist() const
{
return minNearPointDist_;
}
inline Foam::scalar Foam::cv2DControls::minNearPointDist2() const
{
return minNearPointDist2_;
}
inline Foam::scalar Foam::cv2DControls::ppDist() const
{
return ppDist_;
}
// ************************************************************************* //

View File

@ -0,0 +1,251 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2011 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
cv2DMesh
Description
Conformal-Voronoi 2D extruding automatic mesher with grid or read
initial points and point position relaxation with optional
"squarification".
\*---------------------------------------------------------------------------*/
#include "CV2D.H"
#include "argList.H"
#include "MeshedSurfaces.H"
#include "shortEdgeFilter2D.H"
#include "extrude2DMesh.H"
#include "polyMesh.H"
#include "patchToPoly2DMesh.H"
#include "extrudeModel.H"
#include "polyTopoChange.H"
#include "edgeCollapser.H"
#include "relaxationModel.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
int main(int argc, char *argv[])
{
argList::noParallel();
argList::validArgs.clear();
argList::validOptions.insert("pointsFile", "<filename>");
#include "addOverwriteOption.H"
#include "setRootCase.H"
#include "createTime.H"
// Read control dictionary
// ~~~~~~~~~~~~~~~~~~~~~~~
dictionary controlDict(IFstream("system/" + args.executable() + "Dict")());
dictionary shortEdgeFilterDict(controlDict.subDict("shortEdgeFilter"));
dictionary extrusionDict(controlDict.subDict("extrusion"));
Switch extrude(extrusionDict.lookup("extrude"));
const bool overwrite = args.optionFound("overwrite");
autoPtr<relaxationModel> relax
(
relaxationModel::New
(
controlDict.subDict("motionControl"),
runTime
)
);
// Read and triangulation
// ~~~~~~~~~~~~~~~~~~~~~~
CV2D mesh(runTime, controlDict);
if (args.options().found("pointsFile"))
{
fileName pointFileName(IStringStream(args.options()["pointsFile"])());
mesh.insertPoints(pointFileName);
mesh.insertSurfacePointPairs();
mesh.boundaryConform();
}
else
{
mesh.insertGrid();
mesh.insertSurfacePointPairs();
mesh.boundaryConform();
}
while (runTime.loop())
{
Info<< nl << "Time = " << runTime.timeName() << endl;
Info<< "Relaxation = " << relax->relaxation() << endl;
mesh.newPoints(relax->relaxation());
}
mesh.write();
Info<< "Finished Delaunay in = "
<< runTime.cpuTimeIncrement() << " s." << endl;
Info<< "Begin filtering short edges:" << endl;
shortEdgeFilter2D sef(mesh, shortEdgeFilterDict);
sef.filter();
Info<< "Meshed surface after edge filtering :" << endl;
sef.fMesh().writeStats(Info);
Info<< "Write .obj file of the 2D mesh: MeshedSurface.obj" << endl;
sef.fMesh().write("MeshedSurface.obj");
Info<< "Finished filtering in = "
<< runTime.cpuTimeIncrement() << " s." << endl;
Info<< "Begin constructing a polyMesh:" << endl;
patchToPoly2DMesh poly2DMesh
(
sef.fMesh(),
sef.patchNames(),
sef.patchSizes(),
sef.mapEdgesRegion()
);
poly2DMesh.createMesh();
polyMesh pMesh
(
IOobject
(
polyMesh::defaultRegion,
runTime.constant(),
runTime,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
xferMove(poly2DMesh.points()),
xferMove(poly2DMesh.faces()),
xferMove(poly2DMesh.owner()),
xferMove(poly2DMesh.neighbour())
);
Info<< "Constructing patches." << endl;
List<polyPatch*> patches(poly2DMesh.patchNames().size());
forAll(patches, patchI)
{
patches[patchI] = new polyPatch
(
poly2DMesh.patchNames()[patchI],
poly2DMesh.patchSizes()[patchI],
poly2DMesh.patchStarts()[patchI],
patchI,
pMesh.boundaryMesh()
);
}
pMesh.addPatches(patches);
if (extrude)
{
Info<< "Begin extruding the polyMesh:" << endl;
{
// Point generator
autoPtr<extrudeModel> model(extrudeModel::New(extrusionDict));
extrude2DMesh extruder(pMesh, extrusionDict, model());
extruder.addFrontBackPatches();
polyTopoChange meshMod(pMesh.boundaryMesh().size());
extruder.setRefinement(meshMod);
autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(pMesh, false);
pMesh.updateMesh(morphMap);
}
{
edgeCollapser collapser(pMesh);
const edgeList& edges = pMesh.edges();
const pointField& points = pMesh.points();
const boundBox& bb = pMesh.bounds();
const scalar mergeDim = 1E-4 * bb.minDim();
forAll(edges, edgeI)
{
const edge& e = edges[edgeI];
scalar d = e.mag(points);
if (d < mergeDim)
{
Info<< "Merging edge " << e << " since length " << d
<< " << " << mergeDim << endl;
// Collapse edge to e[0]
collapser.collapseEdge(edgeI, e[0]);
}
}
polyTopoChange meshModCollapse(pMesh);
collapser.setRefinement(meshModCollapse);
// Create a mesh from topo changes.
autoPtr<mapPolyMesh> morphMap =
meshModCollapse.changeMesh(pMesh, false);
pMesh.updateMesh(morphMap);
}
if (!overwrite)
{
runTime++;
}
else
{
pMesh.setInstance("constant");
}
}
pMesh.write();
Info<< "Finished extruding in = "
<< runTime.cpuTimeIncrement() << " s." << endl;
Info<< nl << "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,208 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
root "";
case "";
instance "";
local "";
class dictionary;
object cv2DMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
geometry
{
laurence_clean_preciser.stl
{
name laurence_clean_preciser;
type closedTriSurfaceMesh;
//type triSurfaceMesh;
}
// refinementBox
// {
// type searchableBox;
// min (-0.5 0.35 -1000);
// max (-0.5 0.35 1000);
// }
// refinementSphere
// {
// type searchableSphere;
// centre (0.85 0.4 0.0);
// radius 0.01;
// }
}
surfaceConformation
{
locationInMesh (-2.8 0.7 0.5);
pointPairDistanceCoeff 0.005;
minEdgeLenCoeff 0.005;
maxNotchLenCoeff 0.003;
minNearPointDistCoeff 0.0025;
maxQuadAngle 125;
// Insert near-boundary point mirror or point-pairs
insertSurfaceNearestPointPairs yes;
// Mirror near-boundary points rather than insert point-pairs
mirrorPoints no;
// Insert point-pairs vor dual-cell vertices very near the surface
insertSurfaceNearPointPairs yes;
// Maximum number of iterations used in boundaryConform.
maxBoundaryConformingIter 5;
geometryToConformTo
{
laurence_clean_preciser
{
featureMethod extendedFeatureEdgeMesh;
extendedFeatureEdgeMesh "laurence_clean_preciser.extendedFeatureEdgeMesh";
}
}
additionalFeatures
{
}
// Choose if to randomise the initial grid created by insertGrid.
randomiseInitialGrid yes;
// Perturbation fraction, 1 = cell-size.
randomPerturbation 0.1;
}
motionControl
{
defaultCellSize 0.05;
// Assign a priority to all requests for cell sizes, the highest overrules.
defaultPriority 0;
cellSizeControlGeometry
{
laurence_clean_preciser
{
priority 1;
mode bothSides;
cellSizeFunction linearDistance;
linearDistanceCoeffs
{
distanceCellSize 0.05;
surfaceCellSize 0.01;
distance 0.5;
}
uniformCoeffs
{
cellSize 0.01;
}
}
// refinementBox
// {
// priority 1;
// mode outside;
// cellSizeFunction linearDistance;
// linearDistanceCoeffs
// {
// distanceCellSize 0.04;
// surfaceCellSize 0.005;
// distance 0.2;
// }
// }
// refinementSphere
// {
// priority 1;
// mode outside;
// cellSizeFunction linearDistance;
// linearDistanceCoeffs
// {
// distanceCellSize 0.04;
// surfaceCellSize 0.005;
// distance 0.2;
// }
// }
}
relaxationModel adaptiveLinear;
adaptiveLinearCoeffs
{
relaxationStart 0.5;
relaxationEnd 0.0;
}
objOutput no;
// Near-wall region where cells are aligned with the wall specified as a number
// of cell layers
nearWallAlignedDist 3;
}
shortEdgeFilter
{
// Factor to multiply the average of a face's edge lengths by.
// If an edge of that face is smaller than that value then delete it.
shortEdgeFilterFactor 0.2;
// Weighting for the lengths of edges that are attached to the boundaries.
// Used when calculating the length of an edge. Default 2.0.
edgeAttachedToBoundaryFactor 2.0;
}
extrusion
{
extrude on;
extrudeModel linearDirection;
//extrudeModel wedge;
patchInfo
{
//type empty;
//startFace
}
patchType empty;
//patchType wedge;
nLayers 1;
expansionRatio 1.0; //0.9;
linearDirectionCoeffs
{
direction (0 0 1);
thickness 0.1;
}
wedgeCoeffs
{
axisPt (0 0 0);
axis (1 0 0);
angle 10;
}
thickness 0.1;
}

View File

@ -0,0 +1,3 @@
* Displacement limiting
http://en.wikipedia.org/wiki/Geometric_algebra

View File

@ -0,0 +1,147 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2011 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/>.
Class
indexedFace
Description
An indexed form of CGAL::Triangulation_face_base_2<K> used to keep
track of the vertices in the triangulation.
\*---------------------------------------------------------------------------*/
#ifndef indexedFace_H
#define indexedFace_H
#include <CGAL/Triangulation_2.h>
#include "indexedVertex.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace CGAL
{
/*---------------------------------------------------------------------------*\
Class indexedFace Declaration
\*---------------------------------------------------------------------------*/
template<class Gt, class Fb=CGAL::Triangulation_face_base_2<Gt> >
class indexedFace
:
public Fb
{
// Private data
//- The index for this triangle face
// -1: triangle and changed and associated data needs updating
// >=0: index of triangles, set by external update algorithm
int index_;
public:
enum faceTypes
{
UNCHANGED = 0,
CHANGED = -1,
SAVE_CHANGED = -2
};
typedef typename Fb::Vertex_handle Vertex_handle;
typedef typename Fb::Face_handle Face_handle;
template < typename TDS2 >
struct Rebind_TDS
{
typedef typename Fb::template Rebind_TDS<TDS2>::Other Fb2;
typedef indexedFace<Gt, Fb2> Other;
};
indexedFace()
:
Fb(),
index_(CHANGED)
{}
indexedFace(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2)
:
Fb(v0, v1, v2),
index_(CHANGED)
{}
indexedFace
(
Vertex_handle v0,
Vertex_handle v1,
Vertex_handle v2,
Face_handle n0,
Face_handle n1,
Face_handle n2
)
:
Fb(v0, v1, v2, n0, n1, n2),
index_(CHANGED)
{}
void set_vertex(int i, Vertex_handle v)
{
index_ = CHANGED;
Fb::set_vertex(i, v);
}
void set_vertices()
{
index_ = CHANGED;
Fb::set_vertices();
}
void set_vertices(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2)
{
index_ = CHANGED;
Fb::set_vertices(v0, v1, v2);
}
int& faceIndex()
{
return index_;
}
int faceIndex() const
{
return index_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace CGAL
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,268 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2011 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/>.
Class
indexedVertex
Description
An indexed form of CGAL::Triangulation_vertex_base_2<K> used to keep
track of the vertices in the triangulation.
\*---------------------------------------------------------------------------*/
#ifndef indexedVertex_H
#define indexedVertex_H
#include <CGAL/Triangulation_2.h>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace CGAL
{
/*---------------------------------------------------------------------------*\
Class indexedVertex Declaration
\*---------------------------------------------------------------------------*/
template<class Gt, class Vb=CGAL::Triangulation_vertex_base_2<Gt> >
class indexedVertex
:
public Vb
{
// Private data
//- The index for this triangle vertex
int index_;
//- Index of pair-point :
// NEAR_BOUNDARY_POINT : internal near boundary point.
// INTERNAL_POINT : internal point.
// FAR_POINT : far-point.
// >= 0 : part of point-pair. Index of other point.
// Lowest numbered is inside one (master).
int type_;
public:
enum pointTypes
{
NEAR_BOUNDARY_POINT = -4,
INTERNAL_POINT = -3,
MIRROR_POINT = -2,
FAR_POINT = -1
};
typedef typename Vb::Vertex_handle Vertex_handle;
typedef typename Vb::Face_handle Face_handle;
typedef typename Vb::Point Point;
template<typename TDS2>
struct Rebind_TDS
{
typedef typename Vb::template Rebind_TDS<TDS2>::Other Vb2;
typedef indexedVertex<Gt,Vb2> Other;
};
indexedVertex()
:
Vb(),
index_(INTERNAL_POINT),
type_(INTERNAL_POINT)
{}
indexedVertex(const Point& p)
:
Vb(p),
index_(INTERNAL_POINT),
type_(INTERNAL_POINT)
{}
indexedVertex(const Point& p, const int index, const int& type)
:
Vb(p),
index_(index),
type_(type)
{}
indexedVertex(const Point& p, Face_handle f)
:
Vb(f, p),
index_(INTERNAL_POINT),
type_(INTERNAL_POINT)
{}
indexedVertex(Face_handle f)
:
Vb(f),
index_(INTERNAL_POINT),
type_(INTERNAL_POINT)
{}
int& index()
{
return index_;
}
int index() const
{
return index_;
}
int& type()
{
return type_;
}
int type() const
{
return type_;
}
//- Is point a far-point
inline bool farPoint() const
{
return type_ == FAR_POINT;
}
//- Is point internal, i.e. not on boundary
inline bool internalPoint() const
{
return type_ <= INTERNAL_POINT;
}
//- Is point internal and near the boundary
inline bool nearBoundary() const
{
return type_ == NEAR_BOUNDARY_POINT;
}
//- Set the point to be near the boundary
inline void setNearBoundary()
{
type_ = NEAR_BOUNDARY_POINT;
}
//- Is point a mirror point
inline bool mirrorPoint() const
{
return type_ == MIRROR_POINT;
}
//- Either master or slave of pointPair.
inline bool pairPoint() const
{
return type_ >= 0;
}
//- Master of a pointPair is the lowest numbered one.
inline bool ppMaster() const
{
if (type_ > index_)
{
return true;
}
else
{
return false;
}
}
//- Slave of a pointPair is the highest numbered one.
inline bool ppSlave() const
{
if (type_ >= 0 && type_ < index_)
{
return true;
}
else
{
return false;
}
}
//- Either original internal point or master of pointPair.
inline bool internalOrBoundaryPoint() const
{
return internalPoint() || ppMaster();
}
//- Is point near the boundary or part of the boundary definition
inline bool nearOrOnBoundary() const
{
return pairPoint() || mirrorPoint() || nearBoundary();
}
//- Do the two given vertices consitute a boundary point-pair
inline friend bool pointPair
(
const indexedVertex& v0,
const indexedVertex& v1
)
{
return v0.index_ == v1.type_ || v1.index_ == v0.type_;
}
//- Do the three given vertices consitute a boundary triangle
inline friend bool boundaryTriangle
(
const indexedVertex& v0,
const indexedVertex& v1,
const indexedVertex& v2
)
{
return (v0.pairPoint() && pointPair(v1, v2))
|| (v1.pairPoint() && pointPair(v2, v0))
|| (v2.pairPoint() && pointPair(v0, v1));
}
//- Do the three given vertices consitute an outside triangle
inline friend bool outsideTriangle
(
const indexedVertex& v0,
const indexedVertex& v1,
const indexedVertex& v2
)
{
return (v0.farPoint() || v0.ppSlave())
|| (v1.farPoint() || v1.ppSlave())
|| (v2.farPoint() || v2.ppSlave());
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace CGAL
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,323 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2011 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/>.
\*----------------------------------------------------------------------------*/
#include "CV2D.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::CV2D::insertPointPair
(
Triangulation::Finite_vertices_iterator& vit,
const point2D& p,
const label trii,
const label hitSurface
)
{
if
(
!meshControls().mirrorPoints()
|| !insertMirrorPoint(toPoint2D(vit->point()), p)
)
{
pointIndexHit pHit
(
true,
toPoint3D(p),
trii
);
vectorField norm(1);
qSurf_.geometry()[hitSurface].getNormal
(
List<pointIndexHit>(1, pHit),
norm
);
insertPointPair
(
meshControls().ppDist(),
p,
toPoint2D(norm[0])
);
}
vit = Triangulation::Finite_vertices_iterator
(
CGAL::Filter_iterator
<
Triangulation::All_vertices_iterator,
Triangulation::Infinite_tester
>(finite_vertices_end(), vit.predicate(), vit.base())
);
}
bool Foam::CV2D::insertPointPairAtIntersection
(
Triangulation::Finite_vertices_iterator& vit,
const point2D& defVert,
const point2D vertices[],
const scalar maxProtSize2
)
{
bool found = false;
point2D interPoint;
label interTri = -1;
label interHitSurface = -1;
scalar interDist2 = 0;
Face_circulator fcStart = incident_faces(vit);
Face_circulator fc = fcStart;
label vi = 0;
do
{
if (!is_infinite(fc))
{
pointIndexHit pHit;
label hitSurface = -1;
qSurf_.findSurfaceNearestIntersection
(
toPoint3D(defVert),
toPoint3D(vertices[vi]),
pHit,
hitSurface
);
if (pHit.hit())
{
scalar dist2 =
magSqr(toPoint2D(pHit.hitPoint()) - vertices[vi]);
// Check the point is further away than the furthest so far
if (dist2 > interDist2)
{
scalar mps2 = maxProtSize2;
// If this is a boundary triangle reset the tolerance
// to avoid finding a hit point very close to a boundary
// vertex
if (boundaryTriangle(fc))
{
mps2 = meshControls().maxNotchLen2();
}
if (dist2 > mps2)
{
found = true;
interPoint = toPoint2D(pHit.hitPoint());
interTri = pHit.index();
interDist2 = dist2;
interHitSurface = hitSurface;
}
}
}
vi++;
}
} while (++fc != fcStart);
if (found)
{
insertPointPair(vit, interPoint, interTri, interHitSurface);
return true;
}
else
{
return false;
}
}
Foam::label Foam::CV2D::insertBoundaryConformPointPairs
(
const fileName& fName
)
{
label nIntersections = 0;
for
(
Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
vit++
)
{
// Consider only those points part of point-pairs or near boundary
if (!vit->nearOrOnBoundary())
{
continue;
}
// Counter-clockwise circulator
Face_circulator fcStart = incident_faces(vit);
Face_circulator fc = fcStart;
bool infinite = false;
bool changed = false;
do
{
if (is_infinite(fc))
{
infinite = true;
break;
}
else if (fc->faceIndex() < Fb::UNCHANGED)
{
changed = true;
break;
}
} while (++fc != fcStart);
// If the dual-cell is connected to the infinite point or none of the
// faces whose circumcentres it uses have changed ignore
if (infinite || !changed) continue;
fc = fcStart;
label nVerts = 0;
do
{
vertices[nVerts++] = toPoint2D(circumcenter(fc));
if (nVerts == maxNvert)
{
break;
}
} while (++fc != fcStart);
// Check if dual-cell has a large number of faces in which case
// assumed to be in the far-field and reject
if (nVerts == maxNvert) continue;
// Set n+1 vertex to the first vertex for easy circulating
vertices[nVerts] = vertices[0];
// Convert triangle vertex to OpenFOAM point
point2DFromPoint defVert = toPoint2D(vit->point());
scalar maxProtSize2 = meshControls().maxNotchLen2();
if (vit->internalOrBoundaryPoint())
{
// Calculate metrics of the dual-cell
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// The perimeter of the dual-cell
scalar perimeter = 0;
// Twice the area of the dual-cell
scalar areaT2 = 0;
for (int vi=0; vi<nVerts; vi++)
{
vector2D edge(vertices[vi+1] - vertices[vi]);
perimeter += mag(edge);
vector2D otherEdge = defVert - vertices[vi];
areaT2 += mag(edge.x()*otherEdge.y() - edge.y()*otherEdge.x());
}
// If the dual-cell is very small reject refinement
if (areaT2 < meshControls().minEdgeLen2()) continue;
// Estimate the cell width
scalar cellWidth = areaT2/perimeter;
// Check dimensions of dual-cell
/*
// Quick rejection of dual-cell refinement based on it's perimeter
if (perimeter < 2*meshControls().minCellSize()) continue;
// Also check the area of the cell and reject refinement
// if it is less than that allowed
if (areaT2 < meshControls().minCellSize2()) continue;
// Estimate the cell width and reject refinement if it is less than
// that allowed
if (cellWidth < 0.5*meshControls().minEdgeLen()) continue;
*/
if
(
perimeter > 2*meshControls().minCellSize()
&& areaT2 > meshControls().minCellSize2()
&& cellWidth > 0.5*meshControls().minEdgeLen()
)
{
maxProtSize2 = 0.25*meshControls().maxNotchLen2();
}
}
if (insertPointPairAtIntersection(vit, defVert, vertices, maxProtSize2))
{
nIntersections++;
}
}
return nIntersections;
}
void Foam::CV2D::markNearBoundaryPoints()
{
label count = 0;
for
(
Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
vit++
)
{
if (vit->internalPoint())
{
point vert(toPoint3D(vit->point()));
pointIndexHit pHit;
label hitSurface = -1;
qSurf_.findSurfaceNearest
(
vert,
4*meshControls().minCellSize2(),
pHit,
hitSurface
);
if (pHit.hit())
{
vit->setNearBoundary();
++count;
}
}
}
Info<< count << " points marked as being near a boundary" << endl;
}
// ************************************************************************* //

View File

@ -0,0 +1,397 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2011 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/>.
\*----------------------------------------------------------------------------*/
#include "CV2D.H"
#include "plane.H"
#include "triSurfaceTools.H"
#include "unitConversion.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::CV2D::on2DLine(const point2D& p, const linePointRef& line)
{
const point2D& a = toPoint2D(line.start());
const point2D& b = toPoint2D(line.end());
if
(
p.x() < min(a.x(), b.x())
|| p.x() > max(a.x(), b.x())
|| p.y() < min(a.y(), b.y())
|| p.y() > max(a.y(), b.y())
)
{
return false;
}
return true;
}
// Create feature points/edges by creating a triplet in the corner.
// (this triplet will have as its circumcentre the feature)
void Foam::CV2D::insertFeaturePoints()
{
featurePoints_.clear();
label nVert = number_of_vertices();
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
qSurf_.features()
);
if (feMeshes.empty())
{
WarningIn("CV2D::insertFeaturePoints")
<< "Extended Feature Edge Mesh is empty so no feature points will "
<< "be found." << nl
<< " Use: featureMethod extendedFeatureEdgeMesh;" << nl
<< endl;
}
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
const edgeList& edges = feMesh.edges();
const pointField& points = feMesh.points();
if (debug)
{
label nConvex = feMesh.concaveStart() - feMesh.convexStart();
label nConcave = feMesh.mixedStart() - feMesh.concaveStart();
label nMixed = feMesh.nonFeatureStart() - feMesh.mixedStart();
label nExternal = feMesh.internalStart() - feMesh.externalStart();
label nInternal = feMesh.flatStart() - feMesh.internalStart();
label nFlat = feMesh.openStart() - feMesh.flatStart();
label nOpen = feMesh.multipleStart() - feMesh.openStart();
label nMultiple = edges.size() - feMesh.multipleStart();
Info<< "Inserting Feature Points:" << nl
<< " Convex points: " << nConvex << nl
<< " Concave points: " << nConcave << nl
<< " Mixed points: " << nMixed << nl
<< " External edges: " << nExternal << nl
<< " Internal edges: " << nInternal << nl
<< " Flat edges: " << nFlat << nl
<< " Open edges: " << nOpen << nl
<< " Multiple edges: " << nMultiple << endl;
}
// Args: (base point, normal)
// @todo allow user to input this
plane zPlane(vector(0, 0, z_), vector(0, 0, 1));
if (debug)
{
Info<< " plane: " << zPlane << " " << z_ << endl;
}
forAll(edges, edgeI)
{
const edge& e = feMesh.edges()[edgeI];
const point& ep0 = points[e.start()];
const point& ep1 = points[e.end()];
const linePointRef line(ep0, ep1);
scalar intersect = zPlane.lineIntersect(line);
point2D featPoint = toPoint2D(intersect * (ep1 - ep0) + ep0);
if (on2DLine(featPoint, line))
{
vector2DField fpn = toPoint2D(feMesh.edgeNormals(edgeI));
vector2D cornerNormal = sum(fpn);
cornerNormal /= mag(cornerNormal);
if (debug)
{
Info<< nl << " line: " << line << nl
<< " vec: " << line.vec() << nl
<< " featurePoint: " << featPoint << nl
<< " line length: " << line.mag() << nl
<< " intersect: " << intersect << endl;
}
if
(
feMesh.getEdgeStatus(edgeI)
== extendedFeatureEdgeMesh::EXTERNAL
)
{
// Convex Point
Foam::point2D internalPt =
featPoint - meshControls().ppDist()*cornerNormal;
if (debug)
{
Info<< "PREC: " << internalPt << nl
<< " : " << featPoint << nl
<< " : " << meshControls().ppDist() << nl
<< " : " << cornerNormal << endl;
}
featurePoints_.push_back
(
Vb
(
toPoint(internalPt),
nVert,
nVert + 1
)
);
label masterPtIndex = nVert++;
forAll(fpn, nI)
{
const vector n3D(fpn[nI][0], fpn[nI][1], 0.0);
plane planeN = plane(toPoint3D(featPoint), n3D);
Foam::point2D externalPt =
internalPt
+ (
2.0
* planeN.distance(toPoint3D(internalPt))
* fpn[nI]
);
featurePoints_.push_back
(
Vb
(
toPoint(externalPt),
nVert++,
masterPtIndex
)
);
if (debug)
{
Info<< " side point: " << externalPt << endl;
}
}
if (debug)
{
Info<< "Convex Point: " << featPoint << nl
<< " corner norm: " << cornerNormal << nl
<< " reference: " << internalPt << endl;
}
}
else if
(
feMesh.getEdgeStatus(edgeI)
== extendedFeatureEdgeMesh::INTERNAL
)
{
// Concave Point
Foam::point2D externalPt =
featPoint + meshControls().ppDist()*cornerNormal;
Foam::point2D refPt =
featPoint - meshControls().ppDist()*cornerNormal;
label slavePointIndex = 0;
scalar totalAngle =
radToDeg
(
constant::mathematical::pi
+ acos(mag(fpn[0] & fpn[1]))
);
// Number of quadrants the angle should be split into
int nQuads =
int(totalAngle/meshControls().maxQuadAngle()) + 1;
// The number of additional master points needed to
//obtain the required number of quadrants.
int nAddPoints = min(max(nQuads - 2, 0), 2);
// index of reflMaster
label reflectedMaster = nVert + 2 + nAddPoints;
if (debug)
{
Info<< "Concave Point: " << featPoint << nl
<< " corner norm: " << cornerNormal << nl
<< " external: " << externalPt << nl
<< " reference: " << refPt << nl
<< " angle: " << totalAngle << nl
<< " nQuads: " << nQuads << nl
<< " nAddPoints: " << nAddPoints << endl;
}
forAll(fpn, nI)
{
const vector n3D(fpn[nI][0], fpn[nI][1], 0.0);
plane planeN = plane(toPoint3D(featPoint), n3D);
Foam::point2D internalPt =
externalPt
- (
2.0
* planeN.distance(toPoint3D(externalPt))
* fpn[nI]
);
featurePoints_.push_back
(
Vb
(
toPoint(internalPt),
nVert,
reflectedMaster
)
);
slavePointIndex = nVert++;
if (debug)
{
Info<< "Internal Point: " << internalPt << endl;
}
}
if (nAddPoints == 1)
{
// One additional point is the reflection of the slave
// point, i.e., the original reference point
featurePoints_.push_back
(
Vb
(
toPoint(refPt),
nVert++,
reflectedMaster
)
);
if (debug)
{
Info<< "ref Point: " << refPt << endl;
}
}
else if (nAddPoints == 2)
{
point2D reflectedAa =
refPt - ((featPoint - externalPt) & fpn[1])*fpn[1];
featurePoints_.push_back
(
Vb
(
toPoint(reflectedAa),
nVert++,
reflectedMaster
)
);
point2D reflectedBb =
refPt - ((featPoint - externalPt) & fpn[0])*fpn[0];
featurePoints_.push_back
(
Vb
(
toPoint(reflectedBb),
nVert++,
reflectedMaster
)
);
if (debug)
{
Info<< "refA Point: " << reflectedAa << nl
<< "refb Point: " << reflectedBb << endl;
}
}
featurePoints_.push_back
(
Vb
(
toPoint(externalPt),
nVert++,
slavePointIndex
)
);
}
else
{
WarningIn("void Foam::CV2D::insertFeaturePoints()")
<< "Feature Edge " << edges[edgeI] << nl
<< " points(" << points[edges[edgeI].start()]
<< ", " << points[edges[edgeI].end()] << ")" << nl
<< " is not labelled as either concave or convex, it"
<< " is labelled as (#2 = flat): "
<< feMesh.getEdgeStatus(edgeI) << endl;
}
}
else
{
WarningIn("void Foam::CV2D::insertFeaturePoints()")
<< "Point " << featPoint << " is not on the line "
<< line << endl;
}
}
}
// Insert the feature points.
reinsertFeaturePoints();
if (meshControls().objOutput())
{
writePoints("feat_allPoints.obj", false);
writeFaces("feat_allFaces.obj", false);
writeFaces("feat_faces.obj", true);
writeTriangles("feat_triangles.obj", true);
}
}
void Foam::CV2D::reinsertFeaturePoints()
{
for
(
std::list<Vb>::iterator vit=featurePoints_.begin();
vit != featurePoints_.end();
++vit
)
{
insertPoint
(
toPoint2D(vit->point()),
vit->index(),
vit->type()
);
}
}
// ************************************************************************* //

View File

@ -0,0 +1,114 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2011 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/>.
\*----------------------------------------------------------------------------*/
#include "CV2D.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::CV2D::insertSurfaceNearPointPairs()
{
Info<< "insertSurfaceNearPointPairs: ";
label nNearPoints = 0;
for
(
Triangulation::Finite_edges_iterator eit = finite_edges_begin();
eit != finite_edges_end();
eit++
)
{
Vertex_handle v0h = eit->first->vertex(cw(eit->second));
Vertex_handle v1h = eit->first->vertex(ccw(eit->second));
if (v0h->ppMaster() && v1h->ppMaster())
{
point2DFromPoint v0(toPoint2D(v0h->point()));
point2DFromPoint v1(toPoint2D(v1h->point()));
// Check that the two triangle vertices are further apart than the
// minimum cell size
if (magSqr(v1 - v0) > meshControls().minCellSize2())
{
point2D e0(toPoint2D(circumcenter(eit->first)));
point2D e1
(
toPoint2D(circumcenter(eit->first->neighbor(eit->second)))
);
// Calculate the length^2 of the edge normal to the surface
scalar edgeLen2 = magSqr(e0 - e1);
if (edgeLen2 < meshControls().minNearPointDist2())
{
pointIndexHit pHit;
label hitSurface = -1;
qSurf_.findSurfaceNearest
(
toPoint3D(e0),
meshControls().minEdgeLen2(),
pHit,
hitSurface
);
if (pHit.hit())
{
vectorField norm(1);
qSurf_.geometry()[hitSurface].getNormal
(
List<pointIndexHit>(1, pHit),
norm
);
insertPointPair
(
meshControls().ppDist(),
toPoint2D(pHit.hitPoint()),
toPoint2D(norm[0])
);
nNearPoints++;
// Correct the edge iterator for the change in the
// number of edges following the point-pair insertion
eit = Finite_edges_iterator
(
finite_edges_end().base(),
eit.predicate(),
eit.base()
);
}
}
}
}
}
Info<< nNearPoints << " point-pairs inserted" << endl;
}
// ************************************************************************* //

View File

@ -0,0 +1,244 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2011 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/>.
\*----------------------------------------------------------------------------*/
#include "CV2D.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::CV2D::dualCellSurfaceIntersection
(
const Triangulation::Finite_vertices_iterator& vit
) const
{
Triangulation::Edge_circulator ecStart = incident_edges(vit);
Triangulation::Edge_circulator ec = ecStart;
do
{
if (!is_infinite(ec))
{
point e0 = toPoint3D(circumcenter(ec->first));
// If edge end is outside bounding box then edge cuts boundary
if (!qSurf_.globalBounds().contains(e0))
{
return true;
}
point e1 = toPoint3D(circumcenter(ec->first->neighbor(ec->second)));
// If other edge end is ouside bounding box then edge cuts boundary
if (!qSurf_.globalBounds().contains(e1))
{
return true;
}
if (magSqr(e1 - e0) > meshControls().minEdgeLen2())
{
if (qSurf_.findSurfaceAnyIntersection(e0, e1))
{
return true;
}
}
}
} while (++ec != ecStart);
return false;
}
void Foam::CV2D::insertPointPairs
(
const DynamicList<point2D>& nearSurfacePoints,
const DynamicList<point2D>& surfacePoints,
const DynamicList<label>& surfaceTris,
const DynamicList<label>& surfaceHits,
const fileName fName
)
{
if (meshControls().mirrorPoints())
{
forAll(surfacePoints, ppi)
{
insertMirrorPoint
(
nearSurfacePoints[ppi],
surfacePoints[ppi]
);
}
}
else
{
forAll(surfacePoints, ppi)
{
pointIndexHit pHit
(
true,
toPoint3D(surfacePoints[ppi]),
surfaceTris[ppi]
);
vectorField norm(1);
qSurf_.geometry()[surfaceHits[ppi]].getNormal
(
List<pointIndexHit>(1, pHit),
norm
);
insertPointPair
(
meshControls().ppDist(),
surfacePoints[ppi],
toPoint2D(norm[0])
);
}
}
Info<< surfacePoints.size() << " point-pairs inserted" << endl;
if (meshControls().objOutput())
{
OFstream str(fName);
label vertI = 0;
forAll(surfacePoints, ppi)
{
meshTools::writeOBJ(str, toPoint3D(surfacePoints[ppi]));
vertI++;
}
Info<< "insertPointPairs: Written " << surfacePoints.size()
<< " inserted point-pair locations to file "
<< str.name() << endl;
}
}
void Foam::CV2D::insertSurfaceNearestPointPairs()
{
Info<< "insertSurfaceNearestPointPairs: ";
label nSurfacePointsEst =
min
(
number_of_vertices(),
size_t(10*sqrt(scalar(number_of_vertices())))
);
DynamicList<point2D> nearSurfacePoints(nSurfacePointsEst);
DynamicList<point2D> surfacePoints(nSurfacePointsEst);
DynamicList<label> surfaceTris(nSurfacePointsEst);
DynamicList<label> surfaceHits(nSurfacePointsEst);
// Local references to surface mesh addressing
// const pointField& localPoints = qSurf_.localPoints();
// const labelListList& edgeFaces = qSurf_.edgeFaces();
// const vectorField& faceNormals = qSurf_.faceNormals();
// const labelListList& faceEdges = qSurf_.faceEdges();
for
(
Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
vit++
)
{
if (vit->internalPoint())
{
point2DFromPoint vert(toPoint2D(vit->point()));
pointIndexHit pHit;
label hitSurface = -1;
qSurf_.findSurfaceNearest
(
toPoint3D(vert),
4*meshControls().minCellSize2(),
pHit,
hitSurface
);
if (pHit.hit())
{
vit->setNearBoundary();
// Reference to the nearest triangle
// const labelledTri& f = qSurf_[hitSurface];
// // Find where point is on triangle.
// // Note tolerance needed is relative one
// // (used in comparing normalized [0..1] triangle coordinates).
// label nearType, nearLabel;
// triPointRef
// (
// localPoints[f[0]],
// localPoints[f[1]],
// localPoints[f[2]]
// ).classify(pHit.hitPoint(), nearType, nearLabel);
// // If point is on a edge check if it is an internal feature
// bool internalFeatureEdge = false;
// if (nearType == triPointRef::EDGE)
// {
// label edgeI = faceEdges[hitSurface][nearLabel];
// const labelList& eFaces = edgeFaces[edgeI];
// if
// (
// eFaces.size() == 2
// && (faceNormals[eFaces[0]] & faceNormals[eFaces[1]])
// < -0.2
// )
// {
// internalFeatureEdge = true;
// }
// }
if (dualCellSurfaceIntersection(vit)) //&& !internalFeatureEdge)
{
nearSurfacePoints.append(vert);
surfacePoints.append(toPoint2D(pHit.hitPoint()));
surfaceTris.append(pHit.index());
surfaceHits.append(hitSurface);
}
}
}
}
insertPointPairs
(
nearSurfacePoints,
surfacePoints,
surfaceTris,
surfaceHits,
"surfaceNearestIntersections.obj"
);
}
// ************************************************************************* //

View File

@ -0,0 +1,505 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
\*---------------------------------------------------------------------------*/
#include "shortEdgeFilter2D.H"
namespace Foam
{
defineTypeNameAndDebug(shortEdgeFilter2D, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::shortEdgeFilter2D::shortEdgeFilter2D
(
const Foam::CV2D& cv2Dmesh,
const dictionary& dict
)
:
cv2Dmesh_(cv2Dmesh),
shortEdgeFilterFactor_(readScalar(dict.lookup("shortEdgeFilterFactor"))),
edgeAttachedToBoundaryFactor_
(
dict.lookupOrDefault<scalar>("edgeAttachedToBoundaryFactor", 2.0)
),
patchNames_(wordList()),
patchSizes_(labelList()),
mapEdgesRegion_()
{
point2DField points2D;
faceList faces;
cv2Dmesh.calcDual
(
points2D,
faces,
patchNames_,
patchSizes_,
mapEdgesRegion_
);
pointField points(points2D.size());
forAll(points, ip)
{
points[ip] = cv2Dmesh.toPoint3D(points2D[ip]);
}
points2D.clear();
ms_ = MeshedSurface<face>(xferMove(points), xferMove(faces));
Info<< "Meshed surface stats before edge filtering :" << endl;
ms_.writeStats(Info);
if (debug)
{
writeInfo(Info);
ms_.write("MeshedSurface_preFilter.obj");
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::shortEdgeFilter2D::~shortEdgeFilter2D()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void
Foam::shortEdgeFilter2D::filter()
{
// These are global indices.
const pointField& points = ms_.points();
const edgeList& edges = ms_.edges();
const faceList& faces = ms_.faces();
const labelList& meshPoints = ms_.meshPoints();
const labelList& boundaryPoints = ms_.boundaryPoints();
label maxChain = 0;
label nPointsToRemove = 0;
labelList pointsToRemove(ms_.points().size(), -1);
// List of number of vertices in a face.
labelList newFaceVertexCount(faces.size(), -1);
forAll(faces, faceI)
{
newFaceVertexCount[faceI] = faces[faceI].size();
}
// Check if the point is a boundary point. Flag if it is so that
// it will not be deleted.
boolList boundaryPointFlags(points.size(), false);
// This has been removed, otherwise small edges on the boundary are not
// removed.
/* forAll(boundaryPointFlags, pointI)
{
forAll(boundaryPoints, bPoint)
{
if (meshPoints[boundaryPoints[bPoint]] == pointI)
{
boundaryPointFlags[pointI] = true;
}
}
}*/
// Check if an edge has a boundary point. It it does the edge length
// will be doubled when working out its length.
Info<< " Marking edges attached to boundaries." << endl;
boolList edgeAttachedToBoundary(edges.size(), false);
forAll(edges, edgeI)
{
const edge& e = edges[edgeI];
const label startVertex = e.start();
const label endVertex = e.end();
forAll(boundaryPoints, bPoint)
{
if
(
boundaryPoints[bPoint] == startVertex
|| boundaryPoints[bPoint] == endVertex
)
{
edgeAttachedToBoundary[edgeI] = true;
}
}
}
forAll(edges, edgeI)
{
const edge& e = edges[edgeI];
// get the vertices of that edge.
const label startVertex = e.start();
const label endVertex = e.end();
scalar edgeLength =
mag
(
points[meshPoints[e.start()]]
- points[meshPoints[e.end()]]
);
if (edgeAttachedToBoundary[edgeI])
{
edgeLength *= edgeAttachedToBoundaryFactor_;
}
scalar shortEdgeFilterValue = 0.0;
const labelList& psEdges = ms_.pointEdges()[startVertex];
const labelList& peEdges = ms_.pointEdges()[endVertex];
forAll(psEdges, psEdgeI)
{
const edge& psE = edges[psEdges[psEdgeI]];
if (edgeI != psEdges[psEdgeI])
{
shortEdgeFilterValue +=
mag
(
points[meshPoints[psE.start()]]
-points[meshPoints[psE.end()]]
);
}
}
forAll(peEdges, peEdgeI)
{
const edge& peE = edges[peEdges[peEdgeI]];
if (edgeI != peEdges[peEdgeI])
{
shortEdgeFilterValue +=
mag
(
points[meshPoints[peE.start()]]
-points[meshPoints[peE.end()]]
);
}
}
shortEdgeFilterValue *=
shortEdgeFilterFactor_
/(psEdges.size() + peEdges.size() - 2);
if (edgeLength < shortEdgeFilterValue)
{
bool flagDegenerateFace = false;
const labelList& pFaces = ms_.pointFaces()[startVertex];
forAll(pFaces, pFaceI)
{
const face& f = ms_.localFaces()[pFaces[pFaceI]];
forAll(f, fp)
{
// If the edge is part of this face...
if (f[fp] == endVertex)
{
// If deleting vertex would create a triangle, don't!
if (newFaceVertexCount[pFaces[pFaceI]] < 4)
{
flagDegenerateFace = true;
}
else
{
newFaceVertexCount[pFaces[pFaceI]]--;
}
}
// If the edge is not part of this face...
else
{
// Deleting vertex of a triangle...
if (newFaceVertexCount[pFaces[pFaceI]] < 3)
{
flagDegenerateFace = true;
}
}
}
}
// This if statement determines whether a point should be deleted.
if
(
pointsToRemove[meshPoints[startVertex]] == -1
&& pointsToRemove[meshPoints[endVertex]] == -1
&& !boundaryPointFlags[meshPoints[startVertex]]
&& !flagDegenerateFace
)
{
pointsToRemove[meshPoints[startVertex]] =
meshPoints[endVertex];
++nPointsToRemove;
}
}
}
label totalNewPoints = points.size() - nPointsToRemove;
pointField newPoints(totalNewPoints, vector(0, 0, 0));
labelList newPointNumbers(points.size(), -1);
label numberRemoved=0;
forAll(points, pointI)
{
// If the point is NOT going to be removed.
if (pointsToRemove[pointI] == -1)
{
newPoints[pointI-numberRemoved] = points[pointI];
newPointNumbers[pointI] = pointI-numberRemoved;
}
else
{
numberRemoved++;
}
}
// Need a new faceList
faceList newFaces(faces.size());
label newFaceI = 0;
labelList newFace;
label newFaceSize = 0;
// Now need to iterate over the faces and remove points. Global index.
forAll(faces, faceI)
{
const face& f = faces[faceI];
newFace.clear();
newFace.setSize(f.size());
newFaceSize = 0;
forAll(f, fp)
{
label pointI = f[fp];
// If not removing the point, then add it to the new face.
if (pointsToRemove[pointI] == -1)
{
newFace[newFaceSize++] = newPointNumbers[pointI];
}
else
{
label newPointI = pointsToRemove[pointI];
// Replace deleted point with point that it is being
// collapsed to.
if
(
f.nextLabel(fp) != newPointI
&& f.prevLabel(fp) != newPointI
)
{
label pChain = newPointI;
label totalChain = 0;
for (label nChain = 0; nChain <= totalChain; ++nChain)
{
if (newPointNumbers[pChain] != -1)
{
newFace[newFaceSize++] = newPointNumbers[pChain];
newPointNumbers[pointI]
= newPointNumbers[pChain];
maxChain = max(totalChain, maxChain);
}
else
{
WarningIn("shortEdgeFilter")
<< "Point " << pChain
<< " marked for deletion as well as point "
<< pointI << nl
<< " Incrementing maxChain by 1 from "
<< totalChain << " to " << totalChain + 1
<< endl;
totalChain++;
}
pChain = pointsToRemove[pChain];
}
}
else
{
if (newPointNumbers[newPointI] != -1)
{
newPointNumbers[pointI] = newPointNumbers[newPointI];
}
}
}
}
newFace.setSize(newFaceSize);
if (newFace.size() > 2)
{
newFaces[newFaceI++] = face(newFace);
}
else
{
FatalErrorIn("shortEdgeFilter")
<< "Only " << newFace.size() << " in face " << faceI
<< exit(FatalError);
}
}
newFaces.setSize(newFaceI);
MeshedSurface<face> fMesh
(
xferMove(newPoints),
xferMove(newFaces),
xferCopy(List<surfZone>())
);
const Map<int>& fMeshPointMap = fMesh.meshPointMap();
// Reset patchSizes_
patchSizes_.clear();
patchSizes_.setSize(patchNames_.size(), 0);
label equalEdges = 0;
label notFound = 0;
label matches = 0;
label negativeLabels = 0;
forAll(newPointNumbers, pointI)
{
if (newPointNumbers[pointI] == -1)
{
WarningIn("shortEdgeFilter")
<< pointI << " will be deleted and " << newPointNumbers[pointI]
<< ", so it will not be replaced. "
<< "This will cause edges to be deleted." << endl;
}
}
// Create new EdgeMap.
Info<< "Creating new EdgeMap." << endl;
EdgeMap<label> newMapEdgesRegion(mapEdgesRegion_.size());
for
(
label bEdgeI = ms_.nInternalEdges();
bEdgeI < edges.size();
++bEdgeI
)
{
label p1 = meshPoints[edges[bEdgeI][0]];
label p2 = meshPoints[edges[bEdgeI][1]];
edge e(p1, p2);
if (mapEdgesRegion_.found(e))
{
if
(
newPointNumbers[p1] != -1
&& newPointNumbers[p2] != -1
)
{
if (newPointNumbers[p1] != newPointNumbers[p2])
{
label region = mapEdgesRegion_.find(e)();
newMapEdgesRegion.insert
(
edge
(
fMeshPointMap[newPointNumbers[p1]],
fMeshPointMap[newPointNumbers[p2]]
),
region
);
patchSizes_[region]++;
matches++;
}
else
{
equalEdges++;
}
}
else
{
negativeLabels++;
}
}
else
{
notFound++;
}
}
if (debug)
{
Info<< "EdgeMapping :" << nl
<< " Matches : " << matches << nl
<< " Equal : " << equalEdges << nl
<< " Negative : " << negativeLabels << nl
<< " Not Found: " << notFound << endl;
}
mapEdgesRegion_.transfer(newMapEdgesRegion);
ms_.transfer(fMesh);
Info<< " Maximum number of chained collapses = " << maxChain << endl;
if (debug)
{
writeInfo(Info);
}
}
void Foam::shortEdgeFilter2D::writeInfo(Ostream& os)
{
os << "Short Edge Filtering Information:" << nl
<< " shortEdgeFilterFactor: " << shortEdgeFilterFactor_ << nl
<< " edgeAttachedToBoundaryFactor: " << edgeAttachedToBoundaryFactor_
<< endl;
forAll(patchNames_, patchI)
{
os << " Patch " << patchNames_[patchI]
<< ", size " << patchSizes_[patchI] << endl;
}
os << " There are " << mapEdgesRegion_.size()
<< " boundary edges." << endl;
os << " Mesh Info:" << nl
<< " Points: " << ms_.nPoints() << nl
<< " Faces: " << ms_.size() << nl
<< " Edges: " << ms_.nEdges() << nl
<< " Internal: " << ms_.nInternalEdges() << nl
<< " External: " << ms_.nEdges() - ms_.nInternalEdges()
<< endl;
}
// ************************************************************************* //

View File

@ -0,0 +1,136 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 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/>.
Class
Foam::shortEdgeFilter2D
Description
This class filters short edges generated by the CV2D mesher.
SourceFiles
shortEdgeFilter2D.C
\*---------------------------------------------------------------------------*/
#ifndef shortEdgeFilter2D_H
#define shortEdgeFilter2D_H
#include "MeshedSurfaces.H"
#include "CV2D.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class shortEdgeFilter2D Declaration
\*---------------------------------------------------------------------------*/
class shortEdgeFilter2D
{
// Private data
//- Description of data_
const CV2D& cv2Dmesh_;
MeshedSurface<face> ms_;
const scalar shortEdgeFilterFactor_;
const scalar edgeAttachedToBoundaryFactor_;
wordList patchNames_;
labelList patchSizes_;
EdgeMap<label> mapEdgesRegion_;
// Private Member Functions
//- Disallow default bitwise copy construct
shortEdgeFilter2D(const shortEdgeFilter2D&);
//- Disallow default bitwise assignment
void operator=(const shortEdgeFilter2D&);
public:
//- Runtime type information
ClassName("shortEdgeFilter2D");
// Constructors
shortEdgeFilter2D
(
const CV2D& cv2Dmesh,
const dictionary& dict
);
//- Destructor
~shortEdgeFilter2D();
// Access Functions
const wordList& patchNames() const
{
return patchNames_;
}
const labelList& patchSizes() const
{
return patchSizes_;
}
const EdgeMap<label>& mapEdgesRegion() const
{
return mapEdgesRegion_;
}
const MeshedSurface<face>& fMesh() const
{
return ms_;
}
// Member Functions
void filter();
void writeInfo(Ostream& os);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,5 @@
wmake/rules/General/CGAL
-lboost_thread
-lboost_thread-mt

View File

@ -152,12 +152,10 @@ bool Foam::cellSizeControlSurfaces::evalCellSizeFunctions
Foam::cellSizeControlSurfaces::cellSizeControlSurfaces Foam::cellSizeControlSurfaces::cellSizeControlSurfaces
( (
const conformalVoronoiMesh& cvMesh,
const searchableSurfaces& allGeometry, const searchableSurfaces& allGeometry,
const dictionary& motionControlDict const dictionary& motionControlDict
) )
: :
cvMesh_(cvMesh),
allGeometry_(allGeometry), allGeometry_(allGeometry),
surfaces_(), surfaces_(),
cellSizeFunctions_(), cellSizeFunctions_(),
@ -223,7 +221,6 @@ Foam::cellSizeControlSurfaces::cellSizeControlSurfaces
cellSizeFunction::New cellSizeFunction::New
( (
surfaceSubDict, surfaceSubDict,
cvMesh,
surface surface
) )
); );
@ -286,38 +283,38 @@ Foam::scalar Foam::cellSizeControlSurfaces::cellSize
bool anyFunctionFound = evalCellSizeFunctions(pt, size); bool anyFunctionFound = evalCellSizeFunctions(pt, size);
if (!anyFunctionFound) //if (!anyFunctionFound)
{ //{
// Check if the point in question was actually inside the domain, if // // Check if the point in question was actually inside the domain, if
// not, then it may be falling back to an inappropriate default size. // // not, then it may be falling back to an inappropriate default size.
if (cvMesh_.geometryToConformTo().outside(pt)) // if (cvMesh_.geometryToConformTo().outside(pt))
{ // {
pointIndexHit surfHit; // pointIndexHit surfHit;
label hitSurface; // label hitSurface;
cvMesh_.geometryToConformTo().findSurfaceNearest // cvMesh_.geometryToConformTo().findSurfaceNearest
( // (
pt, // pt,
sqr(GREAT), // sqr(GREAT),
surfHit, // surfHit,
hitSurface // hitSurface
); // );
if (!surfHit.hit()) // if (!surfHit.hit())
{ // {
FatalErrorIn // FatalErrorIn
( // (
"Foam::scalar Foam::cellSizeControlSurfaces::cellSize" // "Foam::scalar Foam::cellSizeControlSurfaces::cellSize"
"(" // "("
"const point& pt" // "const point& pt"
") const" // ") const"
) // )
<< "Point " << pt << " did not find a nearest surface point" // << "Point " << pt << " did not find a nearest surface point"
<< nl << exit(FatalError) << endl; // << nl << exit(FatalError) << endl;
} // }
} // }
} //}
return size; return size;
} }

View File

@ -57,9 +57,6 @@ class cellSizeControlSurfaces
{ {
// Private data // Private data
//- Reference to the conformalVoronoiMesh holding this object
const conformalVoronoiMesh& cvMesh_;
//- Reference to the searchableSurfaces object holding all geometry data //- Reference to the searchableSurfaces object holding all geometry data
const searchableSurfaces& allGeometry_; const searchableSurfaces& allGeometry_;
@ -106,7 +103,6 @@ public:
// searchableSurfaces // searchableSurfaces
cellSizeControlSurfaces cellSizeControlSurfaces
( (
const conformalVoronoiMesh& cvMesh,
const searchableSurfaces& allGeometry, const searchableSurfaces& allGeometry,
const dictionary& motionControlDict const dictionary& motionControlDict
); );

View File

@ -45,12 +45,10 @@ cellSizeFunction::cellSizeFunction
( (
const word& type, const word& type,
const dictionary& cellSizeFunctionDict, const dictionary& cellSizeFunctionDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
) )
: :
dictionary(cellSizeFunctionDict), dictionary(cellSizeFunctionDict),
cvMesh_(cvMesh),
surface_(surface), surface_(surface),
coeffsDict_(subDict(type + "Coeffs")), coeffsDict_(subDict(type + "Coeffs")),
sideMode_(), sideMode_(),
@ -99,7 +97,6 @@ cellSizeFunction::cellSizeFunction
autoPtr<cellSizeFunction> cellSizeFunction::New autoPtr<cellSizeFunction> cellSizeFunction::New
( (
const dictionary& cellSizeFunctionDict, const dictionary& cellSizeFunctionDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
) )
{ {
@ -130,7 +127,7 @@ autoPtr<cellSizeFunction> cellSizeFunction::New
return autoPtr<cellSizeFunction> return autoPtr<cellSizeFunction>
( (
cstrIter()(cellSizeFunctionDict, cvMesh, surface) cstrIter()(cellSizeFunctionDict, surface)
); );
} }

View File

@ -81,9 +81,6 @@ protected:
// Protected data // Protected data
//- Reference to the conformalVoronoiMesh holding this cvs object
const conformalVoronoiMesh& cvMesh_;
//- Reference to the searchableSurface that cellSizeFunction //- Reference to the searchableSurface that cellSizeFunction
// relates to // relates to
const searchableSurface& surface_; const searchableSurface& surface_;
@ -120,10 +117,9 @@ public:
dictionary, dictionary,
( (
const dictionary& cellSizeFunctionDict, const dictionary& cellSizeFunctionDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
), ),
(cellSizeFunctionDict, cvMesh, surface) (cellSizeFunctionDict, surface)
); );
@ -134,7 +130,6 @@ public:
( (
const word& type, const word& type,
const dictionary& cellSizeFunctionDict, const dictionary& cellSizeFunctionDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
); );
@ -145,7 +140,6 @@ public:
static autoPtr<cellSizeFunction> New static autoPtr<cellSizeFunction> New
( (
const dictionary& cellSizeFunctionDict, const dictionary& cellSizeFunctionDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
); );

View File

@ -42,11 +42,10 @@ addToRunTimeSelectionTable(cellSizeFunction, linearDistance, dictionary);
linearDistance::linearDistance linearDistance::linearDistance
( (
const dictionary& initialPointsDict, const dictionary& initialPointsDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
) )
: :
cellSizeFunction(typeName, initialPointsDict, cvMesh, surface), cellSizeFunction(typeName, initialPointsDict, surface),
surfaceCellSize_(readScalar(coeffsDict().lookup("surfaceCellSize"))), surfaceCellSize_(readScalar(coeffsDict().lookup("surfaceCellSize"))),
distanceCellSize_(readScalar(coeffsDict().lookup("distanceCellSize"))), distanceCellSize_(readScalar(coeffsDict().lookup("distanceCellSize"))),
distance_(readScalar(coeffsDict().lookup("distance"))), distance_(readScalar(coeffsDict().lookup("distance"))),

View File

@ -87,7 +87,6 @@ public:
linearDistance linearDistance
( (
const dictionary& initialPointsDict, const dictionary& initialPointsDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
); );

View File

@ -42,11 +42,10 @@ addToRunTimeSelectionTable(cellSizeFunction, linearSpatial, dictionary);
linearSpatial::linearSpatial linearSpatial::linearSpatial
( (
const dictionary& initialPointsDict, const dictionary& initialPointsDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
) )
: :
cellSizeFunction(typeName, initialPointsDict, cvMesh, surface), cellSizeFunction(typeName, initialPointsDict, surface),
referencePoint_(coeffsDict().lookup("referencePoint")), referencePoint_(coeffsDict().lookup("referencePoint")),
referenceCellSize_(readScalar(coeffsDict().lookup("referenceCellSize"))), referenceCellSize_(readScalar(coeffsDict().lookup("referenceCellSize"))),
direction_(coeffsDict().lookup("direction")), direction_(coeffsDict().lookup("direction")),

View File

@ -85,7 +85,6 @@ public:
linearSpatial linearSpatial
( (
const dictionary& initialPointsDict, const dictionary& initialPointsDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
); );

View File

@ -46,11 +46,10 @@ addToRunTimeSelectionTable
surfaceOffsetLinearDistance::surfaceOffsetLinearDistance surfaceOffsetLinearDistance::surfaceOffsetLinearDistance
( (
const dictionary& initialPointsDict, const dictionary& initialPointsDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
) )
: :
cellSizeFunction(typeName, initialPointsDict, cvMesh, surface), cellSizeFunction(typeName, initialPointsDict, surface),
surfaceCellSize_(readScalar(coeffsDict().lookup("surfaceCellSize"))), surfaceCellSize_(readScalar(coeffsDict().lookup("surfaceCellSize"))),
distanceCellSize_(readScalar(coeffsDict().lookup("distanceCellSize"))), distanceCellSize_(readScalar(coeffsDict().lookup("distanceCellSize"))),
surfaceOffset_(readScalar(coeffsDict().lookup("surfaceOffset"))), surfaceOffset_(readScalar(coeffsDict().lookup("surfaceOffset"))),

View File

@ -94,7 +94,6 @@ public:
surfaceOffsetLinearDistance surfaceOffsetLinearDistance
( (
const dictionary& initialPointsDict, const dictionary& initialPointsDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
); );

View File

@ -41,11 +41,10 @@ addToRunTimeSelectionTable(cellSizeFunction, uniform, dictionary);
uniform::uniform uniform::uniform
( (
const dictionary& initialPointsDict, const dictionary& initialPointsDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
) )
: :
cellSizeFunction(typeName, initialPointsDict, cvMesh, surface), cellSizeFunction(typeName, initialPointsDict, surface),
cellSize_(readScalar(coeffsDict().lookup("cellSize"))) cellSize_(readScalar(coeffsDict().lookup("cellSize")))
{} {}

View File

@ -69,7 +69,6 @@ public:
uniform uniform
( (
const dictionary& initialPointsDict, const dictionary& initialPointsDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
); );

View File

@ -41,11 +41,10 @@ addToRunTimeSelectionTable(cellSizeFunction, uniformDistance, dictionary);
uniformDistance::uniformDistance uniformDistance::uniformDistance
( (
const dictionary& initialPointsDict, const dictionary& initialPointsDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
) )
: :
cellSizeFunction(typeName, initialPointsDict, cvMesh, surface), cellSizeFunction(typeName, initialPointsDict, surface),
cellSize_(readScalar(coeffsDict().lookup("cellSize"))), cellSize_(readScalar(coeffsDict().lookup("cellSize"))),
distance_(readScalar(coeffsDict().lookup("distance"))), distance_(readScalar(coeffsDict().lookup("distance"))),
distanceSqr_(sqr(distance_)) distanceSqr_(sqr(distance_))

View File

@ -75,7 +75,6 @@ public:
uniformDistance uniformDistance
( (
const dictionary& initialPointsDict, const dictionary& initialPointsDict,
const conformalVoronoiMesh& cvMesh,
const searchableSurface& surface const searchableSurface& surface
); );

View File

@ -1925,13 +1925,13 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
), ),
geometryToConformTo_ geometryToConformTo_
( (
*this, runTime_,
rndGen_,
allGeometry_, allGeometry_,
cvMeshDict.subDict("surfaceConformation") cvMeshDict.subDict("surfaceConformation")
), ),
cellSizeControl_ cellSizeControl_
( (
*this,
allGeometry_, allGeometry_,
cvMeshDict.subDict("motionControl") cvMeshDict.subDict("motionControl")
), ),
@ -1958,7 +1958,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
relaxationModel::New relaxationModel::New
( (
cvMeshDict.subDict("motionControl"), cvMeshDict.subDict("motionControl"),
*this runTime_
) )
), ),
faceAreaWeightModel_ faceAreaWeightModel_

View File

@ -30,12 +30,14 @@ License
Foam::conformationSurfaces::conformationSurfaces Foam::conformationSurfaces::conformationSurfaces
( (
const conformalVoronoiMesh& cvMesh, const Time& runTime,
Random& rndGen,
const searchableSurfaces& allGeometry, const searchableSurfaces& allGeometry,
const dictionary& surfaceConformationDict const dictionary& surfaceConformationDict
) )
: :
cvMesh_(cvMesh), runTime_(runTime),
rndGen_(rndGen),
allGeometry_(allGeometry), allGeometry_(allGeometry),
features_(), features_(),
locationInMesh_(surfaceConformationDict.lookup("locationInMesh")), locationInMesh_(surfaceConformationDict.lookup("locationInMesh")),
@ -139,9 +141,9 @@ Foam::conformationSurfaces::conformationSurfaces
IOobject IOobject
( (
feMeshName, feMeshName,
cvMesh_.time().constant(), runTime_.time().constant(),
"extendedFeatureEdgeMesh", "extendedFeatureEdgeMesh",
cvMesh_.time(), runTime_.time(),
IOobject::MUST_READ, IOobject::MUST_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
) )
@ -211,9 +213,9 @@ Foam::conformationSurfaces::conformationSurfaces
IOobject IOobject
( (
feMeshName, feMeshName,
cvMesh_.time().constant(), runTime_.time().constant(),
"extendedFeatureEdgeMesh", "extendedFeatureEdgeMesh",
cvMesh_.time(), runTime_.time(),
IOobject::MUST_READ, IOobject::MUST_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
) )
@ -232,7 +234,7 @@ Foam::conformationSurfaces::conformationSurfaces
// Extend the global bounds to stop the bound box sitting on the surfaces // Extend the global bounds to stop the bound box sitting on the surfaces
// to be conformed to // to be conformed to
globalBounds_ = globalBounds_.extend(cvMesh_.rndGen(), 1e-4); globalBounds_ = globalBounds_.extend(rndGen_, 1e-4);
// Look at all surfaces at determine whether the locationInMesh point is // Look at all surfaces at determine whether the locationInMesh point is
// inside or outside each, to establish a signature for the domain to be // inside or outside each, to establish a signature for the domain to be
@ -730,7 +732,7 @@ void Foam::conformationSurfaces::findEdgeNearestByType
void Foam::conformationSurfaces::writeFeatureObj(const fileName& prefix) const void Foam::conformationSurfaces::writeFeatureObj(const fileName& prefix) const
{ {
OFstream ftStr(cvMesh_.time().path()/prefix + "_allFeatures.obj"); OFstream ftStr(runTime_.time().path()/prefix + "_allFeatures.obj");
Pout<< nl << "Writing all features to " << ftStr.name() << endl; Pout<< nl << "Writing all features to " << ftStr.name() << endl;

View File

@ -56,8 +56,9 @@ class conformationSurfaces
{ {
// Private data // Private data
//- Reference to the conformalVoronoiMesh holding this object const Time& runTime_;
const conformalVoronoiMesh& cvMesh_;
Random& rndGen_;
//- Reference to the searchableSurfaces object holding all geometry data //- Reference to the searchableSurfaces object holding all geometry data
const searchableSurfaces& allGeometry_; const searchableSurfaces& allGeometry_;
@ -114,7 +115,8 @@ public:
// searchableSurfaces // searchableSurfaces
conformationSurfaces conformationSurfaces
( (
const conformalVoronoiMesh& cvMesh, const Time& runTime,
Random& rndGen,
const searchableSurfaces& allGeometry, const searchableSurfaces& allGeometry,
const dictionary& surfaceConformationDict const dictionary& surfaceConformationDict
); );

View File

@ -41,13 +41,13 @@ addToRunTimeSelectionTable(relaxationModel, adaptiveLinear, dictionary);
adaptiveLinear::adaptiveLinear adaptiveLinear::adaptiveLinear
( (
const dictionary& relaxationDict, const dictionary& relaxationDict,
const conformalVoronoiMesh& cvMesh const Time& runTime
) )
: :
relaxationModel(typeName, relaxationDict, cvMesh), relaxationModel(typeName, relaxationDict, runTime),
relaxationStart_(readScalar(coeffDict().lookup("relaxationStart"))), relaxationStart_(readScalar(coeffDict().lookup("relaxationStart"))),
relaxationEnd_(readScalar(coeffDict().lookup("relaxationEnd"))), relaxationEnd_(readScalar(coeffDict().lookup("relaxationEnd"))),
lastTimeValue_(cvMesh_.time().timeOutputValue()), lastTimeValue_(runTime_.time().timeOutputValue()),
relaxation_(relaxationStart_) relaxation_(relaxationStart_)
{} {}
@ -56,24 +56,24 @@ adaptiveLinear::adaptiveLinear
scalar adaptiveLinear::relaxation() scalar adaptiveLinear::relaxation()
{ {
if (cvMesh_.time().timeOutputValue() > lastTimeValue_) if (runTime_.time().timeOutputValue() > lastTimeValue_)
{ {
scalar currentRelxation = relaxation_; scalar currentRelaxation = relaxation_;
relaxation_ -= relaxation_ -=
(relaxation_ - relaxationEnd_) (relaxation_ - relaxationEnd_)
/( /(
( (
cvMesh_.time().endTime().value() runTime_.time().endTime().value()
- cvMesh_.time().timeOutputValue() - runTime_.time().timeOutputValue()
) )
/(cvMesh_.time().timeOutputValue() - lastTimeValue_) /(runTime_.time().timeOutputValue() - lastTimeValue_)
+ 1 + 1
); );
lastTimeValue_ = cvMesh_.time().timeOutputValue(); lastTimeValue_ = runTime_.time().timeOutputValue();
return currentRelxation; return currentRelaxation;
} }
return relaxation_; return relaxation_;

View File

@ -83,7 +83,7 @@ public:
adaptiveLinear adaptiveLinear
( (
const dictionary& relaxationDict, const dictionary& relaxationDict,
const conformalVoronoiMesh& cvMesh const Time& runTime
); );

View File

@ -41,10 +41,10 @@ addToRunTimeSelectionTable(relaxationModel, rampHoldFall, dictionary);
rampHoldFall::rampHoldFall rampHoldFall::rampHoldFall
( (
const dictionary& relaxationDict, const dictionary& relaxationDict,
const conformalVoronoiMesh& cvMesh const Time& runTime
) )
: :
relaxationModel(typeName, relaxationDict, cvMesh), relaxationModel(typeName, relaxationDict, runTime),
rampStartRelaxation_(readScalar(coeffDict().lookup("rampStartRelaxation"))), rampStartRelaxation_(readScalar(coeffDict().lookup("rampStartRelaxation"))),
holdRelaxation_(readScalar(coeffDict().lookup("holdRelaxation"))), holdRelaxation_(readScalar(coeffDict().lookup("holdRelaxation"))),
fallEndRelaxation_(readScalar(coeffDict().lookup("fallEndRelaxation"))), fallEndRelaxation_(readScalar(coeffDict().lookup("fallEndRelaxation"))),
@ -62,10 +62,10 @@ rampHoldFall::rampHoldFall
scalar rampHoldFall::relaxation() scalar rampHoldFall::relaxation()
{ {
scalar t = cvMesh_.time().timeOutputValue(); scalar t = runTime_.time().timeOutputValue();
scalar tStart = cvMesh_.time().startTime().value(); scalar tStart = runTime_.time().startTime().value();
scalar tEnd = cvMesh_.time().endTime().value(); scalar tEnd = runTime_.time().endTime().value();
scalar tSpan = tEnd - tStart; scalar tSpan = tEnd - tStart;
if (tSpan < VSMALL) if (tSpan < VSMALL)

View File

@ -88,7 +88,7 @@ public:
rampHoldFall rampHoldFall
( (
const dictionary& relaxationDict, const dictionary& relaxationDict,
const conformalVoronoiMesh& cvMesh const Time& runTime
); );

View File

@ -43,11 +43,11 @@ relaxationModel::relaxationModel
( (
const word& type, const word& type,
const dictionary& relaxationDict, const dictionary& relaxationDict,
const conformalVoronoiMesh& cvMesh const Time& runTime
) )
: :
dictionary(relaxationDict), dictionary(relaxationDict),
cvMesh_(cvMesh), runTime_(runTime),
coeffDict_(subDict(type + "Coeffs")) coeffDict_(subDict(type + "Coeffs"))
{} {}
@ -57,7 +57,7 @@ relaxationModel::relaxationModel
autoPtr<relaxationModel> relaxationModel::New autoPtr<relaxationModel> relaxationModel::New
( (
const dictionary& relaxationDict, const dictionary& relaxationDict,
const conformalVoronoiMesh& cvMesh const Time& runTime
) )
{ {
word relaxationModelTypeName word relaxationModelTypeName
@ -85,7 +85,7 @@ autoPtr<relaxationModel> relaxationModel::New
<< exit(FatalError); << exit(FatalError);
} }
return autoPtr<relaxationModel>(cstrIter()(relaxationDict, cvMesh)); return autoPtr<relaxationModel>(cstrIter()(relaxationDict, runTime));
} }

View File

@ -37,7 +37,7 @@ SourceFiles
#define relaxationModel_H #define relaxationModel_H
#include "point.H" #include "point.H"
#include "conformalVoronoiMesh.H" #include "Time.H"
#include "dictionary.H" #include "dictionary.H"
#include "autoPtr.H" #include "autoPtr.H"
#include "runTimeSelectionTables.H" #include "runTimeSelectionTables.H"
@ -61,7 +61,7 @@ protected:
// Protected data // Protected data
//- Reference to the conformalVoronoiMesh holding this cvControls object //- Reference to the conformalVoronoiMesh holding this cvControls object
const conformalVoronoiMesh& cvMesh_; const Time& runTime_;
//- Method coeffs dictionary //- Method coeffs dictionary
dictionary coeffDict_; dictionary coeffDict_;
@ -92,9 +92,9 @@ public:
dictionary, dictionary,
( (
const dictionary& relaxationDict, const dictionary& relaxationDict,
const conformalVoronoiMesh& cvMesh const Time& runTime
), ),
(relaxationDict, cvMesh) (relaxationDict, runTime)
); );
@ -105,7 +105,7 @@ public:
( (
const word& type, const word& type,
const dictionary& relaxationDict, const dictionary& relaxationDict,
const conformalVoronoiMesh& cvMesh const Time& runTime
); );
@ -115,7 +115,7 @@ public:
static autoPtr<relaxationModel> New static autoPtr<relaxationModel> New
( (
const dictionary& relaxationDict, const dictionary& relaxationDict,
const conformalVoronoiMesh& cvMesh const Time& runTime
); );

View File

@ -4,6 +4,8 @@ set -x
wmake libso extrudeModel wmake libso extrudeModel
wmake extrudeMesh wmake extrudeMesh
wmake libso extrudeToRegionMesh/createShellMesh
wmake extrudeToRegionMesh wmake extrudeToRegionMesh

View File

@ -1,4 +1,3 @@
createShellMesh.C
extrudeToRegionMesh.C extrudeToRegionMesh.C
EXE = $(FOAM_APPBIN)/extrudeToRegionMesh EXE = $(FOAM_APPBIN)/extrudeToRegionMesh

View File

@ -1,11 +1,13 @@
EXE_INC = \ EXE_INC = \
-I../extrudeModel/lnInclude \ -I../extrudeModel/lnInclude \
-IcreateShellMesh/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude -I$(LIB_SRC)/dynamicMesh/lnInclude
EXE_LIBS = \ EXE_LIBS = \
-lextrudeModel \ -lextrudeModel \
-lcreateShellMesh \
-lfiniteVolume \ -lfiniteVolume \
-lmeshTools \ -lmeshTools \
-ldynamicMesh -ldynamicMesh

View File

@ -0,0 +1,3 @@
createShellMesh.C
LIB = $(FOAM_LIBBIN)/libcreateShellMesh

View File

@ -0,0 +1,9 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude -DFULLDEBUG -g -O0
LIB_LIBS = \
-lfiniteVolume \
-lmeshTools \
-ldynamicMesh

View File

@ -509,7 +509,6 @@ void Foam::createShellMesh::setRefinement
} }
} }
// Introduce original points // Introduce original points
// ~~~~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~~~~~
@ -526,13 +525,12 @@ void Foam::createShellMesh::setRefinement
); );
pointToPointMap.append(pointI); pointToPointMap.append(pointI);
//Pout<< "Added bottom point " << addedPointI // Pout<< "Added bottom point " << pointToPointMap[pointI]
// << " at " << patch_.localPoints()[pointI] // << " at " << patch_.localPoints()[pointI]
// << " from point " << pointI // << " from point " << pointI
// << endl; // << endl;
} }
// Introduce new points (one for every region) // Introduce new points (one for every region)
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -543,6 +541,7 @@ void Foam::createShellMesh::setRefinement
point pt = patch_.localPoints()[pointI]; point pt = patch_.localPoints()[pointI];
point disp = firstLayerDisp[regionI]; point disp = firstLayerDisp[regionI];
for (label layerI = 0; layerI < nLayers; layerI++) for (label layerI = 0; layerI < nLayers; layerI++)
{ {
pt += disp; pt += disp;
@ -676,7 +675,7 @@ void Foam::createShellMesh::setRefinement
{ {
FatalErrorIn("createShellMesh::setRefinement(..)") FatalErrorIn("createShellMesh::setRefinement(..)")
<< "external/feature edge:" << edgeI << "external/feature edge:" << edgeI
<< " has " << eFaces.size() << " connected extruded faces " << " has " << eFaces.size() << " connected extruded faces"
<< " but only " << ePatches.size() << " but only " << ePatches.size()
<< " boundary faces defined." << exit(FatalError); << " boundary faces defined." << exit(FatalError);
} }

View File

@ -0,0 +1,8 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wclean libso extrude2DMesh
wclean
# ----------------------------------------------------------------- end-of-file

View File

@ -0,0 +1,8 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wmake libso extrude2DMesh
wmake
# ----------------------------------------------------------------- end-of-file

View File

@ -1,4 +1,3 @@
extrude2DMesh.C
extrude2DMeshApp.C extrude2DMeshApp.C
EXE = $(FOAM_APPBIN)/extrude2DMesh EXE = $(FOAM_APPBIN)/extrude2DMesh

View File

@ -1,10 +1,12 @@
EXE_INC = \ EXE_INC = \
/* -DFULLDEBUG -g -O0 */ \ -I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \ -I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude -Iextrude2DMesh/lnInclude \
-I../extrude/extrudeModel/lnInclude
EXE_LIBS = \ EXE_LIBS = \
-lmeshTools \ -lsurfMesh \
-ldynamicMesh \ -ldynamicMesh \
-lfiniteVolume -lextrude2DMesh \
-lextrudeModel

View File

@ -1,263 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
\*---------------------------------------------------------------------------*/
#include "extrude2DMesh.H"
#include "polyMesh.H"
#include "polyTopoChange.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(extrude2DMesh, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from mesh
Foam::extrude2DMesh::extrude2DMesh(const polyMesh& mesh)
:
mesh_(mesh)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::extrude2DMesh::setRefinement
(
const direction extrudeDir,
const scalar thickness,
const label frontPatchI,
polyTopoChange& meshMod
) const
{
for (label cellI = 0; cellI < mesh_.nCells(); cellI++)
{
meshMod.addCell
(
-1, //masterPointID,
-1, //masterEdgeID,
-1, //masterFaceID,
cellI, //masterCellID,
mesh_.cellZones().whichZone(cellI) //zoneID
);
}
// Generate points
// ~~~~~~~~~~~~~~~
forAll(mesh_.points(), pointI)
{
meshMod.addPoint
(
mesh_.points()[pointI],
pointI,
-1, // zoneID
true // inCell
);
}
//Info<< "Adding offsetted points." << nl << endl;
forAll(mesh_.points(), pointI)
{
point newPoint(mesh_.points()[pointI]);
newPoint[extrudeDir] += thickness;
meshMod.addPoint
(
newPoint,
pointI,
-1, // zoneID
true // inCell
);
}
// Generate faces
// ~~~~~~~~~~~~~~
const faceList& faces = mesh_.faces();
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
label zoneID = mesh_.faceZones().whichZone(faceI);
bool zoneFlip = false;
if (zoneID != -1)
{
const faceZone& fZone = mesh_.faceZones()[zoneID];
zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
}
face newFace(4);
const face& f = faces[faceI];
newFace[0] = f[0];
newFace[1] = f[1];
newFace[2] = f[1]+mesh_.nPoints();
newFace[3] = f[0]+mesh_.nPoints();
meshMod.addFace
(
newFace,
mesh_.faceOwner()[faceI], // own
mesh_.faceNeighbour()[faceI], // nei
-1, // masterPointID
-1, // masterEdgeID
faceI, // masterFaceID
false, // flipFaceFlux
-1, // patchID
zoneID, // zoneID
zoneFlip // zoneFlip
);
}
forAll(patches, patchI)
{
label startFaceI = patches[patchI].start();
label endFaceI = startFaceI + patches[patchI].size();
for (label faceI = startFaceI; faceI < endFaceI; faceI++)
{
label zoneID = mesh_.faceZones().whichZone(faceI);
bool zoneFlip = false;
if (zoneID != -1)
{
const faceZone& fZone = mesh_.faceZones()[zoneID];
zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
}
face newFace(4);
const face& f = faces[faceI];
newFace[0] = f[0];
newFace[1] = f[1];
newFace[2] = f[1]+mesh_.nPoints();
newFace[3] = f[0]+mesh_.nPoints();
meshMod.addFace
(
newFace,
mesh_.faceOwner()[faceI], // own
-1, // nei
-1, // masterPointID
-1, // masterEdgeID
faceI, // masterFaceID
false, // flipFaceFlux
patchI, // patchID
zoneID, // zoneID
zoneFlip // zoneFlip
);
}
}
// Generate front and back faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
forAll(mesh_.cells(), cellI)
{
const cell& cFaces = mesh_.cells()[cellI];
// Make a loop out of faces.
const face& f = faces[cFaces[0]];
face frontFace(cFaces.size());
frontFace[0] = f[0];
label nextPointI = f[1];
label nextFaceI = cFaces[0];
for (label i = 1; i < frontFace.size(); i++)
{
frontFace[i] = nextPointI;
// Find face containing pointI
forAll(cFaces, cFaceI)
{
label faceI = cFaces[cFaceI];
if (faceI != nextFaceI)
{
const face& f = faces[faceI];
if (f[0] == nextPointI)
{
nextPointI = f[1];
nextFaceI = faceI;
break;
}
else if (f[1] == nextPointI)
{
nextPointI = f[0];
nextFaceI = faceI;
break;
}
}
}
}
// Add back face.
meshMod.addFace
(
frontFace.reverseFace(),
cellI, // own
-1, // nei
-1, // masterPointID
-1, // masterEdgeID
cFaces[0], // masterFaceID
false, // flipFaceFlux
frontPatchI, // patchID
-1, // zoneID
false // zoneFlip
);
// Offset to create front face.
forAll(frontFace, fp)
{
frontFace[fp] += mesh_.nPoints();
}
meshMod.addFace
(
frontFace,
cellI, // own
-1, // nei
-1, // masterPointID
-1, // masterEdgeID
cFaces[0], // masterFaceID
false, // flipFaceFlux
frontPatchI, // patchID
-1, // zoneID
false // zoneFlip
);
}
}
// ************************************************************************* //

View File

@ -0,0 +1,5 @@
extrude2DMesh/extrude2DMesh.C
patchToPoly2DMesh/patchToPoly2DMesh.C
LIB = $(FOAM_LIBBIN)/libextrude2DMesh

View File

@ -0,0 +1,13 @@
EXE_INC = \
/* -DFULLDEBUG -g -O0 */ \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(FOAM_APP)/utilities/mesh/generation/extrude/extrudeModel/lnInclude
LIB_LIBS = \
-lmeshTools \
-ldynamicMesh \
-lsurfMesh \
-lfiniteVolume

View File

@ -0,0 +1,588 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
\*---------------------------------------------------------------------------*/
#include "extrude2DMesh.H"
#include "polyMesh.H"
#include "polyTopoChange.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(extrude2DMesh, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::extrude2DMesh::check2D() const
{
const faceList& faces = mesh_.faces();
forAll(faces, faceI)
{
if (faces[faceI].size() != 2)
{
FatalErrorIn("extrude2DMesh.C")
<< "Face " << faceI << " size " << faces[faceI].size()
<< " is not of size 2 so mesh is not proper two-dimensional."
<< exit(FatalError);
}
}
}
//void Foam::extrude2DMesh::findExtrudeDirection()
//{
// scalar minRange = GREAT;
// for (direction dir = 0; dir < 3; dir++)
// {
// scalarField cmpts(mesh_.points().component(dir));
// scalar range = max(cmpts)-min(cmpts);
// Info<< "Direction:" << dir << " range:" << range << endl;
// if (range < minRange)
// {
// minRange = range;
// extrudeDir_ = dir;
// }
// }
// Info<< "Extruding in direction " << extrudeDir_
// << " with thickness " << thickness_ << nl
// << endl;
//}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::extrude2DMesh::extrude2DMesh
(
polyMesh& mesh,
const dictionary& dict,
const extrudeModel& model
)
:
mesh_(mesh),
dict_(dict),
//patchDict_(dict.subDict("patchInfo")),
model_(model),
modelType_(dict.lookup("extrudeModel")),
patchType_(dict.lookup("patchType")),
frontPatchI_(-1),
backPatchI_(-1)
{
check2D();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::extrude2DMesh::~extrude2DMesh()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::extrude2DMesh::addFrontBackPatches()
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
frontPatchI_ = patches.findPatchID("front");
backPatchI_ = patches.findPatchID("back");
// Add patch.
List<polyPatch*> newPatches(patches.size() + 2);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
newPatches[patchI] =
pp.clone
(
patches,
newPatches.size(),
pp.size(),
pp.start()
).ptr();
}
if (frontPatchI_ == -1)
{
frontPatchI_ = patches.size();
newPatches[frontPatchI_] =
polyPatch::New
(
patchType_,
"front",
0,
mesh_.nFaces(),
frontPatchI_,
patches
).ptr();
// newPatches[frontPatchI_] = polyPatch::New
// (
// "front",
// patchDict_,
// frontPatchI_,
// patches
// ).ptr();
Info<< "Adding patch " << newPatches[frontPatchI_]->name()
<< " at index " << frontPatchI_
<< " for front faces." << nl << endl;
}
if (backPatchI_ == -1)
{
backPatchI_ = patches.size() + 1;
newPatches[backPatchI_] =
polyPatch::New
(
patchType_,
"back",
0,
mesh_.nFaces(),
backPatchI_,
patches
).ptr();
// newPatches[frontPatchI_] = polyPatch::New
// (
// "back",
// patchDict_,
// backPatchI_,
// patches
// ).ptr();
Info<< "Adding patch " << newPatches[backPatchI_]->name()
<< " at index " << backPatchI_
<< " for back faces." << nl << endl;
}
mesh_.removeBoundary();
mesh_.addPatches(newPatches);
}
void Foam::extrude2DMesh::setRefinement
(
polyTopoChange& meshMod
)
{
const label nLayers = model_.nLayers();
const pointField& points = mesh_.points();
label nFaces = 0;
for (label layer = 0; layer < nLayers; ++layer)
{
label offset = layer * mesh_.nCells();
forAll(mesh_.cells(), cellI)
{
meshMod.addCell
(
-1, //masterPointID,
-1, //masterEdgeID,
-1, //masterFaceID,
cellI + offset, //masterCellID,
mesh_.cellZones().whichZone(cellI) //zoneID
);
}
}
// Generate points
// ~~~~~~~~~~~~~~~
for (label layer = 0; layer <= nLayers; ++layer)
{
label offset = layer * points.size();
forAll(points, pointI)
{
// Don't need the surface normal for either linearDirection or
// wedge. Will need to add to be able to use others.
point newPoint = model_
(
points[pointI],
vector(),
layer
);
meshMod.addPoint
(
newPoint,
pointI + offset,
-1, // zoneID
true // inCell
);
}
Pout<< "Added " << points.size() << " points to layer "
<< layer << endl;
}
// Generate faces
// ~~~~~~~~~~~~~~
const faceList& faces = mesh_.faces();
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
for (label layer = 0; layer < nLayers; ++layer)
{
label currentLayerOffset = layer * mesh_.nPoints();
label nextLayerOffset = currentLayerOffset + mesh_.nPoints();
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
label zoneID = mesh_.faceZones().whichZone(faceI);
bool zoneFlip = false;
if (zoneID != -1)
{
const faceZone& fZone = mesh_.faceZones()[zoneID];
zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
}
face newFace(4);
const face& f = faces[faceI];
newFace[0] = f[0] + currentLayerOffset;
newFace[1] = f[1] + currentLayerOffset;
newFace[2] = f[1] + nextLayerOffset;
newFace[3] = f[0] + nextLayerOffset;
//{
// vector n = newFace.normal(pointField(meshMod.points()));
// label own = mesh_.faceOwner()[faceI];
// const labelList& ownPoints = mesh_.cellPoints()[own];
// point ownCc = sum(pointField(mesh_.points(), ownPoints))/ownPoints.size();
// label nei = mesh_.faceNeighbour()[faceI];
// const labelList& neiPoints = mesh_.cellPoints()[nei];
// point neiCc = sum(pointField(mesh_.points(), neiPoints))/neiPoints.size();
// vector d = neiCc - ownCc;
// Pout<< "face:" << faceI << " at:" << f.centre(mesh_.points()) << endl
// << " own:" << own << " at:" << ownCc << endl
// << " nei:" << nei << " at:" << neiCc << endl
// << " sign:" << (n & d) << endl
// << endl;
//}
label offset = layer * mesh_.nCells();
meshMod.addFace
(
newFace,
mesh_.faceOwner()[faceI] + offset, // own
mesh_.faceNeighbour()[faceI] + offset, // nei
-1, // masterPointID
-1, // masterEdgeID
nFaces++, // masterFaceID
false, // flipFaceFlux
-1, // patchID
zoneID, // zoneID
zoneFlip // zoneFlip
);
if (debug)
{
Info<< newFace << " "
<< mesh_.faceOwner()[faceI] + offset << " "
<< mesh_.faceNeighbour()[faceI] + offset << " "
<< nFaces - 1
<< endl;
}
}
}
forAll(patches, patchI)
{
for (label layer=0; layer < nLayers; layer++)
{
label currentLayerOffset = layer*mesh_.nPoints();
label nextLayerOffset = currentLayerOffset + mesh_.nPoints();
label startFaceI = patches[patchI].start();
label endFaceI = startFaceI + patches[patchI].size();
for (label faceI = startFaceI; faceI < endFaceI; faceI++)
{
label zoneID = mesh_.faceZones().whichZone(faceI);
bool zoneFlip = false;
if (zoneID != -1)
{
const faceZone& fZone = mesh_.faceZones()[zoneID];
zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
}
face newFace(4);
const face& f = faces[faceI];
newFace[0] = f[0] + currentLayerOffset;
newFace[1] = f[1] + currentLayerOffset;
newFace[2] = f[1] + nextLayerOffset;
newFace[3] = f[0] + nextLayerOffset;
label offset = layer * mesh_.nCells();
meshMod.addFace
(
newFace,
mesh_.faceOwner()[faceI] + offset, // own
-1, // nei
-1, // masterPointID
-1, // masterEdgeID
nFaces++, // masterFaceID
false, // flipFaceFlux
patchI, // patchID
zoneID, // zoneID
zoneFlip // zoneFlip
);
if (debug)
{
Info<< newFace << " "
<< mesh_.faceOwner()[faceI] + offset << " "
<< nFaces - 1
<< endl;
}
}
}
}
// Add extra internal faces that need special treatment for owners and
// neighbours.
forAll(mesh_.cells(), cellI)
{
const cell& cFaces = mesh_.cells()[cellI];
face frontFace(cFaces.size());
// Make a loop out of faces.
label nextFaceI = cFaces[0];
const face& f = faces[nextFaceI];
label nextPointI;
if (mesh_.faceOwner()[nextFaceI] == cellI)
{
frontFace[0] = f[0];
nextPointI = f[1];
}
else
{
frontFace[0] = f[1];
nextPointI = f[0];
}
for (label i = 1; i < frontFace.size(); i++)
{
frontFace[i] = nextPointI;
// Find face containing pointI
forAll(cFaces, cFaceI)
{
label faceI = cFaces[cFaceI];
if (faceI != nextFaceI)
{
const face& f = faces[faceI];
if (f[0] == nextPointI)
{
nextPointI = f[1];
nextFaceI = faceI;
break;
}
else if (f[1] == nextPointI)
{
nextPointI = f[0];
nextFaceI = faceI;
break;
}
}
}
}
for (label layer = 0; layer < nLayers - 1; ++layer)
{
// Offset to create front face.
forAll(frontFace, fp)
{
frontFace[fp] += mesh_.nPoints();
}
label offset = layer * mesh_.nCells();
label nei = -1;
if (layer != nLayers - 1)
{
nei = cellI + offset + mesh_.nCells();
}
meshMod.addFace
(
frontFace,
cellI + offset, // own
nei, // nei
-1, // masterPointID
-1, // masterEdgeID
nFaces++, // masterFaceID
false, // flipFaceFlux
-1, // patchID
-1, // zoneID
false // zoneFlip
);
if (debug)
{
Info<< frontFace << " "
<< cellI + offset << " "
<< nei << " "
<< nFaces - 1
<< endl;
}
}
}
// Generate front and back faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
forAll(mesh_.cells(), cellI)
{
const cell& cFaces = mesh_.cells()[cellI];
face frontFace(cFaces.size());
// Make a loop out of faces.
label nextFaceI = cFaces[0];
const face& f = faces[nextFaceI];
label nextPointI;
if (mesh_.faceOwner()[nextFaceI] == cellI)
{
frontFace[0] = f[0];
nextPointI = f[1];
}
else
{
frontFace[0] = f[1];
nextPointI = f[0];
}
for (label i = 1; i < frontFace.size(); i++)
{
frontFace[i] = nextPointI;
// Find face containing pointI
forAll(cFaces, cFaceI)
{
label faceI = cFaces[cFaceI];
if (faceI != nextFaceI)
{
const face& f = faces[faceI];
if (f[0] == nextPointI)
{
nextPointI = f[1];
nextFaceI = faceI;
break;
}
else if (f[1] == nextPointI)
{
nextPointI = f[0];
nextFaceI = faceI;
break;
}
}
}
}
// Add back face.
meshMod.addFace
(
frontFace.reverseFace(),
cellI, // own
-1, // nei
-1, // masterPointID
-1, // masterEdgeID
nFaces++, // masterFaceID
false, // flipFaceFlux
backPatchI_, // patchID
-1, // zoneID
false // zoneFlip
);
if (debug)
{
Info<< nl<<frontFace.reverseFace() << " "
<< cellI << " "
<< nFaces - 1
<< endl;
}
// Offset to create front face.
forAll(frontFace, fp)
{
frontFace[fp] += mesh_.nPoints()* (nLayers);
}
label offset = (nLayers - 1) * mesh_.nCells();
meshMod.addFace
(
frontFace,
cellI + offset, // own
-1, // nei
-1, // masterPointID
-1, // masterEdgeID
nFaces++, // masterFaceID
false, // flipFaceFlux
frontPatchI_, // patchID
-1, // zoneID
false // zoneFlip
);
if (debug)
{
Info<< frontFace << " "
<< cellI + offset << " "
<< nFaces - 1
<< endl;
}
}
}
// ************************************************************************* //

View File

@ -25,8 +25,8 @@ Class
Foam::extrude2DMesh Foam::extrude2DMesh
Description Description
Given 2D mesh insert all the topology changes to extrude. Does not work Given a 2D mesh insert all the topology changes to extrude. Does not work
in parallel in parallel.
SourceFiles SourceFiles
extrude2DMesh.C extrude2DMesh.C
@ -36,9 +36,12 @@ SourceFiles
#ifndef extrude2DMesh_H #ifndef extrude2DMesh_H
#define extrude2DMesh_H #define extrude2DMesh_H
#include "typeInfo.H"
#include "label.H" #include "label.H"
#include "scalar.H" #include "scalar.H"
#include "typeInfo.H" #include "labelList.H"
#include "dictionary.H"
#include "extrudeModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -50,6 +53,7 @@ class polyMesh;
class polyTopoChange; class polyTopoChange;
class mapPolyMesh; class mapPolyMesh;
class mapDistributePolyMesh; class mapDistributePolyMesh;
class polyBoundaryMesh;
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class extrude2DMesh Declaration Class extrude2DMesh Declaration
@ -60,8 +64,34 @@ class extrude2DMesh
// Private data // Private data
//- Reference to 2D mesh //- Reference to 2D mesh
const polyMesh& mesh_; polyMesh& mesh_;
const dictionary dict_;
//const dictionary patchDict_;
const extrudeModel& model_;
const word modelType_;
const word patchType_;
label frontPatchI_;
label backPatchI_;
// Private Member Functions
//- Check the mesh is 2D
void check2D() const;
//- Find extrusion direction
//void findExtrudeDirection();
//- Disallow default bitwise copy construct
extrude2DMesh(const extrude2DMesh&);
//- Disallow default bitwise assignment
void operator=(const extrude2DMesh&);
public: public:
@ -70,21 +100,25 @@ public:
// Constructors // Constructors
extrude2DMesh
(
polyMesh&,
const dictionary& dict,
const extrudeModel& model
);
//- Construct from mesh
extrude2DMesh(const polyMesh&); //- Destructor
~extrude2DMesh();
// Member Functions // Member Functions
//- Add front and back patches
void addFrontBackPatches();
//- Play commands into polyTopoChange to extrude mesh. //- Play commands into polyTopoChange to extrude mesh.
void setRefinement void setRefinement(polyTopoChange&);
(
const direction extrudeDir,
const scalar thickness,
const label frontPatchI,
polyTopoChange&
) const;
//- Force recalculation of locally stored data on topological change //- Force recalculation of locally stored data on topological change
void updateMesh(const mapPolyMesh&) void updateMesh(const mapPolyMesh&)
@ -93,6 +127,17 @@ public:
//- Force recalculation of locally stored data for mesh distribution //- Force recalculation of locally stored data for mesh distribution
void distribute(const mapDistributePolyMesh&) void distribute(const mapDistributePolyMesh&)
{} {}
label frontPatchI() const
{
return frontPatchI_;
}
label backPatchI() const
{
return backPatchI_;
}
}; };

View File

@ -0,0 +1,333 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
\*---------------------------------------------------------------------------*/
#include "patchToPoly2DMesh.H"
#include "PatchTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::patchToPoly2DMesh::flipFaceOrder()
{
const edgeList& edges = patch_.edges();
const faceList& localFaces = patch_.localFaces();
const labelList& meshPoints = patch_.meshPoints();
Info<< "Flipping face order if necessary." << endl;
forAll(edges, edgeI)
{
const edge& e = edges[edgeI];
faces_[edgeI].setSize(2);
label edgeOwner = owner_[edgeI];
const face& f = localFaces[edgeOwner];
label fp = findIndex(f, e[0]);
if (f.nextLabel(fp) != e[1])
{
Info<< "Flipping face " << faces_[edgeI] << endl;
faces_[edgeI][0] = meshPoints[e[1]];
faces_[edgeI][1] = meshPoints[e[0]];
}
else
{
faces_[edgeI][0] = meshPoints[e[0]];
faces_[edgeI][1] = meshPoints[e[1]];
}
}
}
void Foam::patchToPoly2DMesh::createNeighbours()
{
const edgeList& edges = patch_.edges();
const labelListList& edgeFaces = patch_.edgeFaces();
Info<< "Calculating neighbours." << endl;
forAll(edges, edgeI)
{
const labelList& eFaces = edgeFaces[edgeI];
if (eFaces.size() == 2)
{
if (owner_[edgeI] == eFaces[0])
{
neighbour_[edgeI] = eFaces[1];
}
else
{
neighbour_[edgeI] = eFaces[0];
}
}
else if (eFaces.size() == 1)
{
continue;
}
else
{
FatalErrorIn("polyMesh neighbour construction")
<< abort(FatalError);
}
}
}
Foam::labelList Foam::patchToPoly2DMesh::internalFaceOrder()
{
const labelListList& cellFaces = patch_.faceEdges();
labelList oldToNew(owner_.size(), -1);
label newFaceI = 0;
forAll(cellFaces, cellI)
{
const labelList& cFaces = cellFaces[cellI];
// Neighbouring cells
SortableList<label> nbr(cFaces.size(), -1);
forAll(cFaces, cfI)
{
if (cFaces[cfI] < neighbour_.size())
{
label nbrFaceI = neighbour_[cFaces[cfI]];
// Internal face. Get cell on other side.
if (nbrFaceI == cellI)
{
nbrFaceI = owner_[cellI];
}
if (cellI < nbrFaceI)
{
// CellI is master
nbr[cfI] = nbrFaceI;
}
}
}
nbr.sort();
forAll(nbr, i)
{
if (nbr[i] != -1)
{
oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++;
}
}
}
return oldToNew;
}
void Foam::patchToPoly2DMesh::addPatchFacesToFaces()
{
const labelList& meshPoints = patch_.meshPoints();
label offset = patch_.nInternalEdges();
face f(2);
forAll(patchNames_, patchI)
{
forAllConstIter(EdgeMap<label>, mapEdgesRegion_, eIter)
{
if (eIter() == patchI)
{
f[0] = meshPoints[eIter.key().start()];
f[1] = meshPoints[eIter.key().end()];
faces_[offset++] = f;
}
}
}
f.clear();
}
void Foam::patchToPoly2DMesh::addPatchFacesToOwner()
{
const label nInternalEdges = patch_.nInternalEdges();
const faceList& faces = patch_.faces();
const label nExternalEdges = patch_.edges().size() - nInternalEdges;
const labelList& meshPoints = patch_.meshPoints();
// Reorder patch faces on owner list.
labelList newOwner = owner_;
label nMatched = 0;
for
(
label bFaceI = nInternalEdges;
bFaceI < faces_.size();
++bFaceI
)
{
const face& e = faces_[bFaceI];
bool matched = false;
for
(
label bEdgeI = nInternalEdges;
bEdgeI < faces_.size();
++bEdgeI
)
{
if
(
e[0] == meshPoints[patch_.edges()[bEdgeI][0]]
&& e[1] == meshPoints[patch_.edges()[bEdgeI][1]]
)
{
const face& f = faces[owner_[bEdgeI]];
label fp = findIndex(f, e[0]);
newOwner[bFaceI] = owner_[bEdgeI];
if (f.nextLabel(fp) != e[1])
{
Info<< "Flipping" << endl;
faces_[bFaceI][0] = e[1];
faces_[bFaceI][1] = e[0];
}
nMatched++;
matched = true;
}
else if
(
e[0] == meshPoints[patch_.edges()[bEdgeI][1]]
&& e[1] == meshPoints[patch_.edges()[bEdgeI][0]]
)
{
Info<< "Warning: Wrong orientation." << endl;
nMatched++;
matched = true;
}
}
if (!matched)
{
Info<< "No match for edge." << endl;
}
}
if (nMatched != nExternalEdges)
{
Info<< "Number of matched edges, " << nMatched
<< ", does not match number of external edges, "
<< nExternalEdges << endl;
}
owner_ = newOwner.xfer();
}
void Foam::patchToPoly2DMesh::createPolyMeshComponents()
{
flipFaceOrder();
createNeighbours();
// New function for returning a map of old faces to new faces.
labelList oldToNew = internalFaceOrder();
inplaceReorder(oldToNew, faces_);
inplaceReorder(oldToNew, owner_);
inplaceReorder(oldToNew, neighbour_);
// Add patches.
addPatchFacesToFaces();
addPatchFacesToOwner();
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
//- Construct from a primitivePatch
Foam::patchToPoly2DMesh::patchToPoly2DMesh
(
const MeshedSurface<face>& patch,
const wordList& patchNames,
const labelList& patchSizes,
const EdgeMap<label>& mapEdgesRegion
)
:
patch_(patch),
patchNames_(patchNames),
patchSizes_(patchSizes),
patchStarts_(patchNames.size(), 0),
mapEdgesRegion_(mapEdgesRegion),
points_(patch.points()),
faces_(patch.nEdges()),
owner_(PatchTools::edgeOwner(patch)),
neighbour_(patch.nInternalEdges())
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::patchToPoly2DMesh::~patchToPoly2DMesh()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::patchToPoly2DMesh::createMesh()
{
createPolyMeshComponents();
label startFace = patch_.nInternalEdges();
forAll(patchNames_, patchI)
{
patchStarts_[patchI] = startFace;
startFace += patchSizes_[patchI];
}
}
// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -0,0 +1,165 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 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 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/>.
Class
Foam::patchToPoly2DMesh
Description
Convert a primitivePatch into a 2D polyMesh.
SourceFiles
patchToPoly2DMesh.C
\*---------------------------------------------------------------------------*/
#ifndef patchToPoly2DMesh_H
#define patchToPoly2DMesh_H
#include "EdgeMap.H"
#include "MeshedSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class patchToPoly2DMesh Declaration
\*---------------------------------------------------------------------------*/
class patchToPoly2DMesh
{
// Private data
const MeshedSurface<face>& patch_;
const wordList& patchNames_;
const labelList& patchSizes_;
labelList patchStarts_;
const EdgeMap<label>& mapEdgesRegion_;
pointField points_;
faceList faces_;
labelList owner_;
labelList neighbour_;
//- Description of data_
void flipFaceOrder();
void createNeighbours();
labelList internalFaceOrder();
void addPatchFacesToFaces();
void addPatchFacesToOwner();
void createPolyMeshComponents();
// Private Member Functions
//- Disallow default bitwise copy construct
patchToPoly2DMesh(const patchToPoly2DMesh&);
//- Disallow default bitwise assignment
void operator=(const patchToPoly2DMesh&);
public:
// Constructors
//- Construct from a primitivePatch
patchToPoly2DMesh
(
const MeshedSurface<face>& patch,
const wordList& patchNames,
const labelList& patchSizes,
const EdgeMap<label>& mapEdgesRegion
);
//- Destructor
~patchToPoly2DMesh();
// Member Functions
// Access
pointField& points()
{
return points_;
}
faceList& faces()
{
return faces_;
}
labelList& owner()
{
return owner_;
}
labelList& neighbour()
{
return neighbour_;
}
const wordList& patchNames() const
{
return patchNames_;
}
const labelList& patchSizes() const
{
return patchSizes_;
}
const labelList& patchStarts() const
{
return patchStarts_;
}
// Check
// Edit
void createMesh();
// Write
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -28,13 +28,6 @@ Description
Takes 2D mesh (all faces 2 points only, no front and back faces) and Takes 2D mesh (all faces 2 points only, no front and back faces) and
creates a 3D mesh by extruding with specified thickness. creates a 3D mesh by extruding with specified thickness.
Usage
- extrude2DMesh thickness
\param thickness \n
Thickness (in metre) of slab.
Note Note
Not sure about the walking of the faces to create the front and back faces. Not sure about the walking of the faces to create the front and back faces.
@ -43,155 +36,267 @@ Note
#include "argList.H" #include "argList.H"
#include "Time.H" #include "Time.H"
#include "polyMesh.H" #include "polyMesh.H"
#include "polyTopoChange.H"
#include "extrude2DMesh.H" #include "extrude2DMesh.H"
#include "emptyPolyPatch.H" #include "extrudeModel.H"
#include "polyTopoChange.H"
#include "MeshedSurface.H"
#include "edgeCollapser.H"
#include "addPatchCellLayer.H"
#include "patchToPoly2DMesh.H"
using namespace Foam; using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
enum ExtrudeMode
{
POLYMESH2D,
MESHEDSURFACE
};
namespace Foam
{
template<>
const char* NamedEnum<ExtrudeMode, 2>::names[] =
{
"polyMesh2D",
"MeshedSurface"
};
}
static const NamedEnum<ExtrudeMode, 2> ExtrudeModeNames;
//pointField moveInitialPoints
//(
// primitiveFacePatch& fMesh,
// const extrudeModel& model
//)
//{
// pointField layer0Points(fMesh.nPoints());
// pointField layer1Points(fMesh.nPoints());
// pointField displacement(fMesh.nPoints());
// forAll(layer0Points, pointI)
// {
// const labelList& meshPoints = fMesh.meshPoints();
// label meshPointI = meshPoints[pointI];
// layer0Points[meshPointI] = model
// (
// fMesh.points()[meshPointI],
// fMesh.pointNormals()[pointI],
// 0
// );
// layer1Points[meshPointI] = model
// (
// fMesh.points()[meshPointI],
// fMesh.pointNormals()[pointI],
// 1
// );
// displacement[pointI] =
// layer1Points[meshPointI]
// - layer0Points[meshPointI];
// }
// fMesh.movePoints(layer0Points);
// return displacement;
//}
// Main program: // Main program:
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
# include "addOverwriteOption.H" #include "addOverwriteOption.H"
argList::validArgs.append("thickness");
# include "setRootCase.H" argList::validArgs.append("surfaceFormat");
# include "createTime.H"
runTime.functionObjects().off(); #include "setRootCase.H"
# include "createPolyMesh.H"
const word oldInstance = mesh.pointsInstance(); Info<< "Create time\n" << endl;
Time runTimeExtruded
(
Time::controlDictName,
args.rootPath(),
args.caseName()
);
runTimeExtruded.functionObjects().off();
const scalar thickness = args.argRead<scalar>(1);
const bool overwrite = args.optionFound("overwrite"); const bool overwrite = args.optionFound("overwrite");
const ExtrudeMode surfaceFormat = ExtrudeModeNames[args[1]];
// Check that mesh is 2D Info<< "Extruding from " << ExtrudeModeNames[surfaceFormat] << endl;
// ~~~~~~~~~~~~~~~~~~~~~
const faceList& faces = mesh.faces(); IOdictionary extrude2DMeshDict
forAll(faces, faceI) (
IOobject
(
"extrude2DMeshDict",
runTimeExtruded.system(),
runTimeExtruded,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
);
// Point generator
autoPtr<extrudeModel> model(extrudeModel::New(extrude2DMeshDict));
autoPtr<MeshedSurface<face> > fMesh;
autoPtr<polyMesh> mesh;
autoPtr<polyTopoChange> meshMod;
labelListList extrudeEdgePatches;
if (surfaceFormat == MESHEDSURFACE)
{ {
if (faces[faceI].size() != 2) fMesh.set(new MeshedSurface<face>("MeshedSurface.obj"));
EdgeMap<label> edgeRegionMap;
wordList patchNames(1, "default");
labelList patchSizes(1, fMesh->nEdges() - fMesh->nInternalEdges());
const edgeList& edges = fMesh->edges();
forAll(edges, edgeI)
{ {
FatalErrorIn(args.executable()) if (!fMesh->isInternalEdge(edgeI))
<< "Face " << faceI << " size " << faces[faceI].size() {
<< " is not of size 2 so mesh is not proper two-dimensional." edgeRegionMap.insert(edges[edgeI], 0);
<< exit(FatalError);
} }
} }
patchToPoly2DMesh poly2DMesh
(
fMesh(),
patchNames,
patchSizes,
edgeRegionMap
);
// Find extrude direction poly2DMesh.createMesh();
// ~~~~~~~~~~~~~~~~~~~~~~
scalar minRange = GREAT; mesh.set
direction extrudeDir = 4; //illegal value. (
new polyMesh
(
IOobject
(
polyMesh::defaultRegion,
runTimeExtruded.constant(),
runTimeExtruded,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
xferMove(poly2DMesh.points()),
xferMove(poly2DMesh.faces()),
xferMove(poly2DMesh.owner()),
xferMove(poly2DMesh.neighbour())
)
);
for (direction dir = 0; dir < 3; dir++) Info<< "Constructing patches." << endl;
{ List<polyPatch*> patches(poly2DMesh.patchNames().size());
scalarField cmpts(mesh.points().component(dir));
scalar range = max(cmpts)-min(cmpts);
Info<< "Direction:" << dir << " range:" << range << endl;
if (range < minRange)
{
minRange = range;
extrudeDir = dir;
}
}
Info<< "Extruding in direction " << extrudeDir
<< " with thickness " << thickness << nl
<< endl;
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Add front and back patch
// ~~~~~~~~~~~~~~~~~~~~~~~~
label frontPatchI = patches.findPatchID("frontAndBack");
if (frontPatchI == -1)
{
// Add patch.
List<polyPatch*> newPatches(patches.size()+1);
forAll(patches, patchI) forAll(patches, patchI)
{ {
const polyPatch& pp = patches[patchI]; patches[patchI] = new polyPatch
newPatches[patchI] = pp.clone
( (
patches, poly2DMesh.patchNames()[patchI],
newPatches.size(), poly2DMesh.patchSizes()[patchI],
pp.size(), poly2DMesh.patchStarts()[patchI],
pp.start() patchI,
).ptr(); mesh->boundaryMesh()
}
frontPatchI = patches.size();
newPatches[frontPatchI] = new emptyPolyPatch
(
"frontAndBack",
0,
mesh.nFaces(),
frontPatchI,
patches
); );
Info<< "Adding empty patch " << newPatches[frontPatchI]->name()
<< " at index " << frontPatchI
<< " for front and back faces." << nl << endl;
mesh.removeBoundary();
mesh.addPatches(newPatches);
} }
mesh->addPatches(patches);
}
// Topo changes container. Initialise with number of patches. else if (surfaceFormat == POLYMESH2D)
polyTopoChange meshMod(mesh.boundaryMesh().size()); {
mesh.set
(
new polyMesh
(
Foam::IOobject
(
Foam::polyMesh::defaultRegion,
runTimeExtruded.timeName(),
runTimeExtruded,
Foam::IOobject::MUST_READ
)
)
);
}
// Engine to extrude mesh // Engine to extrude mesh
extrude2DMesh extruder(mesh); extrude2DMesh extruder(mesh(), extrude2DMeshDict, model());
// Insert changes into meshMod extruder.addFrontBackPatches();
extruder.setRefinement
( meshMod.set(new polyTopoChange(mesh().boundaryMesh().size()));
extrudeDir,
thickness, extruder.setRefinement(meshMod());
frontPatchI,
meshMod
);
// Create a mesh from topo changes. // Create a mesh from topo changes.
autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false); autoPtr<mapPolyMesh> morphMap = meshMod->changeMesh(mesh(), false);
mesh.updateMesh(morphMap); mesh->updateMesh(morphMap);
if (!overwrite)
{ {
runTime++; edgeCollapser collapser(mesh());
}
else const edgeList& edges = mesh->edges();
const pointField& points = mesh->points();
const boundBox& bb = mesh->bounds();
const scalar mergeDim = 1E-4 * bb.minDim();
forAll(edges, edgeI)
{ {
mesh.setInstance(oldInstance); const edge& e = edges[edgeI];
scalar d = e.mag(points);
if (d < mergeDim)
{
Info<< "Merging edge " << e << " since length " << d
<< " << " << mergeDim << nl;
// Collapse edge to e[0]
collapser.collapseEdge(edgeI, e[0]);
} }
}
polyTopoChange meshModCollapse(mesh());
collapser.setRefinement(meshModCollapse);
// Create a mesh from topo changes.
autoPtr<mapPolyMesh> morphMap
= meshModCollapse.changeMesh(mesh(), false);
mesh->updateMesh(morphMap);
}
mesh->setInstance(runTimeExtruded.constant());
// Take over refinement levels and write to new time directory. // Take over refinement levels and write to new time directory.
Pout<< "Writing extruded mesh to time " << runTime.timeName() << nl Pout<< "Writing extruded mesh to time = " << runTimeExtruded.timeName()
<< endl; << nl << endl;
mesh.write(); mesh().write();
Pout<< "End\n" << endl; Pout<< "End\n" << endl;

View File

@ -0,0 +1,46 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
root "";
case "";
instance "";
local "";
class dictionary;
object extrude2DMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
extrudeModel linearDirection;
//extrudeModel wedge;
patchType empty;
//patchType wedge;
nLayers 1;
expansionRatio 1.0;
linearDirectionCoeffs
{
direction (0 0 1);
thickness 0.1;
}
wedgeCoeffs
{
axisPt (0 0 0);
axis (1 0 0);
angle 10;
}

View File

@ -49,6 +49,7 @@ Note
#include "triSurface.H" #include "triSurface.H"
#include "OFstream.H" #include "OFstream.H"
#include "OSspecific.H" #include "OSspecific.H"
#include "Time.H"
using namespace Foam; using namespace Foam;
@ -83,9 +84,25 @@ int main(int argc, char *argv[])
"factor", "factor",
"geometry scaling factor - default is 1" "geometry scaling factor - default is 1"
); );
argList::addOption
(
"writePrecision",
"label",
"write to output with the specified precision"
);
argList args(argc, argv); argList args(argc, argv);
if (args.optionFound("writePrecision"))
{
label writePrecision = args.optionRead<label>("writePrecision");
IOstream::defaultPrecision(writePrecision);
Sout.precision(writePrecision);
Info<< "Output write precision set to " << writePrecision << endl;
}
const fileName importName = args[1]; const fileName importName = args[1];
const fileName exportName = args[2]; const fileName exportName = args[2];

View File

@ -44,6 +44,7 @@ Description
#include "indexedOctree.H" #include "indexedOctree.H"
#include "treeDataEdge.H" #include "treeDataEdge.H"
#include "unitConversion.H" #include "unitConversion.H"
#include "plane.H"
using namespace Foam; using namespace Foam;
@ -94,6 +95,56 @@ void deleteBox
} }
bool onLine(const point& p, const linePointRef& line)
{
const point& a = line.start();
const point& b = line.end();
if
(
( p.x() < min(a.x(), b.x()) || p.x() > max(a.x(), b.x()) )
|| ( p.y() < min(a.y(), b.y()) || p.y() > max(a.y(), b.y()) )
|| ( p.z() < min(a.z(), b.z()) || p.z() > max(a.z(), b.z()) )
)
{
return false;
}
return true;
}
// Deletes all edges inside/outside bounding box from set.
void deleteEdges
(
const triSurface& surf,
const plane& cutPlane,
List<surfaceFeatures::edgeStatus>& edgeStat
)
{
const pointField& points = surf.points();
const labelList& meshPoints = surf.meshPoints();
forAll(edgeStat, edgeI)
{
const edge& e = surf.edges()[edgeI];
const point& p0 = points[meshPoints[e.start()]];
const point& p1 = points[meshPoints[e.end()]];
const linePointRef line(p0, p1);
// If edge does not intersect the plane, delete.
scalar intersect = cutPlane.lineIntersect(line);
point featPoint = intersect * (p1 - p0) + p0;
if (!onLine(featPoint, line))
{
edgeStat[edgeI] = surfaceFeatures::NONE;
}
}
}
void drawHitProblem void drawHitProblem
( (
label fI, label fI,
@ -270,6 +321,12 @@ int main(int argc, char *argv[])
"manifoldEdgesOnly", "manifoldEdgesOnly",
"remove any non-manifold (open or more than two connected faces) edges" "remove any non-manifold (open or more than two connected faces) edges"
); );
argList::addOption
(
"plane",
"(nx ny nz)(z0 y0 z0)",
"use a plane to create feature edges for 2D meshing"
);
# ifdef ENABLE_CURVATURE # ifdef ENABLE_CURVATURE
argList::addBoolOption argList::addBoolOption
@ -454,6 +511,21 @@ int main(int argc, char *argv[])
} }
} }
if (args.optionFound("plane"))
{
plane cutPlane(args.optionLookup("plane")());
deleteEdges
(
surf,
cutPlane,
edgeStat
);
Info<< "Only edges that intersect the plane with normal "
<< cutPlane.normal() << " and base point " << cutPlane.refPoint()
<< " will be included as feature edges."<< endl;
}
surfaceFeatures newSet(surf); surfaceFeatures newSet(surf);
newSet.setFromStatus(edgeStat); newSet.setFromStatus(edgeStat);
@ -475,7 +547,6 @@ int main(int argc, char *argv[])
<< endl; << endl;
// Extracting and writing a extendedFeatureEdgeMesh // Extracting and writing a extendedFeatureEdgeMesh
extendedFeatureEdgeMesh feMesh extendedFeatureEdgeMesh feMesh
( (
newSet, newSet,

View File

@ -9,4 +9,5 @@ CGAL_LIBS = \
-L$(MPFR_ARCH_PATH)/lib \ -L$(MPFR_ARCH_PATH)/lib \
-L$(BOOST_ARCH_PATH)/lib \ -L$(BOOST_ARCH_PATH)/lib \
-lmpfr \ -lmpfr \
-lboost_thread -I/usr/lib64 \
-lboost_thread-mt