Skip to content
Snippets Groups Projects
Commit f6f7b3d7 authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

optimize hybrid meshes

parent b12589a1
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Mesh/meshGRegionDelaunayInsertion.cpp
+ 55
− 378
View file @ f6f7b3d7
...@@ -9,7 +9,6 @@ ...@@ -9,7 +9,6 @@
#include "GmshMessage.h" #include "GmshMessage.h"
#include "robustPredicates.h" #include "robustPredicates.h"
#include "OS.h" #include "OS.h"
//#include "BackgroundMesh.h"
#include "meshGRegion.h" #include "meshGRegion.h"
#include "meshGRegionLocalMeshMod.h" #include "meshGRegionLocalMeshMod.h"
#include "meshGRegionDelaunayInsertion.h" #include "meshGRegionDelaunayInsertion.h"
...@@ -18,7 +17,7 @@ ...@@ -18,7 +17,7 @@
#include "MTriangle.h" #include "MTriangle.h"
#include "Numeric.h" #include "Numeric.h"
#include "Context.h" #include "Context.h"
#include "HilbertCurve.h" #include "delaunay3d.h"
int MTet4::radiusNorm = 2; int MTet4::radiusNorm = 2;
static double LIMIT_ = 1; static double LIMIT_ = 1;
...@@ -45,6 +44,7 @@ static bool isActive(MTet4 *t, double limit_, int &active) ...@@ -45,6 +44,7 @@ static bool isActive(MTet4 *t, double limit_, int &active)
return false; return false;
} }
int MTet4::inCircumSphere(const double *p) const int MTet4::inCircumSphere(const double *p) const
{ {
double pa[3] = {base->getVertex(0)->x(), double pa[3] = {base->getVertex(0)->x(),
...@@ -64,6 +64,7 @@ int MTet4::inCircumSphere(const double *p) const ...@@ -64,6 +64,7 @@ int MTet4::inCircumSphere(const double *p) const
return (result > 0) ? 1 : 0; return (result > 0) ? 1 : 0;
} }
static int faces[4][3] = {{0,1,2}, {0,2,3}, {0,3,1}, {1,3,2}}; static int faces[4][3] = {{0,1,2}, {0,2,3}, {0,3,1}, {1,3,2}};
struct faceXtet{ struct faceXtet{
...@@ -143,30 +144,6 @@ void connectTets_vector2(std::vector<MTet4*> &t, std::vector<faceXtet> &conn) ...@@ -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> template <class ITER>
void connectTets(ITER beg, ITER end, std::set<MFace, Less_Face> *allEmbeddedFaces = 0) void connectTets(ITER beg, ITER end, std::set<MFace, Less_Face> *allEmbeddedFaces = 0)
{ {
...@@ -371,40 +348,6 @@ void recurFindCavity(std::vector<faceXtet> & shell, ...@@ -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 ) void printTets (const char *fn, std::list<MTet4*> &cavity, bool force = false )
{ {
FILE *f = Fopen (fn,"w"); FILE *f = Fopen (fn,"w");
...@@ -626,52 +569,12 @@ static void setLcs(MTetrahedron *t, std::map<MVertex*, double> &vSizes, ...@@ -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, GRegion *getRegionFromBoundingFaces(GModel *model,
std::set<GFace *> &faces_bound) std::set<GFace *> &faces_bound)
{ {
GModel::riter git = model->firstRegion(); 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()){ while (git != model->lastRegion()){
std::list <GFace *> _faces = (*git)->faces(); 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()){ if(_faces.size() == faces_bound.size()){
bool ok = true; bool ok = true;
for (std::list<GFace *>::iterator it = _faces.begin(); it != _faces.end(); ++it){ for (std::list<GFace *>::iterator it = _faces.begin(); it != _faces.end(); ++it){
...@@ -684,38 +587,6 @@ GRegion *getRegionFromBoundingFaces(GModel *model, ...@@ -684,38 +587,6 @@ GRegion *getRegionFromBoundingFaces(GModel *model,
return 0; 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, void non_recursive_classify(MTet4 *t, std::list<MTet4*> &theRegion,
std::set<GFace*> &faces_bound, GRegion *bidon, std::set<GFace*> &faces_bound, GRegion *bidon,
...@@ -947,6 +818,12 @@ void adaptMeshGRegion::operator () (GRegion *gr) ...@@ -947,6 +818,12 @@ void adaptMeshGRegion::operator () (GRegion *gr)
//template <class CONTAINER, class DATA> //template <class CONTAINER, class DATA>
void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm) 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 ; typedef std::list<MTet4 *> CONTAINER ;
CONTAINER allTets; CONTAINER allTets;
...@@ -1053,18 +930,6 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm) ...@@ -1053,18 +930,6 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
} }
// printf("coucou\n"); // 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()){ if (!newTets.size()){
break; break;
} }
...@@ -1235,6 +1100,7 @@ static void memoryCleanup(MTet4Factory &myFactory, std::set<MTet4*, compareTet4P ...@@ -1235,6 +1100,7 @@ static void memoryCleanup(MTet4Factory &myFactory, std::set<MTet4*, compareTet4P
void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify) void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
{ {
//printf("sizeof MTet4 = %d sizeof MTetrahedron %d sizeof(MVertex) %d\n", //printf("sizeof MTet4 = %d sizeof MTetrahedron %d sizeof(MVertex) %d\n",
// sizeof(MTet4), sizeof(MTetrahedron), sizeof(MVertex)); // sizeof(MTet4), sizeof(MTetrahedron), sizeof(MVertex));
...@@ -1245,6 +1111,7 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify) ...@@ -1245,6 +1111,7 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
std::set<MTet4*, compareTet4Ptr> &allTets = myFactory.getAllTets(); std::set<MTet4*, compareTet4Ptr> &allTets = myFactory.getAllTets();
int NUM = 0; int NUM = 0;
// printTets ("before.pos", cavity, true);
{ // leave this in a block so the map gets deallocated directly { // leave this in a block so the map gets deallocated directly
std::map<MVertex*, double> vSizesMap; std::map<MVertex*, double> vSizesMap;
...@@ -1332,8 +1199,8 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify) ...@@ -1332,8 +1199,8 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
double t1 = Cpu(); double t1 = Cpu();
while(1){ while(1){
// break; if (COUNT_MISS_2 > 100000)break;
if (ITER > maxVert)break; if (ITER >= maxVert)break;
if(allTets.empty()){ if(allTets.empty()){
Msg::Error("No tetrahedra in region %d %d", gr->tag(), allTets.size()); Msg::Error("No tetrahedra in region %d %d", gr->tag(), allTets.size());
break; break;
...@@ -1429,12 +1296,30 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify) ...@@ -1429,12 +1296,30 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
} }
} }
else{ else{
// printf("coucou 2 % cavity size %d\n",ITER,cavity.size()); // printf("%d %d %d %d\n",worst->tet()->getVertex(0)->getNum()
// for (std::list<MTet4*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc){ // ,worst->tet()->getVertex(1)->getNum()
// MTetrahedron *toto = (*itc)->tet(); // ,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); // 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]); 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); myFactory.changeTetRadius(allTets.begin(), 0.0);
COUNT_MISS_2++; COUNT_MISS_2++;
for (std::list<MTet4*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc) (*itc)->setDeleted(false); 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) ...@@ -1708,233 +1593,25 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
///// do a 3D delaunay mesh assuming a set of vertices ///// do a 3D delaunay mesh assuming a set of vertices
static void initialCube (std::vector<MVertex*> &v, // void insertVerticesInRegion (GRegion *gr)
MVertex *box[8], // {
std::vector<MTet4*> &t){ // // compute edges that should not be
SBoundingBox3d bbox ; // std::set<MEdge,Less_Edge> bnd;
for (size_t i=0;i<v.size();i++){ // std::list<GFace*> f_list = gr->faces();
MVertex *pv = v[i]; // for (std::list<GFace*>::iterator it = f_list.begin(); it != f_list.end(); ++it){
bbox += SPoint3(pv->x(),pv->y(),pv->z()); // GFace *gf = *it;
} // for (i = 0;i< gf->triangles.size(); i++) {
bbox *= 1.3; // for (j = 0; j < 3; j++) {
box[0] = new MVertex (bbox.min().x(),bbox.min().y(),bbox.min().z()); // bnd.insert(gf->triangles[i]->getEdge(j));
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 delaunayMeshIn3D(std::vector<MVertex*> &v, std::vector<MTetrahedron*> &result, bool removeBox) 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);
double t1 = Cpu(); double t1 = Cpu();
delaunayTriangulation (1, 1, v, result);
/// 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;
}
double t2 = Cpu(); 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()); Msg::Info("Tetrahedrization of %d points in %g seconds",v.size(),t2-t1);
// 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);
}
Mesh/meshGRegionLocalMeshMod.cpp
+ 3
− 372
View file @ f6f7b3d7
...@@ -8,6 +8,9 @@ ...@@ -8,6 +8,9 @@
#include "GRegion.h" #include "GRegion.h"
#include "GmshMessage.h" #include "GmshMessage.h"
#include "Numeric.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 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}}; 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) ...@@ -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) // another one (of course, not the other one on the unswappable edge)
// after that crap, the sliver is trashed // 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, bool collapseVertex(std::vector<MTet4 *> &newTets,
MTet4 *t, MTet4 *t,
int iVertex, int iVertex,
...@@ -836,303 +767,3 @@ bool smoothVertexOptimize(MTet4 *t, int iVertex, const qmTetrahedron::Measures & ...@@ -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
Mesh/meshGRegionLocalMeshMod.h
+ 2
− 3
View file @ f6f7b3d7
...@@ -15,6 +15,8 @@ ...@@ -15,6 +15,8 @@
enum localMeshModAction {GMSH_DOIT, GMSH_EVALONLY}; enum localMeshModAction {GMSH_DOIT, GMSH_EVALONLY};
void LaplaceSmoothing (GRegion *gr);
bool edgeSwap(std::vector<MTet4*> &newTets, MTet4 *tet, bool edgeSwap(std::vector<MTet4*> &newTets, MTet4 *tet,
int iLocalEdge, const qmTetrahedron::Measures &cr); int iLocalEdge, const qmTetrahedron::Measures &cr);
...@@ -37,7 +39,4 @@ bool egeSplit(std::vector<MTet4*> &newTets, MTet4 *tet, ...@@ -37,7 +39,4 @@ bool egeSplit(std::vector<MTet4*> &newTets, MTet4 *tet,
MVertex *newVertex, int iLocalEdge, MVertex *newVertex, int iLocalEdge,
const qmTetrahedron::Measures &cr); const qmTetrahedron::Measures &cr);
bool sliverRemoval(std::vector<MTet4*> &newTets, MTet4 *t,
const qmTetrahedron::Measures &cr);
#endif #endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment