diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index c55eacd260275a58e16de5969204c3dc9ebaa3f7..d616de0380df27c46f9d5d998e74cfbf2cd004fa 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -365,6 +365,16 @@ void Get_Options(int argc, char *argv[])
else
Msg::Fatal("Missing number");
}
+ else if(!strcmp(argv[i] + 1, "edgelmin")) {
+ i++;
+ if(argv[i] != NULL) {
+ CTX.mesh.tolerance_edge_length = atof(argv[i++]);
+ if( CTX.mesh.tolerance_edge_length <= 0.0)
+ Msg::Fatal("Tolerance for model edge length must be > 0 (here %g)", CTX.mesh.tolerance_edge_length);
+ }
+ else
+ Msg::Fatal("Missing number");
+ }
else if(!strcmp(argv[i] + 1, "epslc1d")) {
i++;
if(argv[i] != NULL) {
diff --git a/Common/Context.h b/Common/Context.h
index 3e4161a615c5dfdec903055efee20f8743d32a8d..76fe3a90d58e099f899299d869fe75df0cbcd049 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -158,7 +158,7 @@ class Context_T {
int optimize, optimize_netgen, refine_steps;
int quality_type, label_type;
double quality_inf, quality_sup, radius_inf, radius_sup;
- double scaling_factor, lc_factor, rand_factor, lc_integration_precision, lc_min, lc_max;
+ double scaling_factor, lc_factor, rand_factor, lc_integration_precision, lc_min, lc_max, tolerance_edge_length;
int lc_from_points, lc_from_curvature, lc_extend_from_boundary;
int dual, voronoi, draw_skin_only;
int light, light_two_side, light_lines;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 9f854006e7203b3ea4757f6365be6a818cdd7b62..db33263c66e1a6445c4fe8ae9cb1538780ef04d9 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -938,7 +938,7 @@ StringXNumber MeshOptions_Number[] = {
"Factor applied to all characteristic lengths" },
{ F|O, "CharacteristicLengthMin" , opt_mesh_lc_min, 0.0 ,
"Minimum characteristic length" },
- { F|O, "CharacteristicLengthMax" , opt_mesh_lc_max, 1.e22 ,
+ { F|O, "CharacteristicLengthMax" , opt_mesh_lc_max, 1.e22,
"Maximum characteristic length" },
{ F|O, "CharacteristicLengthFromCurvature" , opt_mesh_lc_from_curvature , 0. ,
"Compute characteristic lengths from curvature" },
@@ -1093,6 +1093,8 @@ StringXNumber MeshOptions_Number[] = {
"Display size of tangent vectors (in pixels)" },
{ F|O, "Tetrahedra" , opt_mesh_tetrahedra , 1. ,
"Display mesh tetrahedra?" },
+ { F|O, "ToleranceEdgeLength" , opt_mesh_tolerance_edge_length, 0.0,
+ "Skip a model edge in mesh generation if its length is less than user's defined tolerance" },
{ F|O, "Triangles" , opt_mesh_triangles , 1. ,
"Display mesh triangles?" },
diff --git a/Common/Options.cpp b/Common/Options.cpp
index e251ff772726ee55cf8be0230f3c59a345bc4deb..b6423408ed7ef0167a3974394a15cdaa06dc323c 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -4497,6 +4497,14 @@ double opt_mesh_lc_max(OPT_ARGS_NUM)
return CTX.mesh.lc_max;
}
+double opt_mesh_tolerance_edge_length(OPT_ARGS_NUM)
+{
+ if(action & GMSH_SET)
+ CTX.mesh.tolerance_edge_length = val;
+ return CTX.mesh.tolerance_edge_length;
+}
+
+
double opt_mesh_lc_from_curvature(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index d46301269a01616dfe2702230bf402b02d319643..e2b4408be53144862f2d57f357648416a9e69e5b 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -429,6 +429,7 @@ double opt_mesh_explode(OPT_ARGS_NUM);
double opt_mesh_scaling_factor(OPT_ARGS_NUM);
double opt_mesh_lc_min(OPT_ARGS_NUM);
double opt_mesh_lc_max(OPT_ARGS_NUM);
+double opt_mesh_tolerance_edge_length(OPT_ARGS_NUM);
double opt_mesh_lc_factor(OPT_ARGS_NUM);
double opt_mesh_lc_from_curvature(OPT_ARGS_NUM);
double opt_mesh_lc_from_points(OPT_ARGS_NUM);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index a6f2aef38072865375f45e607a64572ab39c78bf..20a5accb0ae39325b976b89ca360d54c0c71a5f4 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1109,7 +1109,6 @@ void mesh_options_ok_cb(CALLBACK_ARGS)
opt_mesh_light_two_side(0, GMSH_SET, WID->mesh_butt[18]->value());
opt_mesh_smooth_normals(0, GMSH_SET, WID->mesh_butt[19]->value());
opt_mesh_light_lines(0, GMSH_SET, WID->mesh_butt[20]->value());
-
opt_mesh_nb_smoothing(0, GMSH_SET, WID->mesh_value[0]->value());
opt_mesh_lc_factor(0, GMSH_SET, WID->mesh_value[2]->value());
opt_mesh_lc_min(0, GMSH_SET, WID->mesh_value[25]->value());
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 1d3bcc087f501a54a250001bc38ae191b82b4be5..78ba2ade13c6b3415a0a36e4674f867e8fc1b8a9 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -19,7 +19,7 @@
#endif
GEdge::GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1)
- : GEntity(model, tag), v0(_v0), v1(_v1)
+ : GEntity(model, tag), _tooSmall(false), v0(_v0), v1(_v1)
{
if(v0) v0->addEdge(this);
if(v1 && v1 != v0) v1->addEdge(this);
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index d9a01a4577a3d667db86a4ac833902f6fe10edab..5eff3e278394cfebc0d7ad9a08a8c7941567bb71 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -20,6 +20,7 @@ class ExtrudeParams;
class GEdge : public GEntity {
private:
double _length;
+ bool _tooSmall;
protected:
GVertex *v0, *v1;
@@ -99,7 +100,8 @@ class GEdge : public GEntity {
virtual double prescribedMeshSizeAtVertex() const { return meshAttributes.meshSize; }
// true if start == end and no more than 2 segments
- bool isMeshDegenerated() const{ return (v0 == v1 && mesh_vertices.size() < 2); }
+ void setTooSmall ( bool b) {_tooSmall = b;}
+ bool isMeshDegenerated() const{ return _tooSmall || (v0 == v1 && mesh_vertices.size() < 2); }
// number of types of elements
int getNumElementTypes() const { return 1; }
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index eacdb7a9602c3c61b3fa76f5d499062b27ded796..ec185692cc0d197de5d4048119ae2fd89d718301 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -201,6 +201,7 @@ SPoint2 gmshEdge::reparamOnFace(GFace *face, double epar,int dir) const
}
}
+
if(s->Typ == MSH_SURF_REGL){
Curve *C[4];
for(int i = 0; i < 4; i++)
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 08e101ea84f69139e2ae88c6f9798ec5e236a21a..95802691d3549d1da14cc9f0fd4378240cb06169 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -186,6 +186,53 @@ GPoint gmshFace::point(double par1, double par2) const
GPoint gmshFace::closestPoint(const SPoint3 & qp, const double initialGuess[2]) const
{
+
+ if (s->Typ == MSH_SURF_PLAN && !s->geometry){
+ double XP = qp.x();
+ double YP = qp.y();
+ double ZP = qp.z();
+ double a = meanPlane.a;
+ double b = meanPlane.b;
+ double c = meanPlane.c;
+ double d = meanPlane.d;
+ // X = P + H N, find y such that
+ // a X_1 + b X_2 + c X_3 + d = 0
+ // a ( XP + y XN) + b ( YP + y YN) + c ( ZP + y ZN) + d = 0
+ // H ( a XN + b YN + c ZN ) = - (a XP + b YP + c ZP + d)
+ const double H = -(a*XP + b*YP + c *ZP + d)/(a*a + b*b + c*c);
+ const double X = XP + H * a;
+ const double Y = YP + H * b;
+ const double Z = ZP + H * c;
+ // now compute parametric coordinates
+ double x,y,z,VX[3], VY[3];
+ getMeanPlaneData(VX, VY, x, y, z);
+ // XP = X + u VX + v VY
+ // We are sure to be on the plane
+ double M[3][2] = {{VX[0],VY[0]},{VX[1],VY[1]},{VX[2],VY[2]}};
+ double MN[2][2];
+ double B[3] = {XP-x,YP-y,ZP-z};
+ double BN[2],UV[2];
+ for (int i=0;i<2;i++){
+ BN[i] = 0;
+ for (int k=0;k<3;k++){
+ BN[i] += B[k] * M[k][i];
+ }
+ }
+ for (int i=0;i<2;i++){
+ for (int j=0;j<2;j++){
+ MN[i][j] = 0;
+ for (int k=0;k<3;k++){
+ MN[i][j] += M[k][i] * M[k][j];
+ }
+ }
+ }
+ sys2x2(MN,BN,UV);
+
+ // GPoint test = point (UV[0],UV[1]);
+ // printf("%g %g %g vs %g %g %g\n",XP,YP,ZP,test.x(),test.y(),test.z());
+ return GPoint(XP, YP, ZP, this, UV);
+ }
+
Vertex v;
double u[2] = {initialGuess[0],initialGuess[1]};
v.Pos.X = qp.x();
@@ -193,7 +240,7 @@ GPoint gmshFace::closestPoint(const SPoint3 & qp, const double initialGuess[2])
v.Pos.Z = qp.z();
bool result = ProjectPointOnSurface(s, v, u);
if (!result)
- printf("Project Point on surface %d \n",result);
+ printf("Project Point on surface %d\n",result);
return GPoint(v.Pos.X, v.Pos.Y, v.Pos.Z, this, u);
}
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index b7e2f3ea4b197733ef369a4573e50953f0c38a25..ee68065662d849b8f05b7f43c557d0a822d84ded 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -163,7 +163,9 @@ double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double
if(lc <= 0.){
Msg::Error("Wrong characteristic length lc = %g", lc);
- lc = l1;
+ printf("%g %g \n", CTX.mesh.lc_min, CTX.mesh.lc_max);
+
+ lc = l1;
}
return lc * CTX.mesh.lc_factor;
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 4e382bc428819a197928f442183a3773f3e6e2e2..3cfcb07d1d89ff2e8a69c3d20e9030e8970d0495 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -25,6 +25,107 @@
extern Context_T CTX;
+
+static MVertex* isEquivalentTo ( std::multimap<MVertex *, MVertex *> & m, MVertex *v )
+{
+ std::multimap<MVertex *, MVertex *>::iterator it = m.lower_bound(v);
+ std::multimap<MVertex *, MVertex *>::iterator ite = m.upper_bound(v);
+ if (it == ite) return v;
+ MVertex *res = it->second; ++it;
+ while (it !=ite){
+ res = std::min(res,it->second);++it;
+ }
+ if (res < v) return isEquivalentTo ( m, res) ;
+ return res;
+}
+
+static void buildASetOfEquivalentMeshVertices ( GFace *gf , std::multimap<MVertex *, MVertex *> & equivalent , std::map<GVertex*,MVertex*> &bm)
+{
+ // an edge is degenerated when is length is considered to be
+ // zero. In some cases, a model edge can be considered as too
+ // small an is ignored.
+
+ // for taking that into account, we loop over the edges
+ // and create pairs of MVertices that are considered as
+ // equal.
+
+ std::list<GEdge*> edges = gf->edges();
+ std::list<GEdge*> emb_edges = gf->embeddedEdges();
+ std::list<GEdge*>::iterator it = edges.begin();
+
+ while(it != edges.end()){
+ if((*it)->isMeshDegenerated()){
+ MVertex *va = *((*it)->getBeginVertex()->mesh_vertices.begin());
+ MVertex *vb = *((*it)->getEndVertex()->mesh_vertices.begin());
+ if (va != vb){
+ equivalent.insert(std::make_pair (va,vb));
+ equivalent.insert(std::make_pair (vb,va));
+ bm[(*it)->getBeginVertex()] = va;
+ bm[(*it)->getEndVertex()] = vb;
+ printf("%d equivalent to %d\n",va->getNum(),vb->getNum());
+
+ }
+ }
+ ++it;
+ }
+
+ it = emb_edges.begin();
+ while(it != emb_edges.end()){
+ if((*it)->isMeshDegenerated()){
+ MVertex *va = *((*it)->getBeginVertex()->mesh_vertices.begin());
+ MVertex *vb = *((*it)->getEndVertex()->mesh_vertices.begin());
+ if (va != vb){
+ equivalent.insert(std::make_pair (va,vb));
+ equivalent.insert(std::make_pair (vb,va));
+ bm[(*it)->getBeginVertex()] = va;
+ bm[(*it)->getEndVertex()] = vb;
+ }
+ }
+ ++it;
+ }
+}
+
+struct geomTresholdVertexEquivalence
+{
+ // Initial MVertex associated to one given MVertex
+ std::map<GVertex *, MVertex *> backward_map;
+ // initiate the forward and backward maps
+ geomTresholdVertexEquivalence (GModel *g);
+ // restores the initial state
+ ~geomTresholdVertexEquivalence ();
+};
+
+
+geomTresholdVertexEquivalence :: geomTresholdVertexEquivalence (GModel *g)
+{
+ std::multimap<MVertex *, MVertex *> equivalenceMap;
+ for (GModel::fiter it = g->firstFace(); it != g->lastFace(); ++it)
+ buildASetOfEquivalentMeshVertices ( *it, equivalenceMap, backward_map );
+ // build the structure that identifiate geometrically equivalent
+ // mesh vertices.
+ for (std::map<GVertex*,MVertex *>::iterator it = backward_map.begin() ; it != backward_map.end() ; ++it){
+ GVertex *g = it->first;
+ MVertex *v = it->second;
+ MVertex *other = isEquivalentTo (equivalenceMap,v);
+ if (v != other){
+ printf("Finally : %d equivalent to %d\n",v->getNum(),other->getNum());
+ g->mesh_vertices.clear();
+ g->mesh_vertices.push_back(other);
+ }
+ }
+}
+
+geomTresholdVertexEquivalence :: ~geomTresholdVertexEquivalence ()
+{
+ // restore the initial data
+ for (std::map<GVertex*,MVertex *>::iterator it = backward_map.begin() ; it != backward_map.end() ; ++it){
+ GVertex *g = it->first;
+ MVertex *v = it->second;
+ g->mesh_vertices.clear();
+ g->mesh_vertices.push_back(v);
+ }
+}
+
template<class T>
static void GetQualityMeasure(std::vector<T*> &ele,
double &gamma, double &gammaMin, double &gammaMax,
@@ -282,6 +383,9 @@ static void Mesh2D(GModel *m)
Msg::StatusBar(1, true, "Meshing 2D...");
double t1 = Cpu();
+ // skip short mesh edges
+ geomTresholdVertexEquivalence inst (m);
+
// boundary layers are special: their generation (including vertices
// and curve meshes) is global as it depends on a smooth normal
// field generated from the surface mesh of the source surfaces
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 4e5a3c78da071685e959fa65e06c150cb279cf9f..638601d5242bd007c403e28780ba801239288705 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -436,24 +436,31 @@ void getFaceVertices(GFace *gf,
MVertex *v;
const double t1 = points(k, 0);
const double t2 = points(k, 1);
- if(!reparamOK || linear || gf->geomType() == GEntity::DiscreteSurface){
+ if(linear || gf->geomType() == GEntity::DiscreteSurface){
SPoint3 pc = face.interpolate(t1, t2);
v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
}
else{
double X(0),Y(0),Z(0),GUESS[2]={0,0};
+
for (int j=0; j<incomplete->getNumVertices(); j++){
double sf ; incomplete->getShapeFunction(j,t1,t2,0,sf);
MVertex *vt = incomplete->getVertex(j);
X += sf * vt->x();
Y += sf * vt->y();
Z += sf * vt->z();
- GUESS[0] += sf * pts[j][0];
- GUESS[1] += sf * pts[j][1];
+ if (reparamOK){
+ GUESS[0] += sf * pts[j][0];
+ GUESS[1] += sf * pts[j][1];
+ }
+ }
+ if (reparamOK){
+ GPoint gp = gf->closestPoint(SPoint3(X,Y,Z),GUESS);
+ v = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
+ }
+ else{
+ v = new MVertex(X,Y,Z, gf);
}
- GPoint gp = gf->closestPoint(SPoint3(X,Y,Z),GUESS);
- // printf("%g %g %g -- %g %g %g\n",X,Y,Z,gp.x(),gp.y(),gp.z());
- v = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
}
faceVertices[p].push_back(v);
gf->mesh_vertices.push_back(v);
@@ -718,13 +725,13 @@ void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
ve, nPts + 1));
}
else{
- if (gf->geomType() == GEntity::Plane){
+ if (0 && gf->geomType() == GEntity::Plane){
getFaceVertices(gf, t, vf, faceVertices, linear, nPts);
}
else{
- MTriangleN incpl (t->getVertex(0), t->getVertex(1), t->getVertex(2),
- ve, nPts + 1);
- getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts);
+ MTriangleN incpl (t->getVertex(0), t->getVertex(1), t->getVertex(2),
+ ve, nPts + 1);
+ getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts);
}
ve.insert(ve.end(), vf.begin(), vf.end());
triangles2.push_back
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 67f7ca73b23ef96fd2cd7d9b9f375797966b117c..6302508a5aabf07f72f15f9738662200b79a7e90 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -292,6 +292,10 @@ void meshGEdge::operator() (GEdge *ge)
if(length == 0.0)
Msg::Debug("Curve %d has a zero length", ge->tag());
+ // TEST
+ if (length < CTX.mesh.tolerance_edge_length)ge->setTooSmall(true);
+
+
// Integrate detJ/lc du
double a;
int N;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 985c21489a52f60a080439b595b5667499b507cf..c50a9662369ea965c905630a82c261abc29f7bed 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -316,6 +316,7 @@ static bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
return true;
}
+
// Builds An initial triangular mesh that respects the boundaries of
// the domain, including embedded points and surfaces
@@ -346,18 +347,21 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
it = emb_edges.begin();
while(it != emb_edges.end()){
- all_vertices.insert((*it)->mesh_vertices.begin(),
- (*it)->mesh_vertices.end() );
- all_vertices.insert((*it)->getBeginVertex()->mesh_vertices.begin(),
- (*it)->getBeginVertex()->mesh_vertices.end());
- all_vertices.insert((*it)->getEndVertex()->mesh_vertices.begin(),
- (*it)->getEndVertex()->mesh_vertices.end());
+ if(!(*it)->isMeshDegenerated()){
+ all_vertices.insert((*it)->mesh_vertices.begin(),
+ (*it)->mesh_vertices.end() );
+ all_vertices.insert((*it)->getBeginVertex()->mesh_vertices.begin(),
+ (*it)->getBeginVertex()->mesh_vertices.end());
+ all_vertices.insert((*it)->getEndVertex()->mesh_vertices.begin(),
+ (*it)->getEndVertex()->mesh_vertices.end());
+ }
++it;
}
if (all_vertices.size() < 3){
- Msg::Warning("Cannot triangulate less than 3 vertices");
- return false;
+ Msg::Warning("Mesh Generation of Model Face %d Skipped : Only %d Mesh Vertices on The Contours",gf->tag(),all_vertices.size());
+ gf->meshStatistics.status = GFace::DONE;
+ return true;
}
double *U_ = new double[all_vertices.size()];
@@ -506,7 +510,8 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
}
it = emb_edges.begin();
while(it != emb_edges.end()){
- recover_medge(m, *it, &edgesToRecover, &edgesNotRecovered, 1);
+ if(!(*it)->isMeshDegenerated())
+ recover_medge(m, *it, &edgesToRecover, &edgesNotRecovered, 1);
++it;
}
@@ -581,7 +586,8 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
it = emb_edges.begin();
while(it != emb_edges.end()){
- recover_medge(m, *it, &edgesToRecover, &edgesNotRecovered, 2);
+ if(!(*it)->isMeshDegenerated())
+ recover_medge(m, *it, &edgesToRecover, &edgesNotRecovered, 2);
++it;
}
// compute characteristic lengths at vertices
@@ -1286,6 +1292,7 @@ const int debugSurface = -1;
void meshGFace::operator() (GFace *gf)
{
+
gf->model()->setCurrentMeshEntity(gf);
if (debugSurface >= 0 && gf->tag() != debugSurface){
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index e88a2140bc17cefffb5d5be8b50ba96772e2f9fc..68f649ed82124068b0d3fc5b78c3a038f9550a1b 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -153,8 +153,29 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
int Count = 0;
for(int j = 0; j < 3; j++){
if(v[j]->onWhat()->dim() == 2){
- v[j]->getParameter(0,guess[0]);
- v[j]->getParameter(1,guess[1]);
+ double uu,vv;
+ v[j]->getParameter(0,uu);
+ v[j]->getParameter(1,vv);
+ guess[0] += uu;
+ guess[1] += vv;
+ Count++;
+ }
+ else if (v[j]->onWhat()->dim() == 1){
+ GEdge *ge = (GEdge*)v[j]->onWhat();
+ double UU;
+ v[j]->getParameter(0, UU);
+ SPoint2 param;
+ param = ge->reparamOnFace(gf, UU, 1);
+ guess[0]+=param.x();
+ guess[1]+=param.y();
+ Count++;
+ }
+ else if (v[j]->onWhat()->dim() == 0){
+ SPoint2 param;
+ GVertex *gv = (GVertex*)v[j]->onWhat();
+ param = gv->reparamOnFace(gf,1);
+ guess[0]+=param.x();
+ guess[1]+=param.y();
Count++;
}
}
@@ -171,14 +192,16 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
// PARAMETRIC COORDINATES SHOULD BE SET for the vertex !!!!!!!!!!!!!
// This is not 100 % safe yet, so we reserve that operation for high order
// meshes.
+ Msg::Debug("A new point has been inserted in mesh face %d by the 3D mesher, guess %g %g",gf->tag(),guess[0],guess[1]);
GPoint gp = gf->closestPoint (SPoint3(v[j]->x(), v[j]->y(), v[j]->z()),guess);
-
- Msg::Debug("A new point has been inserted in mesh face %d by the 3D mesher",gf->tag());
- Msg::Debug("The point has been projected back to the surface (%g %g %g) -> (%g %g %g)",
- v[j]->x(), v[j]->y(), v[j]->z(),gp.x(),gp.y(),gp.z());
+ Msg::Debug("The point has been projected back to the surface (%g %g %g) -> (%g %g %g), (%g,%g)",
+ v[j]->x(), v[j]->y(), v[j]->z(),gp.x(),gp.y(),gp.z(),gp.u(),gp.v());
// To be safe, we should ensure that this mesh motion does not lead to an invalid mesh !!!!
+ //v1b = new MVertex(gp.x(),gp.y(),gp.z(),gf);
v1b = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
+ //v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(),gf);
+
}
else{
v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(),gf);