diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 82a10bf31fb33cfe208ee1b64fd5c017f4696c8f..e4603f866c7a4a6c3c34dbfef6e182f59d482f0b 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -952,6 +952,12 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[
Range<double> uu = parBounds(0);
Range<double> vv = parBounds(1);
+ {
+ GPoint pnt = point(initialGuess[0], initialGuess[1]);
+ SPoint3 spnt(pnt.x(), pnt.y(), pnt.z());
+ double dist = queryPoint.distance(spnt);
+ printf("dist = %12.5E\n",dist);
+ }
// FILE *F = fopen ("hop.pos","w");
// fprintf(F,"View \" \" {\n");
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 92d70237321295ba7ba2100c71892751db477221..d4ad6538c7e63e9023a6ac930c6ad5f850f7090a 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -233,6 +233,8 @@ void MQuadrangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
double MQuadrangle::etaShapeMeasure()
{
+ double AR = 1;//(minEdge()/maxEdge());
+
SVector3 v01 (_v[1]->x()-_v[0]->x(),_v[1]->y()-_v[0]->y(),_v[1]->z()-_v[0]->z());
SVector3 v12 (_v[2]->x()-_v[1]->x(),_v[2]->y()-_v[1]->y(),_v[2]->z()-_v[1]->z());
SVector3 v23 (_v[3]->x()-_v[2]->x(),_v[3]->y()-_v[2]->y(),_v[3]->z()-_v[2]->z());
@@ -262,7 +264,7 @@ double MQuadrangle::etaShapeMeasure()
angle = std::max(fabs(90. - a3),angle);
angle = std::max(fabs(90. - a4),angle);
- return sign*(1.-angle/90);
+ return sign*(1.-angle/90) * AR;
}
/// a shape measure for quadrangles
@@ -274,46 +276,6 @@ double MQuadrangle::etaShapeMeasure()
/// 1 + xi , 1 - xi , -(1-xi), -(1+xi)
double MQuadrangle::gammaShapeMeasure(){
return etaShapeMeasure();
- /*
- // xi = -1 eta = -1
- const double dsf_corner1 [2][4] = {{0,0,-.5,.5},{0,0.5,-.5,0}};
- // xi = 1 eta = -1
- const double dsf_corner2 [2][4] = {{0,0,-.5,.5},{0.5,0,0,-.5}};
- // xi = 1 eta = 1
- const double dsf_corner3 [2][4] = {{.5,-.5,0,0},{0.5,0,0,-.5}};
- // xi = -1 eta = 1
- const double dsf_corner4 [2][4] = {{0,0,-.5,.5},{0.5,0,0,-.5}};
- */
- double QT[4] = {qmTriangle(_v[0],_v[1],_v[2],QMTRI_RHO),
- qmTriangle(_v[1],_v[2],_v[3],QMTRI_RHO),
- qmTriangle(_v[2],_v[3],_v[0],QMTRI_RHO),
- qmTriangle(_v[3],_v[0],_v[1],QMTRI_RHO)} ;
- std::sort(QT,QT+4);
- double quality = QT[0]*QT[1] / (QT[2] * QT[3]);
-
- SVector3 v01 (_v[1]->x()-_v[0]->x(),_v[1]->y()-_v[0]->y(),_v[1]->z()-_v[0]->z());
- SVector3 v12 (_v[2]->x()-_v[1]->x(),_v[2]->y()-_v[1]->y(),_v[2]->z()-_v[1]->z());
- SVector3 v23 (_v[3]->x()-_v[2]->x(),_v[3]->y()-_v[2]->y(),_v[3]->z()-_v[2]->z());
- SVector3 v30 (_v[0]->x()-_v[3]->x(),_v[0]->y()-_v[3]->y(),_v[0]->z()-_v[3]->z());
-
- SVector3 a = crossprod(v01,v12);
- SVector3 b = crossprod(v12,v23);
- SVector3 c = crossprod(v23,v30);
- SVector3 d = crossprod(v30,v01);
-
- if (a.z() < 0 || b.z() < 0 || c.z() < 0 || d.z() < 0) return -quality;
-
- if (dot(a,b) < 0 || dot(a,c) < 0 || dot(a,d) < 0 )return -quality;
-
- return quality;
- /*
- double J[3][3];
- double detJ = getJacobian(-1,-1, J);
- double C[2][2] = {{0,0}{0,0}};
- for (int i=0;i<2;i++){
-
- }*/
-
}
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index a4750c55f927294a1ab3377847b1b2bdb862c445..d72e946bd9bd2f4bd9fa32f09e0bb643b38ab679 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -533,8 +533,9 @@ backgroundMesh::backgroundMesh(GFace *_gf, bool cfd)
_octree = new MElementOctree(_triangles);
// compute the mesh sizes at nodes
- if (CTX::instance()->mesh.lcFromPoints)
+ if (CTX::instance()->mesh.lcFromPoints){
propagate1dMesh(_gf);
+ }
else {
std::map<MVertex*, MVertex*>::iterator itv2 = _2Dto3D.begin();
for ( ; itv2 != _2Dto3D.end(); ++itv2){
@@ -623,7 +624,8 @@ static void propagateValuesOnFace(GFace *_gf,
}
// Solve
- _lsys->systemSolve();
+ if (myAssembler.sizeOfR())
+ _lsys->systemSolve();
// save solution
for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 09dc4721b575b9baac02c16297c0773ff8b7f52b..594c21fa3b8ec3e4690c767a03ca398e01ad372e 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -599,7 +599,7 @@ static void Mesh3D(GModel *m)
f.treat_region(gr);
}
//Recombine3D into hex
- if(CTX::instance()->mesh.recombine3DAll || gr->meshAttributes.recombine3D){
+ if(CTX::instance()->mesh.recombine3DAll || gr->meshAttributes.recombine3D){
Recombinator rec;
rec.execute();
Supplementary sup;
diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp
index 0a88eb33496462ab93e0dabe11c01d776399588e..93418da7537eaf8408c2df284514547c27bfeb35 100644
--- a/Mesh/directions3D.cpp
+++ b/Mesh/directions3D.cpp
@@ -573,7 +573,8 @@ void Size_field::solve(GRegion* gr){
//printf("number of tetrahedra = %d\n",count2);
//printf("volume = %f\n",volume);
- system->systemSolve();
+ if (assembler.sizeOfR())
+ system->systemSolve();
for(it=interior.begin();it!=interior.end();it++){
assembler.getDofValue(*it,0,1,val);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index e06b660df5518a9644602f4465d214165f4a184b..fd35db21853044bfc25a82b9229ca64fdf01646e 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1737,6 +1737,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
outputScalarField(m->triangles, name, 1);
}
+
// start mesh generation for periodic face
if(!algoDelaunay2D(gf)){
// need for a BGM for cross field
@@ -1766,32 +1767,46 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
/// This is a structure that we need only for periodic cases
/// We will duplicate the vertices (MVertex) that are on seams
///
- /*
- {
+
+ std::map<MVertex*, MVertex*> equivalence;
+ std::map<MVertex*, SPoint2> parametricCoordinates;
+ if(algoDelaunay2D(gf)){
std::map<MVertex*, BDS_Point*> invertMap;
std::map<BDS_Point*, MVertex*>::iterator it = recoverMap.begin();
- std::map<BDS_Point*, MVertex*>::iterator it = recoverMap.begin();
while(it != recoverMap.end()){
// we have twice vertex MVertex with 2 different coordinates
- std::map<MVertex*, BDS_Point*>::iterator invIt = invertMap.find((*it)->second);
+ MVertex * mv1 = it->second;
+ BDS_Point* bds = it->first;
+ std::map<MVertex*, BDS_Point*>::iterator invIt = invertMap.find(mv1);
if (invIt != invertMap.end()){
- BDS_Point* bds1 = (*invIt)->second;
- BDS_Point* bds2 = (*it)->first;
- MVertex * mv1 = (*it)->second;
// create a new "fake" vertex that will be destroyed afterwards
- double t;
- mv1->getParameter(0,t);
- MVertex * mv2 = new MEdgeVertex (mv1->x(),mv1->y(),mv1->z(),mv1->onWhat(), t);
- (*it)->second = mv2;
- std::map<BDS_Point*, MVertex*>::iterator itTemp = it;
- ++it;
- recoverMap.erase(itTemp);
+ MVertex * mv2 ;
+ if (mv1->onWhat()->dim() == 1) {
+ double t;
+ mv1->getParameter(0,t);
+ mv2 = new MEdgeVertex (mv1->x(),mv1->y(),mv1->z(),mv1->onWhat(), t, ((MEdgeVertex*)mv1)->getLc());
+ }
+ else if (mv1->onWhat()->dim() == 0) {
+ mv2 = new MVertex (mv1->x(),mv1->y(),mv1->z(),mv1->onWhat());
+ }
+ else
+ Msg::Error("error in seam reconstruction");
+
+ it->second = mv2;
+ equivalence[mv2] = mv1;
+ parametricCoordinates[mv2] = SPoint2(bds->u,bds->v);
+ // printf("new vertex created %p that replaces %p %g %g %g\n",mv2,mv1,mv2->x(),mv2->y(),mv2->z());
+ invertMap[mv2] = bds;
+ }
+ else {
+ parametricCoordinates[mv1] = SPoint2(bds->u,bds->v);
+ invertMap[mv1] = bds;
}
- else ++it;
- invertMap[it->second] = it->first;
+ ++it;
}
+ // recoverMap.insert(new_relations.begin(), new_relations.end());
}
- */
+ Msg::Info("%d points that are duplicated for delaunay meshing",equivalence.size());
// fill the small gmsh structures
{
@@ -1847,17 +1862,29 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
outputScalarField(m->triangles, name, 1);
}
+ bool infty = false;
+ if (gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD || gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS)
+ infty = true;
+ if (infty)
+ buildBackGroundMesh (gf, &equivalence, ¶metricCoordinates);
+ // BOUNDARY LAYER
+ modifyInitialMeshForTakingIntoAccountBoundaryLayers(gf);
+
+
if(algoDelaunay2D(gf)){
if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL)
- bowyerWatsonFrontal(gf);
+ bowyerWatsonFrontal(gf, &equivalence, ¶metricCoordinates);
else if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD)
- bowyerWatsonFrontalLayers(gf,true);
+ bowyerWatsonFrontalLayers(gf,true, &equivalence, ¶metricCoordinates);
+ else if(gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS)
+ bowyerWatsonParallelograms(gf,&equivalence, ¶metricCoordinates);
else if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
gf->getMeshingAlgo() == ALGO_2D_AUTO)
- bowyerWatson(gf);
+ bowyerWatson(gf,1000000000, &equivalence, ¶metricCoordinates);
else
meshGFaceBamg(gf);
- laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
+ if (!infty || !(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine))
+ laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing, infty);
}
// delete the mesh
@@ -1872,6 +1899,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
gf->meshStatistics.best_element_shape,
gf->meshStatistics.nbTriangle,
gf->meshStatistics.nbGoodQuality);
+ gf->meshStatistics.status = GFace::DONE;
return true;
}
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 1f066c9b47fa04f46cec22813923d45fd99c2d39..6adf207b67c69b3612265a22159c49481c577899 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -322,6 +322,7 @@ MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric, bidimMeshData * data, GF
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);
+ // printf("%p %p %p %g %g %g %g %g %g %g %g %g circ %g lc %g\n",base->getVertex(0),base->getVertex(1),base->getVertex(2),pa[0],pa[1],pa[2],pb[0], pb[1],pb[2],pc[0],pc[1],pc[2],circum_radius,lc);
circum_radius /= lc;
}
else {
@@ -867,14 +868,16 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
}
}
-void bowyerWatson(GFace *gf, int MAXPNT)
+void bowyerWatson(GFace *gf, int MAXPNT,
+ std::map<MVertex* , MVertex*>* equivalence,
+ std::map<MVertex*, SPoint2> * parametricCoordinates)
{
std::set<MTri3*,compareTri3Ptr> AllTris;
- bidimMeshData DATA;
+ bidimMeshData DATA(equivalence,parametricCoordinates);
buildMeshGenerationDataStructures(gf, AllTris, DATA);
- // _printTris ("before.pos", AllTris, Us,Vs);
+ // if (equivalence)_printTris ("before.pos", AllTris.begin(), AllTris.end(), DATA);
int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, DATA);
// _printTris ("after2.pos", AllTris, Us,Vs);
Msg::Debug("Delaunization of the initial mesh done (%d swaps)", nbSwaps);
@@ -939,7 +942,7 @@ void bowyerWatson(GFace *gf, int MAXPNT)
}
}
#endif
- transferDataStructure(gf, AllTris, DATA);
+ transferDataStructure(gf, AllTris, DATA);
}
/*
@@ -1151,11 +1154,13 @@ void optimalPointFrontalB (GFace *gf,
}
}
-void bowyerWatsonFrontal(GFace *gf)
+void bowyerWatsonFrontal(GFace *gf,
+ std::map<MVertex* , MVertex*>* equivalence,
+ std::map<MVertex*, SPoint2> * parametricCoordinates)
{
std::set<MTri3*,compareTri3Ptr> AllTris;
std::set<MTri3*,compareTri3Ptr> ActiveTris;
- bidimMeshData DATA;
+ bidimMeshData DATA(equivalence,parametricCoordinates);
buildMeshGenerationDataStructures(gf, AllTris, DATA);
@@ -1324,7 +1329,9 @@ void optimalPointFrontalQuadB (GFace *gf,
return;
}
-void buildBackGroundMesh (GFace *gf)
+void buildBackGroundMesh (GFace *gf,
+ std::map<MVertex* , MVertex*>* equivalence,
+ std::map<MVertex*, SPoint2> * parametricCoordinates)
{
// printf("build bak mesh\n");
quadsToTriangles(gf, 100000);
@@ -1351,7 +1358,7 @@ void buildBackGroundMesh (GFace *gf)
int CurvControl = CTX::instance()->mesh.lcFromCurvature;
CTX::instance()->mesh.lcFromCurvature = 0;
// Do a background mesh
- bowyerWatson(gf,4000);
+ bowyerWatson(gf,4000, equivalence,parametricCoordinates);
// Re-enable curv control if asked
CTX::instance()->mesh.lcFromCurvature = CurvControl;
// apply this to the BGM
@@ -1371,12 +1378,14 @@ void buildBackGroundMesh (GFace *gf)
}
-void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
+void bowyerWatsonFrontalLayers(GFace *gf, bool quad,
+ std::map<MVertex* , MVertex*>* equivalence,
+ std::map<MVertex*, SPoint2> * parametricCoordinates)
{
std::set<MTri3*,compareTri3Ptr> AllTris;
std::set<MTri3*,compareTri3Ptr> ActiveTris;
- bidimMeshData DATA;
+ bidimMeshData DATA(equivalence, parametricCoordinates);
if (quad){
LIMIT_ = sqrt(2.) * .99;
@@ -1496,10 +1505,12 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
backgroundMesh::unset();
}
-void bowyerWatsonParallelograms(GFace *gf)
+void bowyerWatsonParallelograms(GFace *gf,
+ std::map<MVertex* , MVertex*>* equivalence,
+ std::map<MVertex*, SPoint2> * parametricCoordinates)
{
std::set<MTri3*,compareTri3Ptr> AllTris;
- bidimMeshData DATA;
+ bidimMeshData DATA(equivalence, parametricCoordinates);
std::vector<MVertex*> packed;
std::vector<SMetric3> metrics;
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index 449cb68678b953c9c8ad58961e8bc0fdc6ba9d86..00ab44a4c4718799592b711a67a67fbfe9bdb5a3 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -24,10 +24,20 @@ struct bidimMeshData
std::map<MVertex*,int> indices;
std::vector<double> Us, Vs, vSizes, vSizesBGM;
std::vector<SMetric3> vMetricsBGM;
+ std::map<MVertex* , MVertex*>* equivalence;
+ std::map<MVertex*, SPoint2> * parametricCoordinates;
inline void addVertex (MVertex* mv, double u, double v, double size, double sizeBGM){
int index = Us.size();
if (mv->onWhat()->dim() == 2)mv->setIndex(index);
else indices[mv] = index;
+ if (parametricCoordinates){
+ std::map<MVertex*, SPoint2>::iterator it = parametricCoordinates->find(mv);
+ if (it != parametricCoordinates->end()){
+ u = it->second.x();
+ v = it->second.y();
+ // printf("%g %g\n",u,v);
+ }
+ }
Us.push_back(u);
Vs.push_back(v);
vSizes.push_back(size);
@@ -37,6 +47,17 @@ struct bidimMeshData
if (mv->onWhat()->dim() == 2)return mv->getIndex();
return indices[mv];
}
+ inline MVertex * equivalent (MVertex *v1) const {
+ if (equivalence){
+ std::map<MVertex* , MVertex*>::iterator it = equivalence->find(v1);
+ if (it == equivalence->end())return 0;
+ return it->second;
+ }
+ return 0;
+ }
+ bidimMeshData (std::map<MVertex* , MVertex*>* e = 0, std::map<MVertex*, SPoint2> *p = 0) : equivalence(e), parametricCoordinates(p)
+ {
+ }
};
@@ -111,11 +132,21 @@ class compareTri3Ptr
void connectTriangles(std::list<MTri3*> &);
void connectTriangles(std::vector<MTri3*> &);
void connectTriangles(std::set<MTri3*,compareTri3Ptr> &AllTris);
-void bowyerWatson(GFace *gf, int MAXPNT= 1000000000);
-void bowyerWatsonFrontal(GFace *gf);
-void bowyerWatsonFrontalLayers(GFace *gf, bool quad);
-void bowyerWatsonParallelograms(GFace *gf);
-void buildBackGroundMesh (GFace *gf);
+void bowyerWatson(GFace *gf, int MAXPNT= 1000000000,
+ std::map<MVertex* , MVertex*>* equivalence= 0,
+ std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
+void bowyerWatsonFrontal(GFace *gf,
+ std::map<MVertex* , MVertex*>* equivalence= 0,
+ std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
+void bowyerWatsonFrontalLayers(GFace *gf, bool quad,
+ std::map<MVertex* , MVertex*>* equivalence= 0,
+ std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
+void bowyerWatsonParallelograms(GFace *gf,
+ std::map<MVertex* , MVertex*>* equivalence= 0,
+ std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
+void buildBackGroundMesh (GFace *gf,
+ std::map<MVertex* , MVertex*>* equivalence= 0,
+ std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
struct edgeXface
{
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 2e8d0a8410646fedd69dfa0a5e0786b2c0fe2d29..f008d10ff51e5f5adc6627cfcab447c884191ae4 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -23,6 +23,7 @@
#include "OS.h"
#include "SVector3.h"
#include "SPoint3.h"
+#include "robustPredicates.h"
#if defined(HAVE_BFGS)
#include "stdafx.h"
@@ -53,7 +54,7 @@ static int diff2 (MElement *q1, MElement *q2)
return 0;
}
-static int SANITY_(GFace *gf,MQuadrangle *q1, MQuadrangle *q2, int count)
+static int SANITY_(GFace *gf,MQuadrangle *q1, MQuadrangle *q2)
{
for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
MQuadrangle *e1 = gf->quadrangles[i];
@@ -118,20 +119,22 @@ static void setLcsInit(MTriangle *t, std::map<MVertex*, double> &vSizes)
}
}
-static void setLcs(MTriangle *t, std::map<MVertex*, double> &vSizes)
+static void setLcs(MTriangle *t, std::map<MVertex*, double> &vSizes, bidimMeshData & data)
{
for(int i = 0; i < 3; i++){
for(int j = i + 1; j < 3; j++){
MVertex *vi = t->getVertex(i);
MVertex *vj = t->getVertex(j);
- double dx = vi->x() - vj->x();
- double dy = vi->y() - vj->y();
- double dz = vi->z() - vj->z();
- double l = sqrt(dx * dx + dy * dy + dz * dz);
- std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
- std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
- if(iti->second < 0 || iti->second > l) iti->second = l;
- if(itj->second < 0 || itj->second > l) itj->second = l;
+ if (vi != data.equivalent(vj) && vj != data.equivalent(vi) ){
+ double dx = vi->x() - vj->x();
+ double dy = vi->y() - vj->y();
+ double dz = vi->z() - vj->z();
+ double l = sqrt(dx * dx + dy * dy + dz * dz);
+ std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
+ std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
+ if(iti->second < 0 || iti->second > l) iti->second = l;
+ if(itj->second < 0 || itj->second > l) itj->second = l;
+ }
}
}
}
@@ -147,7 +150,7 @@ void buildMeshGenerationDataStructures(GFace *gf,
setLcsInit(gf->triangles[i], vSizesMap);
for(unsigned int i = 0;i < gf->triangles.size(); i++)
- setLcs(gf->triangles[i], vSizesMap);
+ setLcs(gf->triangles[i], vSizesMap, data);
// take care of embedded vertices
{
@@ -227,6 +230,26 @@ void transferDataStructure(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris,
if(pp < 0) t->revert();
}
}
+
+ if (data.equivalence){
+ std::vector<MTriangle*> newT;
+ for (unsigned int i=0;i<gf->triangles.size();i++){
+ MTriangle *t = gf->triangles[i];
+ MVertex *v[3];
+ for (int j=0;j<3;j++){
+ v[j] = t->getVertex(j);
+ std::map<MVertex* , MVertex*>::iterator it = data.equivalence->find(v[j]);
+ if (it != data.equivalence->end()){
+ v[j] = it->second;
+ }
+ }
+ newT.push_back(new MTriangle (v[0],v[1],v[2]));
+ delete t;
+ }
+ gf->triangles = newT;
+ }
+
+
}
void buildVertexToTriangle(std::vector<MTriangle*> &eles, v2t_cont &adj)
@@ -303,7 +326,7 @@ void parametricCoordinates(MElement *t, GFace *gf, double u[4], double v[4],
}
}
-double surfaceFaceUV(MElement *t,GFace *gf, bool *concave = 0)
+double surfaceFaceUV(MElement *t,GFace *gf, bool maximal = true)
{
double u[4],v[4];
parametricCoordinates(t,gf,u,v);
@@ -316,8 +339,7 @@ double surfaceFaceUV(MElement *t,GFace *gf, bool *concave = 0)
const double a2 =
0.5*fabs((u[2]-u[1])*(v[3]-v[1])-(u[3]-u[1])*(v[2]-v[1])) +
0.5*fabs((u[0]-u[3])*(v[1]-v[3])-(u[1]-u[3])*(v[0]-v[3])) ;
- if (concave) *concave = fabs(a2-a1) > 1.e-10*(a2 + a1);
- return std::min(a2,a1);
+ return maximal ? std::max(a2,a1) : std::min(a2,a1);
}
}
@@ -469,8 +491,6 @@ static int _removeTwoQuadsNodes(GFace *gf)
if (q1->getNumVertices() == 4 &&
q2->getNumVertices() == 4 &&
touched.find(q1) == touched.end() && touched.find(q2) == touched.end()){
- touched.insert(q1);
- touched.insert(q2);
int comm = 0;
for (int i=0;i<4;i++){
if (q1->getVertex(i) == it->first){
@@ -495,8 +515,18 @@ static int _removeTwoQuadsNodes(GFace *gf)
q2->writePOS(stdout,true,false,false,false,false,false);
return 0;
}
- gf->quadrangles.push_back(new MQuadrangle(v1,v2,v3,v4));
- vtouched.insert(it->first);
+ MQuadrangle *q = new MQuadrangle(v1,v2,v3,v4);
+ double s1 = surfaceFaceUV(q,gf);
+ double s2 = surfaceFaceUV(q1,gf) + surfaceFaceUV(q2,gf);;
+ if (s1 > s2){
+ delete q;
+ }
+ else{
+ touched.insert(q1);
+ touched.insert(q2);
+ gf->quadrangles.push_back(q);
+ vtouched.insert(it->first);
+ }
}
}
it++;
@@ -549,8 +579,8 @@ static bool _isItAGoodIdeaToCollapseThatVertex (GFace *gf,
{
double surface_old = surfaceFaceUV(q,gf);
double surface_new = 0;
- // double worst_quality_old = q->etaShapeMeasure();
- // double worst_quality_new = 1.0;
+ double worst_quality_old = q->etaShapeMeasure();
+ double worst_quality_new = 1.0;
double x = v->x();
double y = v->y();
@@ -568,13 +598,13 @@ static bool _isItAGoodIdeaToCollapseThatVertex (GFace *gf,
for (unsigned int j=0;j<e1.size();++j){
surface_old += surfaceFaceUV(e1[j],gf);
- // worst_quality_old = std::min(worst_quality_old,e1[j]-> etaShapeMeasure());
+ worst_quality_old = std::min(worst_quality_old,e1[j]-> etaShapeMeasure());
for (int k=0;k<e1[j]->getNumVertices();k++){
if (e1[j]->getVertex(k) == v1 && e1[j] != q)
e1[j]->setVertex(k,v);
}
surface_new += surfaceFaceUV(e1[j],gf);
- // worst_quality_new = std::min(worst_quality_new,e1[j]-> etaShapeMeasure());
+ worst_quality_new = std::min(worst_quality_new,e1[j]-> etaShapeMeasure());
for (int k=0;k<e1[j]->getNumVertices();k++){
if (e1[j]->getVertex(k) == v && e1[j] != q)
e1[j]->setVertex(k,v1);
@@ -583,21 +613,21 @@ static bool _isItAGoodIdeaToCollapseThatVertex (GFace *gf,
for (unsigned int j=0;j<e2.size();++j){
surface_old += surfaceFaceUV(e2[j],gf);
- // worst_quality_old = std::min(worst_quality_old,e2[j]-> etaShapeMeasure());
+ worst_quality_old = std::min(worst_quality_old,e2[j]-> etaShapeMeasure());
for (int k=0;k<e2[j]->getNumVertices();k++){
if (e2[j]->getVertex(k) == v2 && e2[j] != q)
e2[j]->setVertex(k,v);
}
surface_new += surfaceFaceUV(e2[j],gf);
- // worst_quality_new = std::min(worst_quality_new,e2[j]-> etaShapeMeasure());
+ worst_quality_new = std::min(worst_quality_new,e2[j]-> etaShapeMeasure());
for (int k=0;k<e2[j]->getNumVertices();k++){
if (e2[j]->getVertex(k) == v && e2[j] != q)
e2[j]->setVertex(k,v2);
}
}
- if ( fabs (surface_old - surface_new ) > 1.e-10 * surface_old) {
- // || FACT*worst_quality_new < worst_quality_old ) {
+ if ( fabs (surface_old - surface_new ) > 1.e-10 * surface_old ){
+ // || FACT*worst_quality_new < worst_quality_old ) {
v->setXYZ(x,y,z);
v->setParameter(0,uu);
@@ -617,35 +647,31 @@ static bool _isItAGoodIdeaToMoveThatVertex (GFace *gf,
double surface_old = 0;
double surface_new = 0;
- double qualityNew[256];
-
GPoint gp = gf->point(after);
if (!gp.succeeded())return false;
SPoint3 pafter (gp.x(),gp.y(),gp.z());
SPoint3 pbefore (v1->x(),v1->y(),v1->z());
- for (unsigned int j=0;j<e1.size();++j)
- surface_old += surfaceFaceUV(e1[j],gf);
-
+ double minq = 1.0;
+ for (unsigned int j=0;j<e1.size();++j){
+ surface_old += surfaceFaceUV(e1[j],gf,false);
+ minq = std::min(e1[j]->etaShapeMeasure(),minq);
+ }
+
v1->setParameter(0,after.x());
v1->setParameter(1,after.y());
v1->setXYZ(pafter.x(),pafter.y(),pafter.z());
- bool badQuality = false;
+ double minq_new = 1.0;
for (unsigned int j=0;j<e1.size();++j){
- surface_new += surfaceFaceUV(e1[j],gf);
- qualityNew[j] = e1[j]->gammaShapeMeasure();
- if (qualityNew[j] < 0.01) {
- badQuality = true;
- break;
- }
+ surface_new += surfaceFaceUV(e1[j],gf,false);
+ minq_new = std::min(e1[j]->etaShapeMeasure(),minq_new);
}
v1->setParameter(0,before.x());
v1->setParameter(1,before.y());
v1->setXYZ(pbefore.x(),pbefore.y(),pbefore.z());
-
- if (badQuality || (surface_new - surface_old) > 1.e-10 * surface_old) {
+ if ((1.+1.e-10)*surface_old < surface_new|| minq_new < minq) {
return false;
}
return true;
@@ -661,8 +687,21 @@ static bool _isItAGoodIdeaToMoveThoseVertices (GFace *gf,
double surface_old = 0;
double surface_new = 0;
+ double umin=1.e22,umax=-1.e22,vmin=1.e22,vmax=-1.e22;
+
+ for (unsigned int i=0; i<v.size();++i){
+ umin = std::min(before[i].x(),umin);
+ vmin = std::min(before[i].y(),vmin);
+ umax = std::max(before[i].x(),umax);
+ vmax = std::max(before[i].y(),vmax);
+ }
+
std::vector<SPoint3> pafter(v.size()), pbefore(v.size());
for (unsigned int i=0; i<v.size();++i){
+ if (after[i].x() > umax)return false;
+ if (after[i].x() < umin)return false;
+ if (after[i].y() > vmax)return false;
+ if (after[i].y() < vmin)return false;
GPoint gp = gf->point(after[i]);
if (!gp.succeeded())return false;
pafter[i] = SPoint3 (gp.x(),gp.y(),gp.z());
@@ -673,23 +712,27 @@ static bool _isItAGoodIdeaToMoveThoseVertices (GFace *gf,
surface_old += surfaceFaceUV(e1[j],gf);
for (unsigned int i=0; i<v.size();++i){
- v[i]->setParameter(0,after[i].x());
- v[i]->setParameter(1,after[i].y());
- v[i]->setXYZ(pafter[i].x(),pafter[i].y(),pafter[i].z());
+ if (v[i]->onWhat()->dim() == 2){
+ v[i]->setParameter(0,after[i].x());
+ v[i]->setParameter(1,after[i].y());
+ v[i]->setXYZ(pafter[i].x(),pafter[i].y(),pafter[i].z());
+ }
}
- bool badQuality = false;
for (unsigned int j=0;j<e1.size();++j){
surface_new += surfaceFaceUV(e1[j],gf);
}
for (unsigned int i=0; i<v.size();++i){
- v[i]->setParameter(0,before[i].x());
- v[i]->setParameter(1,before[i].y());
- v[i]->setXYZ(pbefore[i].x(),pbefore[i].y(),pbefore[i].z());
+ if (v[i]->onWhat()->dim() == 2){
+ v[i]->setParameter(0,before[i].x());
+ v[i]->setParameter(1,before[i].y());
+ v[i]->setXYZ(pbefore[i].x(),pbefore[i].y(),pbefore[i].z());
+ }
}
- if (badQuality || (surface_new - surface_old) > 1.e-10 * surface_old) {
+ // if (fabs(surface_new-surface_old) > 1.e-10 * surface) {
+ if (surface_new > 1.0001*surface_old) {
return false;
}
return true;
@@ -721,26 +764,6 @@ static int _quadWithOneVertexOnBoundary (GFace *gf,
touched.insert(v2);
touched.insert(v3);
touched.insert(v4);
- /*
- std::vector<MVertex*> line;
- for (int j=0;j<e1.size();j++){
- for (int k=0;k<e1[j]->getNumVertices();k++){
- if (e1[j]->getVertex(k) != v1 && e1[j] != q &&
- e1[j]->getVertex(k)->onWhat()->dim() < 2){
- line.push_back(e1[j]->getVertex(k));
- }
- }
- }
- // do not collapse if it's a corner
- if (line.size() == 2){
- if (fabs(angle3Vertices(line[0],v1,line[1]) - M_PI) > 3*M_PI/8.){
- // printf("coucou %g\n",angle3Vertices(line[0],v1,line[1])*180./M_PI);
- return 0;
- }
- }
- */
- // if (line.size() == 2)printf("caca\n");
- // else printf("hohcozbucof\n");
for (unsigned int j=0;j<e2.size();++j){
for (int k=0;k<e2[j]->getNumVertices();k++){
if (e2[j]->getVertex(k) == v2 && e2[j] != q)
@@ -1549,7 +1572,7 @@ static int _defectsRemovalBunin(GFace *gf, int maxCavitySize)
// if a vertex (on boundary) is adjacent to one only quad
// and if that quad is badly shaped, we split this
// quad
-static int _splitFlatQuads(GFace *gf, double minQual)
+static int _splitFlatQuads(GFace *gf, double minQual, std::set<MEdge,Less_Edge> &prioritory)
{
v2t_cont adj;
buildVertexToElement(gf->triangles,adj);
@@ -1576,6 +1599,8 @@ static int _splitFlatQuads(GFace *gf, double minQual)
MVertex *v3 = e->getVertex((k+1)%4);
MVertex *v2 = e->getVertex((k+2)%4);
MVertex *v4 = e->getVertex((k+3)%4);
+ prioritory.insert(MEdge(v2,v3));
+ prioritory.insert(MEdge(v3,v4));
SPoint2 pv1,pv2,pv3,pv4,pb1,pb2,pb3,pb4;
reparamMeshEdgeOnFace (v1,v3,gf,pv1,pv3);
reparamMeshEdgeOnFace (v1,v4,gf,pv1,pv4);
@@ -1778,11 +1803,11 @@ int removeDiamonds(GFace *gf)
return nbRemove;
}
-int splitFlatQuads(GFace *gf, double q)
+int splitFlatQuads(GFace *gf, double q, std::set<MEdge,Less_Edge> &prioritory)
{
int nbRemove = 0;
while(1){
- int x = _splitFlatQuads(gf,q);
+ int x = _splitFlatQuads(gf,q, prioritory);
if (!x)break;
nbRemove += x;
}
@@ -1853,7 +1878,7 @@ struct opti_data_vertex_relocation {
v[2] = _v3;
v[3] = _v4;
for (int i=0;i<4;i++){
- movable[i] = v[0]->onWhat()->dim() == 2;
+ movable[i] = v[i]->onWhat()->dim() == 2;
}
}
@@ -1872,21 +1897,14 @@ struct opti_data_vertex_relocation {
double f() const
{
- /*
- double minq = 1.0;
- for (unsigned int i=0;i<e.size();++i){
- MElement *el = e[i];
- minq = std::min(el->gammaShapeMeasure(), minq);
- }
- if (minq < 0)minq *= 1.1;
- else minq *= 0.;
- */
- double val = 0.0;
+
+ double val = 1.0;
for (unsigned int i=0;i<e.size();++i){
MElement *el = e[i];
// double m = log((el->gammaShapeMeasure()-minq)/(1.-minq));
// val += m*m;//el->gammaShapeMeasure();
- val += el->gammaShapeMeasure();
+ // val =std::min(el->gammaShapeMeasure(),val);
+ val -= el->gammaShapeMeasure();
}
return val;
@@ -1899,14 +1917,10 @@ struct opti_data_vertex_relocation {
for (int j=0;j<nv;j++)dfdu[j] = dfdv[j] = 0;
return;
}
- // printf("coucou2 %22.15E, %22.15E\n", U[i]+eps,V[i]);
const double FU = f();
set_(U[i],V[i]+eps,i);
- // printf("coucou3\n");
const double FV = f();
set_(U[i],V[i],i);
- // printf("coucou4\n");
- // printf("%12.5E %12.5E %12.5E\n",F,FU,FV);
dfdu[i] = movable[i] ? -(F-FU)/eps : 0 ;
dfdv[i] = movable[i] ? -(F-FV)/eps : 0;
}
@@ -1915,6 +1929,7 @@ struct opti_data_vertex_relocation {
{
// if (!movable[iv])return;
// printf("getting point of surface %d (%22.15E,%22.15E)\n",gf->tag(),U,V);
+ if (!movable[iv])return true;
GPoint gp = gf->point(U,V);
if (!gp.succeeded()) return false;//Msg::Error("point not OK");
v[iv]->x() = gp.x();
@@ -1952,22 +1967,25 @@ void bfgs_callback_vertex_relocation (const alglib::real_1d_array& x,
static int _untangleQuad (GFace *gf, MQuadrangle *q,v2t_cont & adj)
{
std::set<MElement*> all;
+
+ int numU=0;
for (int i=0;i<4;i++){
std::vector<MElement*> &adji = adj[q->getVertex(i)];
all.insert(adji.begin(),adji.end());
// FIXME
- if (q->getVertex(i)->onWhat()->dim() != 2)return 0;
+ if (q->getVertex(i)->onWhat()->dim() == 2)numU ++;
}
+ if (numU == 0)return 0;
std::vector<MElement*> lt;
lt.insert(lt.begin(),all.begin(),all.end());
double minq_old = 100.;
for (unsigned int i = 0; i < lt.size(); i++){
- const double q = lt[i]->gammaShapeMeasure();
+ const double q = lt[i]->etaShapeMeasure();
minq_old = std::min(q,minq_old);
}
// printf("-------x--------x------------x-------------------\n");
- // printf ("optimizing quadrangle (minq %12.5E) with BFGS %3d\n",minq_old,OPTI_NUMBER);
+ // printf ("optimizing quadrangle (minq %12.5E) with BFGS\n",minq_old);
alglib::ae_int_t dim = 8;
@@ -1983,8 +2001,10 @@ static int _untangleQuad (GFace *gf, MQuadrangle *q,v2t_cont & adj)
double U[4],V[4];
for (int i=0;i<4;i++){
- q->getVertex(i)->getParameter(0,U[i]);
- q->getVertex(i)->getParameter(1,V[i]);
+ SPoint2 p;
+ reparamMeshVertexOnFace(q->getVertex(i),gf, p);
+ U[i] = p.x();
+ V[i] = p.y();
initial_conditions[2*i] = U[i];
initial_conditions[2*i+1] = V[i];
}
@@ -2007,25 +2027,44 @@ static int _untangleQuad (GFace *gf, MQuadrangle *q,v2t_cont & adj)
// printf("minq = %12.5E\n",minq);
std::vector<SPoint2> before;for(int i=0;i<4;i++)before.push_back(SPoint2(U[i],V[i]));
- std::vector<SPoint2> after;for(int i=0;i<4;i++)after.push_back(SPoint2(x[2*i],x[2*i+1]));
+ std::vector<SPoint2> after;
+ for(int i=0;i<4;i++)
+ if (q->getVertex(i)->onWhat()->dim() == 2)
+ after.push_back(SPoint2(x[2*i],x[2*i+1]));
+ else
+ after.push_back(SPoint2(U[i],V[i]));
std::vector<MVertex*> vs;for(int i=0;i<4;i++)vs.push_back(q->getVertex(i));
bool success = _isItAGoodIdeaToMoveThoseVertices (gf,lt,vs,before,after);
+
+ if (success){
+ for (int i=0;i<4;i++)data.set_(x[2*i],x[2*i+1],i);
+ // sprintf(NNN,"UNTANGLE_cavity_%d_after.pos",OPTI_NUMBER++);
+ // data.print_cavity (NNN);
+ return 1;
+ }
+
+ // double minq_new = 100.;
+ // for (unsigned int i = 0; i < data.e.size(); i++){
+ // const double q = data.e[i]->gammaShapeMeasure();
+ // minq_new = std::min(q,minq_new);
+ // }
+
+ // if (success) printf("YES %g %g\n",minq,minq_new);
- if (!success || minq < 0.0 || minq <= minq_old/2) for (int i=0;i<4;i++)data.set_(U[i],V[i],i);
- else return 1;
+ for (int i=0;i<4;i++)data.set_(U[i],V[i],i);
return 0;
+ // return 0;
// else printf("OKIDOKI-UNTANGLE\n");
- // sprintf(NNN,"UNTANGLE_cavity_%d_after.pos",OPTI_NUMBER++);
- // data.print_cavity (NNN);
+ return 1;
}
// use optimization
-static void _relocateVertexOpti(GFace *gf, MVertex *ver,
- const std::vector<MElement*> <)
+static int _relocateVertexOpti(GFace *gf, MVertex *ver,
+ const std::vector<MElement*> <)
{
// return;
- if(ver->onWhat()->dim() != 2)return;
+ if(ver->onWhat()->dim() != 2)return 0;
double minq_old = 100.;
// printf("------------------------------------------------\n");
@@ -2075,15 +2114,12 @@ static void _relocateVertexOpti(GFace *gf, MVertex *ver,
bool success = _isItAGoodIdeaToMoveThatVertex (gf, lt, ver,SPoint2(U,V),SPoint2(x[0],x[1]));
if (!success || minq < 0 || minq <= minq_old/2) data.set_(U,V,0);
- else {
-//printf("OKIDOKI\n");
- }
// if (rep.terminationtype != 4){
// data.set_(U,V);
// }
// sprintf(NNN,"cavity_%d_after.pos",OPTI_NUMBER++);
// data.print_cavity (NNN);
-
+ return success;
}
#endif
@@ -2146,11 +2182,6 @@ void _relocateVertex(GFace *gf, MVertex *ver,
ver->z() = pt.z();
}
}
- else{
-#if defined(HAVE_BFGS)
- // _relocateVertexOpti(gf, ver, lt);
-#endif
- }
}
}
@@ -2178,6 +2209,30 @@ void laplaceSmoothing(GFace *gf, int niter, bool infinity_norm)
}
}
+int optiSmoothing(GFace *gf, int niter, bool infinity_norm)
+{
+ v2t_cont adj;
+ buildVertexToElement(gf->triangles, adj);
+ buildVertexToElement(gf->quadrangles, adj);
+ int N = 0;
+ for(int i = 0; i < niter; i++){
+ v2t_cont::iterator it = adj.begin();
+ while (it != adj.end()){
+ bool doit = false;
+ for (unsigned int j=0;j<it->second.size();j++){
+ if (it->second[j]->gammaShapeMeasure() < .05){
+ doit = true;
+ }
+ }
+ if (doit){
+ N += _relocateVertexOpti(gf, it->first, it->second);
+ }
+ ++it;
+ }
+ }
+ return N;
+}
+
int untangleInvalidQuads(GFace *gf, int niter)
@@ -2190,7 +2245,7 @@ int untangleInvalidQuads(GFace *gf, int niter)
buildVertexToElement(gf->quadrangles, adj);
for(int i = 0; i < niter; i++){
for (unsigned int j=0;j<gf->quadrangles.size();j++){
- if (gf->quadrangles[j]->gammaShapeMeasure() < 0.1){
+ if (gf->quadrangles[j]->etaShapeMeasure() < 0.1){
N += _untangleQuad (gf, gf->quadrangles[j], adj);
}
}
@@ -2199,8 +2254,34 @@ int untangleInvalidQuads(GFace *gf, int niter)
return N;
}
+static int orientationOK (GFace *gf, MVertex *v1, MVertex *v2, MVertex *v3){
+ SPoint2 p1,p2,p3;
+ reparamMeshVertexOnFace(v1,gf, p1);
+ reparamMeshVertexOnFace(v2,gf, p2);
+ reparamMeshVertexOnFace(v3,gf, p3);
+ if (robustPredicates::orient2d (p1,p2,p3) < 0)return true;
+ return false;
+}
-int _edgeSwapQuadsForBetterQuality(GFace *gf)
+static int allowSwap (GFace *gf, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4){
+ SPoint2 p1,p2,p3,p4;
+ reparamMeshVertexOnFace(v1,gf, p1);
+ reparamMeshVertexOnFace(v2,gf, p2);
+ reparamMeshVertexOnFace(v3,gf, p3);
+ reparamMeshVertexOnFace(v4,gf, p4);
+ if (robustPredicates::orient2d (p1,p2,p3) *
+ robustPredicates::orient2d (p1,p2,p4) < 0 &&
+ robustPredicates::orient2d (p3,p4,p1) *
+ robustPredicates::orient2d (p3,p4,p2) > 0)return true;
+ return false;
+}
+
+static double myShapeMeasure(MElement *e){
+ return e->etaShapeMeasure();
+}
+
+
+int _edgeSwapQuadsForBetterQuality(GFace *gf, double eps, std::set<MEdge,Less_Edge> &prioritory)
{
e2t_cont adj;
//buildEdgeToElement(gf->triangles, adj);
@@ -2219,9 +2300,12 @@ int _edgeSwapQuadsForBetterQuality(GFace *gf)
MVertex *v2 = it->first.getVertex(1);
MElement *e1 = it->second.first;
MElement *e2 = it->second.second;
- double worst_quality_old = std::min(e1->gammaShapeMeasure(),e2->gammaShapeMeasure());
+ double worst_quality_old = std::min(myShapeMeasure(e1),myShapeMeasure(e2));
+ bool P = prioritory.find(MEdge(v1,v2)) != prioritory.end();
+ if (P) worst_quality_old = 0.0;
- if (worst_quality_old < .3 && (
+ if ((v1->onWhat()->dim() != 2 || v2->onWhat()->dim() != 2) &&
+ worst_quality_old < eps && (
deleted.find(e1) == deleted.end() ||
deleted.find(e2) == deleted.end())){
MVertex *v12=0,*v11=0,*v22=0,*v21=0;
@@ -2248,6 +2332,7 @@ int _edgeSwapQuadsForBetterQuality(GFace *gf)
return 0;
}
+ // edges have to intersect
MQuadrangle *q1A, *q2A, *q1B, *q2B;
if (rot1 > 0.0 && rot2 > 0.0){
q1A = new MQuadrangle (v11,v22,v2,v12);
@@ -2261,40 +2346,49 @@ int _edgeSwapQuadsForBetterQuality(GFace *gf)
q1B = new MQuadrangle (v11,v12,v21,v1);
q2B = new MQuadrangle (v21,v12,v2,v22);
}
- double worst_quality_A = std::min(q1A->gammaShapeMeasure(),q2A->gammaShapeMeasure());
- double worst_quality_B = std::min(q1B->gammaShapeMeasure(),q2B->gammaShapeMeasure());
-
- bool ca1,ca2,cb1,cb2;
- double old_surface = surfaceFaceUV(e1,gf) + surfaceFaceUV(e2,gf) ;
- double new_surface_A = surfaceFaceUV(q1A,gf,&ca1) + surfaceFaceUV(q2A,gf,&ca2) ;
- double new_surface_B = surfaceFaceUV(q1B,gf,&cb1) + surfaceFaceUV(q2B,gf,&cb2) ;
-
- if (!ca1 && !ca2 &&
- fabs(old_surface - new_surface_A) < 1.e-10 * old_surface &&
- 0.8*worst_quality_A > worst_quality_old &&
- 0.8*worst_quality_A > worst_quality_B &&
- SANITY_(gf,q1A,q2A,12121)){
+ double worst_quality_A = std::min(myShapeMeasure(q1A),myShapeMeasure(q2A));
+ double worst_quality_B = std::min(myShapeMeasure(q1B),myShapeMeasure(q2B));
+
+ bool PA = prioritory.find(MEdge(v11,v22)) == prioritory.end();
+ bool PB = prioritory.find(MEdge(v12,v21)) == prioritory.end();
+
+ double old_surface = surfaceFaceUV(e1,gf) + surfaceFaceUV(e2,gf) ;
+ double new_surface_A = surfaceFaceUV(q1A,gf) + surfaceFaceUV(q2A,gf) ;
+ double new_surface_B = surfaceFaceUV(q1B,gf) + surfaceFaceUV(q2B,gf) ;
+
+ bool doA=false, doB=false;
+ if (old_surface > new_surface_A)doA = true;
+ if (old_surface > new_surface_B)doB = true;
+ if ( worst_quality_A > worst_quality_B){
+ if ( worst_quality_A > worst_quality_old && !doB) doA = true;
+ }
+ else{
+ if ( worst_quality_B > worst_quality_old && !doA) doB = true;
+ }
+
+ if(!allowSwap(gf,v1,v2,v11,v22)) doA = false;
+ if(!allowSwap(gf,v1,v2,v21,v12)) doB = false;
+
+ if (!PA)doA = false;
+ if (!PB)doB = false;
+
+
+ if (doA && SANITY_(gf,q1A,q2A)){
deleted.insert(e1);
deleted.insert(e2);
created.push_back(q1A);
created.push_back(q2A);
delete q1B;
delete q2B;
- //printf("edge swap performed -- 1\n");
COUNT++;
}
- else if (!cb1 && !cb2 &&
- fabs(old_surface - new_surface_B) < 1.e-10 * old_surface &&
- 0.8*worst_quality_B > worst_quality_old &&
- 0.8*worst_quality_B > worst_quality_A &&
- SANITY_(gf,q1B,q2B,12121)){
+ else if (doB && SANITY_(gf,q1B,q2B)){
deleted.insert(e1);
deleted.insert(e2);
created.push_back(q1B);
created.push_back(q2B);
delete q1A;
delete q2A;
- //printf("edge swap performed -- 2\n");
COUNT++;
}
else {
@@ -2320,13 +2414,13 @@ int _edgeSwapQuadsForBetterQuality(GFace *gf)
return COUNT;
}
-static int edgeSwapQuadsForBetterQuality (GFace *gf)
+static int edgeSwapQuadsForBetterQuality (GFace *gf, double eps, std::set<MEdge,Less_Edge> &prioritory)
{
int COUNT = 0;
while(1){
- int k = _edgeSwapQuadsForBetterQuality (gf);
+ int k = _edgeSwapQuadsForBetterQuality (gf, eps, prioritory);
COUNT += k;
- if (!k || COUNT > 10)break;
+ if (!k || COUNT > 40)break;
}
Msg::Debug("%i swap quads performed",COUNT);
return COUNT;
@@ -2397,6 +2491,41 @@ static std::vector<MVertex*> computeBoundingPoints (const std::vector<MElement*>
return result;
}
+static MQuadrangle* buildNewQuad(MVertex *first,
+ MVertex *newV,
+ MElement *e,
+ const std::vector<MElement*> & E){
+ int found[3] = {0,0,0};
+ for (unsigned int i=0;i<E.size();i++){
+ if (E[i] != e){
+ for (int j=0;j<E[i]->getNumVertices();j++){
+ MVertex *v = E[i]->getVertex(j);
+ if (v == e->getVertex(0))found[0]=1;
+ else if (v == e->getVertex(1))found[1]=1;
+ else if (v == e->getVertex(2))found[2]=1;
+ }
+ }
+ }
+ int start = -1;
+ int nextVertex=-1;
+ for (int i=0;i<3;i++){
+ if (e->getVertex(i) == first) start = i;
+ if (!found[i])nextVertex = i ;
+ }
+ if (nextVertex == (start+1)%3){
+ return new MQuadrangle (e->getVertex(start),
+ e->getVertex((start+1)%3),
+ e->getVertex((start+2)%3),
+ newV);
+ }
+ else{
+ return new MQuadrangle (e->getVertex(start),
+ newV,
+ e->getVertex((start+1)%3),
+ e->getVertex((start+2)%3));
+ }
+}
+
int postProcessExtraEdges (GFace *gf, std::vector<std::pair<MElement*,MElement*> > &toProcess)
{
// some pairs of elements that should be recombined are not neighbor elements (see our paper)
@@ -2454,47 +2583,17 @@ int postProcessExtraEdges (GFace *gf, std::vector<std::pair<MElement*,MElement*>
p.y());
gf->mesh_vertices.push_back(newVertex);
int start1 = 0,start2 = 0;
- int orientation1=1, orientation2=1;
for (int k=0;k<3;k++){
if (t1->getVertex(k) == it->first){
- // if (t1->getVertex((k+1)%3)->onWhat()->dim
start1 = k;
}
- if (t2->getVertex(k) == it->first)start2 = k;
+ if (t2->getVertex(k) == it->first){
+ start2 = k;
+ }
}
- /*
- FIX THAT !!!!!
- x--------x--------x
- \ t1 / \ t2 /
- \ / \ /
- \ / \ /
- \/ \/
-
- */
- MQuadrangle *q1;
- if (orientation1)
- q1 = new MQuadrangle(t1->getVertex(start1),
- newVertex,
- t1->getVertex((start1+1)%3),
- t1->getVertex((start1+2)%3));
- else
- q1 = new MQuadrangle(newVertex,
- t1->getVertex((start1+2)%3),
- t1->getVertex((start1+1)%3),
- t1->getVertex(start1));
-
- MQuadrangle *q2;
- if (orientation2)
- q2 = new MQuadrangle(t2->getVertex(start2),
- newVertex,
- t2->getVertex((start2+1)%3),
- t2->getVertex((start2+2)%3));
- else
- q2 = new MQuadrangle(newVertex,
- t2->getVertex((start2+2)%3),
- t2->getVertex((start2+1)%3),
- t2->getVertex(start2));
+ MQuadrangle *q1 =buildNewQuad(t1->getVertex(start1),newVertex,t1,it->second);
+ MQuadrangle *q2 =buildNewQuad(t2->getVertex(start2),newVertex,t2,it->second);
std::vector<MElement*> newAdj;
newAdj.push_back(q1);
@@ -2915,7 +3014,7 @@ int recombineWithBlossom(GFace *gf, double dx, double dy,
elist[2*i] = t2n[pairs[i].t1];
elist[2*i+1] = t2n[pairs[i].t2];
//elen [i] = (int) pairs[i].angle;
- elen [i] = (int) pairs[i].total_cost; //addition for class Temporary
+ // elen [i] = (int) pairs[i].total_cost; //addition for class Temporary
double angle = atan2(pairs[i].n1->y()-pairs[i].n2->y(),
pairs[i].n1->x()-pairs[i].n2->x());
@@ -2927,14 +3026,14 @@ int recombineWithBlossom(GFace *gf, double dx, double dy,
if (angle < 0) angle += M_PI/2;
if (angle > M_PI/2) angle -= M_PI/2;
}
- //elen [i] = (int) 180. * fabs(angle-M_PI/4)/M_PI;
+ elen [i] = (int) 180. * fabs(angle-M_PI/4)/M_PI;
int NB = 0;
if (pairs[i].n1->onWhat()->dim() < 2)NB++;
if (pairs[i].n2->onWhat()->dim() < 2)NB++;
if (pairs[i].n3->onWhat()->dim() < 2)NB++;
if (pairs[i].n4->onWhat()->dim() < 2)NB++;
- if (elen[i] > gf->meshAttributes.recombineAngle && NB > 2) {elen[i] = 1000;}
+ if (elen[i] > 180 && NB > 2) {printf("angle %dE\n",elen[i]);elen[i] = 1000;}
}
if (cubicGraph){
@@ -3121,19 +3220,19 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
for (unsigned int i = 0; i < pairs.size(); ++i){
elist[2*i] = t2n[pairs[i].t1];
elist[2*i+1] = t2n[pairs[i].t2];
- //elen [i] = (int) pairs[i].angle;
- elen [i] = (int) pairs[i].total_cost; //addition for class Temporary
- double angle = atan2(pairs[i].n1->y()-pairs[i].n2->y(),
- pairs[i].n1->x()-pairs[i].n2->x());
+ elen [i] = (int) 1000*exp(-pairs[i].angle);
+ //elen [i] = (int) pairs[i].total_cost; //addition for class Temporary
+ // double angle = atan2(pairs[i].n1->y()-pairs[i].n2->y(),
+ // pairs[i].n1->x()-pairs[i].n2->x());
//double x = .5*(pairs[i].n1->x()+pairs[i].n2->x());
//double y = .5*(pairs[i].n1->y()+pairs[i].n2->y());
//double opti = atan2(y,x);
//angle -= opti;
- while (angle < 0 || angle > M_PI/2){
- if (angle < 0) angle += M_PI/2;
- if (angle > M_PI/2) angle -= M_PI/2;
- }
+ // while (angle < 0 || angle > M_PI/2){
+ // if (angle < 0) angle += M_PI/2;
+ // if (angle > M_PI/2) angle -= M_PI/2;
+ // }
//elen [i] = (int) 180. * fabs(angle-M_PI/4)/M_PI;
int NB = 0;
@@ -3141,7 +3240,9 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
if (pairs[i].n2->onWhat()->dim() < 2)NB++;
if (pairs[i].n3->onWhat()->dim() < 2)NB++;
if (pairs[i].n4->onWhat()->dim() < 2)NB++;
- if (elen[i] > gf->meshAttributes.recombineAngle && NB > 2) {elen[i] = 1000;}
+ //if (NB > 2)printf("%d %d\n",elen[i],NB);
+ if (elen[i] > (int)1000*exp(.1) && NB > 2) {elen[i] = 5000;}
+ else if (elen[i] >= 1000 && NB > 2) {elen[i] = 10000;}
}
if (cubicGraph){
@@ -3228,7 +3329,7 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
// FIXME !!
if (an == 100000 /*|| an == 1000*/){
toProcess.push_back(std::make_pair(n2t[i1],n2t[i2]));
- Msg::Debug("Extra edge found in blossom algorithm, optimization will be required");
+ Msg::Warning("Extra edge found in blossom algorithm, optimization will be required");
}
else{
MElement *t1 = n2t[i1];
@@ -3328,6 +3429,22 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
}
+static double printStats(GFace *gf,const char *message){
+ int nbBad=0;
+ int nbInv=0;
+ double Qav=0;
+ double Qmin=1;
+ for (unsigned int i=0;i<gf->quadrangles.size();i++){
+ double Q = gf->quadrangles[i]->etaShapeMeasure() ;
+ if (Q <= 0.0)nbInv ++;
+ if (Q <= 0.1)nbBad ++;
+ Qav += Q;
+ Qmin = std::min(Q,Qmin);
+ }
+ Msg::Info("%s : %5d quads %4d invalid quads %4d quads with Q < 0.1 Avg Q = %12.5E Min %12.5E",message, gf->quadrangles.size(), nbInv, nbBad,Qav/gf->quadrangles.size(),Qmin);
+ return Qmin;
+}
+
void recombineIntoQuads(GFace *gf,
bool topologicalOpti,
bool nodeRepositioning)
@@ -3348,81 +3465,62 @@ void recombineIntoQuads(GFace *gf,
int success = _recombineIntoQuads(gf, 0);
if (saveAll) gf->model()->writeMSH("raw.msh");
- if(haveParam && nodeRepositioning)
+ printStats (gf, "BEFORE OPTIMIZATION");
+ if(haveParam && nodeRepositioning){
laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing);
-
+ }
// blossom-quad algo
if(success && CTX::instance()->mesh.algoRecombine == 1){
if(topologicalOpti){
if(haveParam){
if (saveAll) gf->model()->writeMSH("smoothed.msh");
- int COUNT = 0, ITER=0;
- char NAME[256];
-
- double tFQ=0,tTQ=0,tDI=0,tLA=0,tSW=0,tBU=0;
- double oldoptistatus[5] = {0,0,0,0,0};
- double optistatus[5] = {0,0,0,0,0};
+ int ITER=0;
+ int ITERB = 0;
+ int optistatus[6] = {0,0,0,0,0,0};
+ std::set<MEdge,Less_Edge> prioritory;
while(1){
- double t6 = Cpu();
int maxCavitySize = CTX::instance()->mesh.bunin;
- optistatus[4] = _defectsRemovalBunin(gf,maxCavitySize);
- if(optistatus[4] && saveAll){
- sprintf(NAME,"iter%dB.msh",COUNT++); gf->model()->writeMSH(NAME);
- }
-
- double t7 = Cpu();
- double t1 = Cpu();
- optistatus[0] = splitFlatQuads(gf, .3);
- if(optistatus[0] && saveAll){
- sprintf(NAME,"iter%dSF.msh",COUNT++); gf->model()->writeMSH(NAME);
- }
- double t2 = Cpu();
+ optistatus[0] = (ITERB == 1) ?splitFlatQuads(gf, .01, prioritory) : 0;
optistatus[1] = removeTwoQuadsNodes(gf);
- if(optistatus[1] && saveAll){
- sprintf(NAME,"iter%dTQ.msh",COUNT++); gf->model()->writeMSH(NAME);
- }
- double t3 = Cpu();
- optistatus[2] = removeDiamonds(gf);
- if(optistatus[2] && saveAll){
- sprintf(NAME,"iter%dD.msh",COUNT++); gf->model()->writeMSH(NAME);
- }
-
- double t4 = Cpu();
- untangleInvalidQuads(gf,CTX::instance()->mesh.nbSmoothing);
- laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing,true);
- double t5 = Cpu();
- optistatus[3] = edgeSwapQuadsForBetterQuality(gf);
-
- if(optistatus[3] && saveAll){
- sprintf(NAME,"iter%dS.msh",COUNT++); gf->model()->writeMSH(NAME);
- }
- double t5b = Cpu();
- tFQ += (t2-t1);
- tTQ += (t3-t2);
- tDI += (t4-t3);
- tLA += (t5-t4);
- tSW += (t5b-t5);
- tBU += (t7-t6);
- if (!(optistatus[0]+optistatus[1]+optistatus[2]+optistatus[3]+optistatus[4]))
- break;
- bool theSame = true;
- for (int i=0;i<5;i++){
- if (optistatus[i] != oldoptistatus[i])theSame = false;
- oldoptistatus[i] = optistatus[i];
- }
- if(theSame)break;
- if (ITER++ >= 20)break;
+ optistatus[4] = _defectsRemovalBunin(gf,maxCavitySize);
+ optistatus[2] = removeDiamonds(gf) ;
+ laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing,true);
+ optistatus[3] = edgeSwapQuadsForBetterQuality(gf,.01, prioritory);
+ optistatus[5] = optiSmoothing(gf,CTX::instance()->mesh.nbSmoothing,true);
+ optistatus[5] = (ITERB == 1) ?untangleInvalidQuads(gf,CTX::instance()->mesh.nbSmoothing) : 0;
+
+ double bad = printStats (gf, "IN OPTIMIZATION");
+ if (bad > .1)break;
+ if (ITER == 10){
+ ITERB = 1;
+ }
+ if (ITER > 20)break;
+ int nb = 0;
+ for (int i=0;i<6;i++)nb += optistatus[i];
+ if (!nb && ITERB == 0)ITERB = 1;
+ else if (!nb )break;
+ // printf("%d %d %d %d %d %d\n",optistatus[0],optistatus[1],optistatus[2],optistatus[3],optistatus[4],optistatus[5]);
+ ITER ++;
}
- Msg::Debug("Opti times FQ(%4.3f) tQ(%4.3f) DI(%4.3f) LA(%4.3f) "
- "SW(%4.3f) BUN(%4.3f)",tFQ,tTQ,tDI,tLA,tSW,tBU);
}
- if(haveParam) laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing, true);
+ if(haveParam) {
+ laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing, true);
+ optiSmoothing(gf,CTX::instance()->mesh.nbSmoothing,true);
+ }
+ // untangleInvalidQuads(gf,CTX::instance()->mesh.nbSmoothing);
}
double t2 = Cpu();
Msg::Info("Blossom recombination algorithm completed (%g s)", t2 - t1);
- quadsToTriangles(gf, 0.01);
- // if(haveParam) laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
+ quadsToTriangles(gf, .01);
+ if(haveParam) {
+ laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
+ // optiSmoothing(gf,CTX::instance()->mesh.nbSmoothing,true);
+ }
+ // removeDiamonds(gf);
+ // removeTwoQuadsNodes(gf);
+ printStats (gf, "AFTER OPTIMIZATION");
+
return;
}
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index 010607dddb2856a1e029bf057f3ed8e0cde4c89a..04987bff9447d1f0374874cd8d26f1fbdf1efd2b 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -158,6 +158,10 @@ struct RecombineTriangle
else if(t2->getVertex(1) != n1 && t2->getVertex(1) != n2) n4 = t2->getVertex(1);
else if(t2->getVertex(2) != n1 && t2->getVertex(2) != n2) n4 = t2->getVertex(2);
+ MQuadrangle q (n1,n3,n2,n4);
+ angle = q.etaShapeMeasure();
+
+ /*
double a1 = 180 * angle3Vertices(n1, n4, n2) / M_PI;
double a2 = 180 * angle3Vertices(n4, n2, n3) / M_PI;
double a3 = 180 * angle3Vertices(n2, n3, n1) / M_PI;
@@ -166,6 +170,7 @@ struct RecombineTriangle
angle = std::max(fabs(90. - a2),angle);
angle = std::max(fabs(90. - a3),angle);
angle = std::max(fabs(90. - a4),angle);
+ */
cost_quality = 1.0 - std::max(1.0 - angle/90.0,0.0); //addition for class Temporary
cost_alignment = Temporary::compute_alignment(me,_t1,_t2); //addition for class Temporary
total_cost = Temporary::compute_total_cost(cost_quality,cost_alignment); //addition for class Temporary
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 603923906ceb168c3bf0b365afd702b27e2bf4d6..2b33eef5629c330d1c48b6e336af17e77ebf3129 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -305,6 +305,8 @@ bool insertVertexB(std::list<faceXtet> &shell,
++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 (fabs(oldVolume - newVolume) < 1.e-10 * oldVolume &&
!onePointIsTooClose){
connectTets_vector(new_cavity.begin(), new_cavity.end());
diff --git a/Mesh/simple3D.cpp b/Mesh/simple3D.cpp
index 6a764903c08aa652c5c2cfc05c440c0f3740f0db..485c2c3ea6f581187c0a719963d0604f1937d2f3 100644
--- a/Mesh/simple3D.cpp
+++ b/Mesh/simple3D.cpp
@@ -345,7 +345,7 @@ void Filler::treat_region(GRegion* gr){
std::set<MVertex*> temp;
std::set<MVertex*>::iterator it;
RTree<Node*,double,3,double> rtree;
-
+
Frame_field::init_region(gr);
Size_field::init_region(gr);
Size_field::solve(gr);
@@ -704,4 +704,4 @@ void Filler::print_node(Node* node,std::ofstream& file){
/*********static declarations*********/
-std::vector<MVertex*> Filler::new_vertices;
\ No newline at end of file
+std::vector<MVertex*> Filler::new_vertices;
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index f20fcf562d49bd9d83fdfe96b9e1ab18503ff028..4bd812b21ade39184b7419df2bbc0bc92e789597 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -18,7 +18,7 @@
#include "BackgroundMesh.h"
#include "intersectCurveSurface.h"
-static const double FACTOR = .81;
+static const double FACTOR = .71;
static const int NUMDIR = 3;
//static const double DIRS [NUMDIR] = {0.0};
static const double DIRS [NUMDIR] = {0.0, M_PI/20.,-M_PI/20.};
@@ -133,7 +133,7 @@ bool inExclusionZone (SPoint2 &p,
if (!backgroundMesh::current()->inDomain(p.x(),p.y(),0)) return true;
my_wrapper w (p);
- double _min[2] = {p.x()-1.e-8, p.y()-1.e-8},_max[2] = {p.x()+1.e-8,p.y()+1.e-8};
+ double _min[2] = {p.x()-1.e-3, p.y()-1.e-3},_max[2] = {p.x()+1.e-3,p.y()+1.e-3};
rtree.Search(_min,_max,rtree_callback,&w);
return w._tooclose;
@@ -259,14 +259,15 @@ bool compute4neighbors (GFace *gf, // the surface
,covar1[0],covar1[1],covar2[0],covar2[1],l1,l2,size_1,size_2,size_param_1,size_param_2,M,N,E,s1.x(),s1.y(),s2.x(),s2.y());*/
// this is the rectangle in the parameter plane.
- double r1 = 0*1.e-8*(double)rand() / RAND_MAX;
- double r2 = 0*1.e-8*(double)rand() / RAND_MAX;
- double r3 = 0*1.e-8*(double)rand() / RAND_MAX;
- double r4 = 0*1.e-8*(double)rand() / RAND_MAX;
- double r5 = 0*1.e-8*(double)rand() / RAND_MAX;
- double r6 = 0*1.e-8*(double)rand() / RAND_MAX;
- double r7 = 0*1.e-8* (double)rand() / RAND_MAX;
- double r8 = 0*1.e-8*(double)rand() / RAND_MAX;
+ const double EPS = 1.e-6;
+ double r1 = EPS*(double)rand() / RAND_MAX;
+ double r2 = EPS*(double)rand() / RAND_MAX;
+ double r3 = EPS*(double)rand() / RAND_MAX;
+ double r4 = EPS*(double)rand() / RAND_MAX;
+ double r5 = EPS*(double)rand() / RAND_MAX;
+ double r6 = EPS*(double)rand() / RAND_MAX;
+ double r7 = EPS* (double)rand() / RAND_MAX;
+ double r8 = EPS*(double)rand() / RAND_MAX;
double newPoint[4][2] = {{midpoint[0] - covar1[0] * size_param_1 +r1,
midpoint[1] - covar1[1] * size_param_1 +r2},
{midpoint[0] - covar2[0] * size_param_2 +r3,
diff --git a/benchmarks/hex/cube.geo b/benchmarks/hex/cube.geo
index b10cf147f1f6f4c6763913abb91f0f0dcee246aa..5232659f840bedfba0ac988be755ee59919915e6 100644
--- a/benchmarks/hex/cube.geo
+++ b/benchmarks/hex/cube.geo
@@ -4,15 +4,15 @@ Mesh.Recombine3DAll = 1;
Mesh.Smoothing = 0;
c1 = 0.25; //125;
-c2 = 0.25;
-Point(1) = {1, 1, -1, c1};
+c2 = 0.05;
+Point(1) = {1, 1, -1, c2};
Point(2) = {-1, 1, -1, c1};
Point(3) = {-1, -1, -1, c1};
Point(4) = {1, -1, -1, c1};
-Point(5) = {1, 1, 1, c2};
-Point(6) = {-1, 1, 1, c2};
-Point(7) = {1, -1, 1, c2};
-Point(8) = {-1, -1, 1, c2};
+Point(5) = {1, 1, 1, c1};
+Point(6) = {-1, 1, 1, c1};
+Point(7) = {1, -1, 1, c1};
+Point(8) = {-1, -1, 1, c1};
Line(1) = {7, 5};
Line(2) = {5, 6};
Line(3) = {6, 8};
diff --git a/utils/converters/autocad/dxf2geo.cpp b/utils/converters/autocad/dxf2geo.cpp
index 5794f271eb675faca4fdb8715dd3366ba7861739..a7af2d48728aba43889f8cdae5c22ddbc3cc8727 100644
--- a/utils/converters/autocad/dxf2geo.cpp
+++ b/utils/converters/autocad/dxf2geo.cpp
@@ -292,8 +292,10 @@ void addobj(void)
else if(strstr(curobj, "TEXT")) { /* not implemented for now */
}
else if(strstr(curobj, "SHAPE")) { /* these look very hard */
+ printf("SHAPE\n");
}
else if(strstr(curobj, "BLOCK")) { /* these look very hard */
+ printf("BLOCK\n");
}
else if(strstr(curobj, "ENDBLK")) { /* these look very hard */
}
@@ -304,6 +306,7 @@ void addobj(void)
else if(strstr(curobj, "ATTRIB")) { /* not implemented for now */
}
else if(strstr(curobj, "POLYLINE")) { /* these look fairly hard */
+ printf("POLYLINE\n");
}
else if(strstr(curobj, "VERTEX")) { /* these look fairly hard */
if(ints[0] == 192) {
@@ -407,6 +410,9 @@ void addobj(void)
}
else if(strstr(curobj, "DIMENSION")) { /* not implemented for now */
}
+ else {
+ printf("%s\n",curobj);
+ }
}