diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index 137f5e51d90f58c741072e269c7c3c38710cbb1d..cc96dc3bec6c7e9d4455b62b29a9d14406c49801 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -19,9 +19,10 @@ GEdgeCompound::GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound)
int N = _compound.size();
v0 = _orientation[0] ? _compound[0]->getBeginVertex() : _compound[0]->getEndVertex();
v1 = _orientation[N-1] ? _compound[N-1]->getEndVertex() : _compound[N-1]->getBeginVertex();
+ //printf("vertices of compound %d -> %d\n",v0->tag(),v1->tag());
v0->addEdge(this);
v1->addEdge(this);
- //printf("%d -> %d\n",v0->tag(),v1->tag());
+
parametrize();
}
@@ -31,6 +32,7 @@ void GEdgeCompound::orderEdges()
std::list<GEdge*> edges ;
for (int i=0;i<_compound.size();i++){
+ //printf("set compound %d for edge %d \n", tag(), _compound[i]->tag());
_compound[i]->setCompound(this);
edges.push_back(_compound[i]);
}
@@ -40,12 +42,18 @@ void GEdgeCompound::orderEdges()
for (std::list<GEdge*>::iterator it = edges.begin() ; it != edges.end() ; ++it){
GVertex *v1 = (*it)->getBeginVertex();
GVertex *v2 = (*it)->getEndVertex();
- //printf("BEFORE ORDERING segment vert=%d, %d y =%g %g\n", v1->tag(), v2->tag() , v1->y(), v2->y());
+ //printf("BEFORE ORDERING segment vert=%d, %d \n", v1->tag(), v2->tag());
std::map<GVertex*,GEdge*>::iterator it1 = tempv.find(v1);
- if (it1==tempv.end()) tempv.insert(std::make_pair(v1,*it));
+ if (it1==tempv.end()) {
+ //printf("insert v1=%d \n", v1->tag());
+ tempv.insert(std::make_pair(v1,*it));
+ }
else tempv.erase(it1);
std::map<GVertex*,GEdge*>::iterator it2 = tempv.find(v2);
- if (it2==tempv.end()) tempv.insert(std::make_pair(v2,*it));
+ if (it2==tempv.end()){
+ //printf("insert v2=%d \n", v2->tag());
+ tempv.insert(std::make_pair(v2,*it));
+ }
else tempv.erase(it2);
}
diff --git a/Geo/GEdgeCompound.h b/Geo/GEdgeCompound.h
index 127d1dc7e57b772174ab38170d372a2e39d6a3fc..3965e7b2e7fe2c99d98768fb8b66025e2b826f62 100644
--- a/Geo/GEdgeCompound.h
+++ b/Geo/GEdgeCompound.h
@@ -18,11 +18,13 @@ class GEdgeCompound : public GEdge {
std::vector<GEdge*> _compound;
std::vector<int> _orientation;
std::vector<double> _pars;
- void parametrize();
- void orderEdges();
-
- public:
- void getLocalParameter(const double &t, int &iEdge, double & tLoc) const;
+ void parametrize() ;
+ void orderEdges() ;
+public:
+ void getLocalParameter ( const double &t,
+ int &iEdge,
+ double & tLoc) const;
+
GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound);
virtual ~GEdgeCompound();
Range<double> parBounds(int i) const;
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 54b7c170559b81f0cdf9b04905bddc5891c094db..7429517bbc228d5a47c865eb1a9f11ac97168558 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -92,6 +92,8 @@ void GFaceCompound::parametrize() const
void GFaceCompound::getBoundingEdges()
{
+ printf("***** In GFaceCompound: size U0=%d, v0=%d\n ", _U0.size(), _V0.size());
+
if (_U0.size()){
std::list<GEdge*> :: const_iterator it = _U0.begin();
for ( ; it != _U0.end() ; ++it){
@@ -111,9 +113,9 @@ void GFaceCompound::getBoundingEdges()
std::list<GFace*>::iterator it = _compound.begin();
for ( ; it != _compound.end(); ++it){
std::list<GEdge*> ed = (*it)->edges();
- std::list<GEdge*> :: iterator ite = ed.begin();
+ std::list<GEdge*> :: iterator ite = ed.begin();
for ( ; ite != ed.end(); ++ite){
- _touched.insert(*ite);
+ _touched.insert(*ite);
}
}
it = _compound.begin();
@@ -121,13 +123,16 @@ void GFaceCompound::getBoundingEdges()
std::list<GEdge*> ed = (*it)->edges();
std::list<GEdge*> :: iterator ite = ed.begin();
for ( ; ite != ed.end() ; ++ite){
- if (!(*ite)->degenerate(0) && _touched.count(*ite) == 1)_unique.insert(*ite);
+ if (!(*ite)->degenerate(0) && _touched.count(*ite) == 1) {
+ _unique.insert(*ite);
+ }
}
}
std::set<GEdge*>::iterator itf = _unique.begin();
for ( ; itf != _unique.end(); ++itf){
l_edges.push_back(*itf);
+ printf("for edge %d, add face %d \n", (*itf)->tag(), this->tag());
(*itf)->addFace(this);
}
@@ -220,7 +225,6 @@ 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());
@@ -233,11 +237,11 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
//getParameter Point
v->getParameter(0,tGlob);
- //printf("global param for point (%g, %g, %g ) = %g\n", v->x(), v->y(), v->z(), 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());
+ printf("tag compound=%d, beginvertex=%d, endvertex=%d \n", gec->tag(), gec->getBeginVertex()->tag(), gec->getEndVertex()->tag());
if (gec){
@@ -245,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];
@@ -254,7 +258,7 @@ 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();
@@ -262,7 +266,7 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
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);
+ 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;
@@ -291,17 +295,16 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
//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);
+ 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("uloc=%g, vloc=%g \n", uloc,vloc);
- //printf("Vertex coordinates not found in compound face \n");
//exit(0);
return SPoint2(uloc,vloc);
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 24d72df6b1c3664b7cfb76e82b1441dfde962add..f675a006ee1acb8e720373b8934fe58171205bc4 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -94,6 +94,8 @@ class GModel
std::map<int, std::string> physicalNames, elementaryNames;
int partitionSize[2];
+ std::map<int, GEdge*> mesh2Topo;
+
public:
GModel(std::string name="");
virtual ~GModel();
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 40f2c3256f13ed0ef607a7c3586a97f2125c22a4..10d603b5d509394745a6b08a82f03c68eba77ebb 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -122,11 +122,11 @@ int GModel::importGEOInternals()
GEntity *ge = 0;
switch(p->Typ){
case MSH_PHYSICAL_POINT: ge = getVertexByTag(abs(num)); break;
- case MSH_PHYSICAL_LINE:
+ case MSH_PHYSICAL_LINE:
ge = getEdgeByTag(abs(num));
e_compound.push_back(getEdgeByTag(abs(num)));
break;
- case MSH_PHYSICAL_SURFACE:
+ case MSH_PHYSICAL_SURFACE:
ge = getFaceByTag(abs(num));
f_compound.push_back(getFaceByTag(abs(num)));
break;
@@ -181,15 +181,6 @@ int GModel::importGEOInternals()
Msg::Debug("%d Faces", faces.size());
Msg::Debug("%d Regions", regions.size());
- for(viter it = firstVertex(); it != lastVertex(); it++)
- Msg::Debug("Imported GEO vert of Type: %s", (*it)->getTypeString().c_str());
- for(eiter it = firstEdge(); it != lastEdge(); it++)
- Msg::Debug("Imported GEO edge of Type: %s", (*it)->getTypeString().c_str());
- for(fiter it = firstFace(); it != lastFace(); it++)
- Msg::Debug("Imported GEO face of Type: %s", (*it)->getTypeString().c_str());
- for(riter it = firstRegion(); it != lastRegion(); it++)
- Msg::Debug("Imported GEO region of Type: %s", (*it)->getTypeString().c_str());
-
return 1;
}
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 9513dd9dc36131b36961aa3d18b8c62fe32cc771..8836cb01d31d9bff0daebed34e6ead614fed2a94 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -26,6 +26,7 @@
#include "discreteRegion.h"
#include "MElement.h"
#include "GEdgeCompound.h"
+#include "GFaceCompound.h"
#include <iostream> // DBG
@@ -114,7 +115,7 @@ static void createElementMSH(GModel *m, int num, int type, int physical,
void GModel::createTopologyFromMSH(){
- //printf("Dans createTopologyFromMSH \n");
+ printf("***** In createTopologyFromMSH: \n");
std::vector<GEntity*> entities;
getEntities(entities);
@@ -141,79 +142,42 @@ void GModel::createTopologyFromMSH(){
break;
}
}
- //printf("vertices size =%d \n", vertices.size());
- //printf("edges size =%d \n", edges.size());
- //printf("faces size =%d \n", faces.size());
- //printf("regions size =%d \n", regions.size());
+ printf("vertices size =%d \n", vertices.size());
+ printf("edges size =%d \n", edges.size());
+ printf("faces size =%d \n", faces.size());
+ printf("regions size =%d \n", regions.size());
- int tag = 100;
- for (std::vector<discreteEdge*>::iterator edge = edges.begin();
- edge != edges.end(); edge++){
- if (tag < (*edge)->tag() ) tag = (*edge)->tag() + 1;
- }
-
- //For each discreteEdge, build a new GEdgeCompound
- for (std::vector<discreteEdge*>::iterator edge = edges.begin();
- edge != edges.end(); edge++){
-
- //printf("createTopology: %d EDGES, of size=%d\n",(*edge)->tag(), (*edge)->lines.size());
+// for(unsigned int i = 0; i < entities.size(); i++)
+// for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
+// printf("entity =%d dim = %d tag =%d \n", i, entities[i]->dim() , entities[i]->mesh_vertices[j]->getIndex());
+// }
+// exit(1);
- //create a map with the tags of the mesh vertices
- std::map<int, GVertex*> myMap;
- for (std::vector<MLine*>::const_iterator it = (*edge)->lines.begin();
- it != (*edge)->lines.end() ; ++it){
- int tagB = (*it)->getVertex(0)->getNum();
- int tagE = (*it)->getVertex(1)->getNum();
- std::map<int, GVertex*>::iterator it1 = myMap.find(tagB);
- std::map<int, GVertex*>::iterator it2 = myMap.find(tagE);
- if (it1 == myMap.end()){
- GVertex *gvB = new discreteVertex(this,tagB);
- gvB->mesh_vertices.push_back((*it)->getVertex(0));
- gvB->points.push_back(new MPoint(gvB->mesh_vertices.back()));
- myMap.insert(std::make_pair(tagB, gvB));
- }
- if (it2 == myMap.end()){
- GVertex *gvE = new discreteVertex(this,tagE);
- gvE->mesh_vertices.push_back((*it)->getVertex(1));
- gvE->points.push_back(new MPoint(gvE->mesh_vertices.back()));
- myMap.insert(std::make_pair(tagE, gvE));
- }
- }
+ //For each discreteEdge, create Topology
+ //---------------------------------------------------
- // for(std::map<int, GVertex*>::const_iterator it = myMap.begin(); it != myMap.end(); ++it){
-// printf(" tag=%d tagsize=%d\n", it->first, myMap.size());
-// }
+ for (std::vector<discreteEdge*>::iterator edge = edges.begin(); edge != edges.end(); edge++){
- //create a vector composed of plenty of discreteEdges from the Mlines of the original discreteVertex
- std::vector<GEdge*> e_compound;
+ //printf("createTopology: EDGE= %d, of size=%d\n",(*edge)->tag(), (*edge)->lines.size());
- for (std::vector<MLine*>::const_iterator it = (*edge)->lines.begin();
- it != (*edge)->lines.end(); ++it){
- //printf("MLine =%d %d \n", (*it)->getVertex(0)->getNum(), (*it)->getVertex(1)->getNum());
+ (*edge)->orderMLines();
+ (*edge)->parametrize();
+ (*edge)->setBoundVertices(vertices);
- int tagB = (*it)->getVertex(0)->getNum();
- int tagE = (*it)->getVertex(1)->getNum();
- std::map<int,GVertex*>::iterator it1 = myMap.find(tagB);
- std::map<int,GVertex*>::iterator it2 = myMap.find(tagE);
- GVertex *gvB = it1->second;
- GVertex *gvE = it2->second;
- GEdge *temp = new discreteEdge(this, tag, gvB, gvE); //new GEdge corresponding to the MLine
- gvB->addEdge(temp);
- gvE->addEdge(temp);
-
- e_compound.push_back(temp); //add the compound to the GEdge
- tag ++;
+ }
- }
+ //For each discreteFace, create Topology
+ //---------------------------------------------------
- //now, we can create the GEdgeCompound
- GEdge *gec = new GEdgeCompound(this, tag, e_compound);
- add(gec);
+ 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);
+ }
+
return;
}
@@ -657,6 +621,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
entities[i]->mesh_vertices[j]->writeMSH(fp, binary, saveParametric,
scalingFactor);
+
+
if(binary) fprintf(fp, "\n");
if(version >= 2.0){
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 02ece717c664cdcc340d6175698758e7244a15f9..a02eabbe429c708a68081a92d499aace4fc601cf 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -30,8 +30,6 @@ void MVertex::writeMSH(FILE *fp, bool binary, bool saveParametric, double scalin
{
if(_index < 0) return; // negative index vertices are never saved
- // printf("tag = %d index %d pos = %g %g %g\n",onWhat()->tag(),_index,x(),y(),z());
-
int myDim = 0, myTag = 0;
if(saveParametric){
if(onWhat()){
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 1cb0ee5c645804cf54a7e63efa5f302a68c3b2ce..f8f4d7ca588aee9bf56bb5a9c3ea2806b3fb8134 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -8,6 +8,7 @@
#include "discreteEdge.h"
#include "MLine.h"
#include "Numeric.h"
+#include "MPoint.h"
#include <vector>
#include <list>
@@ -27,10 +28,11 @@ discreteEdge::discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1)
}
-void discreteEdge::orderEdge()
+
+void discreteEdge::orderMLines()
{
- printf(" *** ORDERING DISCRETE EDGE %d of size %d \n", this->tag(), lines.size());
+ //printf(" *** ORDERING DISCRETE EDGE %d of size %d \n", this->tag(), lines.size());
std::vector<MLine*> _m ;
std::list<MLine*> segments;
@@ -42,7 +44,6 @@ void discreteEdge::orderEdge()
}
// find a lonly MLine
- std::map<MVertex*,MLine*> boundv;
for (std::list<MLine*>::iterator it = segments.begin() ; it != segments.end() ; ++it){
MVertex *vL = (*it)->getVertex(0);
MVertex *vR = (*it)->getVertex(1);
@@ -87,7 +88,7 @@ void discreteEdge::orderEdge()
for (std::list<MLine*>::iterator it = segments.begin() ; it != segments.end() ; ++it){
MLine *e = *it;
if (e->getVertex(0) == last){
- printf("orientation 1: beginV=%d \n", e->getVertex(0)->getNum());
+ //printf("orientation 1: beginV=%d \n", e->getVertex(0)->getNum());
_m.push_back(e);
segments.erase(it);
_orientation.push_back(1);
@@ -96,7 +97,7 @@ void discreteEdge::orderEdge()
break;
}
else if (e->getVertex(1) == last){
- printf("orientation 0: endV=%d \n", e->getVertex(1)->getNum());
+ //printf("orientation 0: endV=%d \n", e->getVertex(1)->getNum());
_m.push_back(e);
segments.erase(it);
_orientation.push_back(0);
@@ -111,7 +112,7 @@ void discreteEdge::orderEdge()
first = last;
last = temp;
_orientation[0] = 0;
- }
+ }
else {
Msg::Error("Discrete Edge %d is wrong",tag());
return;
@@ -130,13 +131,70 @@ void discreteEdge::orderEdge()
for (int i=0;i<lines.size();i++) _orientation[i] = !_orientation[i] ;
}
- for (int i=0;i<lines.size();i++){
- printf("AFTER ORDERING segment or=%d, vert=%d, %d\n", (int)_orientation[i], lines[i]->getVertex(0)->getNum(), lines[i]->getVertex(1)->getNum() );
+// for (int i=0;i<lines.size();i++){
+// printf("AFTER ORDERING segment or=%d, vert=%d, %d\n", (int)_orientation[i], lines[i]->getVertex(0)->getNum(), lines[i]->getVertex(1)->getNum() );
+// }
+
+
+ return;
+
+}
+
+void discreteEdge::setBoundVertices(std::vector<discreteVertex*> vertices)
+{
+
+ if (boundv.size()==2){
+ //printf("Found a regular open Curve \n");
+ std::vector<GVertex*> bound_vertices;
+ for (std::map<MVertex*,MLine*>::const_iterator iter = boundv.begin(); iter != boundv.end(); iter++){
+ MVertex* vE = (iter)->first;
+ bool existVertex = false;
+ for (std::vector<discreteVertex*>::iterator point = vertices.begin(); point != vertices.end(); point++) {
+ //printf("Discrete point =%d bound vE=%d\n", (*point)->tag(), vE->getNum());
+ if ((*point)->tag() == vE->getNum()){
+ bound_vertices.push_back((*point));
+ existVertex = true;
+ break;
+ }
+ }
+ if(!existVertex){
+ printf(" !!! bound vertex does not exist, create one \n");
+ GVertex *gvB = new discreteVertex(model(),vE->getNum());
+ bound_vertices.push_back(gvB);
+ gvB->mesh_vertices.push_back(vE);
+ gvB->points.push_back(new MPoint(gvB->mesh_vertices.back()));
+ model()->add(gvB);
+ mesh_vertices.erase(std::find(mesh_vertices.begin(), mesh_vertices.end(), vE));
+ }
+ }
+
+ //printf("set bound vertices %d %d \n",bound_vertices[0]->tag(),bound_vertices[1]->tag());
+ v0 = bound_vertices[0];
+ v1 = bound_vertices[1];
+ }
+ else if (boundv.size()==0){
+ //printf("Found a closed Curve \n");
+
+ std::vector<MLine*>::const_iterator it = lines.begin();
+ MVertex* vE = (*it)->getVertex(0);
+ GVertex *gvB = new discreteVertex(model(),vE->getNum());
+ gvB->mesh_vertices.push_back(vE);
+ gvB->points.push_back(new MPoint(gvB->mesh_vertices.back()));
+ model()->add(gvB);
+ mesh_vertices.erase(std::find(mesh_vertices.begin(), mesh_vertices.end(), vE));
+
+ //printf("set bound vertices %d %d , coord = %g %g %g\n",gvB->tag(),gvB->tag(), gvB->x(), gvB->y(), gvB->z());
+ v0 = gvB;
+ v1 = gvB;
+
}
- return;
+ v0->addEdge(this);
+ v1->addEdge(this);
+ return;
}
+
void discreteEdge::parametrize()
{
@@ -144,29 +202,29 @@ void discreteEdge::parametrize()
_pars.push_back(i);
}
-/*
- +---------------+--------------+----------- ... ----------+
- _pars[0]=0 _pars[1]=1 _pars[2]=2 _pars[N+1]=N+1
-*/
+ /*
+ +---------------+--------------+----------- ... ----------+
+ _pars[0]=0 _pars[1]=1 _pars[2]=2 _pars[N+1]=N+1
+ */
-// for (int i=0;i<lines.size()+1;i++){
-// printf("_pars[%d]=%g\n",i, _pars[i] );
-// }
+ // for (int i=0;i<lines.size()+1;i++){
+ // printf("_pars[%d]=%g\n",i, _pars[i] );
+ // }
}
void discreteEdge::getLocalParameter ( const double &t,
- int &iLine,
- double & tLoc) const
+ int &iLine,
+ double & tLoc) const
{
for (iLine=0 ; iLine<lines.size() ;iLine++){
- double tmin = _pars[iLine];
+ double tmin = _pars[iLine];
double tmax = _pars[iLine+1];
if (t >= tmin && t <= tmax){
tLoc = _orientation[iLine] ? (t-tmin)/(tmax-tmin) : 1 - (t-tmin)/(tmax-tmin) ;
//printf("getlocal param t=%g, iLine=%d, tLoc=%g \n", t, iLine, tLoc);
- return;
+ return;
}
}
}
@@ -174,14 +232,19 @@ void discreteEdge::getLocalParameter ( const double &t,
GPoint discreteEdge::point(double par) const
{
- //printf("dans point v0 =%d (%g,%g,%g), v1=%d (%g,%g,%g)\n", v0->tag(), v0->x(),v0->y(), v0->z(), v1->tag(), v1->x(),v1->y(), v1->z());
+ double tLoc;
+ int iEdge;
+ getLocalParameter(par,iEdge,tLoc);
double x, y, z;
-
- x = v0->x() + par * (v1->x()- v0->x());
- y = v0->y() + par * (v1->y()- v0->y());
- z = v0->z() + par * (v1->z()- v0->z());
+ MVertex *vB = lines[iEdge]->getVertex(0);
+ MVertex *vE = lines[iEdge]->getVertex(1);
+
+ x = vB->x() + tLoc * (vE->x()- vB->x());
+ y = vB->y() + tLoc * (vE->y()- vB->y());
+ z = vB->z() + tLoc * (vE->z()- vB->z());
+ //printf("dans point par=%d iEdge =%d, tLoc=%d \n", par, iEdge, tLoc);
//printf("Discrete Edge POint par=%g, x= %g, y = %g, z = %g \n", par, x,y,z);
return GPoint(x,y,z);
}
@@ -189,11 +252,18 @@ GPoint discreteEdge::point(double par) const
SVector3 discreteEdge::firstDer(double par) const
{
+ double tLoc;
+ int iEdge;
+ getLocalParameter(par,iEdge,tLoc);
+
+ MVertex *vB = lines[iEdge]->getVertex(0);
+ MVertex *vE = lines[iEdge]->getVertex(1);
+
double dx,dy,dz;
- double dt = sqrt((v1->x()-v0->x())*(v1->x()-v0->x()) + (v1->y()-v0->y())*(v1->y()-v0->y()) + (v1->z()-v0->z())*(v1->z()-v0->z()));
- dx = (v1->x() - v0->x() ) /dt;
- dy = (v1->y() - v0->y() ) /dt;
- dz = (v1->z() - v0->z() ) /dt;
+ double dt = sqrt((vE->x()-vB->x())*(vE->x()-vB->x()) + (vE->y()-vB->y())*(vE->y()-vB->y()) + (vE->z()-vB->z())*(vE->z()-vB->z()));
+ dx = (vE->x() - vB->x() ) /dt;
+ dy = (vE->y() - vB->y() ) /dt;
+ dz = (vE->z() - vB->z() ) /dt;
//printf("firstDer discreteEdge par=%g, dx=%g, dy=%g dz=%g dt=%g \n", par,dx,dy,dz, dt);
return SVector3(dx,dy,dz);
@@ -202,8 +272,6 @@ SVector3 discreteEdge::firstDer(double par) const
Range<double> discreteEdge::parBounds(int i) const {
- return Range<double>(0, 1);
- //return Range<double>(0, lines.size()-1);
- //return Range<double>(0, _pars[lines.size()-1]);
-
+ return Range<double>(0, lines.size());
+
}
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index c69e9b0564324659d0f60a177e35be6a43ae6132..c6955da8610078c6bfa90338fe46feb337b2671a 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -8,11 +8,13 @@
#include "GModel.h"
#include "GEdge.h"
+#include "discreteVertex.h"
class discreteEdge : public GEdge {
protected:
std::vector<double> _pars;
std::vector<int> _orientation;
+ std::map<MVertex*,MLine*> boundv;
public:
discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1);
virtual ~discreteEdge() {}
@@ -23,9 +25,9 @@ class discreteEdge : public GEdge {
virtual GPoint point(double p) const;
virtual SVector3 firstDer(double par) const;
virtual Range<double> parBounds(int) const;
- void setVertices(GVertex *_v0, GVertex *_v1){ v0 = _v0; v1 = _v1; }
- void orderEdge();
- void parametrize();
+ void parametrize() ;
+ void orderMLines() ;
+ void setBoundVertices( std::vector<discreteVertex*> vertices );
};
#endif
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index d1751560289b70e4c7fa44af31276145088159cb..0ea510ebc20d2b7b00744ad6e09aec3efc822a2e 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -6,6 +6,8 @@
#include "GmshConfig.h"
#include "GmshMessage.h"
#include "discreteFace.h"
+#include "MTriangle.h"
+#include "MEdge.h"
#if !defined(HAVE_GMSH_EMBEDDED)
#include "Geo.h"
@@ -20,6 +22,62 @@ discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
meshStatistics.status = GFace::DONE;
}
+void discreteFace::setBoundEdges(std::vector<discreteEdge*> discr_edges)
+{
+
+ printf("***** In discrete Face: \n");
+
+ printf("bound edges =%d \n", edges().size());
+
+ for (std::vector<discreteEdge*>::iterator it = discr_edges.begin(); it != discr_edges.end(); it++) {
+ l_edges.push_back(*it);
+ l_dirs.push_back(1);
+ (*it)->addFace(this);
+ printf("add Face %d for edge %d\n",this->tag(), (*it)->tag() );
+ }
+
+ printf("bound edges =%d \n", edges().size());
+
+// std::list<MVertex*> mesh_vertices;
+// std::list<MVertex*> temp_vertices;
+// std::list<MEdge> bound_edges;
+// for (int iFace = 0; iFace < getNumMeshElements() ; iFace++) {
+// std::vector<MVertex*> fv;
+// MElement *e = getMeshElement(iFace);
+// e->getFaceVertices(0, fv);
+// for (std::vector<MVertex*>::iterator it_fv = fv.begin(); it_fv != fv.end(); it_fv++) {
+// if (std::find(mesh_vertices.begin(), mesh_vertices.end(), *it_fv) == mesh_vertices.end())
+// mesh_vertices.push_back(*it_fv);
+// }
+// for (int iEdge = 0; iEdge < e->getNumEdges(); iEdge++) {
+// MEdge tmp_edge = e->getEdge(iEdge);
+// if (std::find(bound_edges.begin(),bound_edges.end(),tmp_edge) == bound_edges.end())
+// bound_edges.push_back(tmp_edge);
+// elsex
+// bound_edges.erase(std::find(bound_edges.begin(),bound_edges.end(),tmp_edge));
+// }
+// }
+// printf( "There are %d msh vertices and bound %d msh edges \n ", mesh_vertices.size(), bound_edges.size());
+
+ // for (std::list<MVertex *>::iterator mesh_vertex = mesh_vertices.begin(); mesh_vertex != mesh_vertices.end(); mesh_vertex++) {
+ // printf("mesh vertex on entity with tag %d\n", (*mesh_vertex)->onWhat()->tag());
+ //}
+
+
+// std::list<GEdge*> adj_edges;
+// for (std::list<MEdge>::iterator be = bound_edges.begin(); be != bound_edges.end(); be++) {
+// MVertex *v0 = be->getVertex(0);
+// MVertex *v1 = be->getVertex(1);
+// printf("bound edge on entity with num=%d, tag %d (dim=%d) et num=%d tag %d (dim=%d) \n", v0->getNum(), v0->onWhat()->tag(), v0->onWhat()->dim(),v1->getNum(), v1->onWhat()->tag(), v0->onWhat()->dim());
+// }
+
+
+
+ //exit(1);
+
+}
+
+
GPoint discreteFace::point(double par1, double par2) const
{
Msg::Error("Cannot evaluate point on discrete face");
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index 554bd344b769d90a0bdd8884dd062add28267d68..bca7f3cac8d39b18e8de1a089096cd945ad6f20d 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -8,6 +8,7 @@
#include "GModel.h"
#include "GFace.h"
+#include "discreteEdge.h"
class discreteFace : public GFace {
public:
@@ -18,6 +19,7 @@ class discreteFace : public GFace {
virtual SVector3 normal(const SPoint2 ¶m) const;
virtual GEntity::GeomType geomType() const { return DiscreteSurface; }
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
+ void setBoundEdges( std::vector<discreteEdge*> edges );
};
#endif
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 63ea16f691dc4c32ba73785205f041859b0ea1f4..efd60491fd660f109168ec014531e567fe607b8d 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -378,19 +378,21 @@ 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");
if (gf->geomType() == GEntity::CompoundSurface){
+ printf("replace edges by compound lines \n");
std::set<GEdge*> mySet;
std::list<GEdge*>::iterator it = edges.begin();
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 edge %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 GFaceCompound %d\n",edges.size(),mySet.size(),gf->tag());
edges.clear();
edges.insert(edges.begin(), mySet.begin(), mySet.end());
}
@@ -404,7 +406,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
// build a set with all points of the boundaries
it = edges.begin();
while(it != edges.end()){
- if ((*it)->isSeam(gf)) return false;
+ if ((*it)->isSeam(gf)) return false;
if(!(*it)->isMeshDegenerated()){
for (unsigned int i = 0; i< (*it)->lines.size(); i++){
all_vertices.insert((*it)->lines[i]->getVertex(0));
@@ -446,6 +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());
int count = 0;
SBoundingBox3d bbox;
@@ -455,7 +458,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
reparamMeshVertexOnFace(here, gf, param);
U_[count] = param[0];
V_[count] = param[1];
- //printf("*** meshGFace : %g %g -> u,v = %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);
@@ -480,6 +483,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
// We add to the triangulation a box with 4 points that encloses the
// domain.
{
+
DocRecord doc(all_vertices.size() + 4);
itv = all_vertices.begin();
int j = 0;
diff --git a/benchmarks/2d/square.geo b/benchmarks/2d/square.geo
new file mode 100644
index 0000000000000000000000000000000000000000..7b918bce3e34d653d05aae7ff9b825ad319c4fa6
--- /dev/null
+++ b/benchmarks/2d/square.geo
@@ -0,0 +1,17 @@
+lc=0.2;
+Point(1) = {0, 0, 0,lc};
+Point(2) = {0, 1, 0,lc};
+Point(3) = {1, 1, 0,lc};
+Point(4) = {1, 0, 0,lc};
+Line(1) = {2, 3};
+Line(2) = {3, 4};
+Line(3) = {4, 1};
+Line(4) = {1, 2};
+
+Line Loop(5) = {1, 2, 3, 4};
+Plane Surface(10) = {5};
+
+Compound Line(150)={1,2};
+Compound Line(160)={3,4};
+
+Compound Surface(170)={10} Boundary {{}};
diff --git a/benchmarks/2d/square_CLASS_GEO.geo b/benchmarks/2d/square_CLASS_GEO.geo
new file mode 100644
index 0000000000000000000000000000000000000000..52d4b52fca2ec34d99a57b1b338bdf5ba01ab088
--- /dev/null
+++ b/benchmarks/2d/square_CLASS_GEO.geo
@@ -0,0 +1,14 @@
+Mesh.CharacteristicLengthFactor=6;
+
+Merge "square_CLASS.msh";
+CreateTopology;
+
+//Compound Line(100)={1};
+//Compound Line(200)={2};
+//Compound Line(300)={3};
+//Compound Line(400)={4};
+
+Compound Line(150)={1,2};
+Compound Line(160)={3,4};
+
+Compound Surface(170)={10} Boundary {{}};
\ No newline at end of file