diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 18b57cd21ae4db07594f6d8242ce48caae14d1c7..a703a1587922ae603f892b7cb8dbf35c9e46fe8f 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -9,7 +9,6 @@
#include "GmshMessage.h"
#include "robustPredicates.h"
#include "OS.h"
-//#include "BackgroundMesh.h"
#include "meshGRegion.h"
#include "meshGRegionLocalMeshMod.h"
#include "meshGRegionDelaunayInsertion.h"
@@ -18,7 +17,7 @@
#include "MTriangle.h"
#include "Numeric.h"
#include "Context.h"
-#include "HilbertCurve.h"
+#include "delaunay3d.h"
int MTet4::radiusNorm = 2;
static double LIMIT_ = 1;
@@ -45,6 +44,7 @@ static bool isActive(MTet4 *t, double limit_, int &active)
return false;
}
+
int MTet4::inCircumSphere(const double *p) const
{
double pa[3] = {base->getVertex(0)->x(),
@@ -64,6 +64,7 @@ int MTet4::inCircumSphere(const double *p) const
return (result > 0) ? 1 : 0;
}
+
static int faces[4][3] = {{0,1,2}, {0,2,3}, {0,3,1}, {1,3,2}};
struct faceXtet{
@@ -143,30 +144,6 @@ void connectTets_vector2(std::vector<MTet4*> &t, std::vector<faceXtet> &conn)
}
-// template <class ITER>
-// void connectTets_vector(ITER beg, ITER end)
-// {
-// // std::set<faceXtet> conn;
-// std::vector<faceXtet> conn;
-// while (beg != end){
-// if (!(*beg)->isDeleted()){
-// for (int i = 0; i < 4; i++){
-// faceXtet fxt(*beg, i);
-// std::vector<faceXtet>::iterator found = std::find(conn.begin(), conn.end(), fxt);
-// // std::set<faceXtet>::iterator found = conn.find(fxt);
-// if (found == conn.end())
-// conn.push_back(fxt);
-// // conn.insert(fxt);
-// else if (found->t1 != *beg){
-// found->t1->setNeigh(found->i1, *beg);
-// (*beg)->setNeigh(i, found->t1);
-// }
-// }
-// }
-// ++beg;
-// }
-// }
-
template <class ITER>
void connectTets(ITER beg, ITER end, std::set<MFace, Less_Face> *allEmbeddedFaces = 0)
{
@@ -371,40 +348,6 @@ void recurFindCavity(std::vector<faceXtet> & shell,
}
}
-
-// void nonrecurFindCavity(std::vector<faceXtet> & shell,
-// std::vector<MTet4*> & cavity,
-// MVertex *v ,
-// MTet4 *t,
-// std::stack<MTet4*> &_stack)
-// {
-
-// _stack.push(t);
-// while(!_stack.empty()){
-// t = _stack.top();
-// _stack.pop();
-// if (!t->isDeleted()){
-// t->setDeleted(true);
-// cavity.push_back(t);
-
-// for (int i = 0; i < 4; i++){
-// MTet4 *neigh = t->getNeigh(i) ;
-// faceXtet fxt (t, i);
-// if (!neigh)
-// shell.push_back(fxt);
-// else if (!neigh->isDeleted()){
-// int circ = neigh->inCircumSphere(v);
-// if (circ && (neigh->onWhat() == t->onWhat()))
-// _stack.push(neigh);
-// else
-// shell.push_back(fxt);
-// }
-// }
-// }
-// }
-// // printf("cavity size %d\n",cavity.size());
-// }
-
void printTets (const char *fn, std::list<MTet4*> &cavity, bool force = false )
{
FILE *f = Fopen (fn,"w");
@@ -626,52 +569,12 @@ static void setLcs(MTetrahedron *t, std::map<MVertex*, double> &vSizes,
}
}
-// 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()){
-// }
-// }
-// }
-
-// 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(),
- (*fit)->triangles.end(),
- tri, cmp);
- if(found){
- *gfound = *fit;
- return true;
- }
- ++fit;
- }
- return false;
-}
-
GRegion *getRegionFromBoundingFaces(GModel *model,
std::set<GFace *> &faces_bound)
{
GModel::riter git = model->firstRegion();
- // for (std::set<GFace *>::iterator it = faces_bound.begin();
- // it != faces_bound.end(); ++it){
- // printf(" %d",(*it)->tag());
- // }
- // printf("\n");
while (git != model->lastRegion()){
std::list <GFace *> _faces = (*git)->faces();
- // printf("region %d %d faces\n",(*git)->tag(),_faces.size());
- // for (std::list<GFace *>::iterator it = _faces.begin(); it != _faces.end(); ++it){
- // printf(" %d",(*it)->tag());
- // }
- // printf("\n");
-
if(_faces.size() == faces_bound.size()){
bool ok = true;
for (std::list<GFace *>::iterator it = _faces.begin(); it != _faces.end(); ++it){
@@ -684,38 +587,6 @@ GRegion *getRegionFromBoundingFaces(GModel *model,
return 0;
}
-// 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::Error("a tet is not connected by a boundary face");
-// 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);
-// }
-// }
-// }
-// for (int i = 0; i < 4; i++){
-// if (!FF[i])
-// recur_classify(t->getNeigh(i), theRegion, faces_bound, bidon, model, search );
-// }
-// }
void non_recursive_classify(MTet4 *t, std::list<MTet4*> &theRegion,
std::set<GFace*> &faces_bound, GRegion *bidon,
@@ -947,6 +818,12 @@ void adaptMeshGRegion::operator () (GRegion *gr)
//template <class CONTAINER, class DATA>
void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
{
+ // well, this should not be true !!!
+ // if (gr->hexahedra.size() ||
+ // gr->prisms.size() ||
+ // gr->pyramids.size())return;
+ if (!gr->tetrahedra.size())return;
+
typedef std::list<MTet4 *> CONTAINER ;
CONTAINER allTets;
@@ -1053,18 +930,6 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
}
// printf("coucou\n");
- if (0 && !newTets.size()){
- int nbSlivers = 0;
- int nbSliversWeCanDoSomething = 0;
- for(unsigned int i = 0; i < illegals.size(); i++)
- if(!(illegals[i]->isDeleted())){
- if(sliverRemoval(newTets, illegals[i], qm))
- nbSliversWeCanDoSomething++;
- nbSlivers++;
- }
- Msg::Info("Opti : %d Sliver Removals", nbSliversWeCanDoSomething);
- }
-
if (!newTets.size()){
break;
}
@@ -1235,6 +1100,7 @@ static void memoryCleanup(MTet4Factory &myFactory, std::set<MTet4*, compareTet4P
void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
{
+
//printf("sizeof MTet4 = %d sizeof MTetrahedron %d sizeof(MVertex) %d\n",
// sizeof(MTet4), sizeof(MTetrahedron), sizeof(MVertex));
@@ -1245,6 +1111,7 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
std::set<MTet4*, compareTet4Ptr> &allTets = myFactory.getAllTets();
int NUM = 0;
+ // printTets ("before.pos", cavity, true);
{ // leave this in a block so the map gets deallocated directly
std::map<MVertex*, double> vSizesMap;
@@ -1332,8 +1199,8 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
double t1 = Cpu();
while(1){
- // break;
- if (ITER > maxVert)break;
+ if (COUNT_MISS_2 > 100000)break;
+ if (ITER >= maxVert)break;
if(allTets.empty()){
Msg::Error("No tetrahedra in region %d %d", gr->tag(), allTets.size());
break;
@@ -1429,12 +1296,30 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
}
}
else{
- // printf("coucou 2 % cavity size %d\n",ITER,cavity.size());
- // for (std::list<MTet4*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc){
- // MTetrahedron *toto = (*itc)->tet();
- // toto->xyz2uvw(center,uvw);
- // printf("point outside %12.5E %12.5E %12.5E %12.5E\n",uvw[0], uvw[1], uvw[2],1-uvw[0]-uvw[1]-uvw[2]);
- // }
+ // printf("%d %d %d %d\n",worst->tet()->getVertex(0)->getNum()
+ // ,worst->tet()->getVertex(1)->getNum()
+ // ,worst->tet()->getVertex(2)->getNum()
+ // ,worst->tet()->getVertex(3)->getNum());
+ /* printf("coucou 2 % cavity size %d\n",ITER,cavity.size());
+ for (std::list<MTet4*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc){
+ MTetrahedron *toto = (*itc)->tet();
+ // toto->xyz2uvw(center,uvw);
+ double mat[3][3], b[3], det;
+ toto->getMat(mat);
+ b[0] = center[0] - toto->getVertex(0)->x();
+ b[1] = center[1] - toto->getVertex(0)->y();
+ b[2] = center[2] - toto->getVertex(0)->z();
+ sys3x3(mat, b, uvw, &det);
+ printf("det = %g\n",det);
+ printf("%g %g %g -- \n",toto->getVertex(0)->x(),toto->getVertex(0)->y(),toto->getVertex(0)->z());
+ printf("%g %g %g -- \n",toto->getVertex(1)->x(),toto->getVertex(1)->y(),toto->getVertex(1)->z());
+ printf("%g %g %g -- \n",toto->getVertex(2)->x(),toto->getVertex(2)->y(),toto->getVertex(2)->z());
+ printf("%g %g %g -- \n",toto->getVertex(3)->x(),toto->getVertex(3)->y(),toto->getVertex(3)->z());
+ printf("tet quality %g\n",toto->gammaShapeMeasure());
+ printf("point %g %g %g outside %12.5E %12.5E %12.5E %12.5E\n",center[0],center[1],center[2],uvw[0], uvw[1], uvw[2],1-uvw[0]-uvw[1]-uvw[2]);
+ }
+ getchar();
+ */
myFactory.changeTetRadius(allTets.begin(), 0.0);
COUNT_MISS_2++;
for (std::list<MTet4*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc) (*itc)->setDeleted(false);
@@ -1708,233 +1593,25 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
///// do a 3D delaunay mesh assuming a set of vertices
-static void initialCube (std::vector<MVertex*> &v,
- MVertex *box[8],
- std::vector<MTet4*> &t){
- SBoundingBox3d bbox ;
- for (size_t i=0;i<v.size();i++){
- MVertex *pv = v[i];
- bbox += SPoint3(pv->x(),pv->y(),pv->z());
- }
- bbox *= 1.3;
- box[0] = new MVertex (bbox.min().x(),bbox.min().y(),bbox.min().z());
- box[1] = new MVertex (bbox.max().x(),bbox.min().y(),bbox.min().z());
- box[2] = new MVertex (bbox.max().x(),bbox.max().y(),bbox.min().z());
- box[3] = new MVertex (bbox.min().x(),bbox.max().y(),bbox.min().z());
- box[4] = new MVertex (bbox.min().x(),bbox.min().y(),bbox.max().z());
- box[5] = new MVertex (bbox.max().x(),bbox.min().y(),bbox.max().z());
- box[6] = new MVertex (bbox.max().x(),bbox.max().y(),bbox.max().z());
- box[7] = new MVertex (bbox.min().x(),bbox.max().y(),bbox.max().z());
- std::vector<MTetrahedron*> t_box;
- MTetrahedron *t0 = new MTetrahedron (box[2],box[7],box[3],box[1]);
- MTetrahedron *t1 = new MTetrahedron (box[0],box[7],box[1],box[3]);
- MTetrahedron *t2 = new MTetrahedron (box[6],box[1],box[7],box[2]);
- MTetrahedron *t3 = new MTetrahedron (box[0],box[1],box[7],box[4]);
- MTetrahedron *t4 = new MTetrahedron (box[1],box[4],box[5],box[7]);
- MTetrahedron *t5 = new MTetrahedron (box[1],box[7],box[5],box[6]);
- t.push_back(new MTet4(t0,0.0));
- t.push_back(new MTet4(t1,0.0));
- t.push_back(new MTet4(t2,0.0));
- t.push_back(new MTet4(t3,0.0));
- t.push_back(new MTet4(t4,0.0));
- t.push_back(new MTet4(t5,0.0));
- connectTets(t);
-}
-
-int straddling_segment_intersects_triangle(SPoint3 &p1,SPoint3 &p2,
- SPoint3 &q1,SPoint3 &q2,
- SPoint3 &q3)
-{
- double s1 = robustPredicates::orient3d(p1, p2, q2, q3);
- double s2 = robustPredicates::orient3d(p1, p2, q3, q1);
- double s3 = robustPredicates::orient3d(p1, p2, q1, q2);
-
- if (s1*s2 < 0.0 || s2 * s3 < 0.0) return false;
-
-
- double s4 = robustPredicates::orient3d(q1, q2, q3, p1);
- double s5 = robustPredicates::orient3d(q3, q2, q1, p2);
-
- return (s4*s5 >= 0) ;
-}
-
-static MTet4* search4Tet (MTet4 *t, MVertex *v, int _size,int & ITER) {
- if (t->inCircumSphere(v)) return t;
- SPoint3 p2 (v->x(),v->y(),v->z());
- std::set<MTet4*> path;
- while (1){
- path.insert(t);
- SPoint3 p1 = t->tet()->barycenter();
- int found = -1;
- MTet4 *neighOK = 0;
- for (int i = 0; i < 4; i++){
- MTet4 *neigh = t->getNeigh(i);
- if (neigh && path.find(neigh) == path.end()){
- neighOK = neigh;
- faceXtet fxt (t, i);
-
- SPoint3 q1(fxt.v[0]->x(),fxt.v[0]->y(),fxt.v[0]->z());
- SPoint3 q2(fxt.v[1]->x(),fxt.v[1]->y(),fxt.v[1]->z());
- SPoint3 q3(fxt.v[2]->x(),fxt.v[2]->y(),fxt.v[2]->z());
-
-
- if ( straddling_segment_intersects_triangle (p1,p2,q1,q2,q3)){
- found = i;
- break;
- }
- }
- }
- if (found < 0){
- if (neighOK)t = neighOK;
- else return 0;
- }
- else{
- t = t->getNeigh(found);
- }
- if (t->inCircumSphere(v)) {
- return t;
- }
- if (ITER++ > .5*_size) {
- break;
- }
- }
- return 0;
-}
-
-MTet4 * getTetToBreak (MVertex *v, std::vector<MTet4*> &t, int &NB_GLOBAL_SEARCH, int &ITER){
- // last inserted is used as starting point
- // we know it is not deleted
- unsigned int k = t.size() - 1;
- while(t[k]->isDeleted()){
- k--;
- }
- MTet4 *start = t[k];
- start = search4Tet (start,v,(int)t.size(),ITER);
- if (start)return start;
- // printf("Global Search has to be done\n");
- NB_GLOBAL_SEARCH++;
- for (size_t i = 0;i<t.size();i++){
- if (!t[i]->isDeleted() && t[i]->inCircumSphere(v))return t[i];
- }
- return 0;
-}
-
-bool tetOnBox (MTetrahedron *t, MVertex *box[8]){
- for (size_t i = 0;i<4;i++)
- for (size_t j = 0;j<8;j++)
- if (t->getVertex(i) == box[j])return true;
- return false;
-}
-
-
-void sanityCheck1(MTet4 *t)
-{
-}
-
-
+// void insertVerticesInRegion (GRegion *gr)
+// {
+// // compute edges that should not be
+// std::set<MEdge,Less_Edge> bnd;
+// std::list<GFace*> f_list = gr->faces();
+// for (std::list<GFace*>::iterator it = f_list.begin(); it != f_list.end(); ++it){
+// GFace *gf = *it;
+// for (i = 0;i< gf->triangles.size(); i++) {
+// for (j = 0; j < 3; j++) {
+// bnd.insert(gf->triangles[i]->getEdge(j));
+// }
+// }
+// }
+// }
- void delaunayMeshIn3D(std::vector<MVertex*> &v, std::vector<MTetrahedron*> &result, bool removeBox)
-{
- std::vector<MTet4*> t;
- t.reserve (v.size()*7);
- std::vector<faceXtet> conn;
- std::vector<faceXtet> shell;
- std::vector<MTet4*> cavity;
- MVertex *box[8];
- initialCube (v,box,t);
-
- int NB_GLOBAL_SEARCH = 0;
- int AVG_ITER = 0;
- SortHilbert(v);
+void delaunayMeshIn3D(std::vector<MVertex*> &v, std::vector<MTetrahedron*> &result, bool removeBox) {
double t1 = Cpu();
-
- /// double ta=0,tb=0,tc=0,td=0,T;
-
- for (size_t i=0;i<v.size();i++){
- MVertex *pv = v[i];
-
- int NITER = 0;
- // T = Cpu();
- MTet4 * found = getTetToBreak (pv,t,NB_GLOBAL_SEARCH,NITER);
- // ta += Cpu()-T;
- AVG_ITER += NITER;
- if(!found) {
- Msg::Error("cannot insert a point in 3D Delaunay");
- continue;
- }
- shell.clear();
- cavity.clear();
- // T = Cpu();
- recurFindCavity(shell, cavity, pv, found);
- // tb += Cpu()-T;
- double V = 0.0;
- for (unsigned int k=0;k<cavity.size();k++)V+=fabs(cavity[k]->tet()->getVolume());
-
- std::vector<MTet4*> extended_cavity;
- double Vb = 0.0;
-
- // T = Cpu();
- for (unsigned int count = 0; count < shell.size(); count++){
- const faceXtet &fxt = shell[count];
- MTetrahedron *tr;
- MTet4 *t4;
- MVertex *v0 = fxt.getVertex(0);
- MVertex *v1 = fxt.getVertex(1);
- MVertex *v2 = fxt.getVertex(2);
- MTet4 *otherSide = fxt.t1->getNeigh(fxt.i1);
- if (count < cavity.size()){
- t4 = cavity[count];
- tr = t4->tet() ;
- tr->setVertex(0,v0);
- tr->setVertex(1,v1);
- tr->setVertex(2,v2);
- tr->setVertex(3,pv);
- }
- else{
- tr = new MTetrahedron(v0,v1,v2,pv);
- t4 = new MTet4(tr, 0.0);
- t.push_back(t4);
- }
- Vb+= fabs(tr->getVolume());
- extended_cavity.push_back(t4);
- if (otherSide)
- extended_cavity.push_back(otherSide);
- }
- // tc += Cpu()-T;
-
- if (fabs(Vb-V) > 1.e-8 * (Vb+V))printf("%12.5E %12.5E\n",Vb,V);
-
- // reuse memory --> reinitialize MTet4s
- for (unsigned int k=0;k<std::min(cavity.size(),shell.size());k++){
- cavity[k]->setDeleted(false);
- for (unsigned int l=0;l<4;l++){
- cavity[k]->setNeigh(l,0);
- }
- }
- // T = Cpu();
- connectTets_vector2(extended_cavity,conn);
- // td += Cpu()-T;
- }
-
+ delaunayTriangulation (1, 1, v, result);
double t2 = Cpu();
- Msg::Info("Delaunay 3D done for %d points : CPU = %g, %d global searches, AVG walk size %g",v.size(), t2-t1,NB_GLOBAL_SEARCH,1.+(double)AVG_ITER/v.size());
- // printf("%d tets allocated (to compare with 7 #V = %d)\n",t.size(),7*v.size());
- // printf("%g %g %g %g --> %g(%g)\n",ta,tb,tc,td,t2-t1,ta+tb+tc+td);
-
- // FILE *f = fopen ("tet.pos","w");
- // fprintf(f,"View \"\"{\n");
- for (size_t i = 0;i<t.size();i++){
- if (t[i]->isDeleted() || (removeBox && tetOnBox (t[i]->tet(),box))) delete t[i]->tet();
- else {
- result.push_back(t[i]->tet());
- // t[i]->tet()->writePOS (f, false,false,true,false,false,false);
- }
- delete t[i];
- }
-
- if (removeBox)for (int i=0;i<8;i++)delete box[i];
- else for (int i=0;i<8;i++)v.push_back(box[i]);
-
- // fprintf(f,"};\n");
- // fclose(f);
+ Msg::Info("Tetrahedrization of %d points in %g seconds",v.size(),t2-t1);
}
+
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
index cec6d34b412eeb553ea397097042968678a49f18..79f601add24053d72321df9fb77a763c06685987 100644
--- a/Mesh/meshGRegionLocalMeshMod.cpp
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -8,6 +8,9 @@
#include "GRegion.h"
#include "GmshMessage.h"
#include "Numeric.h"
+#include "MHexahedron.h"
+#include "MPrism.h"
+#include "MPyramid.h"
static int edges[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
static int efaces[6][2] = {{0,2},{0,1},{1,2},{0,3},{2,3},{1,3}};
@@ -480,78 +483,6 @@ void buildVertexCavity_recur(MTet4 *t, MVertex *v, std::vector<MTet4*> &cavity)
// another one (of course, not the other one on the unswappable edge)
// after that crap, the sliver is trashed
-bool sliverRemoval(std::vector<MTet4*> &newTets, MTet4 *t,
- const qmTetrahedron::Measures &cr)
-{
- std::vector<MTet4*> cavity;
- std::vector<MTet4*> outside;
- std::vector<MVertex*> ring;
- MVertex *v1, *v2;
-
- bool isClosed[6];
- int nbSwappable = 0;
- int iSwappable = 0;
- for (int i = 0; i < 6; i++){
- isClosed[i] = buildEdgeCavity(t, i, &v1, &v2, cavity, outside, ring);
- if (isClosed[i]){
- nbSwappable++;
- iSwappable = i;
- }
- }
-
- if (nbSwappable == 0){
- // all edges are on model edges or model faces, which means that
- // nothing can be done
- return false;
- }
- else if (nbSwappable == 1){
- // classical case, the sliver has 5 edges on the boundary
- // try to swap first
- if (edgeSwap(newTets,t,iSwappable,qmTetrahedron::QMTET_ETA))return true;
- // 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());
- t->onWhat()->mesh_vertices.push_back(newv);
-
- if (!edgeSplit(newTets, t, newv, iSwappable, qmTetrahedron::QMTET_ONE)) return false;
- for (int i = 0; i < 4; i++){
- if (newTets[newTets.size() - 1]->tet()->getVertex(i) == newv){
- smoothVertex(newTets[newTets.size() - 1], i, cr);
- smoothVertexOptimize(newTets[newTets.size() - 1], i, cr);
- }
- }
-
- for (unsigned int i = 0; i < newTets.size(); i++){
- MTet4 *new_t = newTets[i];
- if (!(new_t->isDeleted())){
- for (int j = 0; j < 6; j++){
- MVertex *va = new_t->tet()->getVertex(edges[j][0]);
- MVertex *vb = new_t->tet()->getVertex(edges[j][1]);
- if (va == newv &&
- (va != v1 && vb != v1 && va != v2 && vb != v2)){
- collapseVertex(newTets,new_t,edges[j][0],edges[j][1],cr);
- }
- else if (vb == newv &&
- (va != v1 && vb != v1 && va != v2 && vb != v2)){
- collapseVertex(newTets,new_t,edges[j][1],edges[j][0],cr);
- }
- }
- }
- }
- return true;
- }
- else{
- for (int i = 0; i < 4; i++){
- smoothVertex(t, i, cr);
- smoothVertexOptimize(t, i, cr);
- }
- }
- return false;
-}
-
bool collapseVertex(std::vector<MTet4 *> &newTets,
MTet4 *t,
int iVertex,
@@ -836,303 +767,3 @@ bool smoothVertexOptimize(MTet4 *t, int iVertex, const qmTetrahedron::Measures &
}
}
-// Edge split sets ...
-// Here, we only build a list of tets that are a subdivision
-// of the given
-//static int edges[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
-
-// the resulting triangles are 012 and 230 in vector res
-// void splitPrism (std::vector<MTet4 *> &newTets,
-// std::vector<MVertex *> &steinerPoints,
-// MVertex* v1,
-// MVertex* v2,
-// MVertex* v3,
-// MVertex* v4,
-// MVertex* v5,
-// MVertex* v6,
-// const qualityMeasure4Tet &cr,
-// MVertex *res[4],
-// fs_cont &search,
-// GFace *fakeFace){
-// }
-
-// void splitQuadFace (MVertex* newv1,
-// MVertex* newv2,
-// MVertex* v21,
-// MVertex* v11,
-// const qualityMeasure4Tet &cr,
-// MVertex *res[4],
-// fs_cont &search,
-// GFace *fakeFace){
-// GFace* gfound = findInFaceSearchStructure (newv1,newv2,v11,search);
-// if (gfound){
-// res[0] = newv2;
-// res[1] = newv1;
-// res[2] = v11;
-// res[3] = v21;
-// }
-// else{
-// GFace* gfound = findInFaceSearchStructure (newv1,newv2,v21,search);
-// if (gfound){
-// res[0] = newv1;
-// res[1] = newv2;
-// res[2] = v21;
-// res[3] = v11;
-// }
-// // choose the best configuration
-// else{
-// double q1 = std::min(qmTri(newv1,newv2,v11,cr),qmTet(newv2,v11,v21,cr));
-// double q2 = std::min(qmTet(newv1,newv2,v21,cr),qmTet(newv1,v11,v21,cr));
-// if (q1 > q2){
-// res[0] = newv2;
-// res[1] = newv1;
-// res[2] = v11;
-// res[3] = v21;
-// MVertex *p1 = std::min(std::min(newv1,newv2),v11);
-// MVertex *p2 = std::min(std::min(newv2,v11),v21);
-// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p1, std::pair<MTriangle*,GFace*>(new MTriangle(newv1,newv2,v11),fakeFace)));
-// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p2, std::pair<MTriangle*,GFace*>(new MTriangle(newv2,v11,v21),fakeFace)));
-// }
-// else{
-// res[0] = newv1;
-// res[1] = newv2;
-// res[2] = v21;
-// res[3] = v11;
-// MVertex *p1 = std::min(std::min(newv1,newv2),v21);
-// MVertex *p2 = std::min(std::min(newv2,v11),v21);
-// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p1, std::pair<MTriangle*,GFace*>(new MTriangle(newv1,newv2,v21),fakeFace)));
-// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p2, std::pair<MTriangle*,GFace*>(new MTriangle(newv1,v11,v21),fakeFace)));
-// }
-// }
-// }
-// }
-// void splitQuadFace (std::vector<MTet4 *> &newTets,
-// std::vector<MVertex *> &steinerPoints,
-// MVertex* newv1,
-// MVertex* newv2,
-// MVertex* v21,
-// MVertex* v11,
-// MVertex* other,
-// const qualityMeasure4Tet &cr,
-// fs_cont &search,
-// GFace *fakeFace){
-// MTetrahedron *tr3,*tr4;
-// // now we look if there exists a face with the same vertices
-// GFace* gfound = findInFaceSearchStructure (newv1,newv2,v11,search);
-// if (gfound){
-// tr3 = new MTetrahedron ( newv1,newv2,v11,other);
-// tr4 = new MTetrahedron ( newv2,v11,v21,other);
-// }
-// else{
-// GFace* gfound = findInFaceSearchStructure (newv1,newv2,v21,search);
-// if (gfound){
-// tr3 = new MTetrahedron ( newv1,newv2,v21,other);
-// tr4 = new MTetrahedron ( newv1,v11,v21,other);
-// }
-// // choose the best configuration
-// else{
-// double vol;
-// double q1 = std::min(qmTet(newv1,newv2,v11,other,cr,vol),qmTet(newv2,v11,v21,other,cr,vol));
-// double q2 = std::min(qmTet(newv1,newv2,v21,other,cr,vol),qmTet(newv1,v11,v21,other,cr,vol));
-// if (q1 > q2){
-// tr3 = new MTetrahedron ( newv1,newv2,v11,other);
-// tr4 = new MTetrahedron ( newv2,v11,v21,other);
-// MVertex *p1 = std::min(std::min(newv1,newv2),v11);
-// MVertex *p2 = std::min(std::min(newv2,v11),v21);
-// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p1, std::pair<MTriangle*,GFace*>(new MTriangle(newv1,newv2,v11),fakeFace)));
-// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p2, std::pair<MTriangle*,GFace*>(new MTriangle(newv2,v11,v21),fakeFace)));
-// }
-// else{
-// tr3 = new MTetrahedron ( newv1,newv2,v21,other);
-// tr4 = new MTetrahedron ( newv1,v11,v21,other);
-// MVertex *p1 = std::min(std::min(newv1,newv2),v21);
-// MVertex *p2 = std::min(std::min(newv2,v11),v21);
-// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p1, std::pair<MTriangle*,GFace*>(new MTriangle(newv1,newv2,v21),fakeFace)));
-// search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p2, std::pair<MTriangle*,GFace*>(new MTriangle(newv1,v11,v21),fakeFace)));
-// }
-// }
-// }
-// }
-
-// bool splitEdgesOfTet (std::vector<MTet4 *> &newTets,
-// std::vector<MVertex *> &steinerPoints,
-// MTet4 *t1,
-// int nbEdges,
-// int e[6],
-// MVertex* pts[6],
-// const qualityMeasure4Tet &cr,
-// fs_cont &search,
-// GFace *fakeFace){
-// switch(nbEdges){
-// case 1 :
-// {
-// int iE = edges[0];
-// MVertex *newv = pts[0];
-// MVertex *v1 = t1->tet()->getVertex(e[iE][0]);
-// MVertex *v2 = t1->tet()->getVertex(e[iE][1]);
-// int oE = 5-iE;
-// MVertex *o1 = t1->tet()->getVertex(e[oE][0]);
-// MVertex *o2 = t1->tet()->getVertex(e[oE][1]);
-// MTetrahedron *tr1 = new MTetrahedron ( newv,v1,o1,o2);
-// MTetrahedron *tr2 = new MTetrahedron ( newv,v2,o2,o1);
-// MTet4 *t41 = new MTet4 (tr1,cr);
-// MTet4 *t42 = new MTet4 (tr2,cr);
-// newTets.push_back(t41);
-// newTets.push_back(t42);
-// }
-// break;
-// case 2 :
-// {
-// int iE1 = edges[0];
-// int iE2 = edges[1];
-
-// MVertex *newv1 = pts[0];
-// MVertex *newv2 = pts[1];
-
-// MVertex *v11 = t1->tet()->getVertex(e[iE1][0]);
-// MVertex *v12 = t1->tet()->getVertex(e[iE1][1]);
-// MVertex *v21 = t1->tet()->getVertex(e[iE2][0]);
-// MVertex *v22 = t1->tet()->getVertex(e[iE2][1]);
-
-// if (iE1 == 5-iE2){ // two opposite edges
-// MTetrahedron *tr1 = new MTetrahedron ( newv1,newv2,v11,v21);
-// MTetrahedron *tr2 = new MTetrahedron ( newv1,newv2,v21,v12);
-// MTetrahedron *tr3 = new MTetrahedron ( newv1,newv2,v12,v22);
-// MTetrahedron *tr4 = new MTetrahedron ( newv1,newv2,v22,v11);
-// MTet4 *t41 = new MTet4 (tr1,cr);
-// MTet4 *t42 = new MTet4 (tr2,cr);
-// MTet4 *t43 = new MTet4 (tr3,cr);
-// MTet4 *t44 = new MTet4 (tr4,cr);
-// newTets.push_back(t41);
-// newTets.push_back(t42);
-// newTets.push_back(t43);
-// newTets.push_back(t44);
-// }
-// else{ //two adjacend edges
-// MVertex *vsame,*other=0;
-// if (v11 == v21){vsame=v11; v11=v12; v21=v22;}
-// else if (v11 == v22){vsame=v11; v11=v12}
-// else if (v12 == v21){vsame=v12; v21=v22}
-// else if (v12 == v22){vsame=v12;}
-// else throw;
-// for (int i=0;i<4;i++){
-// if (vsame != t1->tet()->getVertex(i) &&
-// v11 != t1->tet()->getVertex(i) &&
-// v21 != t1->tet()->getVertex(i)){
-// other = t1->tet->getVertex(i);
-// break;
-// }
-// }
-// if (!other)throw;
-// MTetrahedron *tr1 = new MTetrahedron ( newv1,newv2,vsame,other);
-// splitQuadFace (newTets,steinerPoints,newv1,newv2,v21,v11,other,cr,search,fakeFace);
-// }
-// }
-// break;
-// case 3 :
-// {
-// int iE1 = edges[0];
-// int iE2 = edges[1];
-// int iE3 = edges[2];
-// MVertex *newv1;
-// MVertex *newv2;
-// MVertex *newv3;
-// std::sort(edges,edges+3);
-// if (iE1 == edges[0])newv1=pts[0];
-// else if (iE2 == edges[0])newv1=pts[1];
-// else if (iE3 == edges[0])newv1=pts[2];
-// if (iE1 == edges[1])newv2=pts[0];
-// else if (iE2 == edges[1])newv2=pts[1];
-// else if (iE3 == edges[1])newv2=pts[2];
-// if (iE1 == edges[2])newv3=pts[0];
-// else if (iE2 == edges[2])newv3=pts[1];
-// else if (iE3 == edges[2])newv3=pts[2];
-// iE1 = edges[0];
-// iE2 = edges[1];
-// iE3 = edges[2];
-// // edges are sorted and points correspond
-
-// mVertex *v[4];
-// // the 3 edges coincide at one vertex
-// int config;
-// if (iE1 == 0 && iE2 == 1 && iE3 == 2){
-// config = 1;
-// v[0] = t1->tet()->getVertex(0);
-// v[1] = t1->tet()->getVertex(1);
-// v[2] = t1->tet()->getVertex(2);
-// v[3] = t1->tet()->getVertex(3);
-// }
-// else if (iE1 == 0 && iE2 == 3 && iE3 == 4){
-// config = 1;
-// v[0] = t1->tet()->getVertex(1);
-// v[1] = t1->tet()->getVertex(0);
-// v[2] = t1->tet()->getVertex(2);
-// v[3] = t1->tet()->getVertex(3);
-// }
-// else if (iE1 == 1 && iE2 == 3 && iE3 == 5){
-// config = 1;
-// v[0] = t1->tet()->getVertex(2);
-// v[1] = t1->tet()->getVertex(0);
-// v[2] = t1->tet()->getVertex(1);
-// v[3] = t1->tet()->getVertex(3);
-// }
-// else if (iE1 == 2 && iE2 == 4 && iE3 == 5){
-// config = 1;
-// v[0] = t1->tet()->getVertex(3);
-// v[1] = t1->tet()->getVertex(0);
-// v[2] = t1->tet()->getVertex(1);
-// v[3] = t1->tet()->getVertex(2);
-// }
-// // three edges of the same face are cut
-// else if (iE1 == 0 && iE2 == 1 && iE3 == 3){
-// config = 2;
-// v[0] = t1->tet()->getVertex(3);
-// v[1] = t1->tet()->getVertex(1);
-// v[2] = t1->tet()->getVertex(0);
-// v[3] = t1->tet()->getVertex(2);
-// }
-// else if (iE1 == 0 && iE2 == 2 && iE3 == 4){
-// config = 2;
-// v[0] = t1->tet()->getVertex(2);
-// v[1] = t1->tet()->getVertex(1);
-// v[2] = t1->tet()->getVertex(0);
-// v[3] = t1->tet()->getVertex(3);
-// }
-// else if (iE1 == 1 && iE2 == 2 && iE3 == 5){
-// config = 2;
-// v[0] = t1->tet()->getVertex(1);
-// v[1] = t1->tet()->getVertex(2);
-// v[2] = t1->tet()->getVertex(0);
-// v[3] = t1->tet()->getVertex(3);
-// }
-// else if (iE1 == 3 && iE2 == 4 && iE3 == 5){
-// config = 2;
-// v[0] = t1->tet()->getVertex(0);
-// v[1] = t1->tet()->getVertex(2);
-// v[2] = t1->tet()->getVertex(1);
-// v[3] = t1->tet()->getVertex(3);
-// }
-// // the three edges for a kind of Z
-// else if (iE1 == 0 && iE2 == 3 && iE3 == 5){
-// config = 3;
-// }
-
-// if (config == 1){
-// }
-// else if (config == 2){
-// }
-// else{
-// throw; // for the moment.
-// }
-// }
-// default :
-// throw;
-// }
-// t1->setDeleted(true);
-
-
-// }
-
-
-
-// collapse
diff --git a/Mesh/meshGRegionLocalMeshMod.h b/Mesh/meshGRegionLocalMeshMod.h
index 9d2433c7085aec04f85512a880494c039dd92359..6504e27feff9a662b4e4392c43d5c34a214a5585 100644
--- a/Mesh/meshGRegionLocalMeshMod.h
+++ b/Mesh/meshGRegionLocalMeshMod.h
@@ -15,6 +15,8 @@
enum localMeshModAction {GMSH_DOIT, GMSH_EVALONLY};
+void LaplaceSmoothing (GRegion *gr);
+
bool edgeSwap(std::vector<MTet4*> &newTets, MTet4 *tet,
int iLocalEdge, const qmTetrahedron::Measures &cr);
@@ -37,7 +39,4 @@ bool egeSplit(std::vector<MTet4*> &newTets, MTet4 *tet,
MVertex *newVertex, int iLocalEdge,
const qmTetrahedron::Measures &cr);
-bool sliverRemoval(std::vector<MTet4*> &newTets, MTet4 *t,
- const qmTetrahedron::Measures &cr);
-
#endif