diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 51d3e37e14ca1bb8327a5160fc480122647e9708..68048b5dcd00d948a0574865e02b90d11e38f3a7 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1198,7 +1198,6 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
GPoint gp = gf->point(U * scalingU, V * scalingV);
if (!gp.succeeded()){
- // printf ("iha\n");
return false;
}
const double oldX = p->X;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 061f171330302aed8e94d08426c68afb4de24b38..a5a8bc0df132c820e392f5b28fc8aa6e603eac0c 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -77,16 +77,16 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
double Crv = 0;
switch(ge->dim()){
case 0:
- Crv = max_edge_curvature((const GVertex *)ge);
- Crv = std::max(max_surf_curvature((const GVertex *)ge), Crv);
- //Crv = max_surf_curvature((const GVertex *)ge);
+ // Crv = max_edge_curvature((const GVertex *)ge);
+ // Crv = std::max(max_surf_curvature((const GVertex *)ge), Crv);
+ Crv = max_surf_curvature((const GVertex *)ge);
break;
case 1:
{
GEdge *ged = (GEdge *)ge;
- Crv = ged->curvature(U);
- Crv = std::max(Crv, max_surf_curvature(ged, U));
- //Crv = max_surf_curvature(ged, U);
+ // Crv = ged->curvature(U)*2;
+ // Crv = std::max(Crv, max_surf_curvature(ged, U));
+ Crv = max_surf_curvature(ged, U);
}
break;
case 2:
diff --git a/Mesh/DivideAndConquer.cpp b/Mesh/DivideAndConquer.cpp
index a7f527ea7aba92faf65d228f222553753e56b4d1..e3e1f53a023c7e4a198dc15a2fdfdda5756fd761 100644
--- a/Mesh/DivideAndConquer.cpp
+++ b/Mesh/DivideAndConquer.cpp
@@ -496,14 +496,15 @@ static int ConvertDListToTriangles(DocRecord *doc)
{
// on suppose que n >= 3. gPointArray est suppose OK.
- STriangle *striangle;
+ // STriangle *striangle;
int n, i, j;
int count = 0, count2 = 0;
PointNumero aa, bb, cc;
PointRecord *gPointArray = doc->points;
n = doc->numPoints;
- striangle = (STriangle *) Malloc(n * sizeof(STriangle));
+ STriangle *striangle = (STriangle *) Malloc(n * sizeof(STriangle));
+
count2 = CountPointsOnHull(n, doc->points);
// nombre de triangles que l'on doit obtenir
@@ -516,8 +517,6 @@ static int ConvertDListToTriangles(DocRecord *doc)
// de points (t_length)
striangle[i].t = ConvertDlistToArray(&gPointArray[i].adjacent,
&striangle[i].t_length);
- striangle[i].info = NULL;
- striangle[i].info_length = 0;
}
// on balaye les noeuds de gauche a droite -> on cree les triangles
@@ -536,11 +535,8 @@ static int ConvertDListToTriangles(DocRecord *doc)
}
}
}
- for(i = 0; i < n; i++)
- Free(striangle[i].t);
- Free(striangle);
doc->numTriangles = count2;
-
+ doc->striangle = striangle;
return 1;
}
@@ -563,7 +559,7 @@ static void RemoveAllDList(int n, PointRecord *pPointArray)
}
DocRecord::DocRecord(int n)
- : numPoints(n), points(NULL), numTriangles(0), triangles(NULL)
+ : numPoints(n), points(NULL), striangle(0),numTriangles(0), triangles(NULL)
{
if(numPoints)
points = (PointRecord*)Malloc(numPoints * sizeof(PointRecord));
@@ -573,6 +569,11 @@ DocRecord::~DocRecord()
{
Free(points);
Free(triangles);
+ if (striangle){
+ for(int i = 0; i < numPoints; i++)
+ Free(striangle[i].t);
+ Free(striangle);
+ }
}
void DocRecord::MakeMeshWithPoints()
@@ -581,4 +582,17 @@ void DocRecord::MakeMeshWithPoints()
BuildDelaunay(this);
ConvertDListToTriangles(this);
RemoveAllDList(numPoints, points);
+ for(int i = 0; i < numPoints; i++)
+ Free(striangle[i].t);
+ Free(striangle);
+ striangle = 0;
}
+
+void DocRecord::MakeVoronoi(){
+ if(numPoints < 3) return;
+ BuildDelaunay(this);
+ // TagPointsOnHull(this);
+ ConvertDListToTriangles(this);
+ RemoveAllDList(numPoints, points);
+}
+
diff --git a/Mesh/DivideAndConquer.h b/Mesh/DivideAndConquer.h
index 7fdd85ac5cabfc8f5381bfc8c6453d2581392c9e..d71705bc02a4d788c0a5c727982857a18b108e0a 100644
--- a/Mesh/DivideAndConquer.h
+++ b/Mesh/DivideAndConquer.h
@@ -26,14 +26,8 @@ struct _CDLIST{
};
typedef struct{
- PointNumero search;
- PointNumero already;
-}HalfTriangle;
-
-typedef struct{
- HalfTriangle *info;
PointNumero *t;
- int t_length, info_length;
+ int t_length;
}STriangle;
typedef struct{
@@ -50,15 +44,22 @@ typedef struct{
PointNumero a, b, c;
}Triangle;
+
class DocRecord{
public:
int numPoints; // number of points
PointRecord *points; // points to triangulate
int numTriangles; // number of triangles
Triangle *triangles; // 2D results
+ STriangle *striangle; // connected points
DocRecord(int n);
+ double & x(int i)
+ {return points[i].where.v;}
+ double & y(int i)
+ {return points[i].where.h;}
~DocRecord();
void MakeMeshWithPoints();
+ void MakeVoronoi();
};
#endif
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 793453aa079ffab397ab82b7d48f6a6340570d4a..5d204dc4895191fea95781a17009fad8d7f11bae 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -149,6 +149,7 @@ gmshSmoothHighOrder${OBJEXT}: gmshSmoothHighOrder.cpp ../Geo/MLine.h \
../Numeric/GmshMatrix.h ../Numeric/gmshElasticity.h \
../Numeric/gmshTermOfFormulation.h ../Numeric/GmshMatrix.h \
../Numeric/gmshLinearSystemGmm.h ../Numeric/gmshLinearSystem.h \
+ ../Numeric/gmshLinearSystemCSR.h ../Numeric/gmshLinearSystem.h \
../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
meshGEdge${OBJEXT}: meshGEdge.cpp ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp
index 1253033e0eb52b93b08d6388007c16fcf2c8f370..2365f7df61cc6d661a89b670e65231d9fcc35f28 100644
--- a/Mesh/gmshSmoothHighOrder.cpp
+++ b/Mesh/gmshSmoothHighOrder.cpp
@@ -21,6 +21,7 @@
#include "gmshLaplace.h"
#include "gmshElasticity.h"
#include "gmshLinearSystemGmm.h"
+#include "gmshLinearSystemCSR.h"
#include "GFace.h"
#include "GRegion.h"
#include "Context.h"
@@ -428,9 +429,9 @@ void gmshHighOrderSmoother::computeMetricVector (GFace *gf,
void gmshHighOrderSmoother::smooth_metric ( std::vector<MElement*> & all, GFace *gf)
{
gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
- // lsys->setNoisy(1);
+ lsys->setNoisy(1);
lsys->setGmres(1);
- lsys->setPrec(1.e-7);
+ lsys->setPrec(1.e-9);
gmshAssembler<double> myAssembler(lsys);
gmshElasticityTerm El(0, 1.0, .333, getTag());
@@ -516,7 +517,7 @@ void gmshHighOrderSmoother::smooth_metric ( std::vector<MElement*> & all, GFace
double dx2 = smooth_metric_ ( v,gf, myAssembler, verticesToMove,El);
printf(" dx2 = %12.5E\n",dx2);
if (fabs(dx2-dx) < 1.e-4 * dx0)break;
- if (iter++ > 10)break;
+ if (iter++ > 2)break;
dx = dx2;
}
@@ -615,10 +616,10 @@ double gmshHighOrderSmoother::smooth_metric_ ( std::vector<MElement*> & v,
void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
{
- gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
+ gmshLinearSystemCSRGmm<double> *lsys = new gmshLinearSystemCSRGmm<double>;
lsys->setNoisy(1);
lsys->setGmres(1);
- lsys->setPrec(1.e-5);
+ lsys->setPrec(5.e-8);
gmshAssembler<double> myAssembler(lsys);
gmshElasticityTerm El(0, 1.0, .333, getTag());
@@ -632,7 +633,7 @@ void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
if (!v.size()) return;
- const int nbLayers = 2;
+ const int nbLayers = 1;
for (int i=0;i<nbLayers;i++){
addOneLayer ( all, v, layer);
v.insert(v.end(),layer.begin(),layer.end());
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 9391abe4e1d06eeb56bfeaaaa5ce6afe4488ed11..5b681f2ae0be998653197e8f9383d59581e5ad12 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -633,7 +633,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
if(RECUR_ITER > 0)
Msg::Warning(":-) Gmsh was able to recover all edges after %d iterations", RECUR_ITER);
- Msg::Info("Boundary Edges recovered for surface %d", gf->tag());
+ Msg::Debug("Boundary Edges recovered for surface %d", gf->tag());
// Look for an edge that is on the boundary for which one of the two
// neighbors has a negative number node. The other triangle is
@@ -1077,8 +1077,8 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
// Buid a BDS_Mesh structure that is convenient for doing the actual
// meshing procedure
BDS_Mesh *m = new BDS_Mesh;
- m->scalingU = fabs(du);
- m->scalingV = fabs(dv);
+ m->scalingU = 1;//fabs(du);
+ m->scalingV = 1;//fabs(dv);
std::vector<std::vector<BDS_Point*> > edgeLoops_BDS;
SBoundingBox3d bbox;
int nbPointsTotal = 0;
@@ -1449,8 +1449,9 @@ void orientMeshGFace::operator()(GFace *gf)
{
gf->model()->setCurrentMeshEntity(gf);
+ return;
if(gf->geomType() == GEntity::ProjectionFace) return;
-
+ if(gf->geomType() == GEntity::CompoundSurface)return;
// in old versions we did not reorient transfinite surface meshes;
// we could add the following to provide backward compatibility:
// if(gf->meshAttributes.Method == MESH_TRANSFINITE) return;
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index b9d0a36bb67eb67e19870185204f4fbbbfae25fd..e2b293f75565d8b26a7a507ffe7fec61a28adcd8 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -100,20 +100,41 @@ double NewGetLc(BDS_Point *p)
std::min(p->lc(), p->lcBGM()) : p->lcBGM();
}
+static double correctLC_ (BDS_Point *p1,BDS_Point *p2, GFace *f,
+ double SCALINGU, double SCALINGV){
+ double l1 = NewGetLc(p1);
+ double l2 = NewGetLc(p2);
+ double l = 0.5*(l1+l2);
+
+ // printf(" %g %g -- %g %g\n",SCALINGU,SCALINGV,l1,l2);
+
+ if(CTX::instance()->mesh.lcFromCurvature)
+ {
+ // GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
+ // 0.5 * (p1->v + p2->v) * SCALINGV));
+ // double l3 = BGM_MeshSize(f,GP.u(),GP.v(),GP.x(),GP.y(),GP.z());
+ double l3=l2;
+ double lcmin = std::min(std::min(l1,l2),l3);
+ l1 = std::min(lcmin*1.2,l1);
+ l2 = std::min(lcmin*1.2,l2);
+ l3 = std::min(lcmin*1.2,l3);
+ l = (l1+l2+l3)/3.0;
+ }
+ return l;
+}
+
double NewGetLc(BDS_Edge *e, GFace *f, double SCALINGU, double SCALINGV)
{
double linearLength = computeEdgeLinearLength(e, f, SCALINGU, SCALINGV);
- double l1 = NewGetLc(e->p1);
- double l2 = NewGetLc(e->p2);
- return 2 * linearLength / (l1 + l2);
+ double l = correctLC_ (e->p1,e->p2,f, SCALINGU, SCALINGV);
+ return linearLength / l;
}
double NewGetLc(BDS_Point *p1,BDS_Point *p2, GFace *f, double su, double sv)
{
double linearLength = computeEdgeLinearLength(p1, p2, f, su, sv);
- double l1 = NewGetLc(p1);
- double l2 = NewGetLc(p2);
- return 2 * linearLength / (l1 + l2);
+ double l = correctLC_ (p1,p2,f, su, sv);
+ return linearLength / l;
}
void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
@@ -415,8 +436,8 @@ void splitEdgePassUnsorted(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
m.scalingU, m.scalingV);
BDS_Point *mid;
- GPoint gpp = gf->point(coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
- coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v);
+ GPoint gpp = gf->point(m.scalingU*(coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u),
+ m.scalingV*(coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v));
if (gpp.succeeded()){
mid = m.add_point(++m.MAXPOINTNUMBER,
coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
@@ -461,7 +482,7 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
double U = coord * e->p1->u + (1 - coord) * e->p2->u;
double V = coord * e->p1->v + (1 - coord) * e->p2->v;
- GPoint gpp = gf->point(U,V);
+ GPoint gpp = gf->point(m.scalingU*U,m.scalingV*V);
if (gpp.succeeded()){
mid = m.add_point(++m.MAXPOINTNUMBER, gpp.x(),gpp.y(),gpp.z());
mid->u = U;
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index a065e5860f658980bf2fe720eeb43e3f50815a44..dc31fb1ea8ae35b537add9a6c9b60f447954a620 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -69,7 +69,9 @@ void circumCenterMetric(double *pa, double *pb, double *pc,
d * (pa[1] * pa[1] - pc[1] * pc[1]) +
2. * b * (pa[0] * pa[1] - pc[0] * pc[1]);
- sys2x2(sys, rhs, x);
+ if (!sys2x2(sys, rhs, x)){
+ // printf("%g %g %g\n",a,b,d);
+ }
Radius2 =
(x[0] - pa[0]) * (x[0] - pa[0]) * a +
@@ -618,15 +620,15 @@ static void insertAPoint(GFace *gf,
uv[0] * vSizes [ptin->tri()->getVertex(1)->getNum()] +
uv[1] * vSizes [ptin->tri()->getVertex(2)->getNum()]);
double lc = BGM_MeshSize(gf,center[0],center[1],p.x(),p.y(),p.z());
- SMetric3 metr = BGM_MeshMetric(gf,center[0],center[1],p.x(),p.y(),p.z());
+ //SMetric3 metr = BGM_MeshMetric(gf,center[0],center[1],p.x(),p.y(),p.z());
// vMetricsBGM.push_back(metr);
vSizesBGM.push_back(lc);
vSizes.push_back(lc1);
Us.push_back(center[0]);
Vs.push_back(center[1]);
- if (!insertVertex(gf, v, center, worst, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
- Us, Vs, metric) || !p.succeeded()) {
+ if (!p.succeeded() || !insertVertex(gf, v, center, worst, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
+ Us, Vs, metric) ) {
// Msg::Debug("2D Delaunay : a cavity is not star shaped");
AllTris.erase(it);
worst->forceRadius(-1);
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 56f77296c3b8f65cbd8f16247c17be55b718b604..51bdb2843c129f0ff58239dad7c19cabd69895d4 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -787,7 +787,6 @@ void insertVerticesInRegion (GRegion *gr)
if(worst->isDeleted()){
myFactory.Free(worst);
allTets.erase(allTets.begin());
- //Msg::Info("Worst tet is deleted");
}
else{
if(ITER++ %5000 == 0)