diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index ad678591bbc191c81139d4c0b5a17ac8d7f14c7a..e94d5707aad8d36433bdf09f029ec931772f10bf 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -88,7 +88,7 @@ class GEdge : public GEntity {
// get second derivative of edge at the given parameter (default
// implentation using central differences)
virtual SVector3 secondDer(double par) const;
-
+
// get the curvature
virtual double curvature(double par) const;
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 7429517bbc228d5a47c865eb1a9f11ac97168558..0f866e70f7390a0b723cb6d8eecc18da487f04b9 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -241,7 +241,7 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
//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());
+ //printf("tag compound=%d, beginvertex=%d, endvertex=%d \n", gec->tag(), gec->getBeginVertex()->tag(), gec->getEndVertex()->tag());
if (gec){
@@ -249,7 +249,7 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
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);
+ //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];
@@ -258,45 +258,47 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
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());
+ //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);
+ double tB = ge->parBounds(0).low();
+ double tE = ge->parBounds(0).high();
+ //printf("tB=%g uv (%g %g) tE=%g (%g %g), tLoc=%g\n", tB, itL->second.x(), itL->second.y(), tE, itR->second.x(), itR->second.y(), 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());
+ tL=tB;
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)){
+ vR->getParameter(0,tR);
+ //tR = j+1;//HACK HERE
+ //if (!vR->getParameter(0,tR)) {
+ //printf("vertex vr %p not medgevertex \n", vR);
+ //exit(1);
+ //}
+ printf("***tLoc=%g tL=%g, tR=%g L=%p (%g,%g,%g) et R= %p (%g,%g,%g)\n", tLoc, tL, tR, vL, vL->x(),vL->y(),vL->z(),vR, vR->x(), vR->y(), vR->z());
+ if (tLoc >= tL && tLoc <= tR){
found = true;
itR = coordinates.find(vR);
- vR->getParameter(0,tR);
+ if (itR == coordinates.end()){
+ printf("vertex %p (%g %g %g) not found\n", vR, vR->x(), vR->y(), vR->z());
+ exit(1);
+ }
break;
}
else{
- xL = xR; yL = yR ; zL = zR;
itL = coordinates.find(vR);
vL = vR;
- vL->getParameter(0,tL);
+ tL = tR;
}
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());
+ if (!found) {
+ vR = v1;
+ tR = tE;
+ }
+ //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
@@ -304,8 +306,7 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
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);
-
- //exit(0);
+ //exit(1);
return SPoint2(uloc,vloc);
}
@@ -435,6 +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());
std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
if (itf == coordinates.end()){
SPoint3 p(0,0,0);
@@ -443,7 +446,9 @@ void GFaceCompound::parametrize(iterationStep step) const
}
else{
itf->second[step]= value;
- }
+ }
+
+
}
}
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 8836cb01d31d9bff0daebed34e6ead614fed2a94..2dbe07da157cffba8bd91af710d64d2439aeacf0 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -153,31 +153,31 @@ void GModel::createTopologyFromMSH(){
// }
// exit(1);
+ //For each discreteFace, create Topology and if needed create discreteEdges
+ //----------------------------------------------------------------------------
+
+ for (std::vector<discreteFace*>::iterator face = faces.begin(); face != faces.end(); face++){
+
+ printf("createTopology: FACE=%d, of size=%d\n",(*face)->tag(), (*face)->getNumMeshElements());
+ (*face)->setBoundEdges(edges);
+
+ }
//For each discreteEdge, create Topology
//---------------------------------------------------
for (std::vector<discreteEdge*>::iterator edge = edges.begin(); edge != edges.end(); edge++){
-
- //printf("createTopology: EDGE= %d, of size=%d\n",(*edge)->tag(), (*edge)->lines.size());
+
+ printf("createTopology: EDGE= %d, of size=%d\n",(*edge)->tag(), (*edge)->lines.size());
(*edge)->orderMLines();
(*edge)->parametrize();
(*edge)->setBoundVertices(vertices);
-
-
+
}
- //For each discreteFace, create Topology
- //---------------------------------------------------
-
- for (std::vector<discreteFace*>::iterator face = faces.begin(); face != faces.end(); face++){
-
- //printf("createTopology: FACE=%d, of size=%d\n",(*face)->tag(), (*face)->getNumMeshElements());
- (*face)->setBoundEdges(edges);
+ //exit(1);
- }
-
return;
}
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 3873f65405abd2d6a0be50c271645584b597c26d..6a17d5e014ad5728ee2f801859a1ae9257dc7fe1 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -133,7 +133,8 @@ class MEdgeVertex : public MVertex{
{
}
virtual ~MEdgeVertex(){}
- virtual bool getParameter(int i, double &par) const { par = _u; return true; }
+ virtual bool getParameter(int i, double &par) const {
+ par = _u; return true; }
virtual bool setParameter(int i, double par){ _u = par; return true; }
double getLc() const { return _lc; }
};
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index f8f4d7ca588aee9bf56bb5a9c3ea2806b3fb8134..64a3689326d958d2d6ebb34fe800ae91caeae695 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -9,6 +9,7 @@
#include "MLine.h"
#include "Numeric.h"
#include "MPoint.h"
+#include "MTriangle.h"
#include <vector>
#include <list>
@@ -211,6 +212,62 @@ void discreteEdge::parametrize()
// printf("_pars[%d]=%g\n",i, _pars[i] );
// }
+ printf("dans discrete edge %d line.size =%d \n", this->tag(), lines.size());
+
+// std::vector<MVertex*> new_mshv;
+// for(int i = 0; i < mesh_vertices.size(); i++){
+// MVertex *v = mesh_vertices[i];
+// int param = i+1;
+// MVertex *newv = new MEdgeVertex(v->x(),v->y(),v->z(), this, param);
+// new_mshv.push_back(newv);
+// newv->setNum(v->getNum());
+// v= newv;
+// }
+// mesh_vertices = new_mshv;
+
+// std::vector<MVertex*> new_mshv;
+// for(int i = 0; i < mesh_vertices.size(); i++){
+// MVertex *vi = mesh_vertices[i];
+// MVertex *mev = new MEdgeVertex(vi->x(),vi->y(),vi->z(), this, i+1);
+// new_mshv.push_back(mev);
+// mev->setNum(vi->getNum());
+
+// newLines.push_back(new MLine(v1, newv));
+// newLines.push_back(new MLine(newv, v2));
+// delete ge->lines[i];
+// for(std::list<GFace*>::iterator it = l_faces.begin(); it != l_faces.end(); ++it){
+// for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+// MTriangle *t = (*it)->triangles[i];
+// for (int j = 0; j < 3; j++){
+// MVertex *v = t->getVertex(j);
+// if (v == vi) v = mev;
+// }
+// }
+// }
+// delete vi ;
+// }
+// mesh_vertices = new_mshv;
+
+
+// std::vector<MLine*> newLines;
+// newLines.push_back(new MLine(v1, newv));
+// delete lines[i];
+// newLines.push_back(new MLine(newv, v2));
+// lines = newLines;
+
+
+// for (int i = 0; i < mesh_vertices.size(); i++){
+// double t1;
+// mesh_vertices[i]->getParameter(0,t1);
+// printf("** AFTER v1=%d t1=%g\n", mesh_vertices[i]->getNum(),t1 );
+// }
+
+// for (int i = 0; i < lines.size(); i++){
+// printf("** AFTER LINES v1=%d v2=%d\n", lines[i]->getVertex(0)->getIndex(), lines[i]->getVertex(1)->getIndex() );
+// }
+
+ //exit(1);
+
}
void discreteEdge::getLocalParameter ( const double &t,