diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp index 65c23277a2ade26b2113d4e78bb9d81f0f8bcff1..12892ffab14d97acf5682fc4e50035ed5442b4be 100644 --- a/Fltk/Callbacks.cpp +++ b/Fltk/Callbacks.cpp @@ -1,4 +1,4 @@ -// $Id: Callbacks.cpp,v 1.480 2006-11-25 23:29:26 geuzaine Exp $ +// $Id: Callbacks.cpp,v 1.481 2006-11-26 01:03:17 geuzaine Exp $ // // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle // @@ -3805,13 +3805,7 @@ void mesh_degree_cb(CALLBACK_ARGS) void mesh_optimize_cb(CALLBACK_ARGS) { - if(CTX.threads_lock) { - Msg(INFO, "I'm busy! Ask me that later..."); - return; - } - CTX.threads_lock = 1; OptimizeMesh(); - CTX.threads_lock = 0; CTX.mesh.changed = ENT_LINE | ENT_SURFACE | ENT_VOLUME; Draw(); Msg(STATUS2N, " "); diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp index 35075149bb73cf3310c206ecbbef707105d080bc..5042c59421dbcad7a0a4e46b84b79af8a348ccdc 100644 --- a/Mesh/meshGFace.cpp +++ b/Mesh/meshGFace.cpp @@ -1,4 +1,4 @@ -// $Id: meshGFace.cpp,v 1.30 2006-11-25 20:37:41 geuzaine Exp $ +// $Id: meshGFace.cpp,v 1.31 2006-11-26 01:03:17 geuzaine Exp $ // // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle // @@ -27,6 +27,7 @@ #include "GFace.h" #include "MVertex.h" #include "MElement.h" +#include "MRep.h" #include "Context.h" #include "GPoint.h" #include "Message.h" @@ -1358,9 +1359,9 @@ void deMeshGFace :: operator() (GFace *gf) gf->triangles.clear(); for (unsigned int i=0;i<gf->quadrangles.size();i++) delete gf->quadrangles[i]; gf->quadrangles.clear(); + if(gf->meshRep) gf->meshRep->destroy(); } - void meshGFace :: operator() (GFace *gf) { if(gf->geomType() == GEntity::DiscreteSurface) return; diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp index 942f4fa15c2937ec222caf0dddcf6973d62865f3..c93ce67dc34ee3644604966e0e25ca18ac3084be 100644 --- a/Mesh/meshGRegion.cpp +++ b/Mesh/meshGRegion.cpp @@ -1,176 +1,174 @@ +// $Id: meshGRegion.cpp,v 1.12 2006-11-26 01:03:17 geuzaine Exp $ +// +// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program 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 this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA. +// +// Please report all bugs and problems to <gmsh@geuz.org>. + #include "meshGRegion.h" #include "GModel.h" #include "GRegion.h" #include "GFace.h" #include "GEdge.h" +#include "MRep.h" #include "BDS.h" #include "Message.h" #include <vector> - - - -void getAllBoundingVertices (GRegion *gr, int importVolumeMesh, std::set<MVertex*> &allBoundingVertices ) +void getAllBoundingVertices(GRegion *gr, std::set<MVertex*> &allBoundingVertices) { std::list<GFace*> faces = gr->faces(); std::list<GFace*>::iterator it = faces.begin(); - - if (importVolumeMesh) - allBoundingVertices.insert ( gr->mesh_vertices.begin() , gr->mesh_vertices.end() ); - 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; + 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) -{ +void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numberedV) +{ std::set<MVertex*> allBoundingVertices; - getAllBoundingVertices (gr, 0, allBoundingVertices ); + getAllBoundingVertices(gr, allBoundingVertices); in.mesh_dim = 3; in.firstnumber = 1; - in.numberofpoints = allBoundingVertices.size(); in.pointlist = new REAL[in.numberofpoints * 3]; in.pointmarkerlist = NULL; std::set<MVertex*>::iterator itv = allBoundingVertices.begin(); - int I=1; - while (itv != allBoundingVertices.end()) - { - in.pointlist[(I-1)*3 + 0] = (*itv)->x(); - in.pointlist[(I-1)*3 + 1] = (*itv)->y(); - in.pointlist[(I-1)*3 + 2] = (*itv)->z(); - (*itv)->setNum(I++); - numberedV.push_back(*itv); - ++itv; - } + int I = 1; + while(itv != allBoundingVertices.end()){ + in.pointlist[(I-1)*3 + 0] = (*itv)->x(); + in.pointlist[(I-1)*3 + 1] = (*itv)->y(); + in.pointlist[(I-1)*3 + 2] = (*itv)->z(); + (*itv)->setNum(I++); + numberedV.push_back(*itv); + ++itv; + } int nbFace = 0; - std::list<GFace*> faces = gr->faces(); std::list<GFace*>::iterator it = faces.begin(); - while (it != faces.end()) - { - GFace *gf = (*it); - nbFace += gf->triangles.size(); - ++it; - } + while(it != faces.end()){ + GFace *gf = (*it); + nbFace += gf->triangles.size(); + ++it; + } in.numberoffacets = nbFace; in.facetlist = new tetgenio::facet[in.numberoffacets]; in.facetmarkerlist = new int[in.numberoffacets]; it = faces.begin(); - I= 0; - while (it != faces.end()) - { - GFace *gf = (*it); - for (unsigned int i = 0; i< gf->triangles.size(); i++) - { - MTriangle *t = gf->triangles[i]; - tetgenio::facet *f = &in.facetlist[I]; - tetgenio::init(f); - f->numberofholes = 0; - f->numberofpolygons = 1; - f->polygonlist = new tetgenio::polygon[f->numberofpolygons]; - tetgenio::polygon *p = &f->polygonlist[0]; - tetgenio::init(p); - p->numberofvertices = 3; - p->vertexlist = new int[p->numberofvertices]; - p->vertexlist[0] = t->getVertex(0)->getNum(); - p->vertexlist[1] = t->getVertex(1)->getNum(); - p->vertexlist[2] = t->getVertex(2)->getNum(); - // Msg(INFO,"Faces %d = %d %d %d",I,p->vertexlist[0] ,p->vertexlist[1] ,p->vertexlist[2] ); - in.facetmarkerlist[I] = gf->tag(); - ++I; - } - ++it; - } - - + I = 0; + while(it != faces.end()){ + GFace *gf = (*it); + for(unsigned int i = 0; i < gf->triangles.size(); i++){ + MTriangle *t = gf->triangles[i]; + tetgenio::facet *f = &in.facetlist[I]; + tetgenio::init(f); + f->numberofholes = 0; + f->numberofpolygons = 1; + f->polygonlist = new tetgenio::polygon[f->numberofpolygons]; + tetgenio::polygon *p = &f->polygonlist[0]; + tetgenio::init(p); + p->numberofvertices = 3; + p->vertexlist = new int[p->numberofvertices]; + p->vertexlist[0] = t->getVertex(0)->getNum(); + p->vertexlist[1] = t->getVertex(1)->getNum(); + p->vertexlist[2] = t->getVertex(2)->getNum(); + in.facetmarkerlist[I] = gf->tag(); + ++I; + } + ++it; + } } -void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, std::vector<MVertex*> & numberedV) +void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, + std::vector<MVertex*> &numberedV) { - for (int i = numberedV.size(); i < out.numberofpoints; i++) - { - MVertex *v = new MVertex (out.pointlist[i * 3 + 0], - out.pointlist[i * 3 + 1], - out.pointlist[i * 3 + 2], gr); - gr->mesh_vertices.push_back(v); - numberedV.push_back(v); - } + for(int i = numberedV.size(); i < out.numberofpoints; i++){ + MVertex *v = new MVertex(out.pointlist[i * 3 + 0], + out.pointlist[i * 3 + 1], + out.pointlist[i * 3 + 2], gr); + gr->mesh_vertices.push_back(v); + numberedV.push_back(v); + } - Msg(INFO,"%d points %d edges and %d faces in the final mesh",out.numberofpoints,out.numberofedges,out.numberoftrifaces); - // Tetgen modifies both surface & edge mesh. So, we recreate all mes + Msg(INFO,"%d points %d edges and %d faces in the final mesh", + out.numberofpoints,out.numberofedges,out.numberoftrifaces); - { - 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++) - { - delete gf->triangles[i]; - } - gf->triangles.clear(); - ++it; - } + // Tetgen modifies both surface & edge mesh. So, we need to + // re-create everything + + 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++) + delete gf->triangles[i]; + gf->triangles.clear(); + ++it; } - // re create 1D mesh - // 2 be done ... - for (int i = 0; i < out.numberofedges; i++) - { - MVertex *v[2]; - v[0] = numberedV[out.edgelist[i * 2 + 0] -1]; - v[1] = numberedV[out.edgelist[i * 2 + 1] -1]; - } + // TODO: re-create 1D mesh + for(int i = 0; i < out.numberofedges; i++){ + MVertex *v[2]; + v[0] = numberedV[out.edgelist[i * 2 + 0] -1]; + v[1] = numberedV[out.edgelist[i * 2 + 1] -1]; + } // re-create the triangular meshes - for (int i = 0; i < out.numberoftrifaces; i++) - { - MVertex *v[3]; - v[0] = numberedV[out.trifacelist[i * 3 + 0] -1]; - v[1] = numberedV[out.trifacelist[i * 3 + 1] -1]; - v[2] = numberedV[out.trifacelist[i * 3 + 2] -1]; - GFace *gf = gr->model()->faceByTag ( out.trifacemarkerlist[i] ); - for (int j=0;j<3;j++) - { - if ( v[j]->onWhat()->dim() == 3 ) - { - v[j]->onWhat()->mesh_vertices.erase(std::find(v[j]->onWhat()->mesh_vertices.begin(), - v[j]->onWhat()->mesh_vertices.end(), - v[j])); - // SPoint2 pp = gf->parFromPoint(SPoint3 (v[j]->x(),v[j]->y(),v[j]->z())); - SPoint2 pp (0,0); - MFaceVertex *v1b = new MFaceVertex (v[j]->x(),v[j]->y(),v[j]->z(),gf,pp.x(),pp.y() ); - gf->mesh_vertices.push_back(v1b); - numberedV[out.trifacelist[i * 3 + j] -1] = v1b; - delete v[j]; - v[j] = v1b; - } - } - MTriangle *t = new MTriangle(v[0],v[1],v[2]); - gf->triangles.push_back(t); + for (int i = 0; i < out.numberoftrifaces; i++){ + MVertex *v[3]; + v[0] = numberedV[out.trifacelist[i * 3 + 0] -1]; + v[1] = numberedV[out.trifacelist[i * 3 + 1] -1]; + v[2] = numberedV[out.trifacelist[i * 3 + 2] -1]; + GFace *gf = gr->model()->faceByTag(out.trifacemarkerlist[i]); + for(int j = 0; j < 3; j++){ + if(v[j]->onWhat()->dim() == 3){ + v[j]->onWhat()->mesh_vertices.erase + (std::find(v[j]->onWhat()->mesh_vertices.begin(), + v[j]->onWhat()->mesh_vertices.end(), + v[j])); + MFaceVertex *v1b = new MFaceVertex(v[j]->x(), v[j]->y(), v[j]->z(), + gf, 0., 0.); + gf->mesh_vertices.push_back(v1b); + numberedV[out.trifacelist[i * 3 + j] -1] = v1b; + delete v[j]; + v[j] = v1b; + } } - for (int i = 0; i < out.numberoftetrahedra; i++) { + MTriangle *t = new MTriangle(v[0], v[1], v[2]); + gf->triangles.push_back(t); + } + for(int i = 0; i < out.numberoftetrahedra; i++){ MVertex *v1 = numberedV[out.tetrahedronlist[i * 4 + 0] -1]; MVertex *v2 = numberedV[out.tetrahedronlist[i * 4 + 1] -1]; MVertex *v3 = numberedV[out.tetrahedronlist[i * 4 + 2] -1]; @@ -178,130 +176,138 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, std::vector<MV MTetrahedron *t = new MTetrahedron(v1,v2,v3,v4); gr->tetrahedra.push_back(t); } - - } - #endif - #if defined(HAVE_NETGEN) + namespace nglib { #include "nglib.h" #include "nglib_addon.h" } using namespace nglib; -Ng_Mesh * buildNetgenStructure (GRegion *gr, int importVolumeMesh, std::vector<MVertex*> & numberedV) +Ng_Mesh *buildNetgenStructure(GRegion *gr, bool importVolumeMesh, + std::vector<MVertex*> &numberedV) { NgAddOn_Init(); - Ng_Mesh *_ngmesh = Ng_NewMesh(); + Ng_Mesh *ngmesh = Ng_NewMesh(); std::set<MVertex*> allBoundingVertices; - getAllBoundingVertices (gr, importVolumeMesh, allBoundingVertices ); + getAllBoundingVertices(gr, allBoundingVertices); - std::set<MVertex*>::iterator itv = allBoundingVertices.begin(); - int I=1; - while (itv != allBoundingVertices.end()) - { + std::set<MVertex*>::iterator itv = allBoundingVertices.begin(); + int I = 1; + while(itv != allBoundingVertices.end()){ + double tmp[3]; + tmp[0] = (*itv)->x(); + tmp[1] = (*itv)->y(); + tmp[2] = (*itv)->z(); + (*itv)->setNum(I++); + numberedV.push_back(*itv); + Ng_AddPoint(ngmesh, tmp); + ++itv; + } + + if(importVolumeMesh){ + for(unsigned int i = 0; i < gr->mesh_vertices.size(); i++){ double tmp[3]; - tmp[0] = (*itv)->x(); - tmp[1] = (*itv)->y(); - tmp[2] = (*itv)->z(); - (*itv)->setNum(I++); - numberedV.push_back(*itv); - Ng_AddPoint(_ngmesh, tmp); - ++itv; + tmp[0] = gr->mesh_vertices[i]->x(); + tmp[1] = gr->mesh_vertices[i]->y(); + tmp[2] = gr->mesh_vertices[i]->z(); + gr->mesh_vertices[i]->setNum(I++); + Ng_AddPoint(ngmesh, tmp); } + } 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]; - int tmp[3]; - tmp[0] = t->getVertex(0)->getNum(); - tmp[1] = t->getVertex(1)->getNum(); - tmp[2] = t->getVertex(2)->getNum(); - Ng_AddSurfaceElement(_ngmesh, NG_TRIG, tmp); - } - ++it; + while(it != faces.end()){ + GFace *gf = (*it); + for(unsigned int i = 0; i< gf->triangles.size(); i++){ + MTriangle *t = gf->triangles[i]; + int tmp[3]; + tmp[0] = t->getVertex(0)->getNum(); + tmp[1] = t->getVertex(1)->getNum(); + tmp[2] = t->getVertex(2)->getNum(); + Ng_AddSurfaceElement(ngmesh, NG_TRIG, tmp); } - - if (importVolumeMesh) - { - for (unsigned int i = 0; i< gr->tetrahedra.size(); i++) - { - MTetrahedron *t = gr->tetrahedra[i]; - int tmp[4]; - tmp[0] = t->getVertex(0)->getNum(); - tmp[1] = t->getVertex(1)->getNum(); - tmp[2] = t->getVertex(2)->getNum(); - tmp[3] = t->getVertex(2)->getNum(); - Ng_AddVolumeElement(_ngmesh, NG_TET, tmp); - } + ++it; + } + + if(importVolumeMesh){ + for(unsigned int i = 0; i< gr->tetrahedra.size(); i++){ + MTetrahedron *t = gr->tetrahedra[i]; + int tmp[4]; + tmp[0] = t->getVertex(0)->getNum(); + tmp[1] = t->getVertex(1)->getNum(); + tmp[2] = t->getVertex(2)->getNum(); + tmp[3] = t->getVertex(3)->getNum(); + Ng_AddVolumeElement(ngmesh, NG_TET, tmp); } - - return _ngmesh; + } + + return ngmesh; } -void TransferVolumeMesh(GRegion *gr, Ng_Mesh * _ngmesh, std::vector<MVertex*> &numberedV) +void TransferVolumeMesh(GRegion *gr, Ng_Mesh *ngmesh, + std::vector<MVertex*> &numberedV) { // Gets total number of vertices of Netgen's mesh - int nbv = Ng_GetNP(_ngmesh); + int nbv = Ng_GetNP(ngmesh); if(!nbv) return; int nbpts = numberedV.size(); // Create new volume vertices - for(int i = nbpts; i < nbv; i++) - { - double tmp[3]; - Ng_GetPoint(_ngmesh, i+1, tmp); - MVertex *v = new MVertex (tmp[0],tmp[1],tmp[2],gr); - numberedV.push_back(v); - gr->mesh_vertices.push_back(v); - } + for(int i = nbpts; i < nbv; i++){ + double tmp[3]; + Ng_GetPoint(ngmesh, i + 1, tmp); + MVertex *v = new MVertex (tmp[0], tmp[1], tmp[2], gr); + numberedV.push_back(v); + gr->mesh_vertices.push_back(v); + } // Get total number of simplices of Netgen's mesh - int nbe = Ng_GetNE(_ngmesh); + int nbe = Ng_GetNE(ngmesh); // Create new volume simplices - for(int i = 0; i < nbe; i++) { + for(int i = 0; i < nbe; i++){ int tmp[4]; - Ng_GetVolumeElement(_ngmesh, i+1, tmp); - MTetrahedron *t = new MTetrahedron ( numberedV [tmp[0]-1], - numberedV [tmp[1]-1], - numberedV [tmp[2]-1], - numberedV [tmp[3]-1]); + Ng_GetVolumeElement(ngmesh, i + 1, tmp); + MTetrahedron *t = new MTetrahedron(numberedV[tmp[0] - 1], + numberedV[tmp[1] - 1], + numberedV[tmp[2] - 1], + numberedV[tmp[3] - 1]); gr->tetrahedra.push_back(t); } } + #endif -void deMeshGRegion :: operator() (GRegion *gr) +void deMeshGRegion::operator() (GRegion *gr) { - for (unsigned int i=0;i<gr->mesh_vertices.size();i++) delete gr->mesh_vertices[i]; + for(unsigned int i = 0; i < gr->mesh_vertices.size(); i++) + delete gr->mesh_vertices[i]; gr->mesh_vertices.clear(); - for (unsigned int i=0;i<gr->tetrahedra.size();i++) delete gr->tetrahedra[i]; + for(unsigned int i = 0; i < gr->tetrahedra.size(); i++) + delete gr->tetrahedra[i]; gr->tetrahedra.clear(); - for (unsigned int i=0;i<gr->hexahedra.size();i++) delete gr->hexahedra[i]; + for(unsigned int i = 0; i < gr->hexahedra.size(); i++) + delete gr->hexahedra[i]; gr->hexahedra.clear(); - for (unsigned int i=0;i<gr->prisms.size();i++) delete gr->prisms[i]; + for(unsigned int i = 0; i < gr->prisms.size(); i++) + delete gr->prisms[i]; gr->prisms.clear(); - for (unsigned int i=0;i<gr->pyramids.size();i++) delete gr->pyramids[i]; + for(unsigned int i = 0; i < gr->pyramids.size(); i++) + delete gr->pyramids[i]; gr->pyramids.clear(); + if(gr->meshRep) gr->meshRep->destroy(); } -int intersect_line_triangle ( double X[3], - double Y[3], - double Z[3] , - double P[3], - double N[3] ) +int intersect_line_triangle(double X[3], double Y[3], double Z[3] , + double P[3], double N[3] ) { double mat[3][3], det; double b[3], res[3]; @@ -328,36 +334,31 @@ int intersect_line_triangle ( double X[3], // Msg(INFO,"going there %g %g %g",res[0],res[1],res[2]); - if(res[0] >= eps_prec && res[0] <= 1.0 - eps_prec && - res[1] >= eps_prec && res[1] <= 1.0 - eps_prec && - 1-res[0]-res[1] >= eps_prec && 1-res[0]-res[1] <= 1.0 - eps_prec ) - { - // the line clearly intersects the triangle - return (res[2] > 0) ? 1 : 0; - } - else if (res[0] < -eps_prec || res[0] > 1.0 + eps_prec || - res[1] < -eps_prec || res[1] > 1.0 + eps_prec || - 1-res[0]-res[1] < -eps_prec || 1-res[0]-res[1] > 1.0 + eps_prec ) - { - // the line clearly does NOT intersect the triangle - return 0; - } - else - { - // the intersection is not robust, try another triangle - return -10000; - } - + if(res[0] >= eps_prec && res[0] <= 1.0 - eps_prec && + res[1] >= eps_prec && res[1] <= 1.0 - eps_prec && + 1 - res[0] - res[1] >= eps_prec && 1 - res[0] - res[1] <= 1.0 - eps_prec){ + // the line clearly intersects the triangle + return (res[2] > 0) ? 1 : 0; + } + else if(res[0] < -eps_prec || res[0] > 1.0 + eps_prec || + res[1] < -eps_prec || res[1] > 1.0 + eps_prec || + 1 - res[0] - res[1] < -eps_prec || 1 - res[0] - res[1] > 1.0 + eps_prec){ + // the line clearly does NOT intersect the triangle + return 0; + } + else{ + // the intersection is not robust, try another triangle + return -10000; + } } -void setRand (double r[6]) +void setRand(double r[6]) { for (int i=0;i<6;i++) - r[i]= 0.0001*((double)rand() / (double)RAND_MAX); + r[i] = 0.0001 * ((double)rand() / (double)RAND_MAX); } - -void meshNormalsPointOutOfTheRegion (GRegion *gr) +void meshNormalsPointOutOfTheRegion(GRegion *gr) { std::list<GFace*> faces = gr->faces(); std::list<GFace*>::iterator it = faces.begin(); @@ -365,72 +366,65 @@ void meshNormalsPointOutOfTheRegion (GRegion *gr) double rrr[6]; setRand(rrr); - while (it != faces.end()) - { - GFace *gf = (*it); - int nb_intersect = 0; - for (unsigned int i = 0; i< gf->triangles.size(); i++) - { - MTriangle *t = gf->triangles[i]; - double X[3] = {t->getVertex(0)->x(),t->getVertex(1)->x(),t->getVertex(2)->x()}; - double Y[3] = {t->getVertex(0)->y(),t->getVertex(1)->y(),t->getVertex(2)->y()}; - double Z[3] = {t->getVertex(0)->z(),t->getVertex(1)->z(),t->getVertex(2)->z()}; - double P[3] = {(X[0]+X[1]+X[2])/3.,(Y[0]+Y[1]+Y[2])/3.,(Z[0]+Z[1]+Z[2])/3.}; - double v1[3] = {X[0]-X[1],Y[0]-Y[1],Z[0]-Z[1]}; - double v2[3] = {X[2]-X[1],Y[2]-Y[1],Z[2]-Z[1]}; - double N[3] ; - prodve ( v1, v2 , N ); - norme (v1); - norme (v2); - norme(N); - N[0]+= rrr[0]*v1[0]+rrr[1]*v2[0]; - N[1]+= rrr[2]*v1[1]+rrr[3]*v2[1]; - N[2]+= rrr[4]*v1[2]+rrr[5]*v2[2]; - norme(N); - std::list<GFace*>::iterator it_b = faces.begin(); - while (it_b != faces.end()) - { - GFace *gf_b = (*it_b); - for (unsigned int i_b = 0; i_b< gf_b->triangles.size(); i_b++) - { - MTriangle *t_b = gf_b->triangles[i_b]; - if (t_b != t) - { - double X_b[3] = {t_b->getVertex(0)->x(),t_b->getVertex(1)->x(),t_b->getVertex(2)->x()}; - double Y_b[3] = {t_b->getVertex(0)->y(),t_b->getVertex(1)->y(),t_b->getVertex(2)->y()}; - double Z_b[3] = {t_b->getVertex(0)->z(),t_b->getVertex(1)->z(),t_b->getVertex(2)->z()}; - int inters = intersect_line_triangle ( X_b,Y_b,Z_b,P,N ); - nb_intersect += inters; - } - } - ++it_b; - } - Msg (INFO,"Region %d Face %d, %d intersect",gr->tag(),gf->tag(),nb_intersect); - if (nb_intersect >= 0) break; // negative value means intersetcion is not "robust" - } - - - if (nb_intersect < 0) - { - setRand(rrr); + while(it != faces.end()){ + GFace *gf = (*it); + int nb_intersect = 0; + for(unsigned int i = 0; i< gf->triangles.size(); i++){ + MTriangle *t = gf->triangles[i]; + double X[3] = {t->getVertex(0)->x(), t->getVertex(1)->x(), t->getVertex(2)->x()}; + double Y[3] = {t->getVertex(0)->y(), t->getVertex(1)->y(), t->getVertex(2)->y()}; + double Z[3] = {t->getVertex(0)->z(), t->getVertex(1)->z(), t->getVertex(2)->z()}; + double P[3] = {(X[0]+X[1]+X[2])/3., (Y[0]+Y[1]+Y[2])/3., (Z[0]+Z[1]+Z[2])/3.}; + double v1[3] = {X[0]-X[1], Y[0]-Y[1], Z[0]-Z[1]}; + double v2[3] = {X[2]-X[1], Y[2]-Y[1], Z[2]-Z[1]}; + double N[3] ; + prodve(v1, v2, N); + norme(v1); + norme(v2); + norme(N); + N[0] += rrr[0] * v1[0] + rrr[1] * v2[0]; + N[1] += rrr[2] * v1[1] + rrr[3] * v2[1]; + N[2] += rrr[4] * v1[2] + rrr[5] * v2[2]; + norme(N); + std::list<GFace*>::iterator it_b = faces.begin(); + while(it_b != faces.end()){ + GFace *gf_b = (*it_b); + for(unsigned int i_b = 0; i_b < gf_b->triangles.size(); i_b++){ + MTriangle *t_b = gf_b->triangles[i_b]; + if(t_b != t){ + double X_b[3] = {t_b->getVertex(0)->x(), t_b->getVertex(1)->x(), + t_b->getVertex(2)->x()}; + double Y_b[3] = {t_b->getVertex(0)->y(), t_b->getVertex(1)->y(), + t_b->getVertex(2)->y()}; + double Z_b[3] = {t_b->getVertex(0)->z(), t_b->getVertex(1)->z(), + t_b->getVertex(2)->z()}; + int inters = intersect_line_triangle( X_b, Y_b, Z_b, P, N); + nb_intersect += inters; + } } - else - { - if (nb_intersect % 2 == 1) // odd nb of intersections: the normal points inside the region - { - for (unsigned int i = 0; i< gf->triangles.size(); i++) - { - gf->triangles[i]->revert(); - } - } - ++it; + ++it_b; + } + Msg (INFO,"Region %d Face %d, %d intersect", + gr->tag(), gf->tag(), nb_intersect); + if (nb_intersect >= 0) break; // negative value means intersection is not "robust" + } + + if (nb_intersect < 0){ + setRand(rrr); + } + else{ + if(nb_intersect % 2 == 1){ + // odd nb of intersections: the normal points inside the region + for (unsigned int i = 0; i< gf->triangles.size(); i++){ + gf->triangles[i]->revert(); } + } + ++it; } + } } - - -void meshGRegion :: operator() (GRegion *gr) +void meshGRegion::operator() (GRegion *gr) { if(gr->geomType() == GEntity::DiscreteVolume) return; @@ -441,8 +435,7 @@ void meshGRegion :: operator() (GRegion *gr) // Send a messsage to the GMSH environment Msg(STATUS2, "Meshing volume %d", gr->tag()); - if(CTX.mesh.algo3d == DELAUNAY_TETGEN || CTX.mesh.algo3d == DELAUNAY_ISO) - { + if(CTX.mesh.algo3d == DELAUNAY_TETGEN || CTX.mesh.algo3d == DELAUNAY_ISO){ #if !defined(HAVE_TETGEN) Msg(GERROR, "Tetgen is not compiled in this version of Gmsh"); #else @@ -452,63 +445,57 @@ void meshGRegion :: operator() (GRegion *gr) std::list<GFace*> faces = gr->faces(); std::list<GFace*> allFaces; GModel::fiter fit = gr->model()->firstFace() ; - while (fit != gr->model()->lastFace()) - { - allFaces.push_back(*fit); - ++fit; - } - gr->set ( allFaces ); + while (fit != gr->model()->lastFace()){ + allFaces.push_back(*fit); + ++fit; + } + gr->set(allFaces); // mesh with tetgen, possibly changing the mesh on boundaries tetgenio in, out; std::vector<MVertex*> numberedV; char opts[128]; - buildTetgenStructure ( gr, in, numberedV); - sprintf(opts, "pe%c", - (CTX.verbosity < 3)? 'Q': (CTX.verbosity > 6)? 'V': '\0'); + buildTetgenStructure(gr, in, numberedV); + sprintf(opts, "pe%c", (CTX.verbosity < 3) ? 'Q': (CTX.verbosity > 6)? 'V': '\0'); tetrahedralize(opts, &in, &out); TransferTetgenMesh(gr, in, out, numberedV); // now do insertion of points - void insertVerticesInRegion (GRegion *gr) ; + void insertVerticesInRegion(GRegion *gr) ; insertVerticesInRegion (gr); // restore the initial set of faces - gr->set ( faces ); - + gr->set(faces); #endif - } - - if(CTX.mesh.algo3d == FRONTAL_NETGEN) - { + } + + if(CTX.mesh.algo3d == FRONTAL_NETGEN){ #if !defined(HAVE_NETGEN) Msg(GERROR, "Netgen is not compiled in this version of Gmsh"); #else // orient the triangles of with respect to this region - meshNormalsPointOutOfTheRegion (gr); + meshNormalsPointOutOfTheRegion(gr); std::vector<MVertex*> numberedV; - Ng_Mesh * _ngmesh = buildNetgenStructure (gr, 0, numberedV); + Ng_Mesh *ngmesh = buildNetgenStructure(gr, false, numberedV); Ng_Meshing_Parameters mp; mp.maxh = 1; mp.fineness = 1; mp.secondorder = 0; - NgAddOn_GenerateVolumeMesh(_ngmesh, &mp); // does not optimize - TransferVolumeMesh(gr, _ngmesh, numberedV); - Ng_DeleteMesh(_ngmesh); + NgAddOn_GenerateVolumeMesh(ngmesh, &mp); // does not optimize + TransferVolumeMesh(gr, ngmesh, numberedV); + Ng_DeleteMesh(ngmesh); Ng_Exit(); #endif - } + } } -void optimizeMeshGRegion :: operator() (GRegion *gr) +void optimizeMeshGRegion::operator() (GRegion *gr) { if(gr->geomType() == GEntity::DiscreteVolume) return; - + #if !defined(HAVE_NETGEN) Msg(GERROR, "Netgen is not compiled in this version of Gmsh"); #else - Msg(GERROR, "FIXME -- DEBUG THIS"); - // import mesh into netgen, including volume tets std::vector<MVertex*> numberedV; - Ng_Mesh * _ngmesh = buildNetgenStructure (gr, 1, numberedV); + Ng_Mesh *ngmesh = buildNetgenStructure(gr, true, numberedV); // delete volume vertices and tets deMeshGRegion dem; dem(gr); @@ -517,9 +504,9 @@ void optimizeMeshGRegion :: operator() (GRegion *gr) mp.maxh = 1; mp.fineness = 1; mp.secondorder = 0; - NgAddOn_OptimizeVolumeMesh(_ngmesh, &mp); - TransferVolumeMesh(gr, _ngmesh, numberedV); - Ng_DeleteMesh(_ngmesh); + NgAddOn_OptimizeVolumeMesh(ngmesh, &mp); + TransferVolumeMesh(gr, ngmesh, numberedV); + Ng_DeleteMesh(ngmesh); Ng_Exit(); #endif } diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp index 901eee1782871498c87b4615daa497ae641f0bc4..66eb8407424ba39b6f854716eaf5d750d161bfcb 100644 --- a/Mesh/meshGRegionDelaunayInsertion.cpp +++ b/Mesh/meshGRegionDelaunayInsertion.cpp @@ -1,3 +1,24 @@ +// $Id: meshGRegionDelaunayInsertion.cpp,v 1.5 2006-11-26 01:03:17 geuzaine Exp $ +// +// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program 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 this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA. +// +// Please report all bugs and problems to <gmsh@geuz.org>. + #include "meshGRegionDelaunayInsertion.h" #include "GRegion.h" #include "Numeric.h"