diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 6013a87bf410620bd35c165aff7b4bf4098a73c4..4daed944bb49c67af9e2787dd588a2f43b26d6af 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -170,4 +170,11 @@ bool reparamMeshVertexOnEdge(MVertex *v, const GEdge *ge, double ¶m);
double angle3Vertices(MVertex *p1, MVertex *p2, MVertex *p3);
+inline double distance (MVertex *v1, MVertex *v2){
+ const double dx = v1->x() - v2->x();
+ const double dy = v1->y() - v2->y();
+ const double dz = v1->z() - v2->z();
+ return sqrt(dx*dx+dy*dy+dz*dz);
+}
+
#endif
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 88ecbba4c8757efd93a3464ac418c50d77da751f..fac7b4106d081f0b628455d6ef1fc14e67e42e75 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -384,14 +384,15 @@ void meshGEdge::operator() (GEdge *ge)
a = smoothPrimitive(ge, sqrt(CTX::instance()->mesh.smoothRatio), Points);
N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.));
}
-
+
// force odd number of points for if blossom is used for recombination
if(ge->meshAttributes.Method != MESH_TRANSFINITE &&
- CTX::instance()->mesh.algoRecombine == 1 && N % 2 == 0){
- if(CTX::instance()->mesh.recombineAll){
+ (CTX::instance()->mesh.algoRecombine == 1 || CTX::instance()->mesh.recombineAll) && N % 2 == 0){
+ // if(CTX::instance()->mesh.recombineAll){
N++;
- }
- else{
+ // }
+ }
+ /* else{
std::list<GFace*> faces = ge->faces();
for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
if((*it)->meshAttributes.recombine){
@@ -401,6 +402,7 @@ void meshGEdge::operator() (GEdge *ge)
}
}
}
+ */
// printFandPrimitive(ge->tag(),Points);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 888d1183559de13ef29bad58ea1ac1bd6897bd1b..db458fe917a77b7b1b1160528555c5a19f426465 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -574,19 +574,16 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
MTri3 *t4;
t4 = new MTri3(t, LL,0,&Us,&Vs,gf);
- double d1 = sqrt((it->v[0]->x() - v->x()) * (it->v[0]->x() - v->x()) +
- (it->v[0]->y() - v->y()) * (it->v[0]->y() - v->y()) +
- (it->v[0]->z() - v->z()) * (it->v[0]->z() - v->z()));
- double d2 = sqrt((it->v[1]->x() - v->x()) * (it->v[1]->x() - v->x()) +
- (it->v[1]->y() - v->y()) * (it->v[1]->y() - v->y()) +
- (it->v[1]->z() - v->z()) * (it->v[1]->z() - v->z()));
- const double MID[3] = {0.5*(it->v[0]->x()+it->v[1]->x()),
- 0.5*(it->v[0]->y()+it->v[1]->y()),
- 0.5*(it->v[0]->z()+it->v[1]->z())};
- double d3 = sqrt((MID[0] - v->x()) * (MID[0] - v->x()) +
- (MID[1] - v->y()) * (MID[1] - v->y()) +
- (MID[2] - v->z()) * (MID[2] - v->z()));
- if ((d1 < LL * .25 || d2 < LL * .25 || d3 < LL * .25) && !force) {
+ // double din = t->getInnerRadius();
+
+ double d1 = distance(it->v[0],v);
+ double d2 = distance(it->v[1],v);
+ double d3 = distance(it->v[0],it->v[1]);
+
+ // avoid angles that are too obtuse
+ double cosv = ((d1*d1+d2*d2-d3*d3)/(2.*d1*d2));
+
+ if ((d1 < LL * .25 || d2 < LL * .25 || cosv < -0.5) && !force) {
onePointIsTooClose = true;
}
@@ -746,35 +743,8 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
MTri3 *ptin = 0;
bool inside = false;
double uv[2];
- // FIXME !!! ----> FIXED (JFR)
- if (0 && MTri3::radiusNorm == 2){
- inside = invMapUV(worst->tri(), center, Us, Vs, uv, 1.e-8);
- if (inside)ptin = worst;
- if (!inside && worst->getNeigh(0)){
- inside |= invMapUV(worst->getNeigh(0)->tri(), center, Us, Vs, uv, 1.e-8);
- if (inside)ptin = worst->getNeigh(0);
- }
- if (!inside && worst->getNeigh(1)){
- inside |= invMapUV(worst->getNeigh(1)->tri(), center, Us, Vs, uv, 1.e-8);
- if (inside)ptin = worst->getNeigh(1);
- }
- if (!inside && worst->getNeigh(2)){
- inside |= invMapUV(worst->getNeigh(2)->tri(), center, Us, Vs, uv, 1.e-8);
- if (inside)ptin = worst->getNeigh(2);
- }
- if (!inside){
- ptin = search4Triangle (worst, center, Us, Vs, AllTris,uv);
- if (ptin)inside = true;
- }
- }
- else {
- ptin = search4Triangle (worst, center, Us, Vs, AllTris,uv);
- // if (!ptin){
- // printf("strange : %g %g seems to be out of the domain for face %d\n",
- // center[0],center[1],gf->tag());
- // }
- if (ptin) inside = true;
- }
+ ptin = search4Triangle (worst, center, Us, Vs, AllTris,uv);
+ if (ptin) inside = true;
if (inside) {
// we use here local coordinates as real coordinates
@@ -1025,24 +995,6 @@ double optimalPointFrontal (GFace *gf,
2 * dir[1] * dir[0] * metric[1] +
dir[1] * dir[1] * metric[2]);
- // const double p = 0.5 * lengthMetric(P, Q, metric); // / RATIO;
- /*
- const double rhoM1 = 0.5 *
- (vSizes[base->getVertex(ip1)->getIndex()] +
- vSizes[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
- const double rhoM2 = 0.5 *
- (vSizesBGM[base->getVertex(ip1)->getIndex()] +
- vSizesBGM[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
- const double rhoM = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
- */
-
- // const double rhoM_hat = std::min(std::max(rhoM, p), (p * p + q * q) / (2 * q));
- // double d = (rhoM_hat + sqrt (rhoM_hat * rhoM_hat - p * p)) / RATIO;
-
- //const double d = 100./RATIO;
- // printf("(%g %g) (%g %g) %g %g %g %g %g %g\n",
- // P[0],P[1],Q[0],Q[1],RATIO,p,q,rhoM,rhoM_hat,d);
- // printf("size %12.5E\n",vSizesBGM[base->getVertex(ip1)->getIndex()]);
const double rhoM1 = 0.5*
(vSizes[base->getVertex(ip1)->getIndex()] +
vSizes[base->getVertex(ip2)->getIndex()] ) ;// * RATIO;
@@ -1053,17 +1005,22 @@ double optimalPointFrontal (GFace *gf,
const double rhoM_hat = rhoM;
const double q = lengthMetric(center, midpoint, metric);
-
const double d = rhoM_hat * sqrt(3.)*0.5;
+ // d is corrected in a way that the mesh size is computed at point newPoint
+
+
// printf("%12.5E %12.5E\n",d,RATIO);
// const double L = d ;
// avoid to go toooooo far
const double L = d > q ? q : d;
+
newPoint[0] = midpoint[0] + L * dir[0]/RATIO;
newPoint[1] = midpoint[1] + L * dir[1]/RATIO;
+
+
return L;// > q ? d : q;
}
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index a210a953556720e0ff4f6172c9c734d4f74f2a7d..f4291844988d6f688d19735e8ecf80320f444dc7 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -3039,6 +3039,8 @@ void recombineIntoQuads(GFace *gf,
}
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);
return;
}
diff --git a/Post/PViewAsSimpleFunction.cpp b/Post/PViewAsSimpleFunction.cpp
index c8102561a3e1b927ba6840c5b1f116e81d29346e..aeaabce2dbb3f0b7ab24538847c2431b57ebdaca 100644
--- a/Post/PViewAsSimpleFunction.cpp
+++ b/Post/PViewAsSimpleFunction.cpp
@@ -12,6 +12,7 @@ double PViewEvaluator::operator() (const double x, const double y, const double
PViewData * pvd = _pv->getData();
double value;
bool found = pvd->searchScalar(x, y, z, &value, _step);
+ printf("found %d %g %g %g %g\n",found,x,y,value,x*x+y*y);
if (found) return value;
return 1.e22;
}
diff --git a/wrappers/gmshpy/gmshCommon.i b/wrappers/gmshpy/gmshCommon.i
index 487740d838a601ef0be63f5098d380083a1a2fcf..b000422999f6cddc93d150298cdcb05d099abe27 100644
--- a/wrappers/gmshpy/gmshCommon.i
+++ b/wrappers/gmshpy/gmshCommon.i
@@ -39,7 +39,6 @@ namespace std {
%template(DoubleVector) std::vector<double, std::allocator<double> >;
%template(DoubleVectorVector) std::vector<std::vector<double, std::allocator<double> > >;
%template(StringVector) std::vector<std::string, std::allocator<std::string> >;
- %template(IntDoubleMap) std::map < int , double >;
}
%pointer_functions(double,doublep)
%pointer_functions(int,intp)
diff --git a/wrappers/gmshpy/gmshPost.i b/wrappers/gmshpy/gmshPost.i
index e5f0a5617452c6605656ac98675184574818320a..c704b32c75379c9c2fce67700c839ff4039cad88 100644
--- a/wrappers/gmshpy/gmshPost.i
+++ b/wrappers/gmshpy/gmshPost.i
@@ -2,6 +2,8 @@
%module gmshPost
%include typemaps.i
%include std_string.i
+%include std_vector.i
+%include std_map.i
%{
#include "GmshConfig.h"
@@ -14,11 +16,20 @@
#endif
%}
+namespace std {
+ %template(DoubleVector) std::vector<double, std::allocator<double> >;
+ %template(IntDoubleVectorMap) std::map < int , std::vector<double, std::allocator<double> > >;
+}
+
+
%include "GmshConfig.h"
#if defined(HAVE_POST)
%include "PView.h"
%include "PViewFactory.h"
%apply double &OUTPUT { double &val}
%include "PViewData.h"
+%include "simpleFunction.h"
+%template(simpleFunctionDouble) simpleFunction<double>;
+%include "PViewAsSimpleFunction.h"
#endif