diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index e1379d1a36eb6590f8d64e5490e323a228b01ddf..6090bd7793181ee9bc6dd2fabc9baebc84fd4b8a 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -271,10 +271,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;
double init[NumInitGuess];
@@ -288,6 +291,8 @@ bool GEdge::XYZToU(const SVector3 &Q, double &u, const double relax) const
err = 1.0;
iter = 1;
+
+
SVector3 dPQ = P - Q;
err2 = dPQ.norm();
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 3cdad0149751a441474ff49106bb10c401e7fdbb..02c3779ff2fb096fa9cbd6ba3100df52b9d7e2fb 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -74,6 +74,10 @@ class GEdge : public GEntity {
GPoint gp = point(p);
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;
diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index 1d6d45be698810a528d694ac602ccc616ab686c7..8ef28d87726faa6b0a474338ef1415060020979b 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -97,7 +97,7 @@ Range<double> GEdgeCompound::parBounds(int i) const {
/*
+--------+-----------+----------- ... ----+
- 0 _par[1] _par[2] _par[N-1]
+ 0 _pars[1] _pars[2] _pars[N-1]
*/
@@ -105,8 +105,9 @@ void GEdgeCompound::getLocalParameter ( const double &t,
int &iEdge,
double & tLoc) const
{
+
for (iEdge=0 ; iEdge<_compound.size() ;iEdge++){
- // printf("iEdge %d par %g\n",iEdge,_pars[iEdge]);
+ //printf("iEdge=%d tmin=%g\n",iEdge,_pars[iEdge]);
double tmin = _pars[iEdge];
double tmax = _pars[iEdge+1];
if (t >= tmin && t <= tmax){
@@ -114,7 +115,7 @@ void GEdgeCompound::getLocalParameter ( const double &t,
tLoc = _orientation[iEdge] ?
b.low() + (t-tmin)/(tmax-tmin) * (b.high()-b.low()) :
b.high() - (t-tmin)/(tmax-tmin) * (b.high()-b.low()) ;
- // printf("%g %g %d\n",t,tLoc,iEdge);
+ //printf("bhigh=%g, blow=%g, global t=%g , tLoc=%g ,iEdge=%d\n",b.high(), b.low(), t,tLoc,iEdge);
return;
}
}
@@ -127,6 +128,11 @@ void GEdgeCompound::parametrize()
Range<double> b = _compound[i]->parBounds(0);
_pars.push_back(_pars[_pars.size()-1]+(b.high() - b.low()));
}
+
+ // for (int i=0;i<_compound.size()+1;i++){
+ // printf("_pars[%d]=%g\n",i, _pars[i] );
+ //}
+
}
@@ -143,7 +149,7 @@ GPoint GEdgeCompound::point(double par) const
double tLoc;
int iEdge;
getLocalParameter(par,iEdge,tLoc);
- // printf("compound %d par %g edge %d loc %g\n",tag(),par,iEdge,tLoc);
+ //printf("compound %d par %g edge %d loc %g\n",tag(),par,iEdge,tLoc);
return _compound[iEdge]->point(tLoc);
}
diff --git a/Geo/GEdgeCompound.h b/Geo/GEdgeCompound.h
index 8171e3ac833ed796fce908d45e3f1e905441b69b..4f9b1e45489570a41bf77aad7612320bc8cadde2 100644
--- a/Geo/GEdgeCompound.h
+++ b/Geo/GEdgeCompound.h
@@ -20,11 +20,11 @@ class GEdgeCompound : public GEdge {
std::vector<double> _pars;
void parametrize() ;
void orderEdges() ;
+
+public:
void getLocalParameter ( const double &t,
int &iEdge,
double & tLoc) const;
-
-public:
GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound);
virtual ~GEdgeCompound();
Range<double> parBounds(int i) const;
@@ -36,7 +36,7 @@ public:
virtual double curvature(double t) const;
virtual int minimumMeshSegments() const;
virtual int minimumDrawSegments() const;
-
+ std::vector<GEdge*> getEdgesOfCompound() const { return _compound; }
};
#endif
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 6268d4a1b28644ce013710bb916876bd6356f562..54b7c170559b81f0cdf9b04905bddc5891c094db 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -16,11 +16,11 @@
class gmshGradientBasedDiffusivity : public gmshFunction<double>
{
- private:
+private:
MElement *_current;
int _iComp;
mutable std::map<MVertex*, SPoint2> _coordinates;
- public:
+public:
gmshGradientBasedDiffusivity (std::map<MVertex*,SPoint2> &coordinates)
: _current (0), _iComp(-1),_coordinates(coordinates){}
void setCurrent (MElement *current){ _current = current; }
@@ -47,10 +47,10 @@ class gmshGradientBasedDiffusivity : public gmshFunction<double>
class gmshDistanceBasedDiffusivity : public gmshFunction<double>
{
- private:
+private:
MElement *_current;
mutable std::map<MVertex*, SPoint3> _coordinates;
- public:
+public:
gmshDistanceBasedDiffusivity (std::map<MVertex*,SPoint3> &coordinates)
: _current (0),_coordinates(coordinates){}
void setCurrent (MElement *current){ _current = current; }
@@ -79,6 +79,7 @@ static void fixEdgeToValue(GEdge *ed, double value, gmshAssembler<double> &myAss
void GFaceCompound::parametrize() const
{
+
if (!oct){
coordinates.clear();
parametrize(ITERD);
@@ -152,6 +153,8 @@ GFaceCompound::~GFaceCompound()
}
}
+
+
static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
std::vector<double> &coord)
{
@@ -213,6 +216,99 @@ static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
return true;
}
+SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
+{
+ parametrize() ;
+ std::map<MVertex*,SPoint3>::iterator it = coordinates.find(v);
+ // return SPoint2(it->second.x(),it->second.y());
+
+ if(it != coordinates.end()){
+ return SPoint2(it->second.x(),it->second.y());
+ }
+ else{
+
+ double tGlob, tLoc;
+ double tL, tR;
+ int iEdge;
+
+ //getParameter Point
+ v->getParameter(0,tGlob);
+ //printf("global param for point (%g, %g, %g ) = %g\n", v->x(), v->y(), v->z(), tGlob);
+
+ //find compound Edge
+ GEdgeCompound *gec = dynamic_cast<GEdgeCompound*>(v->onWhat());
+ //printf("tag compound=%d, beginvertex=%d, endvertex=%d \n", gec->tag(), gec->getBeginVertex()->tag(), gec->getEndVertex()->tag());
+
+ if (gec){
+
+ //compute local parameter on Edge
+ gec->getLocalParameter(tGlob,iEdge,tLoc);
+ std::vector<GEdge*> gev = gec->getEdgesOfCompound();
+ GEdge *ge = gev[iEdge];
+ //printf("iEdge =%d, leftV =%d, rightV=%d, and local param=%g \n", ge->tag(), ge->getBeginVertex()->tag(), ge->getEndVertex()->tag(), tLoc);
+
+ //left and right vertex of the Edge
+ MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
+ MVertex *v1 = ge->getEndVertex()->mesh_vertices[0];
+ std::map<MVertex*,SPoint3>::iterator itL = coordinates.find(v0);
+ std::map<MVertex*,SPoint3>::iterator itR = coordinates.find(v1);
+
+ //for the Edge, find the left and right vertices of the initial 1D mesh and interpolate to find (u,v)
+ //printf("number of vertices =%d for Edge =%d \n", ge->mesh_vertices.size(), ge->tag());
+ MVertex *vL = v0;
+ MVertex *vR = v1;
+ double xL = vL->x();
+ double yL = vL->y();
+ double zL = vL->z();
+ tL = ge->parBounds(0).low();
+ tR = ge->parBounds(0).high();
+ //printf("tLeftEDGE=%g tRightEDGE=%g, tLoc=%g\n", tL, tR, tLoc);
+ int j = 0;
+ //printf("j =%d, veL (%g,%g,%g), -> uv= (%g,%g)\n",j,xL,yL,zL, itL->second.x(), itL->second.y());
+ bool found = false;
+ while (j < ge->mesh_vertices.size()){
+ vR = ge->mesh_vertices[j];
+ double xR=vR->x();
+ double yR=vR->y();
+ double zR=vR->z();
+ //printf("*** L=(%g,%g,%g) et R=(%g,%g,%g)\n",xL,yL,zL, vR->x(), vR->y(), vR->z());
+ //printf("conditions XL %g < %g, yL %g < %g, zL %g < %g \n", fabs(v->x()-xL) , fabs(xL-xR), fabs(v->y()-yL), fabs(yL-yR), fabs(v->z()-zL) , fabs (zL-zR));
+ //printf("conditions XR %g < %g, yR %g < %g, zR %g < %g \n", fabs(v->x()-vR->x()) , fabs(xL-xR),fabs(v->y()-vR->y()), fabs(yL-yR),fabs(v->z()-vR->z()) , fabs (zL-zR));
+ if (fabs(v->x()-xL) <= fabs(xL-xR) && fabs(v->y()-yL) <= fabs(yL-yR) && fabs(v->z()-zL) <= fabs (zL-zR) &&
+ fabs(v->x()-vR->x()) <= fabs(xL-xR) && fabs(v->y()-vR->y()) <= fabs(yL-yR) && fabs(v->z()-vR->z()) <= fabs(zL-zR)){
+ found = true;
+ itR = coordinates.find(vR);
+ vR->getParameter(0,tR);
+ break;
+ }
+ else{
+ xL = xR; yL = yR ; zL = zR;
+ itL = coordinates.find(vR);
+ vL = vR;
+ vL->getParameter(0,tL);
+ }
+ j++;
+ //printf("in while j =%d, vL (%g,%g,%g), -> uv= (%g,%g)\n",j, vL->x(), vL->y(), vL->z(), itL->second.x(), itL->second.y());
+ }
+ if (!found) vR=v1;
+ //printf("vL (%g,%g,%g), -> uv= (%g,%g)\n",vL->x(), vL->y(), vL->z(), itL->second.x(), itL->second.y());
+ //printf("vR (%g,%g,%g), -> uv= (%g,%g)\n",vR->x(), vR->y(), vR->z(), itR->second.x(), itR->second.y());
+ //printf("tL:%g, tR=%g, tLoc=%g \n", tL, tR, tLoc);
+
+ //Linear interpolation between tL and tR
+ double uloc, vloc;
+ uloc = itL->second.x() + (tLoc-tL)/(tR-tL) * (itR->second.x()-itL->second.x());
+ vloc = itL->second.y() + (tLoc-tL)/(tR-tL) * (itR->second.y()-itL->second.y());
+ //printf("uloc=%g, vloc=%g \n", uloc,vloc);
+
+ //printf("Vertex coordinates not found in compound face \n");
+ //exit(0);
+
+ return SPoint2(uloc,vloc);
+ }
+ }
+}
+
void GFaceCompound::parametrize(iterationStep step) const
{
Msg::Info("Parametrizing Surface %d coordinate (ITER %d)", tag(), step);
@@ -383,18 +479,18 @@ double GFaceCompound::curvature(const SPoint2 ¶m) const
return 0.0;
}
-// if (lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
-// SPoint2 pv1, pv2, pv3;
-// bool ok = reparamMeshVertexOnFace(lt->t->getVertex(0), lt->gf, pv1);
-// ok |= reparamMeshVertexOnFace(lt->t->getVertex(1), lt->gf, pv2);
-// ok |= reparamMeshVertexOnFace(lt->t->getVertex(2), lt->gf, pv3);
-// if (ok){
-// SPoint2 pv = pv1*(1.-U-V) + pv2*U + pv3*V;
-// return lt->gf->curvature(pv));
-// }
-// }
-
-// return curvature(lt->t);
+ // if (lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
+ // SPoint2 pv1, pv2, pv3;
+ // bool ok = reparamMeshVertexOnFace(lt->t->getVertex(0), lt->gf, pv1);
+ // ok |= reparamMeshVertexOnFace(lt->t->getVertex(1), lt->gf, pv2);
+ // ok |= reparamMeshVertexOnFace(lt->t->getVertex(2), lt->gf, pv3);
+ // if (ok){
+ // SPoint2 pv = pv1*(1.-U-V) + pv2*U + pv3*V;
+ // return lt->gf->curvature(pv));
+ // }
+ // }
+
+ // return curvature(lt->t);
return 0.;
}
@@ -418,7 +514,7 @@ GPoint GFaceCompound::point(double par1, double par2) const
getTriangle (par1, par2, <, U,V);
SPoint3 p(0, 0, 0);
if (!lt){
- Msg::Warning("Re-Parametrized face %d --> point (%g %g) lies outside the domain", tag(), par1,par2);
+ // Msg::Warning("Re-Parametrized face %d --> point (%g %g) lies outside the domain", tag(), par1,par2);
return GPoint(p.x(),p.y(),p.z(),this);
}
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 669fe7a202ea9aefa5131100265c9d45a892a9ba..44b3fa873429d493cfd95137d89e169e3f5bf024 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -8,6 +8,7 @@
#include "GFace.h"
#include "GEdge.h"
+#include "GEdgeCompound.h"
/*
A GFaceCompound is a model face that is the compound of model faces.
@@ -71,11 +72,9 @@ public:
virtual GEntity::GeomType geomType() const { return CompoundSurface; }
ModelType getNativeType() const { return GmshModel; }
void * getNativePtr() const { return 0; }
- SPoint2 getCoordinates(MVertex *v) const {
- parametrize() ;
- std::map<MVertex*,SPoint3>::iterator it = coordinates.find(v);
- return SPoint2(it->second.x(),it->second.y());
- }
+
+ virtual SPoint2 getCoordinates(MVertex *v) const;
+
virtual bool buildRepresentationCross(){ return false; }
virtual double curvature(const SPoint2 ¶m) const;
private:
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 55dd9cd3024f3250bb7f5810dcf86b5f682a5a5a..ce1ec9bd7f4477a1622bcd2b1c2b11ac79fc75f8 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -234,6 +234,7 @@ void deMeshGEdge::operator() (GEdge *ge)
void meshGEdge::operator() (GEdge *ge)
{
+
ge->model()->setCurrentMeshEntity(ge);
if(ge->geomType() == GEntity::DiscreteCurve) return;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 4d53fb9409fc43e0ff9eeb2c8fb432d1c6dcf792..63ea16f691dc4c32ba73785205f041859b0ea1f4 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -384,13 +384,13 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
while(it != edges.end()){
if ((*it)->getCompound()){
mySet.insert((*it)->getCompound());
- printf("compound edge %d found in %d\n",(*it)->getCompound()->tag(), (*it)->tag());
+ //printf("compound edge %d found in %d\n",(*it)->getCompound()->tag(), (*it)->tag());
}
else
mySet.insert(*it);
++it;
}
- printf("replacing %d edges by %d in the compound %d\n",edges.size(),mySet.size(),gf->tag());
+ //printf("Replacing %d edges by %d in the compound %d\n",edges.size(),mySet.size(),gf->tag());
edges.clear();
edges.insert(edges.begin(), mySet.begin(), mySet.end());
}
@@ -455,7 +455,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
reparamMeshVertexOnFace(here, gf, param);
U_[count] = param[0];
V_[count] = param[1];
- // printf("coucou : %g %g -> %g %g\n",here->x(),here->y(),param.x(),param.y());
+ //printf("*** meshGFace : %g %g -> u,v = %g %g\n",here->x(),here->y(),param.x(),param.y());
(*itv)->setIndex(count);
numbered_vertices[(*itv)->getIndex()] = *itv;
bbox += SPoint3(param.x(), param.y(), 0);