diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index c8b28aff4290bcf31edb0f72c1184fd4c4bf6326..c90c4f2c58b46794b18847dfffc3baff83e2ac22 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -933,13 +933,12 @@ GRegion* OCC_Internals::addRegionToModel(GModel *model, TopoDS_Solid region)
GRegion *gr = getOCCRegionByNativePtr(model, region);
if(gr) return gr;
- // FIXME THE PREVIOUS IMPLEMENTATION WAS BETTER FOR SOME USERS :-)
+ // FIXME THE PREVIOUS IMPLEMENTATION WAS BETTER FOR SOME USERS :-)
buildShapeFromLists(region);
model->destroy();
buildLists();
buildGModel(model);
return getOCCRegionByNativePtr(model, region);
-
// _addShapeToMaps(region);
// buildShapeFromLists(region);
// buildGModel(model);
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 5958aef6b33c5b0abf12806011298e474ffd7930..92e5649694d1a2f6a9dca480bce3265717fc9014 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -179,7 +179,7 @@ void discreteFace::createGeometry()
int order = 1;
int nPart = 2;
- double eta = 5/(2.*3.14);
+ double eta = .5*5/(2.*3.14);
if (!_atlas.empty())return;
double dtSplit = 0.0;
@@ -192,8 +192,8 @@ void discreteFace::createGeometry()
triangulation* init = new triangulation(-1, tem,this);
init->iter = iter++;
allEdg2Tri = init->ed2tri;
-
toSplit.push(init);
+ // printf("%12.5E\n",(toSplit.top())->aspectRatio() );
if((toSplit.top())->genus()!=0 ||
(toSplit.top())->aspectRatio() > eta ||
(toSplit.top())->seamPoint){
diff --git a/Geo/gmshRegion.cpp b/Geo/gmshRegion.cpp
index 2c1de3b2c2a536b979151dcaf505ba633566b946..7583907d0ff5e5fba3dcf694ed0282008acdfa00 100644
--- a/Geo/gmshRegion.cpp
+++ b/Geo/gmshRegion.cpp
@@ -38,17 +38,17 @@ gmshRegion::gmshRegion(GModel *m, ::Volume *volume)
else
Msg::Error("Unknown surface %d", is);
}
- if(v->EmbeddedSurfaces){
- for(int i = 0; i < List_Nbr(v->EmbeddedSurfaces); i++){
- Surface *s;
- List_Read(v->EmbeddedSurfaces, i, &s);
- GFace *gf = m->getFaceByTag(abs(s->Num));
- if(gf)
- addEmbeddedFace(gf);
- else
- Msg::Error("Unknown surface %d", s->Num);
- }
- }
+ // if(v->EmbeddedSurfaces){
+ // for(int i = 0; i < List_Nbr(v->EmbeddedSurfaces); i++){
+ // Surface *s;
+ // List_Read(v->EmbeddedSurfaces, i, &s);
+ // GFace *gf = m->getFaceByTag(abs(s->Num));
+ // if(gf)
+ // addEmbeddedFace(gf);
+ // else
+ // Msg::Error("Unknown surface %d", s->Num);
+ // }
+ // }
resetMeshAttributes();
}
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 55ddb9adfb4321899ef3c7c8f37fbc8bc0abda54..8f4042452deb0364259f1ea37b567eb21fe5430e 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -49,6 +49,74 @@
#include "PViewData.h"
#endif
+class TEST_IF_MESH_IS_COMPATIBLE_WITH_EMBEDDED_ENTITIES {
+public:
+ void operator () (GRegion *gr) {
+ std::list<GEdge*> e = gr->embeddedEdges();
+ std::list<GFace*> f = gr->embeddedFaces();
+ if (e.empty() && f.empty())return;
+ std::map<MEdge,GEdge*,Less_Edge> edges;
+ std::map<MFace,GFace*,Less_Face> faces;
+ std::list<GEdge*>::iterator it = e.begin();
+ std::list<GFace*>::iterator itf = f.begin();
+ for ( ; it != e.end() ; ++it){
+ for (unsigned int i=0;i<(*it)->lines.size(); ++i){
+ if (distance ((*it)->lines[i]->getVertex(0),(*it)->lines[i]->getVertex(1)) > 1.e-12)
+ edges.insert(std::make_pair(MEdge((*it)->lines[i]->getVertex(0),(*it)->lines[i]->getVertex(1)),*it));
+ }
+ }
+ for ( ; itf != f.end() ; ++itf){
+ for (unsigned int i=0;i<(*itf)->triangles.size(); ++i){
+ faces.insert(std::make_pair(MFace((*itf)->triangles[i]->getVertex(0),(*itf)->triangles[i]->getVertex(1),(*itf)->triangles[i]->getVertex(2)),*itf));
+ }
+ }
+ Msg::Info ("Searching for %d embedded mesh edges and %d embedded mesh faces in region %d", edges.size(), faces.size(), gr->tag());
+ for (unsigned int k=0;k<gr->getNumMeshElements();k++){
+ for (int j=0;j<gr->getMeshElement(k)->getNumEdges();j++){
+ edges.erase (gr->getMeshElement(k)->getEdge(j));
+ }
+ for (int j=0;j<gr->getMeshElement(k)->getNumFaces();j++){
+ faces.erase (gr->getMeshElement(k)->getFace(j));
+ }
+ }
+ if (edges.empty() && faces.empty()) {
+ Msg::Info ("All embedded edges and faces are present in the final mesh");
+ }
+ if (edges.size()) {
+ char name[256];
+ sprintf(name,"missingEdgesOnRegion%d.pos",gr->tag());
+ Msg::Error("Region %d : %d mesh edges that should be embedded are missing in the final mesh",gr->tag(), (int)edges.size());
+ Msg::Error("Saving the missing edges in file %s",name);
+ FILE *f = fopen(name,"w");
+ fprintf(f,"View \" \" {\n");
+ for (std::map<MEdge,GEdge*,Less_Edge>::iterator it = edges.begin() ; it != edges.end(); ++it){
+ MVertex *v1 = it->first.getVertex(0);
+ MVertex *v2 = it->first.getVertex(1);
+ fprintf(f,"SL(%g,%g,%g,%g,%g,%g){%d,%d};\n",v1->x(),v1->y(),v1->z(),v2->x(),v2->y(),v2->z(), it->second->tag(),it->second->tag());
+ }
+ fprintf(f,"};\n");
+ fclose(f);
+ }
+ if (faces.size()) {
+ char name[256];
+ sprintf(name,"missingFacesOnRegion%d.pos",gr->tag());
+ Msg::Error("Region %d : %d mesh faces that should be embedded are missing in the final mesh",gr->tag(), (int)faces.size());
+ Msg::Error("Saving the missing faces in file %s",name);
+ FILE *f = fopen(name,"w");
+ fprintf(f,"View \" \" {\n");
+ for (std::map<MFace,GFace*,Less_Face>::iterator it = faces.begin() ; it != faces.end(); ++it){
+ MVertex *v1 = it->first.getVertex(0);
+ MVertex *v2 = it->first.getVertex(1);
+ MVertex *v3 = it->first.getVertex(2);
+ fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",v1->x(),v1->y(),v1->z(),v2->x(),v2->y(),v2->z(),
+ v3->x(),v3->y(),v3->z(),it->second->tag(),it->second->tag(),it->second->tag());
+ }
+ fprintf(f,"};\n");
+ fclose(f);
+ }
+ }
+};
+
template<class T>
static void GetQualityMeasure(std::vector<T*> &ele,
double &gamma, double &gammaMin, double &gammaMax,
@@ -835,6 +903,10 @@ static void Mesh3D(GModel *m)
// ensure that all volume Jacobians are positive
m->setAllVolumesPositive();
+ // std::for_each(m->firstRegion(), m->lastRegion(), optimizeMeshGRegionNetgen());
+ if (Msg::GetVerbosity() > 98)
+ std::for_each(m->firstRegion(), m->lastRegion(), TEST_IF_MESH_IS_COMPATIBLE_WITH_EMBEDDED_ENTITIES ());
+
CTX::instance()->mesh.changed = ENT_ALL;
double t2 = Cpu();
CTX::instance()->meshTimer[2] = t2 - t1;
@@ -850,6 +922,9 @@ void OptimizeMeshNetgen(GModel *m)
// Ensure that all volume Jacobians are positive
m->setAllVolumesPositive();
+ if (Msg::GetVerbosity() > 98)
+ std::for_each(m->firstRegion(), m->lastRegion(), TEST_IF_MESH_IS_COMPATIBLE_WITH_EMBEDDED_ENTITIES ());
+
double t2 = Cpu();
Msg::StatusBar(true, "Done optimizing 3D mesh with Netgen (%g s)", t2 - t1);
}
@@ -863,6 +938,9 @@ void OptimizeMesh(GModel *m)
// Ensure that all volume Jacobians are positive
m->setAllVolumesPositive();
+ if (Msg::GetVerbosity() > 98)
+ std::for_each(m->firstRegion(), m->lastRegion(), TEST_IF_MESH_IS_COMPATIBLE_WITH_EMBEDDED_ENTITIES ());
+
CTX::instance()->mesh.changed = ENT_ALL;
double t2 = Cpu();
Msg::StatusBar(true, "Done optimizing 3D mesh (%g s)", t2 - t1);
diff --git a/Mesh/delaunay3d.cpp b/Mesh/delaunay3d.cpp
index 97d459f7dd93be6c41d3f5f061c4bb02cefdf715..ee88ae6978fc34e4fe0ffd5781ba54556fe73d05 100644
--- a/Mesh/delaunay3d.cpp
+++ b/Mesh/delaunay3d.cpp
@@ -249,7 +249,7 @@ void HilbertSortB::Sort(Vert** vertices, int arraysize, int e, int d,
if (maxDepth > 0) {
if ((depth + 1) == maxDepth) {
- printf("ARGH\n");
+ // printf("ARGH max depth attained\n");
return;
}
}
@@ -544,8 +544,7 @@ static void edgeSwapPass (tetContainer &T)
*/
static void starShapeness (Vert *v, connContainer &bndK,
- std::vector<unsigned int> &_negatives,
- double threshold)
+ std::vector<unsigned int> &_negatives)
{
_negatives.clear();
for (unsigned int i=0; i< bndK.size(); i++) {
@@ -554,7 +553,7 @@ static void starShapeness (Vert *v, connContainer &bndK,
(double*)bndK[i].f.V[1],
(double*)bndK[i].f.V[2],
(double*)v);
- if (val <= threshold ) {
+ if (val <= 0.0 ) {
_negatives.push_back(i);
}
}
@@ -679,17 +678,16 @@ static Tet *tetInsideCavityWithFAce (Face &f, cavityContainer &cavity)
return NULL;
}
-static bool fixDelaunayCavity (double threshold,
- Vert *v,
+static bool fixDelaunayCavity (Vert *v,
cavityContainer &cavity,
connContainer &bndK,
int myThread, int K,
std::vector<unsigned int> & _negatives)
{
- starShapeness (v, bndK, _negatives, threshold);
+ starShapeness (v, bndK, _negatives);
if (_negatives.empty())return false;
-
+ // return true;
// unset all tets of the cavity
for (unsigned int i=0; i< cavity.size(); i++)cavity[i]->unset(myThread,K);
for (unsigned int i=0; i<bndK.size(); i++)if(bndK[i].t)bndK[i].t->unset(myThread,K);
@@ -715,7 +713,7 @@ static bool fixDelaunayCavity (double threshold,
}
removeIsolatedTets(containsV, cavity,bndK,myThread,K);
buildDelaunayBall (cavity,bndK);
- starShapeness (v, bndK, _negatives, threshold);
+ starShapeness (v, bndK, _negatives);
}
return false;
}
@@ -751,6 +749,7 @@ Tet* walk (Tet *t, Vert *v, int maxx, double &totSearch, int thread)
while (1){
totSearch++;
if (t == NULL) {
+ // printf("in an embedded edge\n");
return NULL; // we should NEVER return here
}
// if (t->inSphere(v,thread)) {return t;}
@@ -905,13 +904,39 @@ static Tet* randomTet (int thread, tetContainer &allocator)
}
+int isCavityCompatibleWithEmbeddedEdges(cavityContainer &cavity,
+ connContainer &bndK,
+ edgeContainer &allEmbeddedEdges){
+
+ const unsigned int bSize = bndK.size();
+ std::vector<Edge> ed;
+ for (unsigned int i=0; i<bSize; i++) {
+ for (unsigned int j=0; j<3; j++) {
+ if (bndK[i].f.V[j] > bndK[i].f.V[(j+1)%3]){
+ ed.push_back(Edge(bndK[i].f.V[j], bndK[i].f.V[(j+1)%3]));
+ }
+ }
+ }
+
+ for (unsigned int i=0; i<cavity.size(); i++){
+ for (int j=0;j<6;j++){
+ Edge e = cavity[i]->getEdge(j);
+ if (std::find(ed.begin(), ed.end(), e) == ed.end() && allEmbeddedEdges.find(e)){
+ return 0;
+ }
+ }
+ }
+ return 1;
+}
+
+
//#define _VERBOSE 1
void delaunayTrgl (const unsigned int numThreads,
const unsigned int NPTS_AT_ONCE,
unsigned int Npts,
std::vector<Vert*> assignTo[],
tetContainer &allocator,
- double threshold)
+ edgeContainer *embeddedEdges)
{
#if defined(_VERBOSE)
double totSearchGlob=0;
@@ -983,6 +1008,7 @@ void delaunayTrgl (const unsigned int numThreads,
////////////////////////////////////////////////////////////////////////////////////
for (unsigned int iPGlob=0 ; iPGlob < maxLocSizeK; iPGlob++){
+ // printf("%d vs %d\n",iPGlob,maxLocSizeK);
#if defined(_OPENMP)
#pragma omp barrier
#endif
@@ -995,9 +1021,11 @@ void delaunayTrgl (const unsigned int numThreads,
if(vToAdd[K]){
// In 3D, insertion of a point may lead to deletion of tets !!
if (!Choice[K]->V[0])Choice[K] = randomTet (myThread, allocator);
+ // int nbCoucou=0;
while(1){
t[K] = walk ( Choice[K] , vToAdd[K], Npts, totSearch, myThread);
if (t[K])break;
+ // printf("coucou %d\n",nbCoucou++);
// the domain may not be convex. we then start from a random tet and
// walk from there
Choice[K] = randomTet(rand() % numThreads, allocator);
@@ -1017,7 +1045,11 @@ void delaunayTrgl (const unsigned int numThreads,
cavityK.clear(); bndK.clear();
delaunayCavity2(t[K], NULL, vToAdd[K], cavityK, bndK, myThread, K);
// verify the cavity
- if (fixDelaunayCavity (threshold, vToAdd[K], cavityK, bndK, myThread, K, _negatives)){
+ if (fixDelaunayCavity (vToAdd[K], cavityK, bndK, myThread, K, _negatives)){
+ vToAdd[K] = NULL;
+ invalidCavities [myThread]++;
+ }
+ if (embeddedEdges && !isCavityCompatibleWithEmbeddedEdges(cavityK, bndK, *embeddedEdges)){
vToAdd[K] = NULL;
invalidCavities [myThread]++;
}
@@ -1093,7 +1125,8 @@ void delaunayTrgl (const unsigned int numThreads,
//checkLocalDelaunayness(T,"final");
-
+
+
// printf(" %12.5E %12.5E %12.5E tot %12.5E \n",t2,t3,t4,t2+t3+t4);
#if defined(_VERBOSE)
diff --git a/Mesh/delaunay3d_private.h b/Mesh/delaunay3d_private.h
index fb330ccf28fa6bd59d9c8c115744f06d6cc05ed5..a0360dc244dd27709ab721e111d8d6c3f437cb8a 100644
--- a/Mesh/delaunay3d_private.h
+++ b/Mesh/delaunay3d_private.h
@@ -126,9 +126,10 @@ inline double orientationTestFast (Vert *va, Vert *vb, Vert *vc, Vert *vd)
return orientationTestFast ((double*)va,(double*)vb,(double*)vc,(double*)vd);
}
-struct Edge {
+class Edge {
+public :
Vert *first, *second;
- Edge (Vert *v1, Vert *v2) : first(v1), second(v2)
+ Edge (Vert *v1, Vert *v2) : first(std::min(v1,v2)), second(std::max(v1,v2))
{
}
bool operator == (const Edge &e) const{
@@ -142,6 +143,37 @@ struct Edge {
}
};
+struct edgeContainer
+{
+ std::vector<std::vector<Edge> > _hash;
+ size_t _size;
+ edgeContainer(unsigned int N = 1000000)
+ {
+ _size = 0;
+ _hash.resize(N);
+ }
+
+ inline bool find (const Edge &e) const {
+ size_t h = ((size_t) e.first >> 8) ;
+ const std::vector<Edge> &v = _hash[h %_hash.size()];
+ for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return true;}
+ return false;
+ }
+
+ bool empty () const {return _size == 0;}
+
+ bool addNewEdge (const Edge &e)
+ {
+ size_t h = ((size_t) e.first >> 8) ;
+ std::vector<Edge> &v = _hash[h %_hash.size()];
+ for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return false;}
+ v.push_back(e);
+ _size++;
+ return true;
+ }
+};
+
+
struct Face {
Vert *v[3];
Vert *V[3];
@@ -181,6 +213,7 @@ struct Tet {
T[0] = T[1] = T[2] = T[3] = NULL;
setAllDeleted();
}
+ // inline bool isFace () const {return V[3]==NULL;}
int setVerticesNoTest (Vert *v0, Vert *v1, Vert *v2, Vert *v3)
{
_modified=true;
@@ -358,7 +391,7 @@ void delaunayTrgl(const unsigned int numThreads,
const unsigned int NPTS_AT_ONCE,
unsigned int Npts,
std::vector<Vert*> assignTo[],
- tetContainer &allocator, double threshold = 0.0);
+ tetContainer &allocator, edgeContainer *embedded = 0);
bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread);
#endif
diff --git a/Mesh/delaunay_refinement.cpp b/Mesh/delaunay_refinement.cpp
index a2a1b18123ab1ce54394c6343280b078fa9957e7..4bc2ebf72eea873a16d15d3283d6981bd59f2047 100644
--- a/Mesh/delaunay_refinement.cpp
+++ b/Mesh/delaunay_refinement.cpp
@@ -22,52 +22,47 @@
#include "MTriangle.h"
#include "GRegion.h"
#include "GFace.h"
+#include "MLine.h"
+#include "Context.h"
#include "SBoundingBox3d.h"
+#include "GModel.h"
+#include "Field.h"
-typedef std::set< Edge > edgeContainer2;
-
-long int AVGSEARCH;
-
-struct edgeContainer
-{
- std::set< Edge > _hash2;
- std::vector<std::vector<Edge> > _hash;
- edgeContainer(unsigned int N = 1000000)
- {
- _hash.resize(N);
- }
- bool addNewEdge2 (const Edge &e)
- {
- std::set< Edge >::iterator it = _hash2.find(e);
- if (it != _hash2.end())return false;
- _hash2.insert(e);
- return true;
- }
- bool addNewEdge (const Edge &e)
- {
- size_t h = ((size_t) e.first >> 3) ;
- std::vector<Edge> &v = _hash[h %_hash.size()];
- AVGSEARCH+=v.size();
- for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return false;}
- v.push_back(e);
- return true;
- }
-};
+#define SQR(X) (X)*(X)
+const double EXCLUSION_FACTOR = 0.8;
struct IPT {
- double _x1,_x2,_x3,_x4;
+ double _x1,_x2,_x3,_x4, _x5;
IPT(double x1, double x2, double x3, double x4)
- : _x1(x1), _x2(x2), _x3(x3), _x4(x4)
+ : _x1(x1), _x2(x2), _x3(x3), _x4(x4),_x5(0.0)
{
};
};
+////int C_COUNT = 0;
+
+//double _C1, _C2, _C3, _C4;
+static double GMSHSIZE (const SPoint3 &p, Field *f, double lc_nodal) {
+ double f_field = 1.e22;
+ if (!CTX::instance()->mesh.lcExtendFromBoundary) lc_nodal = 1.e22;
+ // double _a = Cpu();
+ if (f)f_field = (*f)(p.x(), p.y(), p.z());
+ // _C3 += Cpu() - _a;
+
+ double lc = std::min (lc_nodal,f_field);
+ lc = std::max(lc, CTX::instance()->mesh.lcMin);
+ lc = std::min(lc, CTX::instance()->mesh.lcMax);
+ // C_COUNT++;
+ return lc;
+}
+
+
double adaptiveTrapezoidalRule(SPoint3 p1 , SPoint3 p2 ,
double lc1 , double lc2 ,
- double (*f)(const SPoint3 &p, void *),
- void *data, std::vector< IPT > & _result,
+ std::vector< IPT > & _result,
double &dl, std::stack<IPT> &_stack,
- double epsilon = 1.e-5)
+ Field *bgm,
+ double epsilon = 5.e-2)
{
// _stack.clear();
_result.clear();
@@ -78,12 +73,18 @@ double adaptiveTrapezoidalRule(SPoint3 p1 , SPoint3 p2 ,
SPoint3 dp = p2-p1;
// value of f on both sides
- double f1 = lc1; //f(p1 + dp*t1,data);
- double f2 = lc2; //f(p1 + dp*t2,data);
-
+ double f1 = lc1;//GMSHSIZE(p1, bgm, lc1); //f(p1 + dp*t1,data);
+ double f2 = lc2;//GMSHSIZE(p2, bgm, lc2); //f(p1 + dp*t2,data);
+
dl = p1.distance(p2);
- // printf ("edge length %g lc %g %g\n",dl,f1,f2);
+ // adim_lenght of half the edge should be bigger than EXCLUSION_FACTOR
+ // +------x--------+
+
+ if (dl / (2*std::min(lc1,lc2)) <= EXCLUSION_FACTOR){
+ // printf ("edge length %g lc %g %g\n",dl,f1,f2);
+ return EXCLUSION_FACTOR;
+ }
// add one subsegment on the stack
IPT top (t1,t2,f1,f2);
@@ -98,40 +99,68 @@ double adaptiveTrapezoidalRule(SPoint3 p1 , SPoint3 p2 ,
f1 = pp._x3;
f2 = pp._x4;
// mid point
- double t12 = 0.5* (t1+t2);
- SPoint3 pmid = p1 + dp*t12;
- double dt = t2-t1;
+ const double t12 = 0.5* (t1+t2);
+ const SPoint3 pmid = p1 + dp*t12;
+ const double dt = t2-t1;
// average should be compared to mid value
- double f12 = 0.5* (f1+f2);
- if (fabs (f12 - 0.5*(f1+f2)) > epsilon*dt ) {
+ // double f12 = 0.5* (f1+f2);
+ double f12 = GMSHSIZE(pmid, bgm, lc1 + t12 *(lc2-lc1) ); //f(p1 + dp*t2,data);
+ const double AA = dt * dl;
+ const double I1 = AA / f12;
+ const double I2 = 2 * AA / (f1+f2);
+
+ if (fabs (I1 - I2) > epsilon) {
IPT left (t1,t12,f1,f12);
IPT right (t12,t2,f12,f2);
_stack.push(left);
_stack.push(right);
}
else {
- _result.push_back (pp);
// compute the integral using trapezoidal rule on both sides
totalValue += 1./((0.5*f12+0.25*(f1+f2)))*dt;
+ pp._x5 = totalValue*dl;
+ _result.push_back (pp);
}
}
// take into account the real length of the edge
totalValue *= dl;
- // printf("adimensional length %g\n",totalValue);
+ // printf("adimensional length %g %d evals\n",totalValue,_result.size() );
return totalValue;
}
-void saturateEdge(Edge &e, std::vector<Vert*> &S,
- double (*f)(const SPoint3 &p, void *),
- void *data, std::stack<IPT> &temp)
+void saturateEdge(Edge &e, std::vector<Vert*> &S, std::stack<IPT> &temp)
{
std::vector< IPT > _result;
double dl;
SPoint3 p1 = e.first->point();
SPoint3 p2 = e.second->point();
+
+ FieldManager *fields = GModel::current()->getFields();
+ Field *bgm = NULL;
+ if(fields->getBackgroundField() > 0){
+ bgm = fields->get(fields->getBackgroundField());
+ }
+
+ if (1){
+ const double length = p1.distance(p2);
+ const double lc_avg = (e.first->lc() + e.second->lc())*.5;
+ const double l_adim_approx = length / lc_avg;
+ if (l_adim_approx > 10.){
+ SPoint3 pmid = (p1+p2)*0.5;
+ const double lc = GMSHSIZE(pmid, bgm, lc_avg);
+ S.push_back(new Vert(pmid.x(),pmid.y(),pmid.z(),lc));
+ return;
+ }
+ }
+
+ // double _a = Cpu();
+
const double dN = adaptiveTrapezoidalRule(p1,p2,e.first->lc(), e.second->lc(),
- f, data, _result, dl, temp);
+ _result, dl, temp, bgm);
+ // _C1 += Cpu() - _a;
+ // _a = Cpu();
+
const int N = (int) (dN+0.1);
const double interval = dN/N;
double L = 0.0;
@@ -139,33 +168,40 @@ void saturateEdge(Edge &e, std::vector<Vert*> &S,
// printf("edge length %g %d intervals of size %g (%d results)\n",
// dl,N,interval,_result.size());
const unsigned int Nr = _result.size();
+ double _f = interval;
for (unsigned int i=0; i< Nr ; i++) {
const IPT & rr = _result[i];
const double t1 = rr._x1;
const double t2 = rr._x2;
const double f1 = rr._x3;
const double f2 = rr._x4;
+ const double l1 = rr._x3;
+ const double l2 = rr._x4;
const double dL = 2.*(t2-t1) * dl / (f1+f2);
+ const double _f1 = i ? _result[i-1]._x5 : 0.0;
+ const double _f2 = rr._x5;
+
// printf("%g --> %g for %g --> %g\n",L,dL,t1,t2);
- double L0 = L;
while (1) {
- const double t = t1 + (L+interval-L0)*(t2-t1) / dL;
- if (t >= t2*.999) {
+ const double t = t1 + (L+interval-_f1)*(t2-t1) / dL;
+ if (t >= t2*.999) {
break;
}
else {
// printf("%g ",t);
SPoint3 p = p1 * (1.-t) + p2*t;
- double lc = e.first->lc() * (1.-t) + e.second->lc()*t;
- const double dx = 0;//1.e-12 * (double) rand() / RAND_MAX;
- const double dy = 0;//1.e-12 * (double) rand() / RAND_MAX;
- const double dz = 0;//1.e-12 * (double) rand() / RAND_MAX;
+ const double lc = l1 + (_f-_f1)*(l2-l1) / (_f2-_f1);
+ const double dx = 0;//lc * 1.e-12 * (double) rand() / RAND_MAX;
+ const double dy = 0;//lc * 1.e-12 * (double) rand() / RAND_MAX;
+ const double dz = 0;//lc * 1.e-12 * (double) rand() / RAND_MAX;
S.push_back(new Vert(p.x()+dx,p.y()+dy,p.z()+dz,lc));
L += interval;
+ _f += interval;
}
}
}
+ // _C2 += Cpu() - _a;
// printf(" press enter\n");
// getchar();
// printf("%d points added\n",S.size());
@@ -176,11 +212,11 @@ void saturateEdge(Edge &e, std::vector<Vert*> &S,
void saturateEdges(edgeContainer &ec,
tetContainer &T,
int nbThreads,
- std::vector<Vert*> &S,
- double (*f)(const SPoint3 &p, void *), void *data)
+ std::vector<Vert*> &S)
{
+ // _C1 = 0;
+ // _C2 = 0;
std::stack<IPT> temp;
- AVGSEARCH= 0;
// FIXME
const int N = T.size(0);
for (int i=0;i<N;i++){
@@ -191,16 +227,17 @@ void saturateEdges(edgeContainer &ec,
Edge ed = t->getEdge(iEdge);
bool isNew = ec.addNewEdge(ed);
if (isNew){
- saturateEdge (ed, S, f, data, temp);
+ // saturateEdge_A_LA_FACON_PaulLouis(ed, S);
+ saturateEdge (ed, S, temp);
}
}
}
}
+ // printf("timings %12.5E %12.5E\n",_C1,_C2);
}
// F I L T E R I N G
-#define SQR(X) (X)*(X)
class volumePointWithExclusionRegion {
public :
@@ -209,8 +246,7 @@ public :
inline bool inExclusionZone (volumePointWithExclusionRegion *p) const
{
- const double FACTOR = 0.8;
- const double K = FACTOR*p->_v->lc();
+ const double K = EXCLUSION_FACTOR*p->_v->lc();
const double d =
SQR(p->_v->x() - _v->x())+
SQR(p->_v->y() - _v->y())+
@@ -269,9 +305,7 @@ public:
void filterVertices(const int numThreads,
vertexFilter &_filter,
- std::vector<Vert*> &add,
- double (*f)(const SPoint3 &p, void *),
- void *data)
+ std::vector<Vert*> &add)
{
std::vector<int> indices;
SortHilbert(add, indices);
@@ -307,11 +341,6 @@ void filterVertices(const int numThreads,
}
}
-double _fx (const SPoint3 &p, void *)
-{
- return fabs(0.0125 + .02*p.x());
-}
-
/*
static void _print (const char *name, std::vector<Vert*> &T)
{
@@ -338,7 +367,7 @@ void computeAdjacencies (Tet *t, int iFace, connSet &faceToTet)
else{
t->T[iFace] = it->t;
it->t->T[it->i] =t;
- faceToTet.erase(it);
+ // faceToTet.erase(it);
}
}
@@ -348,12 +377,92 @@ bool edgeSwaps(tetContainer &T, int myThread)
return false;
}
+void updateSize (MVertex *v, std::map<MVertex*, double> &s, double d){
+ d = std::max(d, CTX::instance()->mesh.lcMin);
+ d = std::min(d, CTX::instance()->mesh.lcMax);
+ std::map<MVertex*, double>::iterator it = s.find(v);
+ if (it == s.end())s[v]=d;
+ else it->second = std::min(it->second, d);
+}
+
+void createEmbeddedEdges (GRegion *gr,
+ edgeContainer &embedded,
+ std::vector<Vert *> &_vertices){
+ std::list<GEdge*> e = gr->embeddedEdges();
+ for (std::list<GEdge*>::iterator it = e.begin() ; it != e.end(); ++it){
+ for (unsigned int i = 0; i < (*it)->lines.size(); i++){
+ Vert *v0 = _vertices[(*it)->lines[i]->getVertex(0)->getIndex()];
+ Vert *v1 = _vertices[(*it)->lines[i]->getVertex(1)->getIndex()];
+ Edge e (v0,v1);
+ embedded.addNewEdge(e);
+ }
+ }
+}
+
+
+void disconnectEmbeddedFaces (GRegion *gr, connSet &faceToTet, std::vector<Vert *> &_vertices){
+ std::list<GFace*> f = gr->embeddedFaces();
+ for (std::list<GFace*>::iterator it = f.begin() ; it != f.end(); ++it){
+ for (unsigned int i = 0; i < (*it)->triangles.size(); i++){
+ bool ok1 = false;
+ bool ok2 = false;
+ Vert *v0 = _vertices[(*it)->triangles[i]->getVertex(0)->getIndex()];
+ Vert *v1 = _vertices[(*it)->triangles[i]->getVertex(1)->getIndex()];
+ Vert *v2 = _vertices[(*it)->triangles[i]->getVertex(2)->getIndex()];
+ Face f (v0,v1,v2);
+ conn c(f,0,NULL);
+ connSet::iterator itf = faceToTet.find(c);
+ if (itf != faceToTet.end()){
+ Tet *t1 = itf->t;
+ Tet *t2 = itf->t->T[itf->i];
+ // if (!(f == t1->getFace(itf->i)))printf("aie\n");
+ for (int k=0;k<4;k++){
+ if (t1 && t1->T[k] == t2){t1->T[k]=NULL;ok1=true;}
+ if (t2 && t2->T[k] == t1){t2->T[k]=NULL;ok2=true;}
+ }
+ if (!ok1 || !ok2)printf("oops %d %d %p %p\n",ok1,ok2,t1,t2);
+ }
+ else printf("oops2\n");
+ }
+ }
+}
+
+void computeMeshSizes (GRegion *gr, std::map<MVertex*, double> &s){
+ std::list<GEdge*> e = gr->embeddedEdges();
+ for (std::list<GEdge*>::iterator it = e.begin() ; it != e.end(); ++it){
+ for (unsigned int i = 0; i < (*it)->lines.size(); i++){
+ double d = distance((*it)->lines[i]->getVertex(0), (*it)->lines[i]->getVertex(1));
+ updateSize ((*it)->lines[i]->getVertex(0), s, d);
+ updateSize ((*it)->lines[i]->getVertex(1), s, d);
+ }
+ }
+ std::list<GFace*> f ;
+ std::list<GFace*> f1 = gr->faces();
+ std::list<GFace*> f2 = gr->embeddedFaces();
+ f.insert(f.end(),f1.begin(), f1.end());
+ f.insert(f.end(),f2.begin(), f2.end());
+ for (std::list<GFace*>::iterator it = f.begin() ; it != f.end(); ++it){
+ for (unsigned int i = 0; i < (*it)->triangles.size(); i++){
+ for (unsigned int j = 0; j < 3; j++){
+ double d = distance((*it)->triangles[i]->getVertex(j), (*it)->triangles[i]->getVertex((j+1)%3));
+ updateSize ((*it)->triangles[i]->getVertex(j), s, d);
+ updateSize ((*it)->triangles[i]->getVertex((j+1)%3), s, d);
+ }
+ }
+ }
+}
+
void edgeBasedRefinement(const int numThreads,
const int nptsatonce,
GRegion *gr)
{
// fill up old Datastructures
+ edgeContainer embeddedEdges (10000);
+
+ std::map<MVertex*, double> sizes;
+ computeMeshSizes (gr, sizes);
+
tetContainer allocator (numThreads,1000000);
SBoundingBox3d bb;
@@ -383,7 +492,10 @@ void edgeBasedRefinement(const int numThreads,
for (std::set<MVertex*>::iterator it = all.begin();it !=all.end(); ++it){
MVertex *mv = *it;
mv->setIndex(counter);
- Vert *v = new Vert (mv->x(),mv->y(),mv->z(),1.e22, counter);
+ double size = CTX::instance()->mesh.lcMax;
+ std::map<MVertex*,double>::iterator itv = sizes.find(mv);
+ if (itv != sizes.end())size = itv->second;
+ Vert *v = new Vert (mv->x(),mv->y(),mv->z(),size, counter);
_vertices[counter] = v;
bb += SPoint3(v->x(),v->y(),v->z());
_ma[v] = mv;
@@ -408,9 +520,12 @@ void edgeBasedRefinement(const int numThreads,
delete tt;
}
T.clear();
+ disconnectEmbeddedFaces (gr, faceToTet, _vertices);
+ createEmbeddedEdges (gr, embeddedEdges, _vertices);
}
}
+
// do not allow to saturate boundary edges
{
for (unsigned int i=0;i< allocator.size(0);i++) {
@@ -421,36 +536,22 @@ void edgeBasedRefinement(const int numThreads,
for (int k=0;k<3;k++){
Vert *vi = f.V[k];
Vert *vj = f.V[(k+1)%3];
- double l = sqrt ((vi->x()-vj->x())*(vi->x()-vj->x())+
- (vi->y()-vj->y())*(vi->y()-vj->y())+
- (vi->z()-vj->z())*(vi->z()-vj->z()));
ec.addNewEdge(Edge(vi,vj));
- vi->lc() = std::min (l,vi->lc() );
- vj->lc() = std::min (l,vj->lc() );
}
}
}
}
- for (unsigned int i=0;i< allocator.size(0);i++) {
- Tet *tt = allocator (0,i);
- for (int j=0;j<6;j++){
- Edge e = tt->getEdge(j);
- if(e.first->lc() == 1.e22){/*printf("coucou\n");*/e.first->lc() = e.second->lc();}
- else if(e.second->lc() == 1.e22){/*printf("coucou\n");*/e.second->lc() = e.first->lc();}
- }
- }
}
std::vector<Vert*> add_all;
{
- // vertexFilter _filter (bb, 20);
+ // vertexFilter _filter (bb, 20);
vertexFilter _filter;
for (unsigned int i=0;i<_vertices.size();i++){
_filter.insert( _vertices[i] );
}
int iter = 1;
-
Msg::Info("------------------------------------- SATUR FILTR SORTH DELNY TIME TETS");
double __t__ = Cpu();
@@ -458,9 +559,12 @@ void edgeBasedRefinement(const int numThreads,
while(1){
std::vector<Vert*> add;
double t1 = Cpu();
- saturateEdges (ec, allocator, numThreads, add, _fx, NULL);
+ // C_COUNT = 0;
+ saturateEdges (ec, allocator, numThreads, add);
+ // printf("%d calls %d points added",C_COUNT,add.size());
double t2 = Cpu();
- filterVertices (numThreads, _filter, add, _fx, NULL);
+ filterVertices (numThreads, _filter, add);
+ // printf(" %d points remain\n",add.size());
double t3 = Cpu();
if (add.empty())break;
// randomize vertices (EXTREMELY IMPORTANT FOR NOT DETERIORATING PERFORMANCE)
@@ -469,7 +573,7 @@ void edgeBasedRefinement(const int numThreads,
std::vector<int> indices;
SortHilbert(add, indices);
double t4 = Cpu();
- delaunayTrgl (1,1,add.size(), &add,allocator,1.e-28);
+ delaunayTrgl (1,1,add.size(), &add,allocator,embeddedEdges.empty() ? NULL : &embeddedEdges);
double t5 = Cpu();
add_all.insert (add_all.end(), add.begin(), add.end());
Msg::Info("IT %3d %8d points added, timings %5.2f %5.2f %5.2f %5.2f %5.2f %5d",
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index d689a6656879f4504bf93bcaa909e98cf2bc0d48..f11748bdc53e8c4847356b91136fb06b2716b833 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -543,8 +543,14 @@ void meshGEdge::operator() (GEdge *ge)
return;
}
- Msg::Info("Meshing curve %d (%s)", ge->tag(), ge->getTypeString().c_str());
-
+ if (ge->model()->getNumEdges() > 1000){
+ if (ge->tag() % 1000 == 1){
+ Msg::Info("Meshing curve %d/%d (%s)", ge->tag(), ge->model()->getNumEdges(), ge->getTypeString().c_str());
+ }
+ }
+ else {
+ Msg::Info("Meshing curve %d (%s)", ge->tag(), ge->getTypeString().c_str());
+ }
// compute bounds
Range<double> bounds = ge->parBounds(0);
double t_begin = bounds.low();
@@ -576,7 +582,7 @@ void meshGEdge::operator() (GEdge *ge)
int N;
int filterMinimumN = 1;
if(length == 0. && CTX::instance()->mesh.toleranceEdgeLength == 0.){
- Msg::Warning("Curve %d has a zero length", ge->tag());
+ Msg::Debug("Curve %d has a zero length", ge->tag());
a = 0.;
N = 1;
}
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index b89bc6a4563af94d77c693c10b0166be7bc866b1..771dd07a3af152be5bd03c04f978e6e65b8424a7 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1407,8 +1407,9 @@ void bowyerWatsonFrontal(GFace *gf,
gf->mesh_vertices.size(), worst->getRadius());
double newPoint[2], metric[3];
//optimalPointFrontal (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
- if (optimalPointFrontalB (gf,worst,active_edge,DATA,newPoint,metric))
+ if (optimalPointFrontalB (gf,worst,active_edge,DATA,newPoint,metric)){
insertAPoint(gf, AllTris.end(), newPoint, metric, DATA, AllTris, &ActiveTris, worst);
+ }
}
/* if(ITER % 1== 0){
@@ -1419,7 +1420,6 @@ void bowyerWatsonFrontal(GFace *gf,
*/
}
-
nbSwaps = edgeSwapPass(gf, AllTris, SWCR_QUAL, DATA);
/*
char name[245];
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 3b4e541267159e44e5ca61234955c6ef22800a56..eebbf65d53cd63623eb8cb14c7e7eca1fe3b9fd2 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -43,6 +43,8 @@
#include "ANN/ANN.h"
#endif
+
+
// hybrid mesh recovery structure
class splitQuadRecovery {
std::multimap<GEntity*, std::pair<MVertex*,MFace> >_data;
@@ -639,9 +641,9 @@ void MeshDelaunayVolumeTetgen(std::vector<GRegion*> ®ions)
// now do insertion of points
if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_DEL)
- bowyerWatsonFrontalLayers(gr, false);
+ Msg::Error ("Frontal Layers Algo deprecated");
else if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_HEX)
- bowyerWatsonFrontalLayers(gr, true);
+ Msg::Error ("Frontal Layers Algo deprecated");
else if(CTX::instance()->mesh.algo3d == ALGO_3D_MMG3D){
refineMeshMMG(gr);
}
@@ -652,6 +654,7 @@ void MeshDelaunayVolumeTetgen(std::vector<GRegion*> ®ions)
insertVerticesInRegion(gr,2000000000,true);
}
}
+
// crete an initial mesh
if (sqr.buildPyramids (gr->model())){
RelocateVertices(regions, 3);
diff --git a/Mesh/meshGRegionBoundaryRecovery.cpp b/Mesh/meshGRegionBoundaryRecovery.cpp
index 2fd2c789db56a30005f6ba91f52c4d91d897e32f..55ccb56ac340b71d2fdcf76c6028fd95530495f3 100644
--- a/Mesh/meshGRegionBoundaryRecovery.cpp
+++ b/Mesh/meshGRegionBoundaryRecovery.cpp
@@ -425,8 +425,11 @@ bool tetgenmesh::reconstructmesh(void *p)
// Construct a map from points to subfaces.
makepoint2submap(subfaces, idx2shlist, shperverlist);
+ // printf("coucou1\n");
+
// Process the set of PSC edges.
// Remeber that all segments have default marker '-1'.
+ // int COUNTER = 0;
for (std::list<GEdge*>::iterator it = e_list.begin(); it != e_list.end();
++it) {
GEdge *ge = *it;
@@ -480,7 +483,7 @@ bool tetgenmesh::reconstructmesh(void *p)
// It is a dangling segment (not belong to any facets).
// Check if segment [p[0],p[1]] already exists.
// TODO: Change the brute-force search. Slow!
- point *ppt;
+ /* point *ppt;
subsegs->traversalinit();
segloop.sh = shellfacetraverse(subsegs);
while (segloop.sh != NULL) {
@@ -492,7 +495,7 @@ bool tetgenmesh::reconstructmesh(void *p)
break;
}
segloop.sh = shellfacetraverse(subsegs);
- }
+ }*/
if (newseg.sh == NULL) {
makeshellface(subsegs, &newseg);
setshvertices(newseg, p[0], p[1], NULL);
@@ -522,9 +525,12 @@ bool tetgenmesh::reconstructmesh(void *p)
// Boundary recovery.
clock_t t;
+ // printf("coucou2\n");
recoverboundary(t);
+ // printf("coucou3\n");
carveholes();
+ // printf("coucou4\n");
if (subvertstack->objects > 0l) {
suppresssteinerpoints();
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index eedd1c829ff628d9c58eff060831954fb91ce770..e6054d0a4e47112439a3004ca6b5dc3baf057d7f 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -22,7 +22,6 @@
#include "MLine.h"
int MTet4::radiusNorm = 2;
-static double LIMIT_ = 1;
static void createAllEmbeddedEdges (GRegion *gr, std::set<MEdge, Less_Edge> &allEmbeddedEdges)
{
@@ -45,18 +44,6 @@ static void createAllEmbeddedFaces (GRegion *gr, std::set<MFace, Less_Face> &all
}
}
-static bool isActive(MTet4 *t, double limit_, int &active)
-{
- if (t->isDeleted()) return false;
- for (active = 0; active < 4; active++){
- MTet4 *neigh = t->getNeigh(active);
- if (!neigh || (neigh->getRadius() < limit_ && neigh->getRadius() > 0)) {
- return true;
- }
- }
- return false;
-}
-
int MTet4::inCircumSphere(const double *p) const
{
@@ -90,33 +77,29 @@ struct faceXtet{
MVertex *v1 = t1->tet()->getVertex(faces[iFac][1]);
MVertex *v2 = t1->tet()->getVertex(faces[iFac][2]);
- unsorted[0] = v0;
- unsorted[1] = v1;
- unsorted[2] = v2;
-
- v[0] = std::min(std::min(v0,v1),v2);
- v[2] = std::max(std::max(v0,v1),v2);
- v[1] = (v0 != v[0] && v0 != v[2]) ? v0 : (v1 != v[0] && v1 != v[2]) ? v1 : v2;
- //
- // std::sort(v, v + 3);
+ unsorted[0] = v[0] = v0;
+ unsorted[1] = v[1] = v1;
+ unsorted[2] = v[2] = v2;
+ MVertexLessThanNum lll;
+ std::sort(v, v + 3, lll);
}
inline MVertex * getVertex (int i) const { return unsorted[i];}
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]->getNum() < other.v[0]->getNum()) return true;
+ if (v[0]->getNum() > other.v[0]->getNum()) return false;
+ if (v[1]->getNum() < other.v[1]->getNum()) return true;
+ if (v[1]->getNum() > other.v[1]->getNum()) return false;
+ if (v[2]->getNum() < other.v[2]->getNum()) return true;
return false;
}
inline bool operator == (const faceXtet & other) const
{
- return (v[0] == other.v[0] &&
- v[1] == other.v[1] &&
- v[2] == other.v[2] );
+ return (v[0]->getNum() == other.v[0]->getNum() &&
+ v[1]->getNum() == other.v[1]->getNum() &&
+ v[2]->getNum() == other.v[2]->getNum() );
}
bool visible (MVertex *v){
MVertex* v0 = t1->tet()->getVertex(faces[i1][0]);
@@ -184,35 +167,6 @@ void connectTets(ITER beg, ITER end, std::set<MFace, Less_Face> *allEmbeddedFace
}
-static void updateActiveFaces(MTet4 *t, double limit_, std::set<MFace,Less_Face> &front)
-{
- if (t->isDeleted()) return;
- for (int active = 0; active < 4; active++){
- MTet4 *neigh = t->getNeigh(active);
- if (!neigh || (neigh->getRadius() < limit_ && neigh->getRadius() > 0)) {
- faceXtet fxt (t,active);
- MFace me (fxt.v[0],fxt.v[1],fxt.v[2]);
- front.insert(me);
- }
- }
-}
-
-static bool isActive(MTet4 *t, double limit_, int &i, std::set<MFace,Less_Face> *front)
-{
- if (t->isDeleted()) return false;
- for (i = 0; i < 4; i++){
- MTet4 *neigh = t->getNeigh(i);
- if (!neigh || (neigh->getRadius() < limit_ && neigh->getRadius() > 0)) {
- faceXtet fxt (t,i);
- MFace me (fxt.v[0],fxt.v[1],fxt.v[2]);
- if(front->find(me) != front->end()){
- return true;
- }
- }
- }
- return false;
-}
-
void connectTets(std::list<MTet4*> &l) { connectTets(l.begin(), l.end()); }
void connectTets(std::vector<MTet4*> &l) { connectTets(l.begin(), l.end()); }
@@ -367,25 +321,17 @@ void printTets (const char *fn, std::list<MTet4*> &cavity, bool force = false )
bool insertVertexB(std::list<faceXtet> &shell,
std::list<MTet4*> &cavity,
MVertex *v,
+ double lc1,
+ double lc2,
+ std::vector<double> &vSizes,
+ std::vector<double> &vSizesBGM,
MTet4 *t,
MTet4Factory &myFactory,
- std::set<MTet4*,compareTet4Ptr> &allTets,
- std::vector<double> & vSizes,
- std::vector<double> & vSizesBGM,
- std::set<MTet4*,compareTet4Ptr> *activeTets = 0 )
+ std::set<MTet4*,compareTet4Ptr> &allTets)
{
std::vector<faceXtet> conn;
std::vector<MTet4*> new_cavity;
- // check that volume is conserved
- double newVolume = 0;
- double oldVolume = 0;
- std::list<MTet4*>::iterator ittet = cavity.begin();
- std::list<MTet4*>::iterator ittete = cavity.end();
- while (ittet != ittete){
- oldVolume += fabs((*ittet)->getVolume());
- ++ittet;
- }
MTet4** newTets = new MTet4*[shell.size()];
int k = 0;
@@ -395,31 +341,18 @@ bool insertVertexB(std::list<faceXtet> &shell,
bool onePointIsTooClose = false;
while (it != shell.end()){
MTetrahedron *tr = new MTetrahedron(it->getVertex(0), it->getVertex(1), it->getVertex(2), v);
-
- double lc = .25 * (vSizes[tr->getVertex(0)->getIndex()] +
- vSizes[tr->getVertex(1)->getIndex()] +
- vSizes[tr->getVertex(2)->getIndex()] +
- vSizes[tr->getVertex(3)->getIndex()]);
- double lcBGM = .25 * (vSizesBGM[tr->getVertex(0)->getIndex()] +
- vSizesBGM[tr->getVertex(1)->getIndex()] +
- vSizesBGM[tr->getVertex(2)->getIndex()] +
- vSizesBGM[tr->getVertex(3)->getIndex()]);
- double LL = std::min(lc, lcBGM);
-
- MTet4 *t4 = myFactory.Create(tr, vSizes, vSizesBGM);
+ MTet4 *t4 = myFactory.Create(tr, vSizes, vSizesBGM, lc1, lc2);
t4->setOnWhat(t->onWhat());
- double d1 = sqrt((it->v[0]->x() - v->x()) * (it->v[0]->x() - v->x()) +
- (it->v[0]->y() - v->y()) * (it->v[0]->y() - v->y()) +
- (it->v[0]->z() - v->z()) * (it->v[0]->z() - v->z()));
- double d2 = sqrt((it->v[1]->x() - v->x()) * (it->v[1]->x() - v->x()) +
- (it->v[1]->y() - v->y()) * (it->v[1]->y() - v->y()) +
- (it->v[1]->z() - v->z()) * (it->v[1]->z() - v->z()));
- double d3 = sqrt((it->v[2]->x() - v->x()) * (it->v[2]->x() - v->x()) +
- (it->v[2]->y() - v->y()) * (it->v[2]->y() - v->y()) +
- (it->v[2]->z() - v->z()) * (it->v[2]->z() - v->z()));
-
- if (d1 < LL * .05 || d2 < LL * .05 || d3 < LL * .05) onePointIsTooClose = true;
+ double d1 = distance (it->v[0],v);
+ double d2 = distance (it->v[1],v);
+ double d3 = distance (it->v[2],v);
+
+ double lc = Extend2dMeshIn3dVolumes() ? std::min(lc1, lc2) : lc2;
+ if (d1 < lc * .05 || d2 < lc * .05 || d3 < lc * .05) {
+ // printf("coucou %g %g %g %g\n",lc,d1,d2,d3);
+ onePointIsTooClose = true;
+ }
newTets[k++] = t4;
// all new tets are pushed front in order to ba able to destroy
@@ -431,44 +364,19 @@ bool insertVertexB(std::list<faceXtet> &shell,
if (otherSide)
new_cavity.push_back(otherSide);
- // if (!it->t1->isDeleted())throw;
- newVolume += fabs(t4->getVolume());
++it;
}
- // OK, the cavity is star shaped
- // if (onePointIsTooClose)printf("One point is too close\n");
- //if (fabs(oldVolume - newVolume) > 1.e-10 * oldVolume)printf("Volume do not match %22.15E %22.15E %22.15E\n",oldVolume,newVolume,fabs(oldVolume-newVolume)/newVolume);
- // if (!onePointIsTooClose){
- if (fabs(oldVolume - newVolume) < 1.e-10 * oldVolume &&
- !onePointIsTooClose){
- connectTets_vector2(new_cavity,conn);
+ if (!onePointIsTooClose){
+ connectTets_vector2(new_cavity,conn);
allTets.insert(newTets, newTets + shell.size());
-
- if (activeTets){
- for (std::vector<MTet4*>::iterator i = new_cavity.begin(); i != new_cavity.end(); ++i){
- int active_face;
- if(isActive(*i, LIMIT_, active_face) && (*i)->getRadius() > LIMIT_){
- if ((*activeTets).find(*i) == (*activeTets).end())
- (*activeTets).insert(*i);
- }
- }
- }
-
delete [] newTets;
-
return true;
}
- else { // The cavity is NOT star shaped
- // printf("hola %d %g %g %22.15E\n",onePointIsTooClose, oldVolume, newVolume, 100.*fabs(oldVolume - newVolume) / oldVolume);
- // new_cavity.clear();
- // for (unsigned int i = 0; i <shell.size(); i++) new_cavity.push_back(newTets[i]);
- // printTets ("oldCavity.pos",cavity,true);
- // printTets ("newCavity.pos",new_cavity);
- // Msg::Fatal("");
+ else /* one point is too close */ {
for (unsigned int i = 0; i <shell.size(); i++) myFactory.Free(newTets[i]);
delete [] newTets;
- ittet = cavity.begin();
- ittete = cavity.end();
+ std::list<MTet4*>::iterator ittet = cavity.begin();
+ std::list<MTet4*>::iterator ittete = cavity.end();
while(ittet != ittete){
(*ittet)->setDeleted(false);
++ittet;
@@ -477,24 +385,8 @@ bool insertVertexB(std::list<faceXtet> &shell,
}
}
-bool insertVertex(MVertex *v,
- MTet4 *t,
- MTet4Factory &myFactory,
- std::set<MTet4*,compareTet4Ptr> &allTets,
- std::vector<double> & vSizes,
- std::vector<double> & vSizesBGM,
- std::set<MTet4*,compareTet4Ptr> *activeTets = 0 )
-{
- std::list<faceXtet> shell;
- std::list<MTet4*> cavity;
-
- findCavity(shell, cavity, v, t);
-
- return insertVertexB(shell,cavity,v,t,myFactory,allTets,vSizes,vSizesBGM,activeTets);
-}
-
-static void setLcs(MTriangle *t, std::map<MVertex*, double> &vSizes,
- std::set<MVertex*> &bndVertices)
+static void setLcs(MTriangle *t, std::map<MVertex*, double, MVertexLessThanNum> &vSizes,
+ std::set<MVertex*,MVertexLessThanNum> &bndVertices)
{
for(int i = 0; i < 3; i++){
@@ -506,8 +398,8 @@ static void setLcs(MTriangle *t, std::map<MVertex*, double> &vSizes,
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);
+ std::map<MVertex*,double, MVertexLessThanNum>::iterator iti = vSizes.find(vi);
+ std::map<MVertex*,double, MVertexLessThanNum>::iterator itj = vSizes.find(vj);
// use largest edge length
if (iti == vSizes.end() || iti->second < l) vSizes[vi] = l;
if (itj == vSizes.end() || itj->second < l) vSizes[vj] = l;
@@ -539,8 +431,8 @@ static void setLcs(MTriangle *t, std::map<MVertex*, double> &vSizes,
*/
}
-static void setLcs(MTetrahedron *t, std::map<MVertex*, double> &vSizes,
- std::set<MVertex*> &bndVertices)
+static void setLcs(MTetrahedron *t, std::map<MVertex*, double,MVertexLessThanNum> &vSizes,
+ std::set<MVertex*,MVertexLessThanNum> &bndVertices)
{
for (int i = 0; i < 4; i++){
for (int j = i + 1; j < 4; j++){
@@ -550,10 +442,10 @@ static void setLcs(MTetrahedron *t, std::map<MVertex*, double> &vSizes,
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);
- std::set<MVertex*>::iterator itvi = bndVertices.find(vi);
- std::set<MVertex*>::iterator itvj = bndVertices.find(vj);
+ std::map<MVertex*,double,MVertexLessThanNum>::iterator iti = vSizes.find(vi);
+ std::map<MVertex*,double,MVertexLessThanNum>::iterator itj = vSizes.find(vj);
+ std::set<MVertex*,MVertexLessThanNum>::iterator itvi = bndVertices.find(vi);
+ std::set<MVertex*,MVertexLessThanNum>::iterator itvj = bndVertices.find(vj);
// smallest tet edge
if (itvi == bndVertices.end() &&
(iti == vSizes.end() || iti->second > l)) vSizes[vi] = l;
@@ -886,7 +778,7 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
double qq = (*it)->getQuality();
if (qq < qMin){
for (int i = 0; i < 4; i++){
- if (faceSwap(newTets, *it, i, qm)){
+ if (0 && faceSwap(newTets, *it, i, qm)){
nbFSwap++;
break;
}
@@ -1138,8 +1030,8 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
createAllEmbeddedEdges (gr, allEmbeddedEdges);
{ // leave this in a block so the map gets deallocated directly
- std::map<MVertex*, double> vSizesMap;
- std::set<MVertex*> bndVertices;
+ std::map<MVertex*, double,MVertexLessThanNum> vSizesMap;
+ std::set<MVertex*,MVertexLessThanNum> bndVertices;
std::set<MEdge, Less_Edge>::iterator it = allEmbeddedEdges.begin();
while (it != allEmbeddedEdges.end()){
@@ -1149,8 +1041,8 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
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 = vSizesMap.find(vi);
- std::map<MVertex*,double>::iterator itj = vSizesMap.find(vj);
+ std::map<MVertex*,double,MVertexLessThanNum>::iterator iti = vSizesMap.find(vi);
+ std::map<MVertex*,double,MVertexLessThanNum>::iterator itj = vSizesMap.find(vj);
// smallest tet edge
if (iti == vSizesMap.end() || iti->second > l) vSizesMap[vi] = l;
if (itj == vSizesMap.end() || itj->second > l) vSizesMap[vj] = l;
@@ -1165,7 +1057,7 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
}
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
setLcs(gr->tetrahedra[i], vSizesMap, bndVertices);
- for(std::map<MVertex*, double>::iterator it = vSizesMap.begin();
+ for(std::map<MVertex*, double, MVertexLessThanNum>::iterator it = vSizesMap.begin();
it != vSizesMap.end(); ++it){
it->first->setIndex(NUM++);
vSizes.push_back(it->second);
@@ -1245,6 +1137,7 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
double t1 = Cpu();
// MAIN LOOP IN DELAUNAY INSERTION STARTS HERE
+
while(1){
if (COUNT_MISS_2 > 100000)break;
if (ITER >= maxVert)break;
@@ -1281,8 +1174,9 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
base->getVertex(3)->y(),
base->getVertex(3)->z()};
- tetcircumcenter(pa,pb,pc,pd, center,&uvw[0],&uvw[1],&uvw[2] );
-
+ tetcircumcenter(pa,pb,pc,pd, center, NULL, NULL, NULL);
+
+
//// A TEST !!!
std::list<faceXtet> shell;
std::list<MTet4*> cavity;
@@ -1325,14 +1219,18 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
NB_CORRECTION_OF_CAVITY ++;
if(!isCavityCompatibleWithEmbeddedEdges(cavity, shell, allEmbeddedEdges)){
correctedCavityIncompatibleWithEmbeddedEdge = true;
- // printf("ARGH\n");
}
}
- // printTets ("after.pos", cavity, true);
- // compute mesh spacing there
- if(correctedCavityIncompatibleWithEmbeddedEdge || !starShaped || !insertVertexB(shell,cavity,v, worst, myFactory, allTets, vSizes,vSizesBGM)){
+ double lc1 =
+ (1 - uvw[0] - uvw[1] - uvw[2]) * vSizes[worst->tet()->getVertex(0)->getIndex()] +
+ uvw[0] * vSizes[worst->tet()->getVertex(1)->getIndex()] +
+ uvw[1] * vSizes[worst->tet()->getVertex(2)->getIndex()] +
+ uvw[2] * vSizes[worst->tet()->getVertex(3)->getIndex()];
+ double lc2 = BGM_MeshSize(gr, 0, 0, center[0], center[1], center[2]);
+
+
+ if(correctedCavityIncompatibleWithEmbeddedEdge || !starShaped || !insertVertexB(shell,cavity,v, lc1, lc2, vSizes, vSizesBGM, worst, myFactory, allTets)){
COUNT_MISS_1++;
- // printf("coucou 1 %d %d\n",correctedCavityIncompatibleWithEmbeddedEdge, starShaped);
myFactory.changeTetRadius(allTets.begin(), 0.);
for (std::list<MTet4*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc)
(*itc)->setDeleted(false);
@@ -1340,16 +1238,8 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
NUM--;
}
else{
- //}
- double lc1 =
- (1 - uvw[0] - uvw[1] - uvw[2]) * vSizes[worst->tet()->getVertex(0)->getIndex()] +
- uvw[0] * vSizes[worst->tet()->getVertex(1)->getIndex()] +
- uvw[1] * vSizes[worst->tet()->getVertex(2)->getIndex()] +
- uvw[2] * vSizes[worst->tet()->getVertex(3)->getIndex()];
- double lc = BGM_MeshSize(gr, 0, 0, center[0], center[1], center[2]);
- // double lc = std::min(lc1, BGM_MeshSize(gr, 0, 0, center[0], center[1], center[2]));
vSizes.push_back(lc1);
- vSizesBGM.push_back(lc);
+ vSizesBGM.push_back(lc2);
REALCOUNT++;
v->onWhat()->mesh_vertices.push_back(v);
}
@@ -1409,222 +1299,6 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
}
}
-MVertex * optimalPointFrontal(GRegion *gr,
- MTet4 *worst,
- int active_face,
- std::vector<double> &vSizes,
- std::vector<double> &vSizesBGM)
-{
- double centerTet[3], centerFace[3];
- faceXtet fxt ( worst, active_face );
- double pa[3] = {fxt.v[0]->x(),fxt.v[0]->y(),fxt.v[0]->z()};
- double pb[3] = {fxt.v[1]->x(),fxt.v[1]->y(),fxt.v[1]->z()};
- double pc[3] = {fxt.v[2]->x(),fxt.v[2]->y(),fxt.v[2]->z()};
- circumCenterXYZ(pa, pb, pc, centerFace);
- worst->circumcenter(centerTet);
-
- SVector3 dir (centerTet[0] - centerFace[0],
- centerTet[1] - centerFace[1],
- centerTet[2] - centerFace[2]);
- dir.normalize();
-
- SVector3 rDir (pa[0] - centerFace[0],
- pa[1] - centerFace[1],
- pa[2] - centerFace[2]);
-
- // const double a = 2*p *3 /sqrt(3);
- const double a = .025;
- const double tt = a*sqrt(6.)/3;
-
- MVertex *vert = new MVertex(centerFace[0] + tt * dir[0],
- centerFace[1] + tt * dir[1],
- centerFace[2] + tt * dir[2],
- gr);
- return vert;
-}
-
-
-void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
-{
- std::vector<double> vSizes;
- std::vector<double> vSizesBGM;
- MTet4Factory myFactory(1600000);
- std::set<MTet4*, compareTet4Ptr> &allTets = myFactory.getAllTets();
- std::set<MTet4*, compareTet4Ptr> activeTets;
- int NUM = 0;
-
-// if (!backgroundMesh::current()) {
-// // TODO !!!
-// }
-
- if (hex){
- LIMIT_ = sqrt(2.) * .99;
- MTet4::radiusNorm =-1;
- }
-
- { // leave this in a block so the map gets deallocated directly
- std::map<MVertex*, double> vSizesMap;
- std::set<MVertex*> bndVertices;
- for(GModel::fiter it = gr->model()->firstFace(); it != gr->model()->lastFace(); ++it){
- GFace *gf = *it;
- for(unsigned int i = 0; i < gf->triangles.size(); i++){
- setLcs(gf->triangles[i], vSizesMap, bndVertices);
- }
- }
- for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
- setLcs(gr->tetrahedra[i], vSizesMap, bndVertices);
- for(std::map<MVertex*, double>::iterator it = vSizesMap.begin();
- it != vSizesMap.end(); ++it){
- it->first->setIndex(NUM++);
- vSizes.push_back(it->second);
- vSizesBGM.push_back(it->second);
- }
- }
-
- for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
- allTets.insert(myFactory.Create(gr->tetrahedra[i], vSizes,vSizesBGM));
-
- gr->tetrahedra.clear();
- connectTets(allTets.begin(), allTets.end());
-
- fs_cont search;
- buildFaceSearchStructure(gr->model(), search);
-
- for(MTet4Factory::iterator it = allTets.begin(); it != allTets.end(); ++it){
- if(!(*it)->onWhat()){
- std::list<MTet4*> theRegion;
- std::set<GFace *> faces_bound;
- GRegion *bidon = (GRegion*)123;
- double _t1 = Cpu();
- Msg::Debug("start with a non classified tet");
- non_recursive_classify(*it, theRegion, faces_bound, bidon, gr->model(), search);
- double _t2 = Cpu();
- Msg::Debug("found %d tets with %d faces (%g sec for the classification)",
- theRegion.size(), faces_bound.size(), _t2 - _t1);
- GRegion *myGRegion = getRegionFromBoundingFaces(gr->model(), faces_bound);
- // Msg::Info("a region is found %p",myGRegion);
- if(myGRegion) // a geometrical region associated to the list of faces has been found
- for(std::list<MTet4*>::iterator it2 = theRegion.begin();
- it2 != theRegion.end(); ++it2) (*it2)->setOnWhat(myGRegion);
- else // the tets are in the void
- for(std::list<MTet4*>::iterator it2 = theRegion.begin();
- it2 != theRegion.end(); ++it2)(*it2)->setDeleted(true);
- }
- }
- search.clear();
-
- for(MTet4Factory::iterator it = allTets.begin(); it!=allTets.end(); ++it){
- (*it)->setNeigh(0, 0);
- (*it)->setNeigh(1, 0);
- (*it)->setNeigh(2, 0);
- (*it)->setNeigh(3, 0);
- }
- std::set<MFace, Less_Face> allEmbeddedFaces;
- createAllEmbeddedFaces (gr, allEmbeddedFaces);
- connectTets(allTets.begin(),allTets.end(),&allEmbeddedFaces);
-
- int ITER = 0, active_face;
-
- std::set<MFace,Less_Face> _front;
- for(MTet4Factory::iterator it = allTets.begin(); it!=allTets.end(); ++it){
- if(isActive(*it,LIMIT_,active_face) && (*it)->getRadius() > LIMIT_){
- activeTets.insert(*it);
- updateActiveFaces(*it, LIMIT_, _front);
- }
- else if ((*it)->getRadius() < LIMIT_)break;
- }
-
- // insert points
- int ITERATION = 1;
- while (1){
- if (ITERATION == 8)break;
- ITERATION ++;
-
- std::set<MTet4*, compareTet4Ptr> activeTetsNotInFront;
-
- while (1){
-
- if (!activeTets.size())break;
-
- // printf("%d active tets %d tets\n",activeTets.size(),allTets.size());
-
- std::set<MTet4*,compareTet4Ptr>::iterator WORST_ITER = activeTets.begin();
- MTet4 *worst = (*WORST_ITER);
- if(worst->isDeleted()){
- activeTets.erase(WORST_ITER);
- // myFactory.Free(worst);
- }
- else {
- activeTets.erase(WORST_ITER);
- if (isActive(worst, LIMIT_, active_face,&_front) &&
- worst->getRadius() > LIMIT_){
- // printf("worst = %12.5E\n",worst->getRadius());
-
- if(ITER++ % 5000 == 0)
- Msg::Info("%d points created -- Worst tet radius is %g",
- vSizes.size(), worst->getRadius());
-
- MVertex *v = optimalPointFrontal (gr,worst,active_face,vSizes,vSizesBGM);
- v->setIndex(NUM++);
- vSizes.push_back(.025);
- vSizesBGM.push_back(.025);
-
- if(!worst->inCircumSphere(v) ||
- !insertVertex(v, worst, myFactory, allTets, vSizes,vSizesBGM,&activeTets)){
- myFactory.changeTetRadius(allTets.begin(), 0.);
- if(v) delete v;
- }
- else{
- // printf("yeah ! one new vertex \n");
- v->onWhat()->mesh_vertices.push_back(v);
- // if (ITER == 100)break;
- }
- }
- else if (worst->getRadius() > LIMIT_){
- activeTetsNotInFront.insert(worst);
- }
- }
- }
- _front.clear();
- MTet4Factory::iterator it = activeTetsNotInFront.begin();
- for ( ; it!=activeTetsNotInFront.end();++it){
- if((*it)->getRadius() > LIMIT_ && isActive(*it,LIMIT_,active_face)){
- activeTets.insert(*it);
- updateActiveFaces(*it, LIMIT_, _front);
- }
- }
- if (!activeTets.size())break;
- }
-
-
- int nbReloc = 0;
- for (int SM=0;SM<CTX::instance()->mesh.nbSmoothing;SM++){
- for(MTet4Factory::iterator it = allTets.begin(); it != allTets.end(); ++it){
- if (!(*it)->isDeleted()){
- double qq = (*it)->getQuality();
- if (qq < .4)
- for (int i = 0; i < 4; i++){
- if (smoothVertex(*it, i, qmTetrahedron::QMTET_GAMMA)) nbReloc++;
- }
- }
- }
- }
-
-
- while(1){
- if(allTets.begin() == allTets.end()) break;
- MTet4 *worst = *allTets.begin();
- if(!worst->isDeleted()){
- worst->onWhat()->tetrahedra.push_back(worst->tet());
- worst->tet() = 0;
- }
- myFactory.Free(worst);
- allTets.erase(allTets.begin());
- }
- MTet4::radiusNorm = 2;
- LIMIT_ = 1;
-}
-
///// do a 3D delaunay mesh assuming a set of vertices
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index cc75283d2b84dc57ea54ac1bde0551fafea3dad9..a586cd4055e7fae52eae103efe061b2aaa57a7a9 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -89,8 +89,7 @@ class MTet4
double B[4] = {v1->x(), v1->y(), v1->z()};
double C[4] = {v2->x(), v2->y(), v2->z()};
double D[4] = {v3->x(), v3->y(), v3->z()};
- double x,y,z;
- tetcircumcenter (A,B,C,D,res,&x,&y,&z);
+ tetcircumcenter (A,B,C,D,res,NULL,NULL,NULL);
}
void setup(MTetrahedron *t, std::vector<double> &sizes, std::vector<double> &sizesBGM)
@@ -103,6 +102,19 @@ class MTet4
const double dy = base->getVertex(0)->y() - center[1];
const double dz = base->getVertex(0)->z() - center[2];
circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
+ /*
+ if (base->getVertex(0)->getIndex() >= sizes.size() ||
+ base->getVertex(1)->getIndex() >= sizes.size() ||
+ base->getVertex(2)->getIndex() >= sizes.size() ||
+ base->getVertex(3)->getIndex() >= sizes.size()){
+ printf("ERROR %d vs %d %d %d %d\n",sizes.size() ,
+ base->getVertex(0)->getIndex(),
+ base->getVertex(1)->getIndex(),
+ base->getVertex(2)->getIndex(),
+ base->getVertex(3)->getIndex());
+
+ }
+ */
double lc1 = 0.25*(sizes[base->getVertex(0)->getIndex()]+
sizes[base->getVertex(1)->getIndex()]+
sizes[base->getVertex(2)->getIndex()]+
@@ -115,6 +127,43 @@ class MTet4
circum_radius /= lc;
deleted = false;
}
+
+ void setup(MTetrahedron *t, std::vector<double> &sizes, std::vector<double> &sizesBGM, double lcA, double lcB)
+ {
+ base = t;
+ neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
+ double center[3];
+ circumcenter(center);
+ const double dx = base->getVertex(0)->x() - center[0];
+ const double dy = base->getVertex(0)->y() - center[1];
+ const double dz = base->getVertex(0)->z() - center[2];
+ circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
+
+ /*
+ if (base->getVertex(0)->getIndex() >= sizes.size() ||
+ base->getVertex(1)->getIndex() >= sizes.size() ||
+ base->getVertex(2)->getIndex() >= sizes.size()){
+ printf("ERROR %d vs %d %d %d %d\n",sizes.size() ,
+ base->getVertex(0)->getIndex(),
+ base->getVertex(1)->getIndex(),
+ base->getVertex(2)->getIndex(),
+ base->getVertex(3)->getIndex());
+
+ }
+ */
+ double lc1 = 0.25*(sizes[base->getVertex(0)->getIndex()]+
+ sizes[base->getVertex(1)->getIndex()]+
+ sizes[base->getVertex(2)->getIndex()]+
+ lcA);
+ double lcBGM = 0.25*(sizesBGM[base->getVertex(0)->getIndex()]+
+ sizesBGM[base->getVertex(1)->getIndex()]+
+ sizesBGM[base->getVertex(2)->getIndex()]+
+ lcB);
+ double lc = Extend2dMeshIn3dVolumes() ? std::min(lc1, lcBGM) : lcBGM;
+ circum_radius /= lc;
+ deleted = false;
+ }
+
inline GRegion *onWhat() const { return gr; }
inline void setOnWhat(GRegion *g) { gr = g; }
inline bool isDeleted() const { return deleted; }
@@ -181,7 +230,7 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex);
GRegion *getRegionFromBoundingFaces(GModel *model,
std::set<GFace *> &faces_bound);
-class compareTet4Ptr
+class compareTet4Ptr
{
public:
inline bool operator () (const MTet4 *a, const MTet4 *b) const
@@ -245,6 +294,18 @@ class MTet4Factory
t4->setup(t, sizes, sizesBGM);
return t4;
}
+ MTet4 *Create(MTetrahedron * t, std::vector<double> &sizes,
+ std::vector<double> &sizesBGM, double lc1, double lc2)
+ {
+#ifdef _GMSH_PRE_ALLOCATE_STRATEGY_
+ MTet4 *t4 = getAnEmptySlot();
+#else
+ MTet4 *t4 = new MTet4;
+#endif
+ t4->setup(t, sizes, sizesBGM, lc1, lc2);
+ return t4;
+ }
+
void Free(MTet4 *t)
{
if (t->tet()) delete t->tet();
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
index e4cc66e88c27df3bd8e2996db2ed3965c6aa9bda..7091b9c8e3c79d2aba80947e3407d8bbdfb8b323 100644
--- a/Mesh/meshGRegionLocalMeshMod.cpp
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -206,6 +206,12 @@ bool edgeSwap(std::vector<MTet4 *> &newTets,
int iLocalEdge,
const qmTetrahedron::Measures &cr)
{
+
+ //static int edges[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
+ int permut[6] = {0,3,1,2,5,4};
+ iLocalEdge = permut[iLocalEdge];
+
+
std::vector<MTet4*> cavity;
std::vector<MTet4*> outside;
std::vector<MVertex*> ring;
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index ac01ba6ad64998b8e8311b920fab11cd29a139e8..c8cc6c0e9e78a5a0db51b0c5e81b0955af6ad072 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -665,9 +665,13 @@ bool newton_fd(bool (*func)(fullVector<double> &, fullVector<double> &, void *),
fullMatrix<double> J(N, N);
fullVector<double> f(N), feps(N), dx(N);
- for (int iter = 0; iter < MAXIT; iter++){
- if(!func(x, f, data)) return false;
-
+
+ for (int iter = 0; iter < MAXIT; iter++){
+ if (x.norm() > 1.e6)return false;
+ if(!func(x, f, data)) {
+ return false;
+ }
+
bool isZero = false;
for (int k=0; k<N; k++) {
if (f(k) == 0. ) isZero = true;
@@ -680,7 +684,9 @@ bool newton_fd(bool (*func)(fullVector<double> &, fullVector<double> &, void *),
double h = EPS * fabs(x(j));
if(h == 0.) h = EPS;
x(j) += h;
- if(!func(x, feps, data)) return false;
+ if(!func(x, feps, data)) {
+ return false;
+ }
for (int i = 0; i < N; i++){
J(i, j) = (feps(i) - f(i)) / h;
}
@@ -690,13 +696,16 @@ bool newton_fd(bool (*func)(fullVector<double> &, fullVector<double> &, void *),
if (N == 1)
dx(0) = f(0) / J(0, 0);
else
- if (!J.luSolve(f, dx))
+ if (!J.luSolve(f, dx)){
return false;
+ }
for (int i = 0; i < N; i++)
x(i) -= relax * dx(i);
- if(dx.norm() < tolx) return true;
+ if(dx.norm() < tolx) {
+ return true;
+ }
}
return false;
}
diff --git a/benchmarks/boolean/aneurysmBL.py b/benchmarks/boolean/aneurysmBL.py
index c5b15484aef5249d1751ff0f51d1bc34279b4ce1..e1bb89787a96415fc3c435dbbd2cadbe30b1aa66 100644
--- a/benchmarks/boolean/aneurysmBL.py
+++ b/benchmarks/boolean/aneurysmBL.py
@@ -1,4 +1,4 @@
-from dgpy import *
+from gmshpy import *
GmshSetOption('General', 'Verbosity', 99.0)
GmshSetOption('Mesh', 'SaveAll', 1.0)