From 256f01284fae7954e9a726bb47a2a01937bdb817 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Wed, 8 Aug 2012 10:31:19 +0000
Subject: [PATCH] fix pyramids
---
Geo/GRegion.cpp | 51 ++++++++++++-------------------
Geo/GRegion.h | 8 ++---
Geo/OCCRegion.cpp | 13 ++++----
Geo/gmshRegion.cpp | 3 +-
Mesh/meshGRegion.cpp | 73 +++++++++++++++++++++-----------------------
tutorial/t5.geo | 18 +++++------
6 files changed, 75 insertions(+), 91 deletions(-)
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index 0709934519..ba668a3a2e 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -31,17 +31,6 @@ GRegion::~GRegion()
deleteMesh();
}
-void GRegion::set(const std::list<GFace*> f) {
- for (std::list<GFace*>::iterator it = l_faces.begin(); it != l_faces.end() ; ++it){
- (*it)->delRegion(this);
- }
- l_faces = f;
- for (std::list<GFace*>::iterator it = l_faces.begin(); it != l_faces.end() ; ++it){
- (*it)->addRegion(this);
- }
-}
-
-
void GRegion::deleteMesh()
{
for(unsigned int i = 0; i < mesh_vertices.size(); i++) delete mesh_vertices[i];
@@ -62,7 +51,7 @@ void GRegion::deleteMesh()
}
unsigned int GRegion::getNumMeshElements()
-{
+{
return tetrahedra.size() + hexahedra.size() + prisms.size() + pyramids.size() +
polyhedra.size();
}
@@ -108,19 +97,19 @@ MElement *const *GRegion::getStartElementType(int type) const
}
MElement *GRegion::getMeshElement(unsigned int index)
-{
+{
if(index < tetrahedra.size())
return tetrahedra[index];
else if(index < tetrahedra.size() + hexahedra.size())
return hexahedra[index - tetrahedra.size()];
else if(index < tetrahedra.size() + hexahedra.size() + prisms.size())
return prisms[index - tetrahedra.size() - hexahedra.size()];
- else if(index < tetrahedra.size() + hexahedra.size() + prisms.size() +
+ else if(index < tetrahedra.size() + hexahedra.size() + prisms.size() +
pyramids.size())
return pyramids[index - tetrahedra.size() - hexahedra.size() - prisms.size()];
- else if(index < tetrahedra.size() + hexahedra.size() + prisms.size() +
+ else if(index < tetrahedra.size() + hexahedra.size() + prisms.size() +
pyramids.size() + polyhedra.size())
- return polyhedra[index - tetrahedra.size() - hexahedra.size() - prisms.size() -
+ return polyhedra[index - tetrahedra.size() - hexahedra.size() - prisms.size() -
pyramids.size()];
return 0;
}
@@ -152,7 +141,7 @@ SOrientedBoundingBox GRegion::getOBB()
if (!_obb) {
std::vector<SPoint3> vertices;
std::list<GFace*> b_faces = faces();
- for (std::list<GFace*>::iterator b_face = b_faces.begin();
+ for (std::list<GFace*>::iterator b_face = b_faces.begin();
b_face != b_faces.end(); b_face++) {
if((*b_face)->getNumMeshVertices() > 0) {
int N = (*b_face)->getNumMeshVertices();
@@ -177,17 +166,17 @@ SOrientedBoundingBox GRegion::getOBB()
vertices.push_back(pt1);
vertices.push_back(pt2);
}
- }
+ }
else if ((*b_face)->buildSTLTriangulation()) {
for (unsigned int i = 0; i < (*b_face)->stl_vertices.size(); i++){
GPoint p = (*b_face)->point((*b_face)->stl_vertices[i]);
vertices.push_back(SPoint3(p.x(), p.y(), p.z()));
- }
+ }
}
else {
int N = 10;
std::list<GEdge*> b_edges = (*b_face)->edges();
- for (std::list<GEdge*>::iterator b_edge = b_edges.begin();
+ for (std::list<GEdge*>::iterator b_edge = b_edges.begin();
b_edge != b_edges.end(); b_edge++) {
Range<double> tr = (*b_edge)->parBounds(0);
for (int j = 0; j < N; j++) {
@@ -196,7 +185,7 @@ SOrientedBoundingBox GRegion::getOBB()
SPoint3 pt(p.x(), p.y(), p.z());
vertices.push_back(pt);
}
- }
+ }
}
}
_obb = SOrientedBoundingBox::buildOBB(vertices);
@@ -312,7 +301,7 @@ void GRegion::replaceFaces (std::list<GFace*> &new_faces)
std::list<int> newdirs;
for ( ; it != l_faces.end(); ++it, ++it2, ++it3){
(*it)->delRegion(this);
- (*it2)->addRegion(this);
+ (*it2)->addRegion(this);
// if ((*it2)->getBeginVertex() == (*it)->getBeginVertex())
newdirs.push_back(*it3);
// else
@@ -332,13 +321,13 @@ double GRegion::computeSolidProperties (std::vector<double> cg,
double volumez = 0;
double surface = 0;
cg[0] = cg[1] = cg[2] = 0.0;
- for ( ; it != l_faces.end(); ++it,++itdir){
+ for ( ; it != l_faces.end(); ++it,++itdir){
//printf("face %d dir %d %d elements\n",(*it)->tag(),*itdir,(int)(*it)->triangles.size());
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *e = (*it)->triangles[i];
int npt;
IntPt *pts;
- e->getIntegrationPoints (2*(e->getPolynomialOrder()-1)+3, &npt, &pts);
+ e->getIntegrationPoints (2*(e->getPolynomialOrder()-1)+3, &npt, &pts);
for (int j=0;j<npt;j++){
SPoint3 pt;
// compute x,y,z of the integration point
@@ -370,19 +359,19 @@ double GRegion::computeSolidProperties (std::vector<double> cg,
it = l_faces.begin();
itdir = l_dirs.begin();
- inertia[0] =
- inertia[1] =
- inertia[2] =
- inertia[3] =
- inertia[4] =
+ inertia[0] =
+ inertia[1] =
+ inertia[2] =
+ inertia[3] =
+ inertia[4] =
inertia[5] = 0.0;
- for ( ; it != l_faces.end(); ++it,++itdir){
+ for ( ; it != l_faces.end(); ++it,++itdir){
for (unsigned int i = 0; i < (*it)->getNumMeshElements(); ++i){
MElement *e = (*it)->getMeshElement(i);
int npt;
IntPt *pts;
- e->getIntegrationPoints (2*(e->getPolynomialOrder()-1)+3, &npt, &pts);
+ e->getIntegrationPoints (2*(e->getPolynomialOrder()-1)+3, &npt, &pts);
for (int j = 0; j < npt; j++){
SPoint3 pt;
// compute x,y,z of the integration point
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index e6a6733ee6..7f87c3a3fe 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -26,8 +26,8 @@ class GRegion : public GEntity {
protected:
std::list<GFace*> l_faces;
std::list<int> l_dirs;
- GRegionCompound *compound; // this model ede belongs to a compound
-
+ GRegionCompound *compound; // this model ede belongs to a compound
+
// replace faces (for gluing) for specific modelers, we have to
// re-create internal data ...
virtual void replaceFacesInternal (std::list<GFace*> &) {}
@@ -48,7 +48,7 @@ class GRegion : public GEntity {
// get/set faces that bound the region
virtual std::list<GFace*> faces() const{ return l_faces; }
virtual std::list<int> faceOrientations() const{ return l_dirs; }
- void set(const std::list<GFace*> f);
+ inline void set(const std::list<GFace*> f) { l_faces = f; }
// edges that bound the region
virtual std::list<GEdge*> edges() const;
@@ -91,7 +91,7 @@ class GRegion : public GEntity {
// reset the mesh attributes to default values
virtual void resetMeshAttributes();
-
+
// compound
void setCompound(GRegionCompound *grc) { compound = grc; }
GRegionCompound *getCompound() const { return compound; }
diff --git a/Geo/OCCRegion.cpp b/Geo/OCCRegion.cpp
index 9d65e62005..e341df9d35 100644
--- a/Geo/OCCRegion.cpp
+++ b/Geo/OCCRegion.cpp
@@ -26,7 +26,7 @@ void OCCRegion::setup()
for(exp2.Init(s, TopAbs_SHELL); exp2.More(); exp2.Next()){
TopoDS_Shape shell = exp2.Current();
Msg::Debug("OCC Region %d - New Shell",tag());
- for(exp3.Init(shell, TopAbs_FACE); exp3.More(); exp3.Next()){
+ for(exp3.Init(shell, TopAbs_FACE); exp3.More(); exp3.Next()){
TopoDS_Face face = TopoDS::Face(exp3.Current());
GFace *f = getOCCFaceByNativePtr(model(),face);
if(f){
@@ -35,7 +35,7 @@ void OCCRegion::setup()
}
else
Msg::Error("Unknown face in region %d", tag());
- }
+ }
}
Msg::Debug("OCC Region %d with %d faces", tag(), l_faces.size());
}
@@ -45,7 +45,6 @@ GEntity::GeomType OCCRegion::geomType() const
return Unknown;
}
-
bool FaceHaveDifferentOrientations(const TopoDS_Face& aFR,
const TopoDS_Face& aF)
{
@@ -82,19 +81,19 @@ void OCCRegion::replaceFacesInternal(std::list<GFace*> &new_faces)
if (occF){
TopoDS_Face oldf = occF->getTopoDS_Face();
if (oldf.IsSame(_face)){
- _face_replacement = *((TopoDS_Face*)(*it2)->getNativePtr());
+ _face_replacement = *((TopoDS_Face*)(*it2)->getNativePtr());
}
else {
oldf = occF->getTopoDS_FaceOld();
if (oldf.IsSame(_face)){
- _face_replacement = *((TopoDS_Face*)(*it2)->getNativePtr());
+ _face_replacement = *((TopoDS_Face*)(*it2)->getNativePtr());
}
}
}
}
if (_face_replacement.IsNull()){
Msg::Error("cannot find an face for gluing a region");
- }
+ }
if (_face_replacement.IsSame(_face)) {
aBB.Add(_shell_replacement, _face);
@@ -104,7 +103,7 @@ void OCCRegion::replaceFacesInternal(std::list<GFace*> &new_faces)
_face_replacement.Reverse();
aBB.Add(_shell_replacement, _face_replacement);
}
- }
+ }
aBB.Add(_s_replacement, _shell_replacement);
}
s = _s_replacement;
diff --git a/Geo/gmshRegion.cpp b/Geo/gmshRegion.cpp
index 5e0c9889a1..8ce695cebb 100644
--- a/Geo/gmshRegion.cpp
+++ b/Geo/gmshRegion.cpp
@@ -21,6 +21,7 @@ gmshRegion::gmshRegion(GModel *m, ::Volume *volume)
if(f){
l_faces.push_back(f);
l_dirs.push_back(ori);
+ f->addRegion(this);
}
else
Msg::Error("Unknown surface %d", s->Num);
@@ -32,11 +33,11 @@ gmshRegion::gmshRegion(GModel *m, ::Volume *volume)
if(f){
l_faces.push_back(f);
l_dirs.push_back(sign(is));
+ f->addRegion(this);
}
else
Msg::Error("Unknown surface %d", is);
}
-
resetMeshAttributes();
}
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index c75fe74f47..09e5c7295f 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -76,11 +76,10 @@ class splitQuadRecovery {
allFaces.erase(*itf1);
allFaces.erase(*itf2);
allFaces.erase(*itf3);
- printf("some pyramids should be created %d regions\n",(*it)->numRegions());
+ // printf("some pyramids should be created %d regions\n", (*it)->numRegions());
for (int iReg = 0; iReg < (*it)->numRegions(); iReg++){
- printf("one pyramid is created\n");
if (iReg == 1) {
- Msg::Error("cannot build pyramids on non manifold faces ...");
+ Msg::Error("Cannot build pyramids on non manifold faces");
v = new MVertex(v->x(), v->y(), v->z(), (*it)->getRegion(iReg));
}
else
@@ -358,7 +357,7 @@ void getBoundingInfoAndSplitQuads(GRegion *gr,
std::map<MFace, GEntity*, Less_Face>::iterator itx = allBoundingFaces_temp.begin();
for (; itx != allBoundingFaces_temp.end();++itx){
const MFace &f = itx->first;
- // SPLIT that face into 4 triangular faces
+ // split the quad face into 4 triangular faces
if (f.getNumVertices() == 4){
sqr.setEmpty(false);
MVertex *v0 = f.getVertex(0);
@@ -374,7 +373,6 @@ void getBoundingInfoAndSplitQuads(GRegion *gr,
allBoundingFaces[MFace(v2,v3,newv)] = itx->second;
allBoundingFaces[MFace(v3,v0,newv)] = itx->second;
itx->second->mesh_vertices.push_back(newv);
- // myGr->pyramids.push_back(new MPyramid(v0,v1,v2,v3,newv));
allBoundingVertices.insert(v0);
allBoundingVertices.insert(v1);
allBoundingVertices.insert(v2);
@@ -390,26 +388,10 @@ void getBoundingInfoAndSplitQuads(GRegion *gr,
}
}
-void getAllBoundingVertices(GRegion *gr, std::set<MVertex*> &allBoundingVertices)
-{
- std::list<GFace*> faces = gr->faces();
- std::list<GFace*>::iterator it = faces.begin();
-
- while (it != faces.end()){
- GFace *gf = (*it);
- for(unsigned int i = 0; i < gf->triangles.size(); i++){
- MTriangle *t = gf->triangles[i];
- for(int k = 0; k < 3; k++)
- if(allBoundingVertices.find(t->getVertex(k)) == allBoundingVertices.end())
- allBoundingVertices.insert(t->getVertex(k));
- }
- ++it;
- }
-}
-
#if defined(HAVE_TETGEN)
#include "tetgen.h"
-void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numberedV, splitQuadRecovery &sqr)
+void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numberedV,
+ splitQuadRecovery &sqr)
{
std::set<MVertex*> allBoundingVertices;
std::map<MFace,GEntity*,Less_Face> allBoundingFaces;
@@ -479,10 +461,10 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb
void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
std::vector<MVertex*> &numberedV)
{
- // improvement has to be done here : tetgen splits some of the
- // existing edges of the mesh. If those edges are classified on some
- // model faces, new points SHOULD be classified on the model face
- // and get the right set of parametric coordinates.
+ // improvement has to be done here : tetgen splits some of the existing edges
+ // of the mesh. If those edges are classified on some model faces, new points
+ // SHOULD be classified on the model face and get the right set of parametric
+ // coordinates.
for(int i = numberedV.size(); i < out.numberofpoints; i++){
MVertex *v = new MVertex(out.pointlist[i * 3 + 0],
@@ -520,9 +502,8 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
bool needParam = (CTX::instance()->mesh.order > 1 &&
CTX::instance()->mesh.secondOrderExperimental);
- // re-create the triangular meshes FIXME: this can lead to hanging
- // nodes for non manifold geometries (single surface connected to
- // volume)
+ // re-create the triangular meshes FIXME: this can lead to hanging nodes for
+ // non manifold geometries (single surface connected to volume)
for(int i = 0; i < out.numberoftrifaces; i++){
if (out.trifacemarkerlist[i] > 0) {
MVertex *v[3];
@@ -580,9 +561,9 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
v[j]));
MVertex *v1b;
if(needParam){
- // parametric coordinates should be set for the vertex !
- // (this is not 100 % safe yet, so we reserve that operation
- // for high order meshes)
+ // parametric coordinates should be set for the vertex ! (this is
+ // not 100 % safe yet, so we reserve that operation for high order
+ // meshes)
GPoint gp = gf->closestPoint(SPoint3(v[j]->x(), v[j]->y(), v[j]->z()), guess);
if(gp.g()){
v1b = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
@@ -650,8 +631,6 @@ void _relocateVertex(MVertex *ver,
}
}
-
-
void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
{
if(regions.empty()) return;
@@ -686,8 +665,6 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
std::vector<MVertex*> numberedV;
char opts[128];
buildTetgenStructure(gr, in, numberedV, sqr);
- // JFR REMOVED THE -q OPTION BECAUSE MESH SIZES AT VERTICES WERE WRONG ...
- // COULD BE FIXED IF NEEDED
if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_DEL ||
CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_HEX ||
CTX::instance()->mesh.algo3d == ALGO_3D_MMG3D ||
@@ -695,8 +672,9 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
CTX::instance()->mesh.algo2d == ALGO_2D_BAMG){
sprintf(opts, "-pY%c", (Msg::GetVerbosity() < 3) ? 'Q':
(Msg::GetVerbosity() > 6) ? 'V': '\0');
+ // removed -q because mesh sizes at vertices were wrong...
// sprintf(opts, "-q1.5pY%c", (Msg::GetVerbosity() < 3) ? 'Q':
- // (Msg::GetVerbosity() > 6) ? 'V': '\0');
+ // (Msg::GetVerbosity() > 6) ? 'V': '\0');
}
else {
sprintf(opts, "-Ype%c", (Msg::GetVerbosity() < 3) ? 'Q':
@@ -751,7 +729,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
// restore the initial set of faces
gr->set(faces);
- //EMI VORONOI FOR CENRTERLINE
+ // EMI Voronoi for centerlines
if (Msg::GetVerbosity() == 20) {
std::set<SPoint3> candidates;
printVoronoi(gr, candidates);
@@ -796,6 +774,23 @@ namespace nglib {
}
using namespace nglib;
+static void getAllBoundingVertices(GRegion *gr, std::set<MVertex*> &allBoundingVertices)
+{
+ std::list<GFace*> faces = gr->faces();
+ std::list<GFace*>::iterator it = faces.begin();
+
+ while (it != faces.end()){
+ GFace *gf = (*it);
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ MTriangle *t = gf->triangles[i];
+ for(int k = 0; k < 3; k++)
+ if(allBoundingVertices.find(t->getVertex(k)) == allBoundingVertices.end())
+ allBoundingVertices.insert(t->getVertex(k));
+ }
+ ++it;
+ }
+}
+
Ng_Mesh *buildNetgenStructure(GRegion *gr, bool importVolumeMesh,
std::vector<MVertex*> &numberedV)
{
diff --git a/tutorial/t5.geo b/tutorial/t5.geo
index 635b159908..d10f08ae06 100644
--- a/tutorial/t5.geo
+++ b/tutorial/t5.geo
@@ -1,7 +1,7 @@
-/*********************************************************************
+/*********************************************************************
*
* Gmsh tutorial 5
- *
+ *
* Characteristic lengths, arrays of variables, functions, loops
*
*********************************************************************/
@@ -33,7 +33,7 @@ lcar3 = .055;
// truncated cube:
Point(1) = {0.5,0.5,0.5,lcar2}; Point(2) = {0.5,0.5,0,lcar1};
-Point(3) = {0,0.5,0.5,lcar1}; Point(4) = {0,0,0.5,lcar1};
+Point(3) = {0,0.5,0.5,lcar1}; Point(4) = {0,0,0.5,lcar1};
Point(5) = {0.5,0,0.5,lcar1}; Point(6) = {0.5,0,0,lcar1};
Point(7) = {0,0.5,0,lcar1}; Point(8) = {0,1,0,lcar1};
Point(9) = {1,1,0,lcar1}; Point(10) = {0,0,1,lcar1};
@@ -61,7 +61,7 @@ Line Loop(38) = {-14,-13,-12,19}; Plane Surface(39) = {38};
// Instead of using included files, we now use a user-defined function
// in order to carve some holes in the cube:
-Function CheeseHole
+Function CheeseHole
// In the following commands we use the reserved variable name
// `newp', which automatically selects a new point number. This
@@ -102,11 +102,11 @@ Function CheeseHole
// for later reference (we will need these to define the final
// volume):
- theloops[t] = newreg ;
+ theloops[t] = newreg ;
Surface Loop(theloops[t]) = {l8+1,l5+1,l1+1,l2+1,l3+1,l7+1,l6+1,l4+1};
- thehole = newreg ;
+ thehole = newreg ;
Volume(thehole) = theloops[t] ;
Return
@@ -117,8 +117,8 @@ x = 0 ; y = 0.75 ; z = 0 ; r = 0.09 ;
For t In {1:5}
- x += 0.166 ;
- z += 0.166 ;
+ x += 0.166 ;
+ z += 0.166 ;
// We call the `CheeseHole' function:
@@ -127,7 +127,7 @@ For t In {1:5}
// We define a physical volume for each hole:
Physical Volume (t) = thehole ;
-
+
// We also print some variables on the terminal (note that, since
// all variables are treated internally as floating point numbers,
// the format string should only contain valid floating point format
--
GitLab