diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 4cf7f895658079c4b0c724306d5a545574e47bf6..c70d9efcbfb4e440aa86df6cabae346e3d83196a 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.33 2008-01-20 10:10:42 geuzaine Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.34 2008-01-21 21:25:29 geuzaine Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -32,82 +32,82 @@
#include <map>
#include <algorithm>
-int MTet4::inCircumSphere ( const double *p ) const
+int MTet4::inCircumSphere(const double *p) const
{
- double pa[3] = {base->getVertex(0)->x(),base->getVertex(0)->y(),base->getVertex(0)->z()};
- double pb[3] = {base->getVertex(1)->x(),base->getVertex(1)->y(),base->getVertex(1)->z()};
- double pc[3] = {base->getVertex(2)->x(),base->getVertex(2)->y(),base->getVertex(2)->z()};
- double pd[3] = {base->getVertex(3)->x(),base->getVertex(3)->y(),base->getVertex(3)->z()};
- double result = gmsh::insphere(pa, pb, pc, pd, (double*)p) * gmsh::orient3d(pa, pb, pc, pd);
+ double pa[3] = {base->getVertex(0)->x(),
+ base->getVertex(0)->y(),
+ base->getVertex(0)->z()};
+ double pb[3] = {base->getVertex(1)->x(),
+ base->getVertex(1)->y(),
+ base->getVertex(1)->z()};
+ double pc[3] = {base->getVertex(2)->x(),
+ base->getVertex(2)->y(),
+ base->getVertex(2)->z()};
+ double pd[3] = {base->getVertex(3)->x(),
+ base->getVertex(3)->y(),
+ base->getVertex(3)->z()};
+ double result =
+ gmsh::insphere(pa, pb, pc, pd, (double*)p) * gmsh::orient3d(pa, pb, pc, pd);
return (result > 0) ? 1 : 0;
}
-static int faces[4][3] = {{0,1,2},{0,2,3},{0,1,3},{1,2,3}};
+static int faces[4][3] = {{0,1,2}, {0,2,3}, {0,1,3}, {1,2,3}};
-struct faceXtet
-{
+struct faceXtet{
MVertex *v[3];
- MTet4 * t1;
+ MTet4 *t1;
int i1;
- faceXtet ( MTet4 *_t, int iFac)
- : t1(_t),i1(iFac)
+ faceXtet(MTet4 *_t, int iFac) : t1(_t), i1(iFac)
{
- v[0] = t1->tet()->getVertex ( faces [iFac][0] );
- v[1] = t1->tet()->getVertex ( faces [iFac][1] );
- v[2] = t1->tet()->getVertex ( faces [iFac][2] );
- std::sort ( v, v+3 );
+ v[0] = t1->tet()->getVertex(faces[iFac][0]);
+ v[1] = t1->tet()->getVertex(faces[iFac][1]);
+ v[2] = t1->tet()->getVertex(faces[iFac][2]);
+ std::sort(v, v + 3);
}
- inline bool operator < ( const faceXtet & other) const
+ inline bool operator < (const faceXtet & other) const
{
- if (v[0] < other.v[0])return true;
- if (v[0] > other.v[0])return false;
- if (v[1] < other.v[1])return true;
- if (v[1] > other.v[1])return false;
- if (v[2] < other.v[2])return true;
+ if (v[0] < other.v[0]) return true;
+ if (v[0] > other.v[0]) return false;
+ if (v[1] < other.v[1]) return true;
+ if (v[1] > other.v[1]) return false;
+ if (v[2] < other.v[2]) return true;
return false;
}
};
-
template <class ITER>
void connectTets ( ITER beg, ITER end)
{
std::set<faceXtet> conn;
- {
- while (beg != end)
- {
- if (!(*beg)->isDeleted())
- {
- for (int i=0;i<4;i++)
- {
- faceXtet fxt ( *beg , i );
- std::set<faceXtet>::iterator found = conn.find (fxt);
- if (found == conn.end())
- conn.insert ( fxt );
- else if (found->t1 != *beg)
- {
- found->t1->setNeigh(found->i1,*beg);
- (*beg)->setNeigh ( i, found->t1);
- }
- }
- }
- ++beg;
+ while (beg != end){
+ if (!(*beg)->isDeleted()){
+ for (int i = 0; i < 4; i++){
+ faceXtet fxt(*beg, i);
+ std::set<faceXtet>::iterator found = conn.find(fxt);
+ if (found == conn.end())
+ conn.insert(fxt);
+ else if (found->t1 != *beg){
+ found->t1->setNeigh(found->i1, *beg);
+ (*beg)->setNeigh(i, found->t1);
+ }
}
+ }
+ ++beg;
}
}
-void connectTets ( std::list<MTet4*> & l) {connectTets(l.begin(),l.end());}
-void connectTets ( std::vector<MTet4*> &l ) {connectTets(l.begin(),l.end());}
+void connectTets(std::list<MTet4*> &l) { connectTets(l.begin(), l.end()); }
+void connectTets(std::vector<MTet4*> &l) { connectTets(l.begin(), l.end()); }
-void recurFindCavity ( std::list<faceXtet> & shell,
- std::list<MTet4*> & cavity,
- MVertex *v ,
- MTet4 *t)
+void recurFindCavity(std::list<faceXtet> & shell,
+ std::list<MTet4*> & cavity,
+ MVertex *v ,
+ MTet4 *t)
{
-// Msg(INFO,"tet %d %d %d %d",t->tet()->getVertex(0)->getNum(),
-// t->tet()->getVertex(1)->getNum(),
-// t->tet()->getVertex(2)->getNum(),
-// t->tet()->getVertex(3)->getNum());
+ // Msg(INFO,"tet %d %d %d %d",t->tet()->getVertex(0)->getNum(),
+ // t->tet()->getVertex(1)->getNum(),
+ // t->tet()->getVertex(2)->getNum(),
+ // t->tet()->getVertex(3)->getNum());
// invariant : this one has to be inserted in the cavity
// consider this tet deleted
@@ -117,36 +117,34 @@ void recurFindCavity ( std::list<faceXtet> & shell,
// because it violates delaunay criterion
cavity.push_back(t);
- for (int i=0;i<4;i++)
- {
- MTet4 *neigh = t->getNeigh(i) ;
- if (!neigh)
- shell.push_back ( faceXtet ( t, i ) );
- else if (!neigh->isDeleted() )
- {
- int circ = neigh->inCircumSphere ( v );
- if (circ && (neigh->onWhat() == t->onWhat()))
- recurFindCavity ( shell, cavity,v , neigh);
- else
- shell.push_back ( faceXtet ( t, i ) );
- }
+ for (int i = 0; i < 4; i++){
+ MTet4 *neigh = t->getNeigh(i) ;
+ if (!neigh)
+ shell.push_back(faceXtet(t, i));
+ else if (!neigh->isDeleted()){
+ int circ = neigh->inCircumSphere(v);
+ if (circ && (neigh->onWhat() == t->onWhat()))
+ recurFindCavity(shell, cavity, v, neigh);
+ else
+ shell.push_back(faceXtet(t, i));
}
+ }
}
-bool insertVertex (MVertex *v ,
- MTet4 *t ,
- MTet4Factory &myFactory,
- std::set<MTet4*,compareTet4Ptr> &allTets,
- std::vector<double> & vSizes)
+bool insertVertex(MVertex *v,
+ MTet4 *t,
+ MTet4Factory &myFactory,
+ std::set<MTet4*,compareTet4Ptr> &allTets,
+ std::vector<double> & vSizes)
{
- std::list<faceXtet> shell;
- std::list<MTet4*> cavity;
- std::list<MTet4*> new_cavity;
+ std::list<faceXtet> shell;
+ std::list<MTet4*> cavity;
+ std::list<MTet4*> new_cavity;
- recurFindCavity ( shell, cavity, v , t);
+ recurFindCavity(shell, cavity, v, t);
-// Msg(INFO,"%d %d",cavity.size(),NC);
-// if (NC != cavity.size())throw;
+ // Msg(INFO,"%d %d",cavity.size(),NC);
+ // if (NC != cavity.size())throw;
// check that volume is conserved
double newVolume = 0;
@@ -156,12 +154,10 @@ bool insertVertex (MVertex *v ,
// FILE *ff2 = fopen (name2,"w");
// fprintf(ff2,"View\"test\"{\n");
-
std::list<MTet4*>::iterator ittet = cavity.begin();
std::list<MTet4*>::iterator ittete = cavity.end();
- while ( ittet != ittete )
- {
- oldVolume += fabs((*ittet)->getVolume());
+ while (ittet != ittete){
+ oldVolume += fabs((*ittet)->getVolume());
// fprintf(ff2,"SS(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0,0};\n",
// (*ittet)->tet()->getVertex(0)->x(),
// (*ittet)->tet()->getVertex(0)->y(),
@@ -196,16 +192,11 @@ bool insertVertex (MVertex *v ,
MTet4** newTets = new MTet4*[shell.size()];;
int k = 0;
- std::list<faceXtet>::iterator it = shell.begin();
-
+ std::list<faceXtet>::iterator it = shell.begin();
- while (it != shell.end())
- {
- MTetrahedron *tr = new MTetrahedron ( it->v[0],
- it->v[1],
- it->v[2],
- v);
- // Msg(INFO,"shell %d %d %d",it->v[0]->getNum(),it->v[1]->getNum(),it->v[2]->getNum());
+ while (it != shell.end()){
+ MTetrahedron *tr = new MTetrahedron (it->v[0], it->v[1], it->v[2], v);
+ // Msg(INFO,"shell %d %d %d",it->v[0]->getNum(),it->v[1]->getNum(),it->v[2]->getNum());
// fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
// it->v[0]->x(),
// it->v[0]->y(),
@@ -217,20 +208,18 @@ bool insertVertex (MVertex *v ,
// it->v[2]->y(),
// it->v[2]->z());
- MTet4 *t4 = myFactory.Create ( tr , vSizes);
+ MTet4 *t4 = myFactory.Create(tr, vSizes);
t4->setOnWhat(t->onWhat());
- newTets[k++]=t4;
- // all new tets are pushed front in order to
- // ba able to destroy them if the cavity is not
- // star shaped around the new vertex.
- // here, we better use robust perdicates
- // that implies to orient all faces and
- // ensure that the tets are all oriented the same.
- new_cavity.push_back (t4);
+ newTets[k++] = t4;
+ // all new tets are pushed front in order to ba able to destroy
+ // them if the cavity is not star shaped around the new vertex.
+ // here, we better use robust perdicates that implies to orient
+ // all faces and ensure that the tets are all oriented the same.
+ new_cavity.push_back(t4);
MTet4 *otherSide = it->t1->getNeigh(it->i1);
if (otherSide)
- new_cavity.push_back (otherSide);
+ new_cavity.push_back(otherSide);
// if (!it->t1->isDeleted())throw;
newVolume += fabs(t4->getVolume());
++it;
@@ -239,162 +228,138 @@ bool insertVertex (MVertex *v ,
// fclose (ff);
// Msg(INFO,"new cavity of vol %g (%d boundaries)",newVolume,shell.size());
// OK, the cavity is star shaped
- if (fabs(oldVolume - newVolume) < 1.e-10 * oldVolume )
- {
- connectTets ( new_cavity.begin(),new_cavity.end() );
- allTets.insert(newTets,newTets+shell.size());
-
-// ittet = cavity.begin();
-// ittete = cavity.end();
-// while ( ittet != ittete )
-// {
-// myFactory.Free (*ittet);
-// ++ittet;
-// }
- delete [] newTets;
- return true;
- }
- // The cavity is NOT star shaped
- else
- {
- for (unsigned int i=0;i<shell.size();i++)myFactory.Free(newTets[i]);
- delete [] newTets;
- ittet = cavity.begin();
- ittete = cavity.end();
- while ( ittet != ittete )
- {
- (*ittet)->setDeleted(false);
- ++ittet;
- }
- return false;
+ if (fabs(oldVolume - newVolume) < 1.e-10 * oldVolume){
+ connectTets ( new_cavity.begin(),new_cavity.end());
+ allTets.insert(newTets,newTets+shell.size());
+
+// ittet = cavity.begin();
+// ittete = cavity.end();
+// while ( ittet != ittete ){
+// myFactory.Free (*ittet);
+// ++ittet;
+// }
+ delete [] newTets;
+ return true;
+ }
+ else { // The cavity is NOT star shaped
+ for (unsigned int i = 0; i <shell.size(); i++) myFactory.Free(newTets[i]);
+ delete [] newTets;
+ ittet = cavity.begin();
+ ittete = cavity.end();
+ while(ittet != ittete){
+ (*ittet)->setDeleted(false);
+ ++ittet;
}
+ return false;
+ }
}
-static void setLcs ( MTetrahedron *t, std::map<MVertex*,double> &vSizes)
+static void setLcs(MTetrahedron *t, std::map<MVertex*,double> &vSizes)
{
- for (int i=0;i<4;i++)
- {
- for (int j=i+1;j<4;j++)
- {
- MVertex *vi = t->getVertex(i);
- MVertex *vj = t->getVertex(j);
- double dx = vi->x()-vj->x();
- double dy = vi->y()-vj->y();
- double dz = vi->z()-vj->z();
- double l = sqrt(dx*dx+dy*dy+dz*dz);
- std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
- std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
- if (iti==vSizes.end() || iti->second > l)vSizes[vi] = l;
- if (itj==vSizes.end() || itj->second > l)vSizes[vj] = l;
- }
+ for (int i = 0; i < 4; i++){
+ for (int j = i + 1; j < 4; j++){
+ MVertex *vi = t->getVertex(i);
+ MVertex *vj = t->getVertex(j);
+ double dx = vi->x()-vj->x();
+ double dy = vi->y()-vj->y();
+ double dz = vi->z()-vj->z();
+ double l = sqrt(dx * dx + dy * dy + dz * dz);
+ std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
+ std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
+ if (iti == vSizes.end() || iti->second > l) vSizes[vi] = l;
+ if (itj == vSizes.end() || itj->second > l) vSizes[vj] = l;
}
+ }
}
// void recover_volumes( GRegion *gr , std::set<MTet4*,compareTet4Ptr> & allTets )
// {
// std::set<MTet4*,compareTet4Ptr>::iterator it = allTets.begin();
-// for ( ; it != allTets.end() ; ++it )
-// {
-// MTet4 *t = *allTets.begin();
-// if (!t->isDeleted())
-// {
-// }
+// for (; it != allTets.end(); ++it){
+// MTet4 *t = *allTets.begin();
+// if (!t->isDeleted()){
// }
+// }
// }
// 4th argument will disappear when the reclassification of vertices will be done
- bool find_triangle_in_model ( GModel *model , MTriangle *tri, GFace **gfound, bool force)
- {
-
- static compareMTriangleLexicographic cmp;
-
- GModel::fiter fit = model->firstFace() ;
- while (fit != model->lastFace()){
-
- bool found = std::binary_search((*fit)->triangles.begin(),
+bool find_triangle_in_model(GModel *model, MTriangle *tri, GFace **gfound, bool force)
+{
+ static compareMTriangleLexicographic cmp;
+
+ GModel::fiter fit = model->firstFace() ;
+ while (fit != model->lastFace()){
+ bool found = std::binary_search((*fit)->triangles.begin(),
(*fit)->triangles.end(),
- tri,cmp);
- if (found )
- {
- *gfound = *fit;
- return true;
- }
- ++fit;
- }
- return false;
-
- }
-
-
-
-
+ tri, cmp);
+ if(found){
+ *gfound = *fit;
+ return true;
+ }
+ ++fit;
+ }
+ return false;
+}
-GRegion *getRegionFromBoundingFaces (GModel *model,
+GRegion *getRegionFromBoundingFaces(GModel *model,
std::set<GFace *> &faces_bound)
{
- GModel::riter git = model->firstRegion() ;
+ GModel::riter git = model->firstRegion();
while (git != model->lastRegion()){
-
std::list <GFace *> _faces = (*git)->faces();
- if (_faces.size() == faces_bound.size())
- {
- bool ok = true;
- for (std::list <GFace *>::iterator it = _faces.begin();it != _faces.end();++it)
- {
- if (faces_bound.find (*it) == faces_bound.end())ok = false;
- }
- if (ok)return *git;
- }
+ if(_faces.size() == faces_bound.size()){
+ bool ok = true;
+ for (std::list<GFace *>::iterator it = _faces.begin(); it != _faces.end(); ++it){
+ if(faces_bound.find (*it) == faces_bound.end()) ok = false;
+ }
+ if (ok) return *git;
+ }
++git;
}
return 0;
}
-
-void recur_classify ( MTet4 *t ,
- std::list<MTet4*> &theRegion,
- std::set<GFace *> &faces_bound,
- GRegion *bidon ,
- GModel *model,
- const fs_cont &search)
+void recur_classify(MTet4 *t ,
+ std::list<MTet4*> &theRegion,
+ std::set<GFace *> &faces_bound,
+ GRegion *bidon ,
+ GModel *model,
+ const fs_cont &search)
{
if (!t) Msg (GERROR,"a tet is not connected by a boundary face");
- if (t->onWhat())return; // should never return here...
+ if (t->onWhat()) return; // should never return here...
theRegion.push_back(t);
t->setOnWhat(bidon);
bool FF[4] = {0,0,0,0};
- for (int i=0;i<4;i++)
- {
- // if (!t->getNeigh(i) || !t->getNeigh(i)->onWhat())
- {
-
- GFace* gfound = findInFaceSearchStructure ( t->tet()->getVertex ( faces[i][0] ),
- t->tet()->getVertex ( faces[i][1] ),
- t->tet()->getVertex ( faces[i][2] ), search );
- if (gfound)
- {
- FF[i]=true;
- if (faces_bound.find(gfound) == faces_bound.end())
- faces_bound.insert(gfound);
- }
-
-// MTriangle tri ( t->tet()->getVertex ( faces[i][0] ),
-// t->tet()->getVertex ( faces[i][1] ),
-// t->tet()->getVertex ( faces[i][2] ) );
-// GFace *gfound;
-// if (FF[i] = find_triangle_in_model ( model , &tri, &gfound, false))
-// {
-// if (faces_bound.find(gfound) == faces_bound.end())
-// faces_bound.insert(gfound);
-// }
- }
- }
- for (int i=0;i<4;i++)
+ for (int i = 0; i < 4; i++){
+ // if (!t->getNeigh(i) || !t->getNeigh(i)->onWhat())
{
- if (!FF[i]) recur_classify ( t->getNeigh(i) , theRegion, faces_bound, bidon, model,search );
+ GFace* gfound = findInFaceSearchStructure (t->tet()->getVertex(faces[i][0]),
+ t->tet()->getVertex(faces[i][1]),
+ t->tet()->getVertex(faces[i][2]),
+ search);
+ if (gfound){
+ FF[i] = true;
+ if (faces_bound.find(gfound) == faces_bound.end())
+ faces_bound.insert(gfound);
+ }
+
+// MTriangle tri (t->tet()->getVertex (faces[i][0]),
+// t->tet()->getVertex (faces[i][1]),
+// t->tet()->getVertex (faces[i][2]));
+// GFace *gfound;
+// if (FF[i] = find_triangle_in_model(model, &tri, &gfound, false)){
+// if (faces_bound.find(gfound) == faces_bound.end())
+// faces_bound.insert(gfound);
+// }
}
+ }
+ for (int i = 0; i < 4; i++){
+ if (!FF[i])
+ recur_classify(t->getNeigh(i), theRegion, faces_bound, bidon, model, search );
+ }
}
void adaptMeshGRegion::operator () (GRegion *gr)
@@ -419,85 +384,92 @@ void adaptMeshGRegion::operator () (GRegion *gr)
double worst = 1.0;
double avg = 0;
int count=0;
- for (int i=0;i<nbRanges;i++)quality_ranges[i] = 0;
- for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
- {
- if (!(*it)->isDeleted()){
- double vol = fabs((*it)->tet()->getVolume());
- double qual = (*it)->getQuality();
- worst = std::min(qual,worst);
- avg+=qual;
- count++;
- totalVolumeb+=vol;
- for (int i=0;i<nbRanges;i++){
- double low = (double)i/(nbRanges);
- double high = (double)(i+1)/(nbRanges);
- if (qual >= low && qual < high)quality_ranges[i]++;
- }
- }
+ for (int i = 0; i < nbRanges; i++) quality_ranges[i] = 0;
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end(); ++it){
+ if (!(*it)->isDeleted()){
+ double vol = fabs((*it)->tet()->getVolume());
+ double qual = (*it)->getQuality();
+ worst = std::min(qual, worst);
+ avg += qual;
+ count++;
+ totalVolumeb += vol;
+ for (int i = 0; i < nbRanges; i++){
+ double low = (double)i / nbRanges;
+ double high = (double)(i + 1) / nbRanges;
+ if (qual >= low && qual < high) quality_ranges[i]++;
+ }
}
- Msg(INFO,"Adaptation : START with %12.5E QBAD %12.5E QAVG %12.5E",totalVolumeb,worst,avg/count);
- for (int i=0;i<nbRanges;i++){
- double low = (double)i/(nbRanges);
- double high = (double)(i+1)/(nbRanges);
- Msg(INFO,"Opti : %3.2f < QUAL < %3.2f : %9d elements ",low,high,quality_ranges[i]);
+ }
+ Msg(INFO,"Adaptation : START with %12.5E QBAD %12.5E QAVG %12.5E",
+ totalVolumeb, worst, avg / count);
+ for (int i = 0; i < nbRanges; i++){
+ double low = (double)i / nbRanges;
+ double high = (double)(i + 1) / nbRanges;
+ Msg(INFO,"Opti : %3.2f < QUAL < %3.2f : %9d elements ",
+ low, high, quality_ranges[i]);
}
}
double qMin = 0.5;
double sliverLimit = 0.2;
- int nbESwap=0, nbFSwap=0, nbReloc=0, nbCollapse = 0;
+ int nbESwap = 0, nbFSwap = 0, nbReloc = 0, nbCollapse = 0;
while (1){
std::vector<MTet4*> newTets;
-
- for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it){
+ for (CONTAINER::iterator it = allTets.begin(); it!=allTets.end(); ++it){
if (!(*it)->isDeleted()){
- for (int i=0;i<4;i++){
- for (int j=0;j<4;j++){
- if (gmshCollapseVertex(newTets,*it,i,j,QMTET_2)){nbCollapse++;i=j=10;}
+ for (int i = 0; i < 4; i++){
+ for (int j = 0; j < 4; j++){
+ if (gmshCollapseVertex(newTets, *it, i, j, QMTET_2)){
+ nbCollapse++; i = j = 10;
+ }
}
}
}
}
+
+ printf("nbCollapses = %d\n", nbCollapse);
- printf("nbCollapses = %d\n",nbCollapse);
-
- for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it){
+ for (CONTAINER::iterator it = allTets.begin(); it!=allTets.end(); ++it){
if (!(*it)->isDeleted()){
double qq = (*it)->getQuality();
if (qq < qMin){
- for (int i=0;i<4;i++){
- if (gmshFaceSwap(newTets,*it,i,qm)){nbFSwap++;break;}
+ for (int i = 0; i < 4; i++){
+ if (gmshFaceSwap(newTets, *it, i, qm)){
+ nbFSwap++;
+ break;
+ }
}
}
}
}
illegals.clear();
- for (int i=0;i<nbRanges;i++)quality_ranges[i] = 0;
+ for (int i = 0; i < nbRanges; i++) quality_ranges[i] = 0;
- for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
- {
- if (!(*it)->isDeleted()){
- double qq = (*it)->getQuality();
- if (qq < qMin)
- for (int i=0;i<6;i++){
- if (gmshEdgeSwap(newTets,*it,i,qm)) {nbESwap++;break;}
+ for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
+ if (!(*it)->isDeleted()){
+ double qq = (*it)->getQuality();
+ if (qq < qMin)
+ for (int i = 0; i < 6; i++){
+ if (gmshEdgeSwap(newTets, *it, i, qm)) {
+ nbESwap++;
+ break;
}
+ }
if (!(*it)->isDeleted()){
- if (qq < sliverLimit)illegals.push_back(*it);
- for (int i=0;i<nbRanges;i++){
- double low = (double)i/(nbRanges);
- double high = (double)(i+1)/(nbRanges);
- if (qq >= low && qq < high)quality_ranges[i]++;
+ if (qq < sliverLimit) illegals.push_back(*it);
+ for (int i = 0; i < nbRanges; i++){
+ double low = (double)i / nbRanges;
+ double high = (double)(i + 1)/ nbRanges;
+ if (qq >= low && qq < high) quality_ranges[i]++;
}
}
- }
}
-
- if (!newTets.size())break;
+ }
+
+ if (!newTets.size()) break;
// add all the new tets in the container
for(unsigned int i = 0; i < newTets.size(); i++){
@@ -511,65 +483,65 @@ void adaptMeshGRegion::operator () (GRegion *gr)
}
// relocate vertices
- for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
- {
- if (!(*it)->isDeleted()){
- double qq = (*it)->getQuality();
- if (qq < qMin)
- for (int i=0;i<4;i++){
- if (gmshSmoothVertex(*it,i,qm))nbReloc++;
- }
- }
+ for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
+ if (!(*it)->isDeleted()){
+ double qq = (*it)->getQuality();
+ if (qq < qMin)
+ for (int i = 0; i < 4; i++){
+ if (gmshSmoothVertex(*it, i, qm)) nbReloc++;
+ }
}
-
+ }
+
double totalVolumeb = 0.0;
double worst = 1.0;
double avg = 0;
- int count=0;
- for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
- {
- if (!(*it)->isDeleted()){
- double vol = fabs((*it)->tet()->getVolume());
- double qual = (*it)->getQuality();
- worst = std::min(qual,worst);
- avg+=qual;
- count++;
- totalVolumeb+=vol;
- }
+ int count = 0;
+ for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
+ if (!(*it)->isDeleted()){
+ double vol = fabs((*it)->tet()->getVolume());
+ double qual = (*it)->getQuality();
+ worst = std::min(qual, worst);
+ avg += qual;
+ count++;
+ totalVolumeb += vol;
}
+ }
double t2 = Cpu();
- Msg(INFO,"Opti : (%d,%d,%d) = %12.5E QBAD %12.5E QAVG %12.5E (%8.3f sec)",nbESwap,nbFSwap,nbReloc,totalVolumeb,worst,avg/count,t2-t1);
+ Msg(INFO,"Opti : (%d,%d,%d) = %12.5E QBAD %12.5E QAVG %12.5E (%8.3f sec)",
+ nbESwap, nbFSwap, nbReloc, totalVolumeb, worst, avg / count, t2 - t1);
break;
}
int nbSlivers = 0;
for(unsigned int i = 0; i < illegals.size(); i++)
- if(!(illegals[i]->isDeleted())) nbSlivers ++;
+ if(!(illegals[i]->isDeleted())) nbSlivers++;
if (nbSlivers){
- Msg(INFO,"Opti : %d illegal tets are still in the mesh, trying to remove them",nbSlivers);
+ Msg(INFO,"Opti : %d illegal tets are still in the mesh, trying to remove them",
+ nbSlivers);
}
else{
- Msg(INFO,"Opti : no illegal tets in the mesh ;-)",nbSlivers);
+ Msg(INFO,"Opti : no illegal tets in the mesh ;-)", nbSlivers);
}
- for (int i=0;i<nbRanges;i++){
- double low = (double)i/(nbRanges);
- double high = (double)(i+1)/(nbRanges);
- Msg(INFO,"Opti : %3.2f < QUAL < %3.2f : %9d elements ",low,high,quality_ranges[i]);
+ for (int i = 0; i < nbRanges ;i++){
+ double low = (double)i / nbRanges;
+ double high = (double)(i + 1) / nbRanges;
+ Msg(INFO,"Opti : %3.2f < QUAL < %3.2f : %9d elements",
+ low, high, quality_ranges[i]);
}
-
- for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
- {
- if (!(*it)->isDeleted()){
- gr->tetrahedra.push_back((*it)->tet());
- delete *it;
- }
- else{
- delete (*it)->tet();
- delete *it;
- }
+
+ for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
+ if (!(*it)->isDeleted()){
+ gr->tetrahedra.push_back((*it)->tet());
+ delete *it;
+ }
+ else{
+ delete (*it)->tet();
+ delete *it;
}
+ }
}
//template <class CONTAINER, class DATA>
@@ -577,88 +549,94 @@ void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm)
{
typedef std::list<MTet4 *> CONTAINER ;
CONTAINER allTets;
- for(unsigned int i=0;i<gr->tetrahedra.size();i++){
- MTet4 * t = new MTet4(gr->tetrahedra[i],qm);
+ for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
+ MTet4 * t = new MTet4(gr->tetrahedra[i], qm);
t->setOnWhat(gr);
allTets.push_back(t);
}
gr->tetrahedra.clear();
- connectTets(allTets.begin(),allTets.end());
+ connectTets(allTets.begin(), allTets.end());
double t1 = Cpu();
std::vector<MTet4*> illegals;
- const int nbRanges=10;
+ const int nbRanges = 10;
int quality_ranges [nbRanges];
{
double totalVolumeb = 0.0;
double worst = 1.0;
double avg = 0;
- int count=0;
- for (int i=0;i<nbRanges;i++)quality_ranges[i] = 0;
- for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
- {
- if (!(*it)->isDeleted()){
- double vol = fabs((*it)->tet()->getVolume());
- double qual = (*it)->getQuality();
- worst = std::min(qual,worst);
- avg+=qual;
- count++;
- totalVolumeb+=vol;
- for (int i=0;i<nbRanges;i++){
- double low = (double)i/(nbRanges);
- double high = (double)(i+1)/(nbRanges);
- if (qual >= low && qual < high)quality_ranges[i]++;
- }
- }
+ int count = 0;
+ for (int i = 0; i < nbRanges; i++) quality_ranges[i] = 0;
+ for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
+ if (!(*it)->isDeleted()){
+ double vol = fabs((*it)->tet()->getVolume());
+ double qual = (*it)->getQuality();
+ worst = std::min(qual, worst);
+ avg += qual;
+ count++;
+ totalVolumeb += vol;
+ for (int i = 0; i < nbRanges; i++){
+ double low = (double)i / nbRanges;
+ double high = (double)(i + 1) / nbRanges;
+ if (qual >= low && qual < high) quality_ranges[i]++;
+ }
}
- Msg(INFO,"Opti : START with %12.5E QBAD %12.5E QAVG %12.5E",totalVolumeb,worst,avg/count);
- for (int i=0;i<nbRanges;i++){
- double low = (double)i/(nbRanges);
- double high = (double)(i+1)/(nbRanges);
- Msg(INFO,"Opti : %3.2f < QUAL < %3.2f : %9d elements ",low,high,quality_ranges[i]);
+ }
+ Msg(INFO,"Opti : START with %12.5E QBAD %12.5E QAVG %12.5E",
+ totalVolumeb, worst, avg / count);
+ for (int i = 0; i < nbRanges; i++){
+ double low = (double)i / nbRanges;
+ double high = (double)(i + 1) / nbRanges;
+ Msg(INFO,"Opti : %3.2f < QUAL < %3.2f : %9d elements",
+ low, high, quality_ranges[i]);
}
}
-
+
double qMin = 0.5;
double sliverLimit = 0.1;
- int nbESwap=0, nbFSwap=0, nbReloc=0;
+ int nbESwap = 0, nbFSwap = 0, nbReloc = 0;
while (1){
std::vector<MTet4*> newTets;
- for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it){
+ for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
if (!(*it)->isDeleted()){
double qq = (*it)->getQuality();
if (qq < qMin){
- for (int i=0;i<4;i++){
- if (gmshFaceSwap(newTets,*it,i,qm)){nbFSwap++;break;}
+ for (int i = 0; i < 4; i++){
+ if (gmshFaceSwap(newTets, *it, i, qm)){
+ nbFSwap++;
+ break;
+ }
}
}
}
}
-
+
illegals.clear();
- for (int i=0;i<nbRanges;i++)quality_ranges[i] = 0;
+ for (int i = 0; i < nbRanges; i++) quality_ranges[i] = 0;
- for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
- {
- if (!(*it)->isDeleted()){
- double qq = (*it)->getQuality();
- if (qq < qMin)
- for (int i=0;i<6;i++){
- if (gmshEdgeSwap(newTets,*it,i,qm)) {nbESwap++;break;}
+ for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
+ if (!(*it)->isDeleted()){
+ double qq = (*it)->getQuality();
+ if (qq < qMin)
+ for (int i = 0; i < 6; i++){
+ if (gmshEdgeSwap(newTets, *it, i, qm)) {
+ nbESwap++;
+ break;
}
- if (!(*it)->isDeleted()){
- if (qq < sliverLimit)illegals.push_back(*it);
- for (int i=0;i<nbRanges;i++){
- double low = (double)i/(nbRanges);
- double high = (double)(i+1)/(nbRanges);
- if (qq >= low && qq < high)quality_ranges[i]++;
- }
}
+ if (!(*it)->isDeleted()){
+ if (qq < sliverLimit) illegals.push_back(*it);
+ for (int i = 0; i < nbRanges; i++){
+ double low = (double)i / nbRanges;
+ double high = (double)(i + 1) / nbRanges;
+ if (qq >= low && qq < high) quality_ranges[i]++;
+ }
}
}
+ }
if (0 && !newTets.size()){
int nbSlivers = 0;
@@ -667,9 +645,9 @@ void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm)
if(!(illegals[i]->isDeleted())){
if(gmshSliverRemoval(newTets, illegals[i], qm))
nbSliversWeCanDoSomething++;
- nbSlivers ++;
+ nbSlivers++;
}
- Msg(INFO,"Opti : %d Sliver Removals",nbSliversWeCanDoSomething);
+ Msg(INFO, "Opti : %d Sliver Removals", nbSliversWeCanDoSomething);
}
if (!newTets.size()){
@@ -688,65 +666,61 @@ void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm)
}
// relocate vertices
- for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
- {
- if (!(*it)->isDeleted()){
- double qq = (*it)->getQuality();
- if (qq < qMin)
- for (int i=0;i<4;i++){
- if (gmshSmoothVertex(*it,i,qm))nbReloc++;
- }
- }
+ for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it){
+ if (!(*it)->isDeleted()){
+ double qq = (*it)->getQuality();
+ if (qq < qMin)
+ for (int i = 0; i < 4; i++){
+ if (gmshSmoothVertex(*it, i, qm)) nbReloc++;
+ }
}
+ }
double totalVolumeb = 0.0;
double worst = 1.0;
double avg = 0;
- int count=0;
- for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
- {
- if (!(*it)->isDeleted()){
- double vol = fabs((*it)->tet()->getVolume());
- double qual = (*it)->getQuality();
- worst = std::min(qual,worst);
- avg+=qual;
- count++;
- totalVolumeb+=vol;
- }
+ int count = 0;
+ for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
+ if (!(*it)->isDeleted()){
+ double vol = fabs((*it)->tet()->getVolume());
+ double qual = (*it)->getQuality();
+ worst = std::min(qual, worst);
+ avg += qual;
+ count++;
+ totalVolumeb += vol;
}
+ }
double t2 = Cpu();
- Msg(INFO,"Opti : (%d,%d,%d) = %12.5E QBAD %12.5E QAVG %12.5E (%8.3f sec)",nbESwap,nbFSwap,nbReloc,totalVolumeb,worst,avg/count,t2-t1);
+ Msg(INFO, "Opti : (%d,%d,%d) = %12.5E QBAD %12.5E QAVG %12.5E (%8.3f sec)",
+ nbESwap, nbFSwap, nbReloc, totalVolumeb, worst, avg / count, t2 - t1);
}
-
if (illegals.size()){
- Msg(INFO,"Opti : %d illegal tets are still in the mesh",illegals.size());
+ Msg(INFO, "Opti : %d illegal tets are still in the mesh", illegals.size());
}
else{
- Msg(INFO,"Opti : no illegal tets in the mesh ;-)");
+ Msg(INFO, "Opti : no illegal tets in the mesh ;-)");
}
- for (int i=0;i<nbRanges;i++){
- double low = (double)i/(nbRanges);
- double high = (double)(i+1)/(nbRanges);
- Msg(INFO,"Opti : %3.2f < QUAL < %3.2f : %9d elements ",low,high,quality_ranges[i]);
+ for (int i = 0; i < nbRanges; i++){
+ double low = (double)i / nbRanges;
+ double high = (double)(i + 1) / nbRanges;
+ Msg(INFO,"Opti : %3.2f < QUAL < %3.2f : %9d elements",
+ low, high, quality_ranges[i]);
}
- for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it)
- {
- if (!(*it)->isDeleted()){
- gr->tetrahedra.push_back((*it)->tet());
- delete *it;
- }
- else{
- delete (*it)->tet();
- delete *it;
- }
+ for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
+ if (!(*it)->isDeleted()){
+ gr->tetrahedra.push_back((*it)->tet());
+ delete *it;
}
-
+ else{
+ delete (*it)->tet();
+ delete *it;
+ }
+ }
}
-
void insertVerticesInRegion (GRegion *gr)
{
//printf("sizeof MTet4 = %d sizeof MTetrahedron %d sizeof(MVertex) %d\n",
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
index c2f7313ad5589227faa49bb1e2b7cfbf84530572..3c4852efee7e3f046938cc313f1f18774d976a41 100644
--- a/Mesh/meshGRegionLocalMeshMod.cpp
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -50,20 +50,21 @@ bool gmshBuildEdgeCavity(MTet4 *t,
MVertex **v1,MVertex **v2,
std::vector<MTet4*> &cavity,
std::vector<MTet4*> &outside,
- std::vector<MVertex*> &ring){
+ std::vector<MVertex*> &ring)
+{
cavity.clear();
ring.clear();
*v1 = t->tet()->getVertex(edges[iLocalEdge][0]);
*v2 = t->tet()->getVertex(edges[iLocalEdge][1]);
- MVertex *lastinring = t->tet()->getVertex(edges[5-iLocalEdge][0]);
+ MVertex *lastinring = t->tet()->getVertex(edges[5 - iLocalEdge][0]);
ring.push_back(lastinring);
cavity.push_back(t);
while (1){
- MVertex *ov1 = t->tet()->getVertex(edges[5-iLocalEdge][0]);
- MVertex *ov2 = t->tet()->getVertex(edges[5-iLocalEdge][1]);
+ MVertex *ov1 = t->tet()->getVertex(edges[5 - iLocalEdge][0]);
+ MVertex *ov2 = t->tet()->getVertex(edges[5 - iLocalEdge][1]);
int K = ov1 == lastinring ? 1 : 0;
lastinring = ov1 == lastinring ? ov2 : ov1;
// look in the 2 faces sharing this edge the one that has vertex
@@ -77,25 +78,24 @@ bool gmshBuildEdgeCavity(MTet4 *t,
else if (faces[iFace2][0] == edges[5-iLocalEdge][K] ||
faces[iFace2][1] == edges[5-iLocalEdge][K] ||
faces[iFace2][2] == edges[5-iLocalEdge][K] ) iFace = iFace2;
- else {printf("error of connexion\n");throw;}
+ else { Msg(GERROR, "Error of connexion"); throw; }
t=t->getNeigh(iFace);
- if (!t)return false;
- if (t->isDeleted())throw;
- if (t==cavity[0]) break;
+ if (!t) return false;
+ if (t->isDeleted()) throw;
+ if (t == cavity[0]) break;
ring.push_back(lastinring);
cavity.push_back(t);
iLocalEdge = -1;
- for (int i=0;i<6;i++)
- {
- MVertex *a = t->tet()->getVertex(edges[i][0]);
- MVertex *b = t->tet()->getVertex(edges[i][1]);
- if ((a == *v1 && b == *v2) || (a == *v2 && b == *v1)){
- iLocalEdge = i;
- break;
- }
- }
+ for (int i = 0; i < 6; i++){
+ MVertex *a = t->tet()->getVertex(edges[i][0]);
+ MVertex *b = t->tet()->getVertex(edges[i][1]);
+ if ((a == *v1 && b == *v2) || (a == *v2 && b == *v1)){
+ iLocalEdge = i;
+ break;
+ }
+ }
if (iLocalEdge == -1){
- printf("ERROR : loc = %d\n",iLocalEdge);
+ Msg(GERROR, "loc = %d", iLocalEdge);
throw;
}
}
@@ -103,9 +103,6 @@ bool gmshBuildEdgeCavity(MTet4 *t,
return true;
}
-
-
-// return 1 if it is closed and
typedef struct {
int nbr_triangles ; /* number of different triangles */
int (*triangles)[3] ; /* triangles array */
@@ -114,7 +111,6 @@ typedef struct {
int (*trianguls)[5] ; /* retriangulations array */
} SwapPattern ;
-
void BuildSwapPattern3(SwapPattern *sc)
{
static int trgl[][3] = { {0,1,2} };
@@ -204,16 +200,16 @@ void BuildSwapPattern7(SwapPattern *sc)
bool gmshEdgeSwap (std::vector<MTet4 *> &newTets,
MTet4 *tet,
int iLocalEdge,
- const gmshQualityMeasure4Tet &cr){
+ const gmshQualityMeasure4Tet &cr)
+{
std::vector<MTet4*> cavity;
std::vector<MTet4*> outside;
std::vector<MVertex*> ring;
- MVertex *v1,*v2;
- // printf("swapping an edge\n");
+ MVertex *v1, *v2;
- bool closed = gmshBuildEdgeCavity ( tet,iLocalEdge,&v1,&v2,cavity,outside,ring);
+ bool closed = gmshBuildEdgeCavity(tet, iLocalEdge, &v1, &v2, cavity, outside, ring);
- if (!closed)return false;
+ if (!closed) return false;
double volumeRef = 0.0;
double tetQualityRef = 1;
@@ -227,46 +223,46 @@ bool gmshEdgeSwap (std::vector<MTet4 *> &newTets,
SwapPattern sp;
switch (ring.size()){
- case 3 : BuildSwapPattern3 (&sp); break;
- case 4 : BuildSwapPattern4 (&sp); break;
- case 5 : BuildSwapPattern5 (&sp); break;
- case 6 : BuildSwapPattern6 (&sp); break;
- case 7 : BuildSwapPattern7 (&sp); break;
+ case 3 : BuildSwapPattern3(&sp); break;
+ case 4 : BuildSwapPattern4(&sp); break;
+ case 5 : BuildSwapPattern5(&sp); break;
+ case 6 : BuildSwapPattern6(&sp); break;
+ case 7 : BuildSwapPattern7(&sp); break;
default : return false;
}
- // compute qualities of all tets that appear in the patterns
- double tetQuality1[100],tetQuality2[100];
- double volume1[100],volume2[100];
- for (int i=0;i<sp.nbr_triangles;i++){
+ // compute qualities of all tets that appear in the patterns
+ double tetQuality1[100], tetQuality2[100];
+ double volume1[100], volume2[100];
+ for (int i = 0; i < sp.nbr_triangles; i++){
int p1 = sp.triangles[i][0];
int p2 = sp.triangles[i][1];
int p3 = sp.triangles[i][2];
- tetQuality1[i] = qmTet ( ring[p1], ring[p2], ring[p3], v1 , cr , &(volume1[i]));
- tetQuality2[i] = qmTet ( ring[p1], ring[p2], ring[p3], v2 , cr , &(volume2[i]));
+ tetQuality1[i] = qmTet(ring[p1], ring[p2], ring[p3], v1, cr, &(volume1[i]));
+ tetQuality2[i] = qmTet(ring[p1], ring[p2], ring[p3], v2, cr, &(volume2[i]));
}
- // look for the best triangulation, i.e. the one
- // that maximize the minimum element quality
+ // look for the best triangulation, i.e. the one that maximize the
+ // minimum element quality
double minQuality[100];
// for all triangulations
- for (int i=0;i<sp.nbr_trianguls;i++){
+ for (int i = 0; i < sp.nbr_trianguls; i++){
// for all triangles in a triangulation
minQuality[i] = 1;
- double volume=0;
- for (int j=0;j<sp.nbr_triangles_2;j++){
+ double volume = 0;
+ for (int j = 0; j < sp.nbr_triangles_2; j++){
int iT = sp.trianguls[i][j];
- minQuality[i] = std::min (minQuality[i],tetQuality1[iT]);
- minQuality[i] = std::min (minQuality[i],tetQuality2[iT]);
- volume += (volume1[iT] + volume2[iT] );
+ minQuality[i] = std::min(minQuality[i], tetQuality1[iT]);
+ minQuality[i] = std::min(minQuality[i], tetQuality2[iT]);
+ volume += (volume1[iT] + volume2[iT]);
}
- // printf("config %3d qual %12.5E volume %12.5E\n",i,minQuality[i],volume);
- if (fabs(volume-volumeRef) > 1.e-10 * (volume+volumeRef))minQuality[i] = -1;
+ // printf("config %3d qual %12.5E volume %12.5E\n",i,minQuality[i],volume);
+ if (fabs(volume-volumeRef) > 1.e-10 * (volume+volumeRef)) minQuality[i] = -1;
}
int iBest = 0;
double best = -1.0;
- for (int i=0;i<sp.nbr_trianguls;i++){
+ for (int i = 0; i < sp.nbr_trianguls; i++){
if(minQuality[i] > best){
best = minQuality[i];
iBest = i;
@@ -277,11 +273,12 @@ bool gmshEdgeSwap (std::vector<MTet4 *> &newTets,
if (best <= tetQualityRef) return false;
// we have the best configuration, so we swap
- // printf("outside size = %d\n",outside.size());
+ // printf("outside size = %d\n",outside.size());
- // printf("a swap with %d tets reconnect %d tets cavity %d tets\n",ring.size(),outside.size(),cavity.size());
+ // printf("a swap with %d tets reconnect %d tets cavity %d tets\n",
+ // ring.size(),outside.size(),cavity.size());
- for (int j=0;j<sp.nbr_triangles_2;j++){
+ for (int j = 0; j < sp.nbr_triangles_2; j++){
int iT = sp.trianguls[iBest][j];
int p1 = sp.triangles[iT][0];
int p2 = sp.triangles[iT][1];
@@ -289,16 +286,10 @@ bool gmshEdgeSwap (std::vector<MTet4 *> &newTets,
MVertex *pv1 = ring[p1];
MVertex *pv2 = ring[p2];
MVertex *pv3 = ring[p3];
- MTetrahedron *tr1 = new MTetrahedron ( pv1,
- pv2,
- pv3,
- v1);
- MTetrahedron *tr2 = new MTetrahedron ( pv3,
- pv2,
- pv1,
- v2);
- MTet4 *t41 = new MTet4 ( tr1 , tetQuality1[iT]);
- MTet4 *t42 = new MTet4 ( tr2 , tetQuality2[iT]);
+ MTetrahedron *tr1 = new MTetrahedron(pv1, pv2, pv3, v1);
+ MTetrahedron *tr2 = new MTetrahedron (pv3, pv2, pv1, v2);
+ MTet4 *t41 = new MTet4(tr1, tetQuality1[iT]);
+ MTet4 *t42 = new MTet4(tr2, tetQuality2[iT]);
t41->setOnWhat(cavity[0]->onWhat());
t42->setOnWhat(cavity[0]->onWhat());
outside.push_back(t41);
@@ -309,37 +300,32 @@ bool gmshEdgeSwap (std::vector<MTet4 *> &newTets,
for(unsigned int i = 0; i < cavity.size(); i++) cavity[i]->setDeleted(true);
- connectTets ( outside );
+ connectTets(outside);
return true;
}
-bool gmshEdgeSplit (std::vector<MTet4 *> &newTets,
- MTet4 *tet,
- MVertex *newVertex,
- int iLocalEdge,
- const gmshQualityMeasure4Tet &cr){
+bool gmshEdgeSplit(std::vector<MTet4 *> &newTets,
+ MTet4 *tet,
+ MVertex *newVertex,
+ int iLocalEdge,
+ const gmshQualityMeasure4Tet &cr)
+{
std::vector<MTet4*> cavity;
std::vector<MTet4*> outside;
std::vector<MVertex*> ring;
- MVertex *v1,*v2;
+ MVertex *v1, *v2;
- bool closed = gmshBuildEdgeCavity ( tet,iLocalEdge,&v1,&v2,cavity,outside,ring);
- if (!closed)return false;
+ bool closed = gmshBuildEdgeCavity(tet, iLocalEdge, &v1, &v2, cavity, outside, ring);
+ if (!closed) return false;
for(unsigned int j = 0; j < ring.size(); j++){
MVertex *pv1 = ring[j];
- MVertex *pv2 = ring[(j+1)%ring.size()];
- MTetrahedron *tr1 = new MTetrahedron ( pv1,
- pv2,
- newVertex,
- v1);
- MTetrahedron *tr2 = new MTetrahedron ( newVertex,
- pv2,
- pv1,
- v2);
- MTet4 *t41 = new MTet4 ( tr1 , cr) ;
- MTet4 *t42 = new MTet4 ( tr2 , cr);
+ MVertex *pv2 = ring[(j + 1) % ring.size()];
+ MTetrahedron *tr1 = new MTetrahedron(pv1, pv2, newVertex, v1);
+ MTetrahedron *tr2 = new MTetrahedron(newVertex, pv2, pv1, v2);
+ MTet4 *t41 = new MTet4(tr1, cr);
+ MTet4 *t42 = new MTet4(tr2, cr);
t41->setOnWhat(cavity[0]->onWhat());
t42->setOnWhat(cavity[0]->onWhat());
outside.push_back(t41);
@@ -356,30 +342,30 @@ bool gmshEdgeSplit (std::vector<MTet4 *> &newTets,
}
// swap a face i.e. remove a face shared by 2 tets
-bool gmshFaceSwap (std::vector<MTet4 *> &newTets,
- MTet4 *t1,
- int iLocalFace,
- const gmshQualityMeasure4Tet &cr){
-
+bool gmshFaceSwap(std::vector<MTet4 *> &newTets,
+ MTet4 *t1,
+ int iLocalFace,
+ const gmshQualityMeasure4Tet &cr)
+{
MTet4 *t2 = t1->getNeigh(iLocalFace);
- if (!t2)return false;
- if (t1->onWhat() != t2->onWhat())return false;
+ if (!t2) return false;
+ if (t1->onWhat() != t2->onWhat()) return false;
MVertex *v1 = t1->tet()->getVertex(vnofaces[iLocalFace]);
MVertex *f1 = t1->tet()->getVertex(faces[iLocalFace][0]);
MVertex *f2 = t1->tet()->getVertex(faces[iLocalFace][1]);
MVertex *f3 = t1->tet()->getVertex(faces[iLocalFace][2]);
MVertex *v2 = 0;
- for (int i=0;i<4;i++){
+ for (int i = 0; i < 4; i++){
MVertex *v = t2->tet()->getVertex(i);
- if (v != f1 && v != f2 && v != f3 ){
+ if (v != f1 && v != f2 && v != f3){
v2 = v;
break;
}
}
- if (!v2){printf("coucou\n");throw;}
+ if (!v2) throw;
- // printf("%p %p -- %p %p %p\n",v1,v2,f1,f2,f3);
+ // printf("%p %p -- %p %p %p\n",v1,v2,f1,f2,f3);
double vol1 = fabs(t1->tet()->getVolume());
double q1 = t1->getQuality();
@@ -387,17 +373,17 @@ bool gmshFaceSwap (std::vector<MTet4 *> &newTets,
double q2 = t2->getQuality();
double vol3;
- double q3 = qmTet(f1,f2,v1,v2,cr,&vol3);
+ double q3 = qmTet(f1, f2, v1, v2, cr, &vol3);
double vol4;
- double q4 = qmTet(f2,f3,v1,v2,cr,&vol4);
+ double q4 = qmTet(f2, f3, v1, v2, cr, &vol4);
double vol5;
- double q5 = qmTet(f3,f1,v1,v2,cr,&vol5);
+ double q5 = qmTet(f3, f1, v1, v2, cr, &vol5);
- if (fabs ( vol1 + vol2 - vol3 - vol4 - vol5 ) > 1.e-10 * (vol1+vol2)) return false;
- if (std::min(q1,q2) > std::min(std::min(q3,q4),q5)) return false;
- // printf("%g %g %g\n",vol1 + vol2, vol3 + vol4 + vol5,vol1 + vol2 - vol3 - vol4 - vol5);
- // printf("qs = %g %g vs %g %g %g\n",q1,q2,q3,q4,q5);
+ if (fabs(vol1 + vol2 - vol3 - vol4 - vol5) > 1.e-10 * (vol1 + vol2)) return false;
+ if (std::min(q1, q2) > std::min(std::min(q3, q4), q5)) return false;
+ // printf("%g %g %g\n",vol1 + vol2, vol3 + vol4 + vol5,vol1 + vol2 - vol3 - vol4 - vol5);
+ // printf("qs = %g %g vs %g %g %g\n",q1,q2,q3,q4,q5);
std::vector<MTet4*> outside;
for(int i = 0; i < 4; i++){
@@ -424,12 +410,12 @@ bool gmshFaceSwap (std::vector<MTet4 *> &newTets,
t1->setDeleted(true);
t2->setDeleted(true);
- MTetrahedron *tr1 = new MTetrahedron ( f1,f2,v1,v2);
- MTetrahedron *tr2 = new MTetrahedron ( f2,f3,v1,v2);
- MTetrahedron *tr3 = new MTetrahedron ( f3,f1,v1,v2);
- MTet4 *t41 = new MTet4 ( tr1 , q3);
- MTet4 *t42 = new MTet4 ( tr2 , q4);
- MTet4 *t43 = new MTet4 ( tr3 , q5);
+ MTetrahedron *tr1 = new MTetrahedron(f1, f2, v1, v2);
+ MTetrahedron *tr2 = new MTetrahedron(f2, f3, v1, v2);
+ MTetrahedron *tr3 = new MTetrahedron(f3, f1, v1, v2);
+ MTet4 *t41 = new MTet4(tr1, q3);
+ MTet4 *t42 = new MTet4(tr2, q4);
+ MTet4 *t43 = new MTet4(tr3, q5);
t41->setOnWhat(t1->onWhat());
t42->setOnWhat(t1->onWhat());
t43->setOnWhat(t1->onWhat());
@@ -439,29 +425,29 @@ bool gmshFaceSwap (std::vector<MTet4 *> &newTets,
newTets.push_back(t41);
newTets.push_back(t42);
newTets.push_back(t43);
- connectTets ( outside );
+ connectTets(outside);
return true;
}
-
-void gmshBuildVertexCavity_recur ( MTet4 *t,
- MVertex *v,
- std::vector<MTet4*> &cavity){
- // if (recur > 20)printf("oufti %d\n",recur);
+void gmshBuildVertexCavity_recur(MTet4 *t,
+ MVertex *v,
+ std::vector<MTet4*> &cavity)
+{
+ // if (recur > 20)printf("oufti %d\n",recur);
if(t->isDeleted()){
Msg(FATAL,"a deleted triangle is a neighbor of a non deleted triangle");
}
- int iV=-1;
- for (int i=0; i<4; i++){
+ int iV = -1;
+ for (int i = 0; i < 4; i++){
if (t->tet()->getVertex(i) == v){
iV = i;
break;
}
}
- if (iV==-1){
- Msg(FATAL,"trying to build a cavity of tets for a vertex that does not belong to this tet");
+ if (iV == -1){
+ Msg(FATAL, "trying to build a cavity of tets for a vertex that does not belong to this tet");
}
- for (int i=0; i<3; i++){
+ for (int i = 0; i < 3; i++){
MTet4 *neigh = t->getNeigh(vFac[iV][i]);
if (neigh){
bool found = false;
@@ -479,27 +465,25 @@ void gmshBuildVertexCavity_recur ( MTet4 *t,
}
}
-// sliver removal by compound mesh modif
-// postulate : the edge cannot be swopped
-// so we split it, and then collapse the new
-// vertex on another one (of course, not the
-// other one on the unswappable edge)
+// sliver removal by compound mesh modif postulate : the edge cannot
+// be swopped so we split it, and then collapse the new vertex on
+// another one (of course, not the other one on the unswappable edge)
// after that crap, the sliver is trashed
-bool gmshSliverRemoval ( std::vector<MTet4 *> &newTets,
- MTet4 *t,
- const gmshQualityMeasure4Tet &cr){
- // look if
+bool gmshSliverRemoval(std::vector<MTet4*> &newTets,
+ MTet4 *t,
+ const gmshQualityMeasure4Tet &cr)
+{
std::vector<MTet4*> cavity;
std::vector<MTet4*> outside;
std::vector<MVertex*> ring;
- MVertex *v1,*v2;
+ MVertex *v1, *v2;
bool isClosed[6];
int nbSwappable = 0;
int iSwappable = 0;
- for (int i=0;i<6;i++){
- isClosed[i] = gmshBuildEdgeCavity ( t,i,&v1,&v2,cavity,outside,ring);
+ for (int i = 0; i < 6; i++){
+ isClosed[i] = gmshBuildEdgeCavity(t, i, &v1, &v2, cavity, outside, ring);
if (isClosed[i]){
nbSwappable++;
iSwappable = i;
@@ -518,19 +502,18 @@ bool gmshSliverRemoval ( std::vector<MTet4 *> &newTets,
// if it does not work, split, smooth and collapse
MVertex *v1 = t->tet()->getVertex(edges[iSwappable][0]);
MVertex *v2 = t->tet()->getVertex(edges[iSwappable][1]);
- MVertex *newv = new MVertex (0.5*(v1->x()+v2->x()),
- 0.5*(v1->y()+v2->y()),
- 0.5*(v1->z()+v2->z()),t->onWhat());
+ MVertex *newv = new MVertex (0.5 * (v1->x() + v2->x()),
+ 0.5 * (v1->y() + v2->y()),
+ 0.5 * (v1->z() + v2->z()), t->onWhat());
t->onWhat()->mesh_vertices.push_back(newv);
- if (!gmshEdgeSplit(newTets,t,newv,iSwappable,QMTET_ONE))return false;
- for (int i=0;i<4;i++){
- if (newTets[newTets.size()-1]->tet()->getVertex(i) == newv)
- {
- gmshSmoothVertex(newTets[newTets.size()-1], i,cr);
- gmshSmoothVertexOptimize (newTets[newTets.size()-1], i,cr);
- }
- }
+ if (!gmshEdgeSplit(newTets, t, newv, iSwappable, QMTET_ONE)) return false;
+ for (int i = 0; i < 4; i++){
+ if (newTets[newTets.size() - 1]->tet()->getVertex(i) == newv){
+ gmshSmoothVertex(newTets[newTets.size() - 1], i, cr);
+ gmshSmoothVertexOptimize (newTets[newTets.size() - 1], i, cr);
+ }
+ }
for (unsigned int i = 0; i < newTets.size(); i++){
MTet4 *new_t = newTets[i];
@@ -549,38 +532,37 @@ bool gmshSliverRemoval ( std::vector<MTet4 *> &newTets,
}
}
}
-
return true;
}
else{
- for (int i=0;i<4;i++){
- gmshSmoothVertex(t, i,cr);
- gmshSmoothVertexOptimize (t, i,cr);
+ for (int i = 0; i < 4; i++){
+ gmshSmoothVertex(t, i, cr);
+ gmshSmoothVertexOptimize(t, i, cr);
}
}
return false;
}
-bool gmshCollapseVertex ( std::vector<MTet4 *> &newTets,
- MTet4 *t,
- int iVertex,
- int iTarget,
- const gmshQualityMeasure4Tet &cr,
- const gmshLocalMeshModAction action,
- double *minQual){
-
- if(t->isDeleted())throw;
+bool gmshCollapseVertex(std::vector<MTet4 *> &newTets,
+ MTet4 *t,
+ int iVertex,
+ int iTarget,
+ const gmshQualityMeasure4Tet &cr,
+ const gmshLocalMeshModAction action,
+ double *minQual)
+{
+ if(t->isDeleted()) throw;
MVertex *v = t->tet()->getVertex(iVertex);
MVertex *tg = t->tet()->getVertex(iTarget);
- if(v->onWhat()->dim() < 3)return false;
- if(tg->onWhat()->dim() < 3)return false;
+ if(v->onWhat()->dim() < 3) return false;
+ if(tg->onWhat()->dim() < 3) return false;
std::vector<MTet4*> cavity_v;
std::vector<MTet4*> outside;
cavity_v.push_back(t);
- gmshBuildVertexCavity_recur (t,v,cavity_v);
+ gmshBuildVertexCavity_recur(t, v, cavity_v);
std::vector<MTet4*> toDelete;
std::vector<MTet4*> toUpdate;
@@ -605,7 +587,7 @@ bool gmshCollapseVertex ( std::vector<MTet4 *> &newTets,
v->y() = tg->y();
v->z() = tg->z();
- double volume_update=0;
+ double volume_update = 0;
double worstAfter = 1.0;
double newQuals[2000];
@@ -614,24 +596,24 @@ bool gmshCollapseVertex ( std::vector<MTet4 *> &newTets,
double vv;
newQuals[i] = qmTet(toUpdate[i]->tet(),cr,&vv);
worstAfter = std::min(worstAfter,newQuals[i]);
- volume_update+=vv;
+ volume_update += vv;
}
- // printf("%12.5E %12.5E %12.5E %12.5E %d\n",volume,volume_update,worstAfter,worst,toUpdate.size());
+ // printf("%12.5E %12.5E %12.5E %12.5E %d\n",
+ // volume,volume_update,worstAfter,worst,toUpdate.size());
- if (fabs(volume-volume_update) > 1.e-10 * volume || worstAfter < worst)
- {
- v->x() = x;
- v->y() = y;
- v->z() = z;
- return false;
- }
+ if (fabs(volume-volume_update) > 1.e-10 * volume || worstAfter < worst){
+ v->x() = x;
+ v->y() = y;
+ v->z() = z;
+ return false;
+ }
if (action == GMSH_EVALONLY){
*minQual = worstAfter;
return true;
}
// ok we collapse
- computeNeighboringTetsOfACavity (cavity_v,outside);
+ computeNeighboringTetsOfACavity(cavity_v, outside);
for (unsigned int i = 0; i < toUpdate.size(); i++){
MTetrahedron *tr1 = new MTetrahedron
(toUpdate[i]->tet()->getVertex(0) == v ? tg : toUpdate[i]->tet()->getVertex(0),
@@ -651,39 +633,38 @@ bool gmshCollapseVertex ( std::vector<MTet4 *> &newTets,
return true;
}
-bool gmshSmoothVertex( MTet4 *t,
- int iVertex,
- const gmshQualityMeasure4Tet &cr){
-
- if(t->isDeleted())throw;
- if(t->tet()->getVertex(iVertex)->onWhat()->dim() < 3)return false;
+bool gmshSmoothVertex(MTet4 *t,
+ int iVertex,
+ const gmshQualityMeasure4Tet &cr)
+{
+ if(t->isDeleted()) throw;
+ if(t->tet()->getVertex(iVertex)->onWhat()->dim() < 3) return false;
std::vector<MTet4*> cavity;
cavity.push_back(t);
- gmshBuildVertexCavity_recur (t,t->tet()->getVertex(iVertex),cavity);
+ gmshBuildVertexCavity_recur(t, t->tet()->getVertex(iVertex), cavity);
- double xcg=0,ycg=0,zcg=0;
- double vTot=0;
+ double xcg = 0, ycg = 0, zcg = 0;
+ double vTot = 0;
double worst = 1.0;
for(unsigned int i = 0 ; i < cavity.size(); i++){
double volume = fabs(cavity[i]->tet()->getVolume());
double q = cavity[i]->getQuality();
worst = std::min(worst,q);
- xcg += 0.25*(cavity[i]->tet()->getVertex(0)->x()+
- cavity[i]->tet()->getVertex(1)->x()+
- cavity[i]->tet()->getVertex(2)->x()+
- cavity[i]->tet()->getVertex(3)->x())*volume;
- ycg += 0.25*(cavity[i]->tet()->getVertex(0)->y()+
- cavity[i]->tet()->getVertex(1)->y()+
- cavity[i]->tet()->getVertex(2)->y()+
- cavity[i]->tet()->getVertex(3)->y())*volume;
- zcg += 0.25*(cavity[i]->tet()->getVertex(0)->z()+
- cavity[i]->tet()->getVertex(1)->z()+
- cavity[i]->tet()->getVertex(2)->z()+
- cavity[i]->tet()->getVertex(3)->z())*volume;
+ xcg += 0.25 * (cavity[i]->tet()->getVertex(0)->x() +
+ cavity[i]->tet()->getVertex(1)->x() +
+ cavity[i]->tet()->getVertex(2)->x() +
+ cavity[i]->tet()->getVertex(3)->x()) * volume;
+ ycg += 0.25 * (cavity[i]->tet()->getVertex(0)->y() +
+ cavity[i]->tet()->getVertex(1)->y() +
+ cavity[i]->tet()->getVertex(2)->y() +
+ cavity[i]->tet()->getVertex(3)->y()) * volume;
+ zcg += 0.25 * (cavity[i]->tet()->getVertex(0)->z() +
+ cavity[i]->tet()->getVertex(1)->z() +
+ cavity[i]->tet()->getVertex(2)->z() +
+ cavity[i]->tet()->getVertex(3)->z()) * volume;
vTot += volume;
-
}
xcg /= (vTot);
ycg /= (vTot);
@@ -704,14 +685,14 @@ bool gmshSmoothVertex( MTet4 *t,
double volume;
newQuals[i] = qmTet(cavity[i]->tet(),cr,&volume);
volumeAfter += volume;
- worstAfter =std::min(worstAfter,newQuals[i]);
+ worstAfter = std::min(worstAfter,newQuals[i]);
}
if (fabs(volumeAfter-vTot) > 1.e-10 * vTot || worstAfter < worst){
t->tet()->getVertex(iVertex)->x() = x;
t->tet()->getVertex(iVertex)->y() = y;
t->tet()->getVertex(iVertex)->z() = z;
- return false;//gmshSmoothVertexOptimize ( t,iVertex,cr);
+ return false;//gmshSmoothVertexOptimize (t, iVertex, cr);
}
else{
// restore new quality
@@ -724,11 +705,13 @@ bool gmshSmoothVertex( MTet4 *t,
struct smoothVertexData3D{
MVertex *v;
- std::vector < MTet4 * >ts;
+ std::vector<MTet4 *> ts;
double LC;
};
-double smoothing_objective_function_3D(double X, double Y, double Z, MVertex *v, std::vector < MTet4 * >&ts){
+double smoothing_objective_function_3D(double X, double Y, double Z,
+ MVertex *v, std::vector<MTet4 *> &ts)
+{
const double oldX = v->x();
const double oldY = v->y();
const double oldZ = v->z();
@@ -738,9 +721,9 @@ double smoothing_objective_function_3D(double X, double Y, double Z, MVertex *v,
std::vector < MTet4 * >::iterator it = ts.begin();
std::vector < MTet4 * >::iterator ite = ts.end();
- double qMin=1,vol;
- while(it != ite) {
- qMin = std::min(qmTet((*it)->tet(),QMTET_2,&vol),qMin);
+ double qMin = 1, vol;
+ while(it != ite){
+ qMin = std::min(qmTet((*it)->tet(), QMTET_2, &vol), qMin);
++it;
}
v->x() = oldX;
@@ -752,46 +735,48 @@ double smoothing_objective_function_3D(double X, double Y, double Z, MVertex *v,
void deriv_smoothing_objective_function_3D(double X, double Y, double Z,
double &F,
double &dFdX, double &dFdY, double &dFdZ,
- void *data){
+ void *data)
+{
smoothVertexData3D *svd = (smoothVertexData3D*)data;
MVertex *v = svd->v;
const double LARGE = svd->LC*1.e5;
const double SMALL = 1./LARGE;
- F = smoothing_objective_function_3D(X,Y,Z,v,svd->ts);
- double F_X = smoothing_objective_function_3D(X+SMALL,Y,Z,v,svd->ts);
- double F_Y = smoothing_objective_function_3D(X,Y+SMALL,Z,v,svd->ts);
- double F_Z = smoothing_objective_function_3D(X,Y,Z+SMALL,v,svd->ts);
- dFdX = (F_X-F)*LARGE;
- dFdY = (F_Y-F)*LARGE;
- dFdZ = (F_Z-F)*LARGE;
+ F = smoothing_objective_function_3D(X, Y, Z, v, svd->ts);
+ double F_X = smoothing_objective_function_3D(X + SMALL, Y, Z, v, svd->ts);
+ double F_Y = smoothing_objective_function_3D(X, Y + SMALL, Z, v, svd->ts);
+ double F_Z = smoothing_objective_function_3D(X, Y, Z + SMALL, v, svd->ts);
+ dFdX = (F_X - F) * LARGE;
+ dFdY = (F_Y - F) * LARGE;
+ dFdZ = (F_Z - F) * LARGE;
}
-double smooth_obj_3D(double X, double Y, double Z, void *data){
+double smooth_obj_3D(double X, double Y, double Z, void *data)
+{
smoothVertexData3D *svd = (smoothVertexData3D*)data;
- return smoothing_objective_function_3D(X,Y,Z,svd->v,svd->ts);
+ return smoothing_objective_function_3D(X, Y, Z, svd->v, svd->ts);
}
-
-bool gmshSmoothVertexOptimize ( MTet4 *t,
- int iVertex,
- const gmshQualityMeasure4Tet &cr){
+bool gmshSmoothVertexOptimize(MTet4 *t,
+ int iVertex,
+ const gmshQualityMeasure4Tet &cr)
+{
+ if(t->tet()->getVertex(iVertex)->onWhat()->dim() < 3) return false;
- if(t->tet()->getVertex(iVertex)->onWhat()->dim() < 3)return false;
-
smoothVertexData3D vd;
vd.ts.push_back(t);
vd.v = t->tet()->getVertex(iVertex);
vd.LC = 1.0; // WRONG
- gmshBuildVertexCavity_recur (t,t->tet()->getVertex(iVertex),vd.ts);
+ gmshBuildVertexCavity_recur(t, t->tet()->getVertex(iVertex), vd.ts);
- double xopti=vd.v->x();
- double yopti=vd.v->y();
- double zopti=vd.v->z();
+ double xopti = vd.v->x();
+ double yopti = vd.v->y();
+ double zopti = vd.v->z();
double val;
- minimize_3 ( smooth_obj_3D, deriv_smoothing_objective_function_3D, &vd, 4, xopti,yopti,zopti,val);
+ minimize_3(smooth_obj_3D, deriv_smoothing_objective_function_3D, &vd, 4,
+ xopti, yopti, zopti, val);
- double vTot=0;
+ double vTot = 0;
for(unsigned int i = 0; i < vd.ts.size(); i++){
double volume = fabs(vd.ts[i]->tet()->getVolume());
@@ -831,9 +816,6 @@ bool gmshSmoothVertexOptimize ( MTet4 *t,
}
}
-
-
-
// Edge split sets ...
// Here, we only build a list of tets that are a subdivision
// of the given