diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 86fc4b682c423725488897151a572e2a8ec2ae1d..e0e88d0b9e524f48d659fe4f95ae35ef49e645f9 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -186,14 +186,15 @@ void GEdge::writeGEO(FILE *fp)
double umin = bounds.low();
double umax = bounds.high();
fprintf(fp, "p%d = newp;\n", tag());
- for(int i = 1; i < minimumDrawSegments(); i++){
- double u = umin + (double)i / minimumDrawSegments() * (umax - umin);
+ int N = minimumDrawSegments();
+ for(int i = 1; i < N; i++){
+ double u = umin + (double)i / N * (umax - umin);
GPoint p = point(u);
fprintf(fp, "Point(p%d + %d) = {%.16g, %.16g, %.16g};\n",
tag(), i, p.x(), p.y(), p.z());
}
fprintf(fp, "Spline(%d) = {%d", tag(), getBeginVertex()->tag());
- for(int i = 1; i < minimumDrawSegments(); i++)
+ for(int i = 1; i < N; i++)
fprintf(fp, ", p%d + %d", tag(), i);
fprintf(fp, ", %d};\n", getEndVertex()->tag());
}
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 5b2a64f787f6228b821b0be01bdb837184929ed4..cfd2b53a7770c1be3c063dc9e453febc298f6733 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -19,6 +19,7 @@
#include "Context.h"
#include "meshGFaceLloyd.h"
#include "Bindings.h"
+#include "meshGFaceOptimize.h"
#define SQU(a) ((a)*(a))
@@ -1204,6 +1205,30 @@ void GFace::moveToValidRange(SPoint2 &pt) const
}
}
+void GFace::addLayersOfQuads(int nLayers, GVertex *gv, double hmin, double ratio){
+
+ std::list<GEdgeLoop>::iterator it = edgeLoops.begin();
+ for (; it != edgeLoops.end(); ++it){
+ bool found = false;
+ // look if this edge loop has the GVertex as an endpoint
+ for (GEdgeLoop::iter it2 = it->begin(); it2 != it->end(); ++it2){
+ if (it2->ge->getBeginVertex() == gv || it2->ge->getEndVertex() == gv)
+ found = true;
+ }
+ // we found an edge loop with the GVertex that was specified
+ if (found){
+ for (GEdgeLoop::iter it2 = it->begin(); it2 != it->end(); ++it2){
+ GEdge *ge = it2->ge;
+ for (int i=0;i<ge->lines.size();i++){
+ SPoint2 p[2];
+ reparamMeshEdgeOnFace ( ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),gf,p[0],p[1]);
+ }
+ }
+ }
+ }
+}
+
+
#include "Bindings.h"
void GFace::registerBindings(binding *b)
@@ -1223,6 +1248,9 @@ void GFace::registerBindings(binding *b)
mb->setArgNames("quadrangle", NULL);
mb = cb->addMethod("edges", &GFace::edges);
mb->setDescription("return the list of edges bounding this surface");
+ mb = cb->addMethod("addLayersOfQuads", &GFace::addLayersOfQuads);
+ mb->setDescription("insert layers of quads");
+ mb->setArgNames("nLayers","startingVertex", NULL);
/* mb = cb->addMethod("addPolygon", &GFace::addPolygon);
mb->setDescription("insert a polygon mesh element");
mb->setArgNames("polygon", NULL);*/
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 71e76f2cce065b6f61d8b4a041030d130c23ce10..5ff9dff192175a1327a69c540b46007b8a72420b 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -319,6 +319,10 @@ class GFace : public GEntity
// tells if it's a sphere, and if it is, returns parameters
virtual bool isSphere (double &radius, SPoint3 ¢er) const {return false;}
+
+ // add layers of quads
+ void addLayersOfQuads (int nLayers, GVertex *start, double hmin, double factor);
+
};
#endif
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 3030182672711f3645a3c5c2f357449be4377593..f803e7751ed34052e2be6f4143d52a1b1fe5e9ad 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -41,6 +41,7 @@
#include "Generator.h"
#include "meshGFaceOptimize.h"
#include "meshPartition.h"
+#include "HighOrder.h"
#endif
std::vector<GModel*> GModel::list;
@@ -505,6 +506,19 @@ int GModel::mesh(int dimension)
#endif
}
+
+int GModel::setOrderN(int order, int linear, int incomplete)
+{
+#if defined(HAVE_MESH)
+ SetOrderN(this, order, linear, incomplete);
+ return true;
+#else
+ Msg::Error("Mesh module not compiled");
+ return false;
+#endif
+}
+
+
int GModel::getMeshStatus(bool countDiscrete)
{
for(riter it = firstRegion(); it != lastRegion(); ++it)
@@ -2553,6 +2567,9 @@ void GModel::registerBindings(binding *b)
cm = cb->addMethod("mesh", &GModel::mesh);
cm->setArgNames("dim", NULL);
cm->setDescription("Generate a mesh of this model in dimension 'dim'.");
+ cm = cb->addMethod("setOrderN", &GModel::setOrderN);
+ cm->setArgNames("order", "linear", "incomplete", NULL);
+ cm->setDescription(" make the mesh a high order mesh at order N\n linear is 1 if the high order points are not placed on the geometry of the model\n incomplete is 1 if incomplete basis are used.");
cm = cb->addMethod("load", &GModel::load);
cm->setDescription("Merge the file 'filename' in this model, the file can be "
"in any format (guessed from the extension) known by gmsh.");
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 731251c21403419121b2103b64a3caf59d17a0a7..3ead894b4ebeeec469c719d722ceae0506190222 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -370,6 +370,11 @@ class GModel
// mesh the model
int mesh(int dimension);
+ // make the mesh a high order mesh at order N
+ // linear is 1 if the high order points are not placed on the geometry of the model
+ // incomplete is 1 if incomplete basis are used
+ int setOrderN(int order, int linear, int incomplete);
+
// create partition boundaries
void createPartitionBoundaries(int createGhostCells);
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index c600e06584216f6342d584709b4f49bb90e68c9e..b2aeda6599fb73c24e05d591fbe59cccabb0c2d3 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -17,18 +17,24 @@ static Vertex InterpolateCubicSpline(Vertex *v[4], double t, double mat[4][4],
{
Vertex V;
int i, j;
- double vec[4], T[4];
+ double vec[4], T[4];
V.Pos.X = V.Pos.Y = V.Pos.Z = 0.0;
V.lc = (1 - t) * v[1]->lc + t * v[2]->lc;
V.w = (1 - t) * v[1]->w + t * v[2]->w;
- if(derivee) {
+ if(derivee == 1) {
T[3] = 0.;
T[2] = 1.;
T[1] = 2. * t;
T[0] = 3. * t * t;
}
+ else if(derivee == 2) {
+ T[3] = 0.;
+ T[2] = 0.;
+ T[1] = 2.;
+ T[0] = 6. * t;
+ }
else {
T[3] = 1.;
T[2] = t;
@@ -75,28 +81,47 @@ static Vertex InterpolateCubicSpline(Vertex *v[4], double t, double mat[4][4],
vec[j] = 0.0;
}
- if(derivee) {
+ if(derivee == 1) {
V.Pos.X /= ((t2 - t1));
V.Pos.Y /= ((t2 - t1));
V.Pos.Z /= ((t2 - t1));
}
+ else if(derivee == 2) {
+ V.Pos.X /= ((t2 - t1)*(t2 - t1));
+ V.Pos.Y /= ((t2 - t1)*(t2 - t1));
+ V.Pos.Z /= ((t2 - t1)*(t2 - t1));
+ }
return V;
}
// interpolation in the parametric space !
SPoint2 InterpolateCubicSpline(Vertex *v[4], double t, double mat[4][4],
- double t1, double t2, gmshSurface *s)
+ double t1, double t2, gmshSurface *s, int derivee)
{
Vertex V;
int i, j;
double T[4];
- T[3] = 1.;
- T[2] = t;
- T[1] = t * t;
- T[0] = t * t * t;
-
+ if (derivee == 0){
+ T[3] = 1.;
+ T[2] = t;
+ T[1] = t * t;
+ T[0] = t * t * t;
+ }
+ else if (derivee == 1){
+ T[3] = 0.;
+ T[2] = 1;
+ T[1] = 2 * t;
+ T[0] = 3 * t * t;
+ }
+ else if (derivee == 2){
+ T[3] = 0.;
+ T[2] = 0;
+ T[1] = 2;
+ T[0] = 6 * t;
+ }
+
SPoint2 coord [4], p;
for(i = 0; i < 4; i++) {
@@ -131,7 +156,7 @@ static Vertex InterpolateUBS(Curve *Curve, double u, int derivee)
}
if(Curve->geometry){
- SPoint2 pp = InterpolateCubicSpline(v, t, Curve->mat, t1, t2, Curve->geometry);
+ SPoint2 pp = InterpolateCubicSpline(v, t, Curve->mat, t1, t2, Curve->geometry,derivee);
SPoint3 pt = Curve->geometry->point(pp);
Vertex V;
V.Pos.X = pt.x();
@@ -140,7 +165,7 @@ static Vertex InterpolateUBS(Curve *Curve, double u, int derivee)
return V;
}
else
- return InterpolateCubicSpline(v, t, Curve->mat, 0, t1, t2);
+ return InterpolateCubicSpline(v, t, Curve->mat, derivee, t1, t2);
}
// Non Uniform BSplines
@@ -218,17 +243,48 @@ Vertex InterpolateCurve(Curve *c, double u, int derivee)
Vertex V;
- if(derivee) {
- double eps1 = (u == 0) ? 0 : 1.e-5;
- double eps2 = (u == 1) ? 0 : 1.e-5;
- Vertex D[2];
- D[0] = InterpolateCurve(c, u - eps1, 0);
- D[1] = InterpolateCurve(c, u + eps2, 0);
- V.Pos.X = (D[1].Pos.X - D[0].Pos.X) / (eps1 + eps2);
- V.Pos.Y = (D[1].Pos.Y - D[0].Pos.Y) / (eps1 + eps2);
- V.Pos.Z = (D[1].Pos.Z - D[0].Pos.Z) / (eps1 + eps2);
- V.u = u;
- return V;
+ if(derivee==1) {
+ switch (c->Typ) {
+ // case MSH_SEGM_BSPLN:
+ // case MSH_SEGM_BEZIER:
+ // V = InterpolateUBS(c, u, 1);
+ // V.u = u;
+ // break;
+ default :
+ double eps1 = (u == 0) ? 0 : 1.e-5;
+ double eps2 = (u == 1) ? 0 : 1.e-5;
+ Vertex D[2];
+ D[0] = InterpolateCurve(c, u - eps1, 0);
+ D[1] = InterpolateCurve(c, u + eps2, 0);
+ V.Pos.X = (D[1].Pos.X - D[0].Pos.X) / (eps1 + eps2);
+ V.Pos.Y = (D[1].Pos.Y - D[0].Pos.Y) / (eps1 + eps2);
+ V.Pos.Z = (D[1].Pos.Z - D[0].Pos.Z) / (eps1 + eps2);
+ V.u = u;
+ break;
+ }
+ return V;
+ }
+
+ if(derivee==2) {
+ switch (c->Typ) {
+ case MSH_SEGM_BSPLN:
+ case MSH_SEGM_BEZIER:
+ V = InterpolateUBS(c, u, 2);
+ V.u = u;
+ break;
+ default :
+ double eps1 = (u == 0) ? 0 : 1.e-5;
+ double eps2 = (u == 1) ? 0 : 1.e-5;
+ Vertex D[2];
+ D[0] = InterpolateCurve(c, u - eps1, 1);
+ D[1] = InterpolateCurve(c, u + eps2, 1);
+ V.Pos.X = (D[1].Pos.X - D[0].Pos.X) / (eps1 + eps2);
+ V.Pos.Y = (D[1].Pos.Y - D[0].Pos.Y) / (eps1 + eps2);
+ V.Pos.Z = (D[1].Pos.Z - D[0].Pos.Z) / (eps1 + eps2);
+ V.u = u;
+ break;
+ }
+ return V;
}
int N, i;
@@ -344,7 +400,7 @@ Vertex InterpolateCurve(Curve *c, double u, int derivee)
List_Read(c->Control_Points, i + 2, &v[3]);
}
if(c->geometry){
- SPoint2 pp = InterpolateCubicSpline(v, t, c->mat, t1, t2,c->geometry);
+ SPoint2 pp = InterpolateCubicSpline(v, t, c->mat, t1, t2,c->geometry,0);
SPoint3 pt = c->geometry->point(pp);
V.Pos.X = pt.x();
V.Pos.Y = pt.y();
diff --git a/Geo/GeoInterpolation.h b/Geo/GeoInterpolation.h
index 341b751728f1866688849090f9804a1a1815dc0a..87806593fe9c98910a82789dc2091a9bce42234e 100644
--- a/Geo/GeoInterpolation.h
+++ b/Geo/GeoInterpolation.h
@@ -12,5 +12,5 @@ bool iSRuledSurfaceASphere(Surface *s, SPoint3 ¢er, double &radius);
Vertex InterpolateCurve(Curve *Curve, double u, int derivee);
Vertex InterpolateSurface(Surface *s, double u, double v, int derivee, int u_v);
SPoint2 InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],
- double t1, double t2, gmshSurface *s);
+ double t1, double t2, gmshSurface *s, int derivee);
#endif
diff --git a/Geo/MElementOctree.cpp b/Geo/MElementOctree.cpp
index 895e866e30b8c4c73e518b900d3c1f215ba582a4..439bc9d9cc88650dff061e0559fb388d46a8629c 100644
--- a/Geo/MElementOctree.cpp
+++ b/Geo/MElementOctree.cpp
@@ -115,13 +115,13 @@ MElement *MElementOctree::find(double x, double y, double z, int dim)
{
double P[3] = {x, y, z};
MElement *e = (MElement*)Octree_Search(P, _octree);
-
- if (!e){
+
+ if (!e || (dim != -1 && e->getDim() != dim)){
double initialTol = MElement::getTolerance();
double tol = initialTol;
- while (tol < .1){
+ while (tol < 1){
tol *= 10;
- MElement::setTolerance(tol);
+ MElement::setTolerance(tol);
std::vector<GEntity*> entities;
_gm->getEntities(entities);
for(unsigned int i = 0; i < entities.size(); i++){
@@ -129,6 +129,7 @@ MElement *MElementOctree::find(double x, double y, double z, int dim)
e = entities[i]->getMeshElement(j);
if (dim == -1 || e->getDim() == dim){
if (MElementInEle(e, P)){
+ // printf("coucou FOUND\n");
MElement::setTolerance(initialTol);
return e;
}
@@ -136,6 +137,8 @@ MElement *MElementOctree::find(double x, double y, double z, int dim)
}
}
}
+ MElement::setTolerance(initialTol);
+ Msg::Warning("Point %g %g %g not found",x,y,z);
}
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 8ed49986f6e36c3042bdcf3c7ad97036e21249f3..92d0aab1fa3fd62b27f9f6b0f9c2880528226e47 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -61,12 +61,13 @@ double MTriangle::etaShapeMeasure()
double a1 = 180 * angle3Vertices(_v[0], _v[1], _v[2]) / M_PI;
double a2 = 180 * angle3Vertices(_v[1], _v[2], _v[0]) / M_PI;
double a3 = 180 * angle3Vertices(_v[2], _v[0], _v[1]) / M_PI;
- double angle = fabs(60. - a1);
- angle = std::max(fabs(60. - a2),angle);
- angle = std::max(fabs(60. - a3),angle);
+
+ double amin = std::min(std::min(a1,a2),a3);
+ double angle = fabs(60. - amin);
return 1.-angle/60;
}
+
double MTriangle::gammaShapeMeasure()
{
#if defined(HAVE_MESH)
diff --git a/Geo/STensor3.h b/Geo/STensor3.h
index db90bf13a4f417812af6db326c252cd98d5fabeb..c7ddeab713c710ea04bb58ca8a8959228b7d3f61 100644
--- a/Geo/STensor3.h
+++ b/Geo/STensor3.h
@@ -33,7 +33,7 @@ class SMetric3 {
void setMat(const fullMatrix<double> & mat)
{
for (int i = 0; i < 3; i++)
- for (int j = i; j < 3; j++)
+ for (int j = 0; j < 3; j++)
_val[getIndex(i, j)] = mat(i, j);
}
SMetric3(const SMetric3& m)
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index bf13b69555c74b2e392eadac482170c923dcc2e3..413c820c7c75ded899dd8184bfcac00a76936bb7 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -46,6 +46,13 @@ SVector3 gmshEdge::firstDer(double par) const
return SVector3(a.Pos.X, a.Pos.Y, a.Pos.Z);
}
+SVector3 gmshEdge::secondDer(double par) const
+{
+ // printf("coucou mon chou\n");
+ Vertex a = InterpolateCurve(c, par, 2);
+ return SVector3(a.Pos.X, a.Pos.Y, a.Pos.Z);
+}
+
GEntity::GeomType gmshEdge::geomType() const
{
switch (c->Typ){
@@ -144,7 +151,7 @@ SPoint2 gmshEdge::reparamOnFace(const GFace *face, double epar,int dir) const
k = periodic ? k - NbControlPoints + 1: NbControlPoints - 1;
List_Read(c->Control_Points, k, &v[j]);
}
- return InterpolateCubicSpline(v, t, c->mat, t1, t2, c->geometry);
+ return InterpolateCubicSpline(v, t, c->mat, t1, t2, c->geometry,0);
}
case MSH_SEGM_SPLN :
{
@@ -185,7 +192,7 @@ SPoint2 gmshEdge::reparamOnFace(const GFace *face, double epar,int dir) const
else{
List_Read(c->Control_Points, i + 2, &v[3]);
}
- return InterpolateCubicSpline(v, t, c->mat, t1, t2, c->geometry);
+ return InterpolateCubicSpline(v, t, c->mat, t1, t2, c->geometry,0);
}
default:
Msg::Error("Unknown edge type in reparamOnFace");
diff --git a/Geo/gmshEdge.h b/Geo/gmshEdge.h
index 48c220cd4276bd8a25281bdf0c3490ecef8b914f..2915188a28265134b0c374063572067800c4f40a 100644
--- a/Geo/gmshEdge.h
+++ b/Geo/gmshEdge.h
@@ -20,6 +20,7 @@ class gmshEdge : public GEdge {
virtual GeomType geomType() const;
virtual GPoint point(double p) const;
virtual SVector3 firstDer(double par) const;
+ virtual SVector3 secondDer(double par) const;
ModelType getNativeType() const { return GmshModel; }
void * getNativePtr() const { return c; }
virtual std::string getAdditionalInfoString();
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 026f80a21ff84ca2809d5d61f8a6a52e4be2ec35..e05c8bcdc2c5ab8fb90cac09bfc17a06287b3b07 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -32,6 +32,64 @@
// CTX::instance()->mesh.minCircPoints tells the minimum number of points per
// radius of curvature
+SMetric3 buildMetricTangentToCurve (SVector3 &t, double curvature, double &lambda){
+ lambda = 1.e22;
+ if (curvature == 0.0)return SMetric3(1.e-22);
+ SVector3 a;
+ if (t(0) <= t(1) && t(0) <= t(2)){
+ a = SVector3(1,0,0);
+ }
+ else if (t(1) <= t(0) && t(1) <= t(2)){
+ a = SVector3(0,1,0);
+ }
+ else{
+ a = SVector3(0,0,1);
+ }
+ SVector3 b = crossprod (t,a);
+ SVector3 c = crossprod (b,t);
+ b.normalize();
+ c.normalize();
+ t.normalize();
+ lambda = ((2 * M_PI) /( fabs(curvature) * CTX::instance()->mesh.minCircPoints ) );
+
+ SMetric3 curvMetric (1./(lambda*lambda),1.e-10,1.e-10,t,b,c);
+ //SMetric3 curvMetric (1./(lambda*lambda));
+ return curvMetric;
+}
+
+SMetric3 max_edge_curvature_metric(const GVertex *gv, double &length)
+{
+ SMetric3 val (1.e-12);
+ length = 1.e22;
+ std::list<GEdge*> l_edges = gv->edges();
+ for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
+ ite != l_edges.end(); ++ite){
+ GEdge *_myGEdge = *ite;
+ Range<double> range = _myGEdge->parBounds(0);
+ SMetric3 cc;
+ double l;
+ if (gv == _myGEdge->getBeginVertex()) {
+ SVector3 t = _myGEdge->firstDer(range.low());
+ t.normalize();
+ cc = buildMetricTangentToCurve(t,_myGEdge->curvature(range.low()),l);
+ }
+ else {
+ SVector3 t = _myGEdge->firstDer(range.high());
+ t.normalize();
+ cc = buildMetricTangentToCurve(t,_myGEdge->curvature(range.high()),l);
+ }
+ val = intersection(val,cc);
+ length = std::min(length,l);
+ }
+ return val;
+}
+
+SMetric3 max_edge_curvature_metric(const GEdge *ge, double u, double &l){
+ SVector3 t = ge->firstDer(u);
+ t.normalize();
+ return buildMetricTangentToCurve(t,ge->curvature(u),l);
+}
+
static double max_edge_curvature(const GVertex *gv)
{
double val = 0;
@@ -161,16 +219,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)*2;
- // 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/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index d4b818e7878a0bbb8cb90ecf5caa0e5565847ede..91f155181a7911a5bc37855bb6ab94e54af18e5b 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -12,6 +12,8 @@
class MElementOctree;
class GFace;
+class GVertex;
+class GEdge;
class MElement;
class MVertex;
class GEntity;
@@ -42,5 +44,7 @@ double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double
SMetric3 BGM_MeshMetric(GEntity *ge, double U, double V, double X, double Y, double Z);
bool Extend1dMeshIn2dSurfaces();
bool Extend2dMeshIn3dVolumes();
+SMetric3 max_edge_curvature_metric(const GVertex *gv, double &l);
+SMetric3 max_edge_curvature_metric(const GEdge *ge, double u, double &l);
#endif
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index b1288474ffa6a8a6265e881d7bf5625851ba8f24..7f937607c2a555e5e42f46c04677fac33b54c504 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1555,7 +1555,7 @@ public:
const double ll1 = dist*(ratio-1) + hwall_n;
const double lc_n = std::min(ll1,hfar);
const double ll2 = dist*(ratio-1) + hwall_t;
- const double lc_t = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2,hfar));
+ double lc_t = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2,hfar));
SVector3 t1,t2,t3;
double L1,L2,L3;
@@ -1565,8 +1565,22 @@ public:
std::pair<AttractorInfo,SPoint3> pp = cc->getAttractorInfo();
if (pp.first.dim ==1){
GEdge *e = GModel::current()->getEdgeByTag(pp.first.ent);
+
+
+ // the tangent size at this point is the size of the
+ // 1D mesh at this point !
+ // hack : use curvature
+ if(CTX::instance()->mesh.lcFromCurvature){
+ double Crv = e->curvature(pp.first.u);
+ double lc = Crv > 0 ? 2 * M_PI / Crv / CTX::instance()->mesh.minCircPoints : 1.e-22;
+ const double ll2b = dist*(ratio-1) + lc;
+ double lc_tb = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2b,hfar));
+ lc_t = std::min(lc_t,lc_tb);
+ }
+
+
if (dist < hwall_n){
- L1 = lc_t;
+ L1 = lc_t;
L2 = lc_n;
L3 = lc_n;
t1 = e->firstDer(pp.first.u);
@@ -1583,7 +1597,7 @@ public:
// printf("hfar = %g lc = %g dir %g %g \n",hfar,lc,t1.x(),t1.y());
}
else {
- L1 = lc_t;
+ L1 = lc_t;
L2 = lc_n;
L3 = lc_t;
GPoint p = e->point(pp.first.u);
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 6784b09e211022f3b12d5a94f3638a30b622f72a..3882923b207854eda2b973d04cdbb2fc39e9bb4d 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1166,7 +1166,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
if(displ2D) delete displ2D;
if(displ3D) delete displ3D;
- // printJacobians(m, "smoothness.pos");
+ // printJacobians(m, "smoothness.pos");
double t2 = Cpu();
Msg::StatusBar(2, true, "Done meshing order %d (%g s)", order, t2 - t1);
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index 0abbd586a7c32d25780e121548fb4ddc17666641..9dea04cb76c51824508ca0c4bf0b761356aff14c 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -354,8 +354,8 @@ void highOrderSmoother::smooth(GRegion *gr)
void highOrderSmoother::optimize(GFace * gf,
- edgeContainer &edgeVertices,
- faceContainer &faceVertices)
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices)
{
// if (gf->geomType() != GEntity::Plane) return;
@@ -447,7 +447,7 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*> & all, GFace *gf)
lsys->setPrec(5.e-8);
#endif
dofManager<double> myAssembler(lsys);
- elasticityTerm El(0, 1.0, .48, getTag());
+ elasticityTerm El(0, 1.0, .33, getTag());
std::vector<MElement*> layer, v;
double minD;
@@ -517,11 +517,11 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*> & all, GFace *gf)
double dx0 = smooth_metric_(v, gf, myAssembler, verticesToMove, El);
double dx = dx0;
- printf(" dx0 = %12.5E\n", dx0);
+ Msg::Debug(" dx0 = %12.5E", dx0);
int iter = 0;
while(1){
double dx2 = smooth_metric_(v, gf, myAssembler, verticesToMove, El);
- printf(" dx2 = %12.5E\n", dx2);
+ Msg::Debug(" dx2 = %12.5E", dx2);
if (fabs(dx2 - dx) < 1.e-4 * dx0)break;
if (iter++ > 2)break;
dx = dx2;
@@ -560,6 +560,12 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
if (myAssembler.sizeOfR()){
+ // initialize _materialLaw to 1 everywhere
+ for (unsigned int i = 0; i < v.size(); i++){
+ MElement *e = v[i];
+ // _materialLaw[e] = 1.0;
+ }
+ // while convergence -------------------------------------------------------
for (unsigned int i = 0; i < v.size(); i++){
MElement *e = v[i];
int nbNodes = e->getNumVertices();
@@ -575,6 +581,8 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
fullMatrix<double> J23K33(n2, n3);
K33.setAll(0.0);
SElement se(e);
+ // set kappa to the actual kappa of the material law
+ // El.setCompressibility( compressibilityFunction ( _materialLaw[e] ) );
El.elementMatrix(&se, K33);
computeMetricVector(gf, e, J32, J23, D3);
J23K33.gemm(J23, K33, 1, 0);
@@ -590,6 +598,8 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
}
}
myAssembler.systemSolve();
+ // for all element, compute detJ at integration points --> material law
+ // end while convergence -------------------------------------------------------
for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
if ((*it)->onWhat()->dim() == 2){
diff --git a/Mesh/highOrderSmoother.h b/Mesh/highOrderSmoother.h
index 15d7737414e6ac49a9b744fe5644a307657e6848..9ff58bc77043b8d448697d07bfcffe0454472219 100644
--- a/Mesh/highOrderSmoother.h
+++ b/Mesh/highOrderSmoother.h
@@ -26,6 +26,7 @@ class GRegion;
class highOrderSmoother
{
const int _tag;
+ std::map<MElement*,std::vector<double> > _materialLaw;
std::map<MVertex*,SVector3> _straightSidedLocation;
std::map<MVertex*,SVector3> _targetLocation;
int _dim;
@@ -38,6 +39,7 @@ class highOrderSmoother
void moveTo(MVertex *v, const std::map<MVertex*,SVector3> &) const;
public:
highOrderSmoother(int dim) : _tag(111), _dim(dim) {_clean();}
+ highOrderSmoother(GModel *gm);
void add(MVertex * v, const SVector3 &d ) {
_straightSidedLocation[v] = d;
_targetLocation[v] = SPoint3(v->x(),v->y(),v->z());
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 8bf75cd732060990b87b5f21635583fef6f35dd6..8c207127fd9350c47211367fcae2bebec0cef556 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -883,8 +883,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
bowyerWatson(gf);
meshGFaceBamg(gf);
}
- for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++)
- laplaceSmoothing(gf);
+ laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
}
if(debug){
@@ -1465,8 +1464,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
bowyerWatson(gf);
else
meshGFaceBamg(gf);
- for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++)
- laplaceSmoothing(gf);
+ laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
}
// delete the mesh
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index f3be29529ed1485570c517a79203ba73d558bb55..a860fe5c74d72ee1d5ac77a2e862f525d1a3399d 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -12,6 +12,7 @@
#include "MVertex.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
+#include "MLine.h"
#include "BackgroundMesh.h"
#include "Numeric.h"
#include "GmshMessage.h"
@@ -497,10 +498,10 @@ static bool _isItAGoodIdeaToCollapseThatVertex (GFace *gf,
double uu,vv;
v->getParameter(0,uu);
v->getParameter(1,vv);
- SPoint2 p1,p2;
- reparamMeshEdgeOnFace(v1,v2,gf,p1,p2);
- SPoint2 p =(p1+p2)*0.5;
- GPoint gp = gf->point(p);
+ // SPoint2 p1,p2;
+ // reparamMeshEdgeOnFace(v1,v2,gf,p1,p2);
+ // SPoint2 p =(p1+p2)*0.5;
+ // GPoint gp = gf->point(p);
// v->setXYZ(gp.x(),gp.y(),gp.z());
// v->setParameter(0,p.x());
// v->setParameter(1,p.y());
@@ -537,6 +538,7 @@ static bool _isItAGoodIdeaToCollapseThatVertex (GFace *gf,
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);
v->setParameter(1,vv);
@@ -557,22 +559,53 @@ static bool _isItAGoodIdeaToMoveThatVertex (GFace *gf,
double surface_old = 0;
double surface_new = 0;
- for (unsigned int j=0;j<e1.size();++j)
+ double qualityOld[256],qualityNew[256];
+
+ GPoint gp = gf->point(after);
+ 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);
+ qualityOld[j] = (e1[j]->getNumVertices() == 4) ? e1[j]->etaShapeMeasure() : e1[j]->gammaShapeMeasure();
+ }
v1->setParameter(0,after.x());
v1->setParameter(1,after.y());
+ v1->x() = pafter.x();
+ v1->y() = pafter.y();
+ v1->z() = pafter.z();
- for (unsigned int j=0;j<e1.size();++j)
+ for (unsigned int j=0;j<e1.size();++j){
surface_new += surfaceFaceUV(e1[j],gf);
+ qualityNew[j] = (e1[j]->getNumVertices() == 4) ? e1[j]->etaShapeMeasure() : e1[j]->gammaShapeMeasure();
+ }
v1->setParameter(0,before.x());
v1->setParameter(1,before.y());
+ v1->x() = pbefore.x();
+ v1->y() = pbefore.y();
+ v1->z() = pbefore.z();
if ( surface_new - surface_old > 1.e-10 * surface_old) {
return false;
}
return true;
+ int nBetter=0,nWorse=0;
+ int nReallyBadOld=0,nReallyBadNew=0;
+ for (unsigned int j=0;j<e1.size();++j){
+ if (qualityNew[j] >= qualityOld[j])nBetter++;
+ else {
+ nWorse++;
+ }
+ if (qualityNew[j] < 0.2) nReallyBadNew ++;
+ if (qualityOld[j] < 0.2) nReallyBadOld ++;
+ }
+
+ if (nReallyBadNew < nReallyBadOld)return true;
+ if (nReallyBadOld < nReallyBadNew)return false;
+
+ return (nBetter >= nWorse);
}
@@ -589,15 +622,15 @@ static int _quadWithOneVertexOnBoundary (GFace *gf,
MVertex *v2,
MVertex *v3,
MVertex *v4){
- // return 0;
+ return 0;
if (v1->onWhat()->dim() < 2 &&
v2->onWhat()->dim() == 2 &&
v3->onWhat()->dim() == 2 &&
v4->onWhat()->dim() == 2 &&
e2.size() < 5 &&
e4.size() < 5 &&
- ( _isItAGoodIdeaToCollapseThatVertex (gf, e2, e4, q, v2, v4, v4,12.0) ||
- _isItAGoodIdeaToCollapseThatVertex (gf, e2, e4, q, v2, v4, v2,12.0))){
+ _isItAGoodIdeaToCollapseThatVertex (gf, e2, e4, q, v2, v4, v4,12.0)
+ /* || _isItAGoodIdeaToCollapseThatVertex (gf, e2, e4, q, v2, v4, v2,12.0))*/){
touched.insert(v1);
touched.insert(v2);
touched.insert(v3);
@@ -686,8 +719,7 @@ static int _removeDiamonds(GFace *gf)
v1->onWhat()->dim() == 2 &&
v3->onWhat()->dim() == 2 &&
it1->second.size() ==3 && it3->second.size() == 3 &&
- (_isItAGoodIdeaToCollapseThatVertex (gf, it1->second, it3->second, q, v1, v3, v3,10.) ||
- _isItAGoodIdeaToCollapseThatVertex (gf, it1->second, it3->second, q, v1, v3, v1,10.) ) ){
+ _isItAGoodIdeaToCollapseThatVertex (gf, it1->second, it3->second, q, v1, v3, v3,10.)){
touched.insert(v1);
touched.insert(v2);
touched.insert(v3);
@@ -701,13 +733,31 @@ static int _removeDiamonds(GFace *gf)
deleted.insert(v1);
diamonds.insert(q);
}
+ else if (v2->onWhat()->dim() == 2 &&
+ v4->onWhat()->dim() == 2 &&
+ v1->onWhat()->dim() == 2 &&
+ v3->onWhat()->dim() == 2 &&
+ it1->second.size() ==3 && it3->second.size() == 3 &&
+ _isItAGoodIdeaToCollapseThatVertex (gf, it1->second, it3->second, q, v1, v3, v1,10.)){
+ touched.insert(v1);
+ touched.insert(v2);
+ touched.insert(v3);
+ touched.insert(v4);
+ for (unsigned int j=0;j<it3->second.size();++j){
+ for (int k=0;k<it3->second[j]->getNumVertices();k++){
+ if (it3->second[j]->getVertex(k) == v3 && it3->second[j] != q)
+ it3->second[j]->setVertex(k,v1);
+ }
+ }
+ deleted.insert(v3);
+ diamonds.insert(q);
+ }
else if (v1->onWhat()->dim() == 2 &&
v3->onWhat()->dim() == 2 &&
v2->onWhat()->dim() == 2 &&
v4->onWhat()->dim() == 2 &&
it2->second.size() == 3 && it4->second.size() == 3 &&
- (_isItAGoodIdeaToCollapseThatVertex (gf, it2->second, it4->second, q, v2, v4, v4,10.) ||
- _isItAGoodIdeaToCollapseThatVertex (gf, it2->second, it4->second, q, v2, v4, v2,10.) ) ){
+ _isItAGoodIdeaToCollapseThatVertex (gf, it2->second, it4->second, q, v2, v4, v4,10.)){
touched.insert(v1);
touched.insert(v2);
touched.insert(v3);
@@ -721,6 +771,25 @@ static int _removeDiamonds(GFace *gf)
deleted.insert(v2);
diamonds.insert(q);
}
+ else if (v1->onWhat()->dim() == 2 &&
+ v3->onWhat()->dim() == 2 &&
+ v2->onWhat()->dim() == 2 &&
+ v4->onWhat()->dim() == 2 &&
+ it2->second.size() == 3 && it4->second.size() == 3 &&
+ _isItAGoodIdeaToCollapseThatVertex (gf, it2->second, it4->second, q, v2, v4, v2,10.)){
+ touched.insert(v1);
+ touched.insert(v2);
+ touched.insert(v3);
+ touched.insert(v4);
+ for (unsigned int j=0;j<it4->second.size();++j){
+ for (int k=0;k<it4->second[j]->getNumVertices();k++){
+ if (it4->second[j]->getVertex(k) == v4 && it4->second[j] != q)
+ it4->second[j]->setVertex(k,v2);
+ }
+ }
+ deleted.insert(v4);
+ diamonds.insert(q);
+ }
else {
quadrangles2.push_back(q);
}
@@ -755,6 +824,88 @@ int removeDiamonds(GFace *gf){
return nbRemove;
}
+struct p1p2p3 {
+ MVertex *p1,*p2;
+};
+static void _relocateVertexConstrained (GFace *gf,
+ MVertex *ver,
+ const std::vector<MElement*> <){
+ if( ver->onWhat()->dim() == 2){
+
+
+ std::map<MVertex*,p1p2p3> ring;
+ double initu,initv;
+ ver->getParameter(0, initu);
+ ver->getParameter(1, initv);
+ double XX=0,YY=0,ZZ=0;
+ for(unsigned int i = 0; i < lt.size(); i++){
+ for (int j=0;j<lt[i]->getNumEdges();j++){
+ MEdge ed = lt[i]->getEdge(j);
+ if (ed.getVertex(0) != ver && ed.getVertex(1) != ver){
+ std::map<MVertex*,p1p2p3>::iterator it = ring.find(ed.getVertex(0));
+ if (it == ring.end()){
+ p1p2p3 a;
+ a.p1 = ed.getVertex(1);
+ ring[ed.getVertex(0)] = a;
+ }
+ else{
+ it->second.p2 = ed.getVertex(1);
+ }
+ it = ring.find(ed.getVertex(1));
+ if (it == ring.end()){
+ p1p2p3 a;
+ a.p1 = ed.getVertex(0);
+ ring[ed.getVertex(1)] = a;
+ }
+ else{
+ it->second.p2 = ed.getVertex(0);
+ }
+ }
+ }
+ }
+
+ double cu=0,cv=0;
+ std::map<MVertex*,p1p2p3>::iterator it = ring.begin();
+ for (; it!=ring.end();++it){
+ MVertex *center = it->first;
+ MVertex *left = it->second.p1;
+ MVertex *right = it->second.p2;
+ SPoint2 pcenter,pleft,pright;
+ reparamMeshVertexOnFace(center, gf, pcenter);
+ reparamMeshVertexOnFace(left, gf, pleft);
+ reparamMeshVertexOnFace(right, gf, pright);
+ SVector3 vj (pcenter.x()-initu,pcenter.y()-initv,0);
+ SVector3 vjp1 (pleft.x()-initu,pleft.y()-initv,0);
+ SVector3 vjm1 (pright.x()-initu,pright.y()-initv,0);
+ vjp1.normalize();
+ vjm1.normalize();
+ SVector3 bissector = (vjp1+vjm1)*0.5;
+ double dist = vj.norm();
+ cu += (pcenter.x() + dist * bissector.x());
+ cv += (pcenter.y() + dist * bissector.y());
+ }
+ cu/= ring.size();
+ cv/= ring.size();
+ SPoint2 before(initu,initv);
+ SPoint2 after(cu,cv);
+ // if (isSphere){
+ // GPoint gp = gf->closestPoint(SPoint3(XX/fact, YY/fact, ZZ/fact), after);
+ // after = SPoint2(gp.u(),gp.v());
+ // }
+ bool success = _isItAGoodIdeaToMoveThatVertex (gf, lt, ver,before,after);
+ if (success){
+ ver->setParameter(0, after.x());
+ ver->setParameter(1, after.y());
+ GPoint pt = gf->point(after);
+ if(pt.succeeded()){
+ ver->x() = pt.x();
+ ver->y() = pt.y();
+ ver->z() = pt.z();
+ }
+ }
+ }
+}
+
static void _relocateVertex (GFace *gf,
MVertex *ver,
const std::vector<MElement*> <){
@@ -790,7 +941,15 @@ static void _relocateVertex (GFace *gf,
GPoint gp = gf->closestPoint(SPoint3(XX/fact, YY/fact, ZZ/fact), after);
after = SPoint2(gp.u(),gp.v());
}
- bool success = _isItAGoodIdeaToMoveThatVertex (gf, lt, ver,before,after);
+ bool success = false;
+ double factor = 1.0;
+ for (int i=0;i<10;i++){
+ SPoint2 newp = after + before*(1.-factor);
+ success = _isItAGoodIdeaToMoveThatVertex (gf, lt, ver,before,newp);
+ if (success)break;
+ factor *= 0.5;
+ }
+
if (success){
ver->setParameter(0, after.x());
ver->setParameter(1, after.y());
@@ -805,12 +964,12 @@ static void _relocateVertex (GFace *gf,
}
}
-void laplaceSmoothing(GFace *gf)
+void laplaceSmoothing(GFace *gf, int niter)
{
v2t_cont adj;
buildVertexToElement(gf->triangles, adj);
buildVertexToElement(gf->quadrangles, adj);
- for(int i = 0; i < 5; i++){
+ for(int i = 0; i < niter; i++){
v2t_cont :: iterator it = adj.begin();
while (it != adj.end()){
_relocateVertex(gf,it->first,it->second);
@@ -819,9 +978,23 @@ void laplaceSmoothing(GFace *gf)
}
}
+void laplaceSmoothingConstrained(GFace *gf)
+{
+ v2t_cont adj;
+ buildVertexToElement(gf->triangles, adj);
+ buildVertexToElement(gf->quadrangles, adj);
+ for(int i = 0; i < 5; i++){
+ v2t_cont :: iterator it = adj.begin();
+ while (it != adj.end()){
+ _relocateVertexConstrained(gf,it->first,it->second);
+ ++it;
+ }
+ }
+}
+
-int _edgeSwapQuadsForBetterQuality ( GFace *gf ) {
+int _edgeSwapQuadsForBetterQuality ( GFace *gf ) {
e2t_cont adj;
// buildEdgeToElement(gf->triangles, adj);
buildEdgeToElement(gf->quadrangles, adj);
@@ -1566,16 +1739,13 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
}
else{
// TEST
- // FILE *f = fopen (MATCHFILE,"r");
- // int nn,ne;
- // fscanf(f,"%d %d",&nn,&ne);
std::vector<std::pair<MElement*,MElement*> > toProcess;
for (int k=0;k<elist[0];k++){
int i1 = elist[1+3*k],i2 = elist[1+3*k+1],an=elist[1+3*k+2];
// FIXME !!
if (an == 100000){// || an == 1000){
toProcess.push_back(std::make_pair(n2t[i1],n2t[i2]));
- Msg::Debug("Forbidden quad found in blossom quadrangulation, optimization will be required");
+ Msg::Debug("Extra edge found in blossom algorithm, optimization will be required");
}
else{
MElement *t1 = n2t[i1];
@@ -1675,35 +1845,33 @@ void recombineIntoQuads(GFace *gf,
double t1 = Cpu();
- //gf->model()->writeMSH("before.msh");
+ gf->model()->writeMSH("before.msh");
if (topologicalOpti)removeFourTrianglesNodes(gf, false);
int success = _recombineIntoQuads(gf,0);
- //gf->model()->writeMSH("raw.msh");
+ // gf->addLayersOfQuads(1, 0);
+ gf->model()->writeMSH("raw.msh");
if (nodeRepositioning)
- for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++)
- laplaceSmoothing(gf);
+ laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
//blossom-quad algo
if(success && CTX::instance()->mesh.algoRecombine == 1){
if (topologicalOpti){
- //gf->model()->writeMSH("smoothed.msh");
+ gf->model()->writeMSH("smoothed.msh");
int COUNT = 0;
- //char NAME[256];
+ char NAME[256];
while(1){
int x = removeTwoQuadsNodes(gf);
- //if(x){ sprintf(NAME,"iter%dTQ.msh",COUNT++); gf->model()->writeMSH(NAME);}
+ if(x){ sprintf(NAME,"iter%dTQ.msh",COUNT++); gf->model()->writeMSH(NAME);}
int y= removeDiamonds(gf);
- //if(y){ sprintf(NAME,"iter%dD.msh",COUNT++); gf->model()->writeMSH(NAME); }
+ if(y){ sprintf(NAME,"iter%dD.msh",COUNT++); gf->model()->writeMSH(NAME); }
laplaceSmoothing(gf);
int z = 0 ;//edgeSwapQuadsForBetterQuality ( gf );
- //if(z){ sprintf(NAME,"iter%dS.msh",COUNT++); gf->model()->writeMSH(NAME); }
+ if(z){ sprintf(NAME,"iter%dS.msh",COUNT++); gf->model()->writeMSH(NAME); }
if (!(x+y+z))break;
}
edgeSwapQuadsForBetterQuality ( gf );
//if(z){ sprintf(NAME,"iter%dS.msh",COUNT++); gf->model()->writeMSH(NAME); }
- for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++){
- laplaceSmoothing(gf);
- }
+ laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
}
double t2 = Cpu();
Msg::Info("Blossom recombination algorithm completed (%g s)", t2 - t1);
@@ -1712,11 +1880,9 @@ void recombineIntoQuads(GFace *gf,
//simple recombination algo
_recombineIntoQuads(gf,0);
- for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++)
- laplaceSmoothing(gf);
+ laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
_recombineIntoQuads(gf,0);
- for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++)
- laplaceSmoothing(gf);
+ laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
// gf->model()->writeMSH("after.msh");
double t2 = Cpu();
@@ -1724,22 +1890,144 @@ void recombineIntoQuads(GFace *gf,
}
-void quadsToTriangles(GFace *gf){
+
+// give it a try : add one quad layer on the
+void addOneLayerOnContour(GFace *gf, GVertex *gv){
+ /*
+, int nbLayers, double hplus, double factor){
+ // for each vertex
+ std::map<MVertex*,std::vector<MVertex*> >layers;
+ std::vector<MQuadrangle*> newQuads;
+ std::vector<MTriangle*> newTris;
+
+ std::list<GEdgeLoop>::iterator it = gf->edgeLoops.begin();
+ for (; it != gf->edgeLoops.end(); ++it){
+ bool found = false;
+ std::list<GEdge*> ed;
+ for (GEdgeLoop::iter it2 = it->begin(); it2 != it->end(); ++it2){
+ if (it2->ge->getBeginVertex() == gv || it2->ge->getEndVertex() == gv) {
+ found = true;
+ }
+ ed.push_back(it2->ge);
+ }
+ // we found an edge loop with the GVertex that was specified
+ if (found){
+ // compute model vertices that will produce fans
+ for (GEdgeLoop::iter it2 = it->begin(); it2 != it->end(); ++it2){
+ GEdgeLoop::iter it3 = it2; ++it3;
+ GVertex *gv = it2->getEndVertex();
+ GEdgeSigned *before,*after = *it2;
+ if (it3 == it->end()){
+ before = *(it->begin());
+ }
+ else{
+ before = *it2;
+ }
+ }
+
+ for (std::list<GEdge*>::iterator it = ed.begin(); it != ed.end(); ++it){
+ GEdge *ge = *it;
+ for (int i=0;i<ge->lines.size();i++){
+ SPoint2 p[2];
+ reparamMeshEdgeOnFace ( ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),gf,p[0],p[1]);
+ MVertex *vd[2];
+ for (int j=0;j<2;j++){
+ MVertex *v = ge->lines[i]->getVertex(j);
+ std::map<MVertex*,MVertex*>::iterator itv = duplicates.find(v);
+ if (itv == duplicates.end()){
+ vd[j] = new MFaceVertex(v->x(),v->y(),v->z(),gf,p[j].x(),p[j].y());
+ duplicates[v] = vd[j];
+ gf->mesh_vertices.push_back(vd[j]);
+ }
+ else
+ vd[j] = itv->second;
+ }
+ newQuads.push_back(new MQuadrangle(ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),vd[1],vd[0]));
+ }
+ }
+ for (int i=0;i<gf->quadrangles.size();i++){
+ MQuadrangle *q = gf->quadrangles[i];
+ MVertex *vs[4];
+ for (int j=0;j<4;j++){
+ MVertex *v = q->getVertex(j);
+ std::map<MVertex*,MVertex*>::iterator itv = duplicates.find(v);
+ if (itv == duplicates.end()){
+ vs[j] = v;
+ }
+ else{
+ vs[j] = itv->second;
+ }
+ }
+ newQuads.push_back(new MQuadrangle(vs[0],vs[1],vs[2],vs[3]));
+ delete q;
+ }
+ for (int i=0;i<gf->triangles.size();i++){
+ MTriangle *t = gf->triangles[i];
+ MVertex *vs[3];
+ for (int j=0;j<3;j++){
+ MVertex *v = t->getVertex(j);
+ std::map<MVertex*,MVertex*>::iterator itv = duplicates.find(v);
+ if (itv == duplicates.end()){
+ vs[j] = v;
+ }
+ else{
+ vs[j] = itv->second;
+ }
+ }
+ newTris.push_back(new MTriangle(vs[0],vs[1],vs[2]));
+ delete t;
+ }
+
+ gf->triangles = newTris;
+ gf->quadrangles = newQuads;
+ }
+ }
+ */
+}
+
+void quadsToTriangles(GFace *gf, double minqual = -10000.){
+ std::vector<MQuadrangle*> qds;
for (int i=0;i<gf->quadrangles.size();i++){
MQuadrangle *q = gf->quadrangles[i];
- MTriangle *t1 = new MTriangle (q->getVertex(0),q->getVertex(1),q->getVertex(2));
- MTriangle *t2 = new MTriangle (q->getVertex(2),q->getVertex(3),q->getVertex(0));
- gf->triangles.push_back(t1);
- gf->triangles.push_back(t2);
- delete q;
+ if (q->etaShapeMeasure() < minqual){
+ MTriangle *t11 = new MTriangle (q->getVertex(0),q->getVertex(1),q->getVertex(2));
+ MTriangle *t12 = new MTriangle (q->getVertex(2),q->getVertex(3),q->getVertex(0));
+ MTriangle *t21 = new MTriangle (q->getVertex(1),q->getVertex(2),q->getVertex(3));
+ MTriangle *t22 = new MTriangle (q->getVertex(3),q->getVertex(0),q->getVertex(1));
+ double qual1 = std::min(t11->gammaShapeMeasure(),t12->gammaShapeMeasure());
+ double qual2 = std::min(t21->gammaShapeMeasure(),t22->gammaShapeMeasure());
+
+ double surf1 = surfaceFaceUV(t11,gf) + surfaceFaceUV(t12,gf);
+ double surf2 = surfaceFaceUV(t21,gf) + surfaceFaceUV(t22,gf);
+
+ int option = 0;
+ if (surf1 > surf2 * (1.+1.e-8))option = 2;
+ else if (surf2 > surf1 * (1.+1.e-8))option = 1;
+
+ if (option == 1 || (option == 0 && qual1 > qual2)){
+ gf->triangles.push_back(t11);
+ gf->triangles.push_back(t12);
+ delete t21; delete t22;
+ }
+ else {
+ gf->triangles.push_back(t21);
+ gf->triangles.push_back(t22);
+ delete t11; delete t12;
+ }
+ delete q;
+ }
+ else {
+ qds.push_back(q);
+ }
}
- gf->quadrangles.clear();
+ gf->quadrangles = qds;
}
void recombineIntoQuadsIterative(GFace *gf){
recombineIntoQuads(gf);
- // return;
+ quadsToTriangles(gf,0.03);
+ return;
int COUNT = 0;
while (1){
quadsToTriangles(gf);
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index 327d1d49ec1704d7970801914bf717f810824590..200e9c5ab538d1fe16eb159a7bb276335e9f3080 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -14,6 +14,7 @@
#include "STensor3.h"
class GFace;
+class GVertex;
class MVertex;
class edge_angle {
@@ -53,8 +54,9 @@ void buildVertexToTriangle(std::vector<MTriangle*> &, v2t_cont &adj);
void buildEdgeToTriangle(std::vector<MTriangle*> &, e2t_cont &adj);
void buildListOfEdgeAngle(e2t_cont adj, std::vector<edge_angle> &edges_detected,
std::vector<edge_angle> &edges_lonly);
-void laplaceSmoothing(GFace *gf);
+void laplaceSmoothing(GFace *gf, int niter=1);
void edgeSwappingLawson(GFace *gf);
+void addOneLayerOnContour(GFace *gf, GVertex *gv);
enum swapCriterion {SWCR_DEL, SWCR_QUAL, SWCR_NORM, SWCR_CLOSE};
enum splitCriterion {SPCR_CLOSE, SPCR_QUAL, SPCR_ALLWAYS};
diff --git a/benchmarks/2d/naca12_2d.geo b/benchmarks/2d/naca12_2d.geo
index d5d5e3a6b666680d8d67bc61b5582083686c44d4..70ba9b91ff4a8d860daabdb9c452731a6d502651 100644
--- a/benchmarks/2d/naca12_2d.geo
+++ b/benchmarks/2d/naca12_2d.geo
@@ -224,17 +224,6 @@ Line Loop(9) = {6,7,8,5};
Line Loop(10) = {2,3,4,1};
Plane Surface(11) = {9,10};
-Physical Surface(11)={11};
+//Physical Surface(11)={11};
//Point(9999) = {0.6,0,0,1};
-Field[1] = Attractor;
-Field[1].EdgesList = {1, 2, 3, 4};
-Field[1].NNodesByEdge = 2000;
-Field[2] = Laplacian;
-Field[2].Delta = 0.001;
-Field[2].IField = 1;
-Field[3] = MathEval;
-Field[3].F = "20/(F2*3.14)";
-//Background Field = 3;
-Field[4] = Gradient;
-Field[4].Delta = 0.001;
diff --git a/benchmarks/boolean/wikipedia.lua b/benchmarks/boolean/wikipedia.lua
index 1d7a612cd5cc8fca0329e8a61ba5efcac19dc9b7..186fe8bfc4d67de79041fdac79c3f607e2dccbe9 100644
--- a/benchmarks/boolean/wikipedia.lua
+++ b/benchmarks/boolean/wikipedia.lua
@@ -2,8 +2,8 @@
-- from http://en.wikipedia.org/wiki/Constructive_solid_geometry
R = 1.4;
-s = .4;
-t = 1.35;
+s = .7;
+t = 1.25;
myModel = GModel();
myModel:addBlock({-R,-R,-R},{R,R,R});
@@ -24,4 +24,5 @@ myModel2:computeUnion(myTool2,0);
myModel2:computeUnion(myTool3,0);
myModel:computeDifference(myModel2,0);
+
myModel:setAsCurrent();