diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 6090bd7793181ee9bc6dd2fabc9baebc84fd4b8a..a724c9beb94cd85c90619a5008b74550aaa92646 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -255,14 +255,15 @@ GPoint GEdge::closestPoint(const SPoint3 &q, double &t) const
return point(t);
}
-double GEdge::parFromPoint(const SVector3 &Q) const
+double GEdge::parFromPoint(const SPoint3 &P) const
{
double t;
- bool success = XYZToU(Q, t);
+ XYZToU(P.x(), P.y(), P.z(), t);
return t;
}
-bool GEdge::XYZToU(const SVector3 &Q, double &u, const double relax) const
+bool GEdge::XYZToU(const double X, const double Y, const double Z,
+ double &u, const double relax) const
{
const double Precision = 1.e-8;
const int MaxIter = 25;
@@ -271,14 +272,13 @@ bool GEdge::XYZToU(const SVector3 &Q, double &u, const double relax) const
double err, err2;
int iter;
-
Range<double> uu = parBounds(0);
double uMin = uu.low();
double uMax = uu.high();
printf("dans GEdge uMin=%g, uMax=%g \n", uMin, uMax);
- SVector3 P;
+ SVector3 Q(X, Y, Z), P;
double init[NumInitGuess];
@@ -291,8 +291,6 @@ bool GEdge::XYZToU(const SVector3 &Q, double &u, const double relax) const
err = 1.0;
iter = 1;
-
-
SVector3 dPQ = P - Q;
err2 = dPQ.norm();
@@ -316,7 +314,7 @@ bool GEdge::XYZToU(const SVector3 &Q, double &u, const double relax) const
if(relax > 1.e-2) {
Msg::Info("point %g %g %g on edge %d : Relaxation factor = %g",
Q.x(), Q.y(), Q.z(), 0.75 * relax);
- return XYZToU(Q, u, 0.75 * relax);
+ return XYZToU(Q.x(), Q.y(), Q.z(), u, 0.75 * relax);
}
Msg::Error("Could not converge reparametrisation of point (%e,%e,%e) on edge %d",
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index e94d5707aad8d36433bdf09f029ec931772f10bf..d2fa6720a508192ab2032ea2fa9e18fbdbddb06a 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -38,7 +38,6 @@ class GEdge : public GEntity {
GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1);
virtual ~GEdge();
-
// delete mesh data
virtual void deleteMesh();
@@ -78,10 +77,6 @@ class GEdge : public GEntity {
return SVector3(gp.x(), gp.y(), gp.z());
}
- // return the parmater location on the edge given a point in space
- // that is on the edge
- //virtual double parFromPoint(const SPoint3 &) const;
-
// get first derivative of edge at the given parameter
virtual SVector3 firstDer(double par) const = 0;
@@ -143,27 +138,30 @@ class GEdge : public GEntity {
virtual bool periodic(int dim) const { return v0 == v1; }
// true if edge is used in hyperbolic layer on face gf
- virtual bool inHyperbolicLayer(GFace* gf)
+ virtual bool inHyperbolicLayer(GFace *gf)
{
return bl_faces.find(gf) != bl_faces.end();
}
- virtual void flagInHyperbolicLayer(GFace* gf) {bl_faces.insert(gf);}
+ virtual void flagInHyperbolicLayer(GFace *gf) { bl_faces.insert(gf); }
// get bounds of parametric coordinate
virtual Range<double> parBounds(int i) const = 0;
// return the point on the face closest to the given point
- virtual GPoint closestPoint(const SPoint3 & queryPoint,double& param) const;
+ virtual GPoint closestPoint(const SPoint3 &queryPoint, double ¶m) const;
- // try to reparametrise the point on the edge
- virtual double parFromPoint(const SVector3& Q) const;
+ // return the parmater location on the edge given a point in space
+ // that is on the edge
+ virtual double parFromPoint(const SPoint3 &P) const;
- virtual bool XYZToU(const SVector3& Q,double& t,const double relax=0.5) const;
+ // compute the parameter U from a point XYZ
+ virtual bool XYZToU(const double X, const double Y, const double Z,
+ double &U, const double relax=0.5) const;
// compound
- void setCompound (GEdgeCompound *gec) {compound = gec;}
- GEdgeCompound *getCompound () const {return compound;}
+ void setCompound(GEdgeCompound *gec) { compound = gec; }
+ GEdgeCompound *getCompound() const { return compound; }
struct {
char Method;
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 0f866e70f7390a0b723cb6d8eecc18da487f04b9..a98014f468ff7d4cc37e98911e8a792f359e38f4 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -436,8 +436,8 @@ void GFaceCompound::parametrize(iterationStep step) const
for (std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
double value = myAssembler.getDofValue(v,0,1);
- if (step == 1)
- //printf("%p node =%g %g %g , value = %g of size %d \n",v, v->x(), v->y(), v->z(), value, coordinates.size());
+ //if (step == 1)
+ //printf("%p node =%g %g %g , value = %g of size %d \n",v, v->x(), v->y(), v->z(), value, coordinates.size());
std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
if (itf == coordinates.end()){
SPoint3 p(0,0,0);
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index c178903607730177677e83c2c8dca768b1e54dde..827db54120800b47b5672eaa8164a3b7ec1ffcef 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -207,8 +207,6 @@ static double Integration(GEdge *ge, double t1, double t2,
double (*f) (GEdge *e, double X),
std::vector<IntPoint> &Points, double Prec)
{
-
-
IntPoint from, to;
int depth = 0;
@@ -236,12 +234,9 @@ void deMeshGEdge::operator() (GEdge *ge)
void meshGEdge::operator() (GEdge *ge)
{
-
ge->model()->setCurrentMeshEntity(ge);
- if(ge->geomType() == GEntity::DiscreteCurve) {
- return;
- }
+ if(ge->geomType() == GEntity::DiscreteCurve) return;
if(ge->geomType() == GEntity::BoundaryLayerCurve) return;
if(ge->meshAttributes.Method == MESH_NONE) return;
if(CTX::instance()->mesh.meshOnlyVisible && !ge->getVisibility()) return;
@@ -251,7 +246,7 @@ void meshGEdge::operator() (GEdge *ge)
if(MeshExtrudedCurve(ge)) return;
- Msg::Info("** Meshing curve %d (%s)", ge->tag(), ge->getTypeString().c_str());
+ Msg::Info("Meshing curve %d (%s)", ge->tag(), ge->getTypeString().c_str());
// compute bounds
Range<double> bounds = ge->parBounds(0);
@@ -336,7 +331,6 @@ void meshGEdge::operator() (GEdge *ge)
double lc = d/(P1.lc + dlc / dp * (d - P1.p));
GPoint V = ge->point(t);
ge->mesh_vertices[NUMP - 1] = new MEdgeVertex(V.x(), V.y(), V.z(), ge, t, lc);
- // printf("lc = %12.5E %12.5E \n",lc,P1.lc,P2.lc);
NUMP++;
}
else {
@@ -352,7 +346,6 @@ void meshGEdge::operator() (GEdge *ge)
MVertex *v1 = (i == ge->mesh_vertices.size()) ?
ge->getEndVertex()->mesh_vertices[0] : ge->mesh_vertices[i];
ge->lines.push_back(new MLine(v0, v1));
- // printf("New Line v0=%g v1=%g \n", v0->y(), v1->y());
}
if(ge->getBeginVertex() == ge->getEndVertex() &&
@@ -362,14 +355,4 @@ void meshGEdge::operator() (GEdge *ge)
v0->y() = beg_p.y();
v0->z() = beg_p.z();
}
-
- //for (std::vector<MLine*>::iterator it = ge->lines.begin() ; it != ge->lines.end() ; ++it){
- // printf("line with : first = %d last=%d \n", (*it)->getVertex(0)->getIndex(), (*it)->getVertex(1)->getIndex() );
- //}
-
-// //if DiscreteCurve, erase all old MLines
-// if(ge->geomType() == GEntity::DiscreteCurve) {
-// ge->lines.erase(ge->lines.begin(), ge->lines.end()-(N-1));
-// }
-
}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index dc408e069f1ec5772251851b9d5ce7e71d1516f1..f1de893a8d70a8ef8002280e019567122650b2f1 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -378,7 +378,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
std::list<GEdge*> edges = gf->edges();
// here, we will replace edges by their compounds
- printf("***** In meshGFace: \n");
+ //printf("***** In meshGFace: \n");
if (gf->geomType() == GEntity::CompoundSurface){
printf("replace edges by compound lines \n");
std::set<GEdge*> mySet;
@@ -448,7 +448,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
double *V_ = new double[all_vertices.size()];
v_container::iterator itv = all_vertices.begin();
- printf("boundary vertices size = %d \n", all_vertices.size());
+ //printf("boundary vertices size = %d \n", all_vertices.size());
int count = 0;
SBoundingBox3d bbox;
@@ -1089,26 +1089,28 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
it != gf->edgeLoops.end(); it++){
std::vector<BDS_Point* > edgeLoop_BDS;
int nbPointsLocal;
- if(!buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m,
- recover_map, nbPointsLocal, nbPointsTotal,
- 1.e-7*LC2D))
- if(!buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m,
- recover_map, nbPointsLocal, nbPointsTotal,
- 1.e-7 * LC2D, true))
- if(!buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m,
- recover_map, nbPointsLocal, nbPointsTotal,
- 1.e-5 * LC2D))
- if(!buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m,
- recover_map , nbPointsLocal, nbPointsTotal,
- 1.e-5 * LC2D, true))
- if(!buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m,
- recover_map , nbPointsLocal, nbPointsTotal,
- 1.e-3 * LC2D)){
- gf->meshStatistics.status = GFace::FAILED;
- Msg::Error("The 1D Mesh seems not to be forming a closed loop");
- m->scalingU = m->scalingV = 1.0;
- return false;
- }
+ const double fact[4] = {1.e-12, 1.e-7, 1.e-5, 1.e-3};
+ bool ok = false;
+ for(int i = 0; i < 4; i++){
+ if(buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m,
+ recover_map, nbPointsLocal,
+ nbPointsTotal, fact[i] * LC2D)){
+ ok = true;
+ break;
+ }
+ if(buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m,
+ recover_map, nbPointsLocal,
+ nbPointsTotal, fact[i] * LC2D, true)){
+ ok = true;
+ break;
+ }
+ }
+ if(!ok){
+ gf->meshStatistics.status = GFace::FAILED;
+ Msg::Error("The 1D Mesh seems not to be forming a closed loop");
+ m->scalingU = m->scalingV = 1.0;
+ return false;
+ }
nbPointsTotal += nbPointsLocal;
edgeLoops_BDS.push_back(edgeLoop_BDS);
}
diff --git a/benchmarks/misc/adaptive_tet10.msh b/benchmarks/misc/adaptive_tet10.msh
index 570ebbe6ce924bff345d5bf10283693954876196..a08a6fe6dc47e05557b5f51339a3a1f7dc4afc3a 100644
--- a/benchmarks/misc/adaptive_tet10.msh
+++ b/benchmarks/misc/adaptive_tet10.msh
@@ -109,3 +109,40 @@ $NodeData
26 1
27 1
$EndNodeData
+$NodeData
+1
+"2nd order view"
+1
+1.0
+3
+1
+1
+27
+1 12
+2 12
+3 12
+4 12
+5 2
+6 2
+7 2
+8 2
+9 1
+10 1
+11 1
+12 1
+13 1
+14 1
+15 1
+16 1
+17 1
+18 1
+19 1
+20 1
+21 1
+22 1
+23 1
+24 1
+25 1
+26 1
+27 1
+$EndNodeData
diff --git a/doc/TODO.txt b/doc/TODO.txt
index b753e1149e75fcb0061ab6f8e56dc70ce0495f00..4d1c7e23afb19d2d54729d84f618cf3690372d93 100644
--- a/doc/TODO.txt
+++ b/doc/TODO.txt
@@ -1,4 +1,8 @@
-$Id: TODO.txt,v 1.27 2009-05-01 06:37:04 geuzaine Exp $
+$Id: TODO.txt,v 1.28 2009-05-09 07:29:16 geuzaine Exp $
+
+********************************************************************
+
+Add Reverse Cuthill-McKee to reorder the nodes?
********************************************************************