diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index ac25984062206961790560a660b8a7d7384f4a62..d47cf368577e42be37bb9e2e37bfdade73d3dca1 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -92,7 +92,7 @@ void GFaceCompound::parametrize() const
void GFaceCompound::getBoundingEdges()
{
- printf("***** In GFaceCompound: size U0=%d, v0=%d\n ", _U0.size(), _V0.size());
+ //printf("***** In GFaceCompound: size U0=%d, v0=%d\n ", _U0.size(), _V0.size());
if (_U0.size()){
std::list<GEdge*> :: const_iterator it = _U0.begin();
@@ -132,7 +132,7 @@ void GFaceCompound::getBoundingEdges()
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());
+ //printf("for edge %d, add face %d \n", (*itf)->tag(), this->tag());
(*itf)->addFace(this);
}
@@ -223,6 +223,7 @@ static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
{
+
parametrize() ;
std::map<MVertex*,SPoint3>::iterator it = coordinates.find(v);
@@ -237,7 +238,7 @@ 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());
@@ -270,12 +271,11 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
while (j < ge->mesh_vertices.size()){
vR = ge->mesh_vertices[j];
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 (!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);
@@ -299,13 +299,13 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
}
//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("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);
//exit(1);
return SPoint2(uloc,vloc);
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 67051916df6d2308dfb89180444908a7066c41a6..8c161bc9dee772ba2d06472f17756c884f84689c 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -578,7 +578,7 @@ int GModel::indexMeshVertices(bool all)
entities[i]->mesh_vertices[j]->setIndex(-1);
// tag all mesh vertices belonging to elements that need to be saved
- // with 0
+ //with 0
for(unsigned int i = 0; i < entities.size(); i++)
if(all || entities[i]->physicals.size())
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++)
@@ -588,9 +588,21 @@ int GModel::indexMeshVertices(bool all)
// renumber all the mesh vertices tagged with 0
int numVertices = 0;
for(unsigned int i = 0; i < entities.size(); i++)
- for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
- if(!entities[i]->mesh_vertices[j]->getIndex())
- entities[i]->mesh_vertices[j]->setIndex(++numVertices);
+ for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
+ //printf("-->entity with %d elems and index %d and mesh_vertices=%d \n ", entities[i]->getNumMeshElements(), entities[i]->mesh_vertices[j]->getIndex(), entities[i]->mesh_vertices.size());
+ if(!entities[i]->mesh_vertices[j]->getIndex()){
+ int newIndex= ++numVertices;
+ entities[i]->mesh_vertices[j]->setIndex(newIndex);
+ //printf("-->entity with %d elems with index=%d\n ", entities[i]->getNumMeshElements(), newIndex);
+ }
+ }
+
+// printf("numvertices=%d \n", numVertices);
+// for(unsigned int i = 0; i < entities.size(); i++)
+// for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++)
+// for(int k = 0; k < entities[i]->getMeshElement(j)->getNumVertices(); k++)
+// printf("entity with =%d elems, and mesh_vertices =%d index=%d \n", entities[i]->getNumMeshElements(), entities[i]->mesh_vertices.size(), entities[i]->getMeshElement(j)->getVertex(k)->getIndex());
+// exit(1);
return numVertices;
}
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 2dbe07da157cffba8bd91af710d64d2439aeacf0..beef8e7e8f144650abf3141f6695af4fba78fb36 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -120,63 +120,166 @@ void GModel::createTopologyFromMSH(){
std::vector<GEntity*> entities;
getEntities(entities);
- std::vector<discreteVertex*> vertices;
- std::vector<discreteEdge*> edges;
- std::vector<discreteFace*> faces;
- std::vector<discreteRegion*> regions;
+ std::vector<discreteVertex*> Dvertices;
+ std::vector<discreteEdge*> Dedges;
+ std::vector<discreteFace*> Dfaces;
+ std::vector<discreteRegion*> Dregions;
for (std::vector<GEntity*>::iterator entity = entities.begin();
entity != entities.end(); entity++) {
switch ((*entity)->dim()) {
case 0:
- vertices.push_back((discreteVertex*) *entity);
+ Dvertices.push_back((discreteVertex*) *entity);
break;
case 1:
- edges.push_back((discreteEdge*) *entity);
+ Dedges.push_back((discreteEdge*) *entity);
break;
case 2:
- faces.push_back((discreteFace*) *entity);
+ Dfaces.push_back((discreteFace*) *entity);
break;
case 3:
- regions.push_back((discreteRegion*) *entity);
+ Dregions.push_back((discreteRegion*) *entity);
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());
-// 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);
+// printf("vertices size =%d \n", Dvertices.size());
+// printf("edges size =%d \n", Dedges.size());
+// printf("faces size =%d \n", Dfaces.size());
+// printf("regions size =%d \n", Dregions.size());
+
//For each discreteFace, create Topology and if needed create discreteEdges
//----------------------------------------------------------------------------
- for (std::vector<discreteFace*>::iterator face = faces.begin(); face != faces.end(); face++){
+ int initSizeEdges = Dedges.size();
+
+ //find boundary edges of each face and put them in
+ //a map_edges that associates
+ //the MEdges with the tags of the adjacent faces
- printf("createTopology: FACE=%d, of size=%d\n",(*face)->tag(), (*face)->getNumMeshElements());
- (*face)->setBoundEdges(edges);
+ std::map<MEdge, std::vector<int>, Less_Edge > map_edges;
+ for (std::vector<discreteFace*>::iterator face = Dfaces.begin(); face != Dfaces.end(); face++){
+ (*face)->findEdges(map_edges);
}
- //For each discreteEdge, create Topology
- //---------------------------------------------------
+ //create reverse map, for each face find set of MEdges
+ //that are candidate for new discrete Edges
- for (std::vector<discreteEdge*>::iterator edge = edges.begin(); edge != edges.end(); edge++){
+ int num = Dedges.size()+1;
+ std::map<int, std::vector<int> > face2Edges;
+
+ while (!map_edges.empty()){
+
+ //printf("********** new candidate discrete Edge\n");
+ std::vector<MEdge> myEdges;
+ std::vector<int> tagFaces = map_edges.begin()->second;
+ myEdges.push_back(map_edges.begin()->first);
+ map_edges.erase(map_edges.begin());
- printf("createTopology: EDGE= %d, of size=%d\n",(*edge)->tag(), (*edge)->lines.size());
+ std::map<MEdge, std::vector<int>, Less_Edge >::iterator itmap = map_edges.begin();
+ while (itmap != map_edges.end()){
+
+ std::vector<int> tagFaces2 = itmap->second;
+ if (tagFaces2 == tagFaces){
+ myEdges.push_back(itmap->first);
+ map_edges.erase(itmap++);
+ }
+ else
+ itmap++;
+ }
+
+ //if the loaded mesh already contains discrete Edges
+ //check if the candidate discrete Edge does contain any of those
+ //if not, create discreteEdges
+ //create a map face2Edges that associate
+ //for each face the boundary discrete Edges
+
+ if (initSizeEdges != 0 ){
+ //printf(" !!! discrete edges already exist \n");
+ std::vector<int> tagEdges;
+ for(int i = 0; i < myEdges.size(); i++){
+ for (int j=0; j<2; j++) {
+ if (myEdges[i].getVertex(0)->onWhat()->dim() == 1) {
+ int tagEdge = myEdges[i].getVertex(0)->onWhat()->tag();
+ std::vector<int>::iterator itv = std::find(tagEdges.begin(), tagEdges.end(), tagEdge);
+ if (itv == tagEdges.end()) tagEdges.push_back(tagEdge);
+ }
+ }
+ }
+ for (std::vector<int>::iterator itFace = tagFaces.begin(); itFace != tagFaces.end(); itFace++) {
+ std::map<int, std::vector<int> >::iterator it = face2Edges.find(*itFace);
+ if (it == face2Edges.end()) {
+ std::vector<int> allEdges;
+ allEdges.insert(allEdges.begin(), tagEdges.begin(), tagEdges.end());
+ face2Edges.insert(std::make_pair(*itFace,allEdges));
+ }
+ else{
+ std::vector<int> allEdges = it->second;
+ allEdges.insert(allEdges.begin(), tagEdges.begin(), tagEdges.end());
+ it->second = allEdges;
+ }
+ face2Edges.insert(std::make_pair(*itFace, tagEdges));
+ }
+ }
+ else{
+ discreteEdge *e = new discreteEdge(this, num, 0, 0);
+ add(e);
+ Dedges.push_back(e);
+ std::list<MVertex*> all_vertices;
+ for(int i = 0; i < myEdges.size(); i++) {
+ MVertex *v0 = myEdges[i].getVertex(0);
+ MVertex *v1 = myEdges[i].getVertex(1);
+ e->lines.push_back(new MLine( v0, v1));
+ if (std::find(all_vertices.begin(), all_vertices.end(), v0) == all_vertices.end()) all_vertices.push_back(v0);
+ if (std::find(all_vertices.begin(), all_vertices.end(), v1) == all_vertices.end()) all_vertices.push_back(v1);
+ }
+ e->mesh_vertices.insert(e->mesh_vertices.begin(), all_vertices.begin(), all_vertices.end());
+
+ for (std::vector<int>::iterator itFace = tagFaces.begin(); itFace != tagFaces.end(); itFace++) {
+ GFace *dFace = getFaceByTag(abs(*itFace));
+ for (std::list<MVertex*>::iterator itv = all_vertices.begin(); itv != all_vertices.end(); itv++) {
+ std::vector<MVertex*>::iterator itve = std::find(dFace->mesh_vertices.begin(), dFace->mesh_vertices.end(), *itv) ;
+ if (itve != dFace->mesh_vertices.end()) dFace->mesh_vertices.erase(itve);
+ (*itv)->setEntity(e);
+ }
+
+ std::map<int, std::vector<int> >::iterator f2e = face2Edges.find(*itFace);
+ if (f2e == face2Edges.end()){
+ std::vector<int> tagEdges;
+ tagEdges.push_back(num);
+ face2Edges.insert(std::make_pair(*itFace,tagEdges));
+ }
+ else{
+ std::vector<int> tagEdges = f2e->second;
+ tagEdges.push_back(num);
+ f2e->second = tagEdges;
+ }
+ }
+ num++;
+ }
+
+ };
+
+ //set boundary edges for each face
+ for (std::vector<discreteFace*>::iterator face = Dfaces.begin(); face != Dfaces.end(); face++){
+ std::map<int, std::vector<int> >::iterator ite = face2Edges.find((*face)->tag());
+ std::vector<int> myEdges = ite->second;
+ (*face)->setBoundEdges(myEdges);
+ }
+
+
+ //For each discreteEdge, create Topology
+ //---------------------------------------
+ for (std::vector<discreteEdge*>::iterator edge = Dedges.begin(); edge != Dedges.end(); edge++){
+
(*edge)->orderMLines();
+ (*edge)->setBoundVertices();
(*edge)->parametrize();
- (*edge)->setBoundVertices(vertices);
-
- }
- //exit(1);
+ }
return;
@@ -541,6 +644,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
// if there are no physicals we save all the elements
if(noPhysicalGroups()) saveAll = true;
+
+
// get the number of vertices and index the vertices in a continuous
// sequence
int numVertices = indexMeshVertices(saveAll);
@@ -581,6 +686,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
if(n) elements[(*it)->pyramids[0]->getTypeForMSH()] += n;
}
+
+
int numElements = 0;
std::map<int, int>::const_iterator it = elements.begin();
for(; it != elements.end(); ++it)
@@ -621,10 +728,10 @@ 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(binary) fprintf(fp, "\n");
+
if(version >= 2.0){
if(saveParametric)
fprintf(fp, "$EndParametricNodes\n");
@@ -654,6 +761,7 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
for(fiter it = firstFace(); it != lastFace(); ++it)
writeElementsMSH(fp, (*it)->triangles, saveAll, version, binary, num,
(*it)->tag(), (*it)->physicals);
+
writeElementHeaderMSH(binary, fp, elements, MSH_QUA_4, MSH_QUA_9, MSH_QUA_8);
for(fiter it = firstFace(); it != lastFace(); ++it)
writeElementsMSH(fp, (*it)->quadrangles, saveAll, version, binary, num,
diff --git a/Geo/MEdge.h b/Geo/MEdge.h
index e76343c46ba9370daa82fe6e8e51afa67892e19a..8f79e7d83e7ca5b8bc282f5d940b2084a1dac9d2 100644
--- a/Geo/MEdge.h
+++ b/Geo/MEdge.h
@@ -106,6 +106,7 @@ struct Less_Edge : public std::binary_function<MEdge, MEdge, bool> {
if(e1.getMaxVertex() < e2.getMaxVertex()) return true;
return false;
}
+
};
#endif
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index fa48c1dab8886824b97c96df7ff14c7b29ec3ace..d996942034f3d41f77386e6eacda3f88c08357df 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -403,6 +403,7 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
for(int i = 0; i < n; i++)
verts[i] = getVertex(i)->getIndex();
+
if(!binary){
for(int i = 0; i < n; i++)
fprintf(fp, " %d", verts[i]);
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index a02eabbe429c708a68081a92d499aace4fc601cf..1080981bdb7da8c1050cb02c592d6aee72e26326 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -303,6 +303,8 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 ¶m)
{
+
+
if (gf->geomType() == GEntity::CompoundSurface &&
v->onWhat()->dim() < 2){
GFaceCompound *gfc = (GFaceCompound*) gf;
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index c80ed0254d2f6a4b7d171dc5698ea44935499a17..3aadb7f106c7996918fcfe798447c087cbd015ab 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -48,13 +48,12 @@ void discreteEdge::orderMLines()
for (std::list<MLine*>::iterator it = segments.begin() ; it != segments.end() ; ++it){
MVertex *vL = (*it)->getVertex(0);
MVertex *vR = (*it)->getVertex(1);
- //printf("MLIne %d %d = (%g, %g, %g) (%g, %g, %g)\n",vL->getNum(), vR->getNum(), vL->x(),vL->y(),vL->z(), vR->x(),vR->y(),vR->z() );
std::map<MVertex*,MLine*>::iterator it1 = boundv.find(vL);
- if (it1==boundv.end()) boundv.insert(std::make_pair(vL,*it));
- else boundv.erase(it1);
+ if (it1==boundv.end()) boundv.insert(std::make_pair(vL,*it));
+ else boundv.erase(it1);
std::map<MVertex*,MLine*>::iterator it2 = boundv.find(vR);
- if (it2==boundv.end()) boundv.insert(std::make_pair(vR,*it));
- else boundv.erase(it2);
+ if (it2==boundv.end()) boundv.insert(std::make_pair(vR,*it));
+ else boundv.erase(it2);
}
//find the first MLine and erase it from the list segments
@@ -107,7 +106,7 @@ void discreteEdge::orderMLines()
break;
}
}
- if (!found && _orientation[0]){ //reverse orientation of first Line
+ if (!found && _orientation[0]==1){ //reverse orientation of first Line
if (_m.size() == 1 ){
MVertex *temp = first;
first = last;
@@ -141,7 +140,7 @@ void discreteEdge::orderMLines()
}
-void discreteEdge::setBoundVertices(std::vector<discreteVertex*> vertices)
+void discreteEdge::setBoundVertices()
{
if (boundv.size()==2){
@@ -150,7 +149,8 @@ void discreteEdge::setBoundVertices(std::vector<discreteVertex*> 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++) {
+ //for (std::vector<discreteVertex*>::iterator point = vertices.begin(); point != vertices.end(); point++) {
+ for(GModel::viter point = model()->firstVertex(); point != model()->lastVertex(); point++){
//printf("Discrete point =%d bound vE=%d\n", (*point)->tag(), vE->getNum());
if ((*point)->tag() == vE->getNum()){
bound_vertices.push_back((*point));
@@ -159,17 +159,20 @@ void discreteEdge::setBoundVertices(std::vector<discreteVertex*> vertices)
}
}
if(!existVertex){
- printf(" !!! bound vertex does not exist, create one \n");
+ //printf(" !!! bound vertex %d does not exist, create one \n", vE->getNum());
GVertex *gvB = new discreteVertex(model(),vE->getNum());
bound_vertices.push_back(gvB);
+ vE->setEntity(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));
}
+ std::vector<MVertex*>::iterator itve = std::find(mesh_vertices.begin(), mesh_vertices.end(), vE) ;
+ if (itve != mesh_vertices.end()) mesh_vertices.erase(itve);
+
}
- //printf("set bound vertices %d %d \n",bound_vertices[0]->tag(),bound_vertices[1]->tag());
+ //printf("set bound vertices %d %d size mesh-vertices =%d \n",bound_vertices[0]->tag(),bound_vertices[1]->tag(), mesh_vertices.size());
v0 = bound_vertices[0];
v1 = bound_vertices[1];
}
@@ -184,7 +187,7 @@ void discreteEdge::setBoundVertices(std::vector<discreteVertex*> vertices)
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());
+ ///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;
@@ -196,38 +199,77 @@ void discreteEdge::setBoundVertices(std::vector<discreteVertex*> vertices)
return;
}
+
+/*
+ +---------------+--------------+----------- ... ----------+
+ _pars[0]=0 _pars[1]=1 _pars[2]=2 _pars[N+1]=N+1
+*/
void discreteEdge::parametrize()
{
+
for (int i=0;i<lines.size()+1;i++){
_pars.push_back(i);
- }
-
- /*
- +---------------+--------------+----------- ... ----------+
- _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] );
- // }
-
- printf("dans discrete edge %d line.size =%d \n", this->tag(), lines.size());
-
-// create new MEdge Vertices
-// 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;
+ }
+ //Replace MVertex by MedgeVertex
+ //we need to recreate lines and triangles
+ //that contain those new MEdgeVertices
+ std::map<MVertex*, MVertex*> old2new;
-// we should loop over Mlines and MTrinagles to take those new MEdgeVertices into account
+ std::vector<MVertex*> newVertices;
+ std::vector<MLine*> newLines;
+
+ MVertex *vL = getBeginVertex()->mesh_vertices[0];
+ int i = 0;
+ for(i=0 ; i < lines.size()-1; i++){
+ MVertex *vR ;
+ if (_orientation[i] == 1 ) vR = lines[i]->getVertex(1);
+ else vR = lines[i]->getVertex(0);
+ int param = i+1;
+ MVertex *vNEW = new MEdgeVertex(vR->x(),vR->y(),vR->z(), this, param);
+ old2new.insert(std::make_pair(vR,vNEW));
+ newVertices.push_back(vNEW);
+ vNEW->setNum(vR->getNum());
+ newLines.push_back(new MLine(vL, vNEW));
+ _orientation[i] = 1;
+ //printf("i edge =%d , orientation=%d , mesh vertex =%p (%d) newv=%p (%d) \n", i, _orientation[i], vL, vL->getNum(), vNEW, vNEW->getNum());
+ vL = vNEW;
+ }
+ MVertex *vR = getEndVertex()->mesh_vertices[0];
+ newLines.push_back(new MLine(vL, vR));
+ _orientation[i] = 1;
+ //printf("i edge =%d , orientation=%d , mesh vertex =%p (%d) newv=%p (%d) \n", i, _orientation[i], vL, vL->getNum(), vR, vR->getNum());
+
+ mesh_vertices = newVertices;
+ deleteVertexArrays();
+ lines.clear();
+ lines = newLines;
+
+ for(std::list<GFace*>::iterator iFace = l_faces.begin(); iFace != l_faces.end(); ++iFace){
+ std::vector<MTriangle*> newTriangles;
+ for (unsigned int i = 0; i < (*iFace)->triangles.size(); ++i){
+ MTriangle *t = (*iFace)->triangles[i];
+ MVertex *v[3];
+ v[0] = t->getVertex(0);
+ v[1] = t->getVertex(1);
+ v[2] = t->getVertex(2);
+ //printf("old triangle v0=%p (%d) v1=%p (%d) v2=%p (%d) \n",v[0], v[0]->getNum() , v[1],v[1]->getNum() ,v[2], v[2]->getNum());
+ for (int j = 0; j < 3; j++){
+ std::map<MVertex*, MVertex*>::iterator itmap = old2new.find(v[j]);
+ MVertex *vNEW;
+ if (itmap != old2new.end()) {
+ vNEW = itmap->second;
+ v[j]=vNEW;
+ }
+ }
+ //printf(" new triangle v0=%p (%d) v1=%p (%d) v2=%p (%d) \n",v[0], v[0]->getNum() , v[1],v[1]->getNum() ,v[2], v[2]->getNum());
+ newTriangles.push_back(new MTriangle(v[0], v[1], v[2]));
+ }
+ (*iFace)->triangles.clear();
+ (*iFace)->deleteVertexArrays();
+ (*iFace)->triangles = newTriangles;
+ }
// for (int i = 0; i < mesh_vertices.size(); i++){
@@ -235,35 +277,7 @@ void discreteEdge::parametrize()
// 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)
-
-
-// du brol ci-dessous ...
-// std::vector<MLine*> newLines;
-// newLines.push_back(new MLine(v1, newv));
-// delete lines[i];
-// lines = newLines;
-
-// for(int i = 0; i < mesh_vertices.size(); 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;
-// }
-// }
-// }
-// }
-
-
-
-;
+
}
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index c6955da8610078c6bfa90338fe46feb337b2671a..a0cc4eac2b131c520b94200bb6935469b8998d8f 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -27,7 +27,7 @@ class discreteEdge : public GEdge {
virtual Range<double> parBounds(int) const;
void parametrize() ;
void orderMLines() ;
- void setBoundVertices( std::vector<discreteVertex*> vertices );
+ void setBoundVertices();
};
#endif
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index ead70a6610540706030a0e53aecb9ff79d69d9b6..8de241977ef7547f59ed032c9b414a962a1f90df 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -22,57 +22,70 @@ discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
meshStatistics.status = GFace::DONE;
}
-void discreteFace::setBoundEdges(std::vector<discreteEdge*> discr_edges)
-{
+void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge > &map_edges){
+
+ //find the boundary edges
+ 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 (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);
+ else
+ bound_edges.erase(std::find(bound_edges.begin(),bound_edges.end(),tmp_edge));
+ }
+ }
+
+
+ //for the boundary edges, associate the tag of the current discrete face
+ for (std::list<MEdge>::iterator itv = bound_edges.begin() ; itv != bound_edges.end() ; ++itv){
+ std::map<MEdge, std::vector<int> , Less_Edge >::iterator itmap = map_edges.find(*itv);
+ if (itmap == map_edges.end()) {
+ std::vector<int> tagFaces;
+ tagFaces.push_back(this->tag());
+ map_edges.insert(std::make_pair(*itv,tagFaces));
+ }
+ else{
+ std::vector<int> tagFaces = itmap->second;
+ tagFaces.push_back(this->tag());
+ itmap->second = tagFaces;
+ //printf("addface %d %d\n", tagFaces[0], tagFaces[1]);
+ }
+ }
+
+ //printf( "There are %d bound msh edges \n ", map_edges.size());
- 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() );
- }
+}
+
+void discreteFace::setBoundEdges(std::vector<int> tagEdges)
+{
- 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("***** In discrete Face: \n");
+
+
+// for(GModel::eiter it = model()->firstEdge(); it != model()->lastEdge(); ++it){
+// // for (std::vector<discreteEdge*>::iterator it = bound_edges.begin(); it != bound_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( "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());
- //}
+ for (std::vector<int>::iterator it = tagEdges.begin(); it != tagEdges.end(); it++) {
+ GEdge *ge = GModel::current()->getEdgeByTag(abs(*it));
+ l_edges.push_back(ge);
+ l_dirs.push_back(1);
+ ge->addFace(this);
+ //printf("for Face %d add bound edge %d\n",this->tag(), ge->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());
-// }
+ // printf("bound edges =%d \n", edges().size());
-
-
- //exit(1);
}
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index 2e1b4309eafc7da4e46ce3723c8fb86482f144a4..d44f1ca8fc385caedc23644192b469b5940f2b51 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -9,6 +9,7 @@
#include "GModel.h"
#include "GFace.h"
#include "discreteEdge.h"
+#include "MEdge.h"
class discreteFace : public GFace {
public:
@@ -21,7 +22,8 @@ class discreteFace : public GFace {
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
virtual void secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const;
- void setBoundEdges( std::vector<discreteEdge*> edges );
+ void setBoundEdges( std::vector<int> tagEdges );
+ void findEdges( std::map<MEdge, std::vector<int>, Less_Edge > &map_edges);
};
#endif
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index f1de893a8d70a8ef8002280e019567122650b2f1..a61e647fd1482d0bdee43c08d4e1d4a3be0ad4e2 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -380,19 +380,19 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
// here, we will replace edges by their compounds
//printf("***** In meshGFace: \n");
if (gf->geomType() == GEntity::CompoundSurface){
- printf("replace edges by compound lines \n");
+ //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 edge %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 GFaceCompound %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());
}
@@ -458,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);
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index ae9eb85389245d6f984ea173b77bd490ab47bb0e..a4afef5178be07aff05a7572aefaad50dd6de305 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -140,6 +140,8 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
MVertex *v[2];
v[0] = numberedV[out.edgelist[i * 2 + 0] - 1];
v[1] = numberedV[out.edgelist[i * 2 + 1] - 1];
+
+ //implement here the 1D mesh ...
}
// re-create the triangular meshes FIXME: this can lead to hanging
diff --git a/benchmarks/2d/createTopo.msh b/benchmarks/2d/createTopo.msh
new file mode 100644
index 0000000000000000000000000000000000000000..4194b156d74cd74dceabcfd9d7d77a4bca1535a8
--- /dev/null
+++ b/benchmarks/2d/createTopo.msh
@@ -0,0 +1,155 @@
+$MeshFormat
+2 0 8
+$EndMeshFormat
+$Nodes
+58
+1 0 0 0
+2 0 1 0
+3 1 1 0
+4 1 0 0
+5 0 2 0
+6 1 2 0
+7 0.2499999999993998 1 0
+8 0.4999999999990217 1 0
+9 0.7499999999993374 1 0
+10 1 0.7500000000003469 0
+11 1 0.5000000000013878 0
+12 1 0.2500000000010407 0
+13 0.7500000000003469 0 0
+14 0.5000000000013878 0 0
+15 0.2500000000010407 0 0
+16 0 0.2499999999993998 0
+17 0 0.4999999999990217 0
+18 0 0.7499999999993374 0
+19 1 1.2 0
+20 1 1.4 0
+21 1 1.6 0
+22 1 1.8 0
+23 0.7500000000003469 2 0
+24 0.5000000000013878 2 0
+25 0.2500000000010407 2 0
+26 0 1.8 0
+27 0 1.6 0
+28 0 1.4 0
+29 0 1.2 0
+30 0.4838453815992497 0.7001134502703265 0
+31 0.2137219674263763 0.6255840250137896 0
+32 0.7964751095545822 0.3310225115184635 0
+33 0.6380768712362025 0.1712066244048212 0
+34 0.5423577861200293 0.4221404297864205 0
+35 0.2491174661656623 0.3998685906427588 0
+36 0.3919626856685218 0.1855558909268056 0
+37 0.7968308623417641 0.581754889093215 0
+38 0.8369103961582264 0.1504458271848651 0
+39 0.7631804184148582 0.8313434174180454 0
+40 0.178216030367045 0.1670848963137928 0
+41 0.3874034021894878 0.5516464588425142 0
+42 0.2412612248373413 0.845949579213909 0
+43 0.1892871905708922 1.696788567649282 0
+44 0.8107128094292139 1.303211432350718 0
+45 0.7839370157560627 1.556826873049367 0
+46 0.2169163731573559 1.442796305409011 0
+47 0.3493592037515919 1.834522414336105 0
+48 0.6506407962486326 1.165477585663895 0
+49 0.1625751921142504 1.856404308390023 0
+50 0.8374248078858939 1.143595691609977 0
+51 0.3972615131649327 1.21833029569987 0
+52 0.6036364124391213 1.781243294253888 0
+53 0.6042108386940316 1.374669587758423 0
+54 0.3970642534924195 1.625260488761389 0
+55 0.4340625840563876 1.445368591067763 0
+56 0.5700944140718672 1.552093283630357 0
+57 0.8562622380325885 1.789678361217209 0
+58 0.1440296477202814 1.21018776685148 0
+$EndNodes
+$Elements
+88
+37 2 3 0 10 0 31 17 18
+38 2 3 0 10 0 32 11 12
+39 2 3 0 10 0 34 32 33
+40 2 3 0 10 0 36 14 15
+41 2 3 0 10 0 14 36 33
+42 2 3 0 10 0 33 13 14
+43 2 3 0 10 0 37 10 11
+44 2 3 0 10 0 34 33 36
+45 2 3 0 10 0 35 17 31
+46 2 3 0 10 0 32 37 11
+47 2 3 0 10 0 32 34 37
+48 2 3 0 10 0 37 34 30
+49 2 3 0 10 0 13 38 4
+50 2 3 0 10 0 12 38 32
+51 2 3 0 10 0 38 12 4
+52 2 3 0 10 0 39 9 3
+53 2 3 0 10 0 10 39 3
+54 2 3 0 10 0 33 38 13
+55 2 3 0 10 0 33 32 38
+56 2 3 0 10 0 34 36 35
+57 2 3 0 10 0 39 37 30
+58 2 3 0 10 0 37 39 10
+59 2 3 0 10 0 39 8 9
+60 2 3 0 10 0 30 8 39
+61 2 3 0 10 0 16 40 1
+62 2 3 0 10 0 40 15 1
+63 2 3 0 10 0 36 40 35
+64 2 3 0 10 0 40 36 15
+65 2 3 0 10 0 30 41 31
+66 2 3 0 10 0 35 41 34
+67 2 3 0 10 0 41 35 31
+68 2 3 0 10 0 41 30 34
+69 2 3 0 10 0 35 16 17
+70 2 3 0 10 0 16 35 40
+71 2 3 0 10 0 7 42 2
+72 2 3 0 10 0 18 42 31
+73 2 3 0 10 0 42 18 2
+74 2 3 0 10 0 42 7 8
+75 2 3 0 10 0 30 42 8
+76 2 3 0 10 0 42 30 31
+77 2 3 0 20 0 24 25 47
+78 2 3 0 20 0 8 9 48
+79 2 3 0 20 0 46 43 27
+80 2 3 0 20 0 45 44 20
+81 2 3 0 20 0 20 44 19
+82 2 3 0 20 0 27 43 26
+83 2 3 0 20 0 21 45 20
+84 2 3 0 20 0 28 46 27
+85 2 3 0 20 0 49 5 26
+86 2 3 0 20 0 25 5 49
+87 2 3 0 20 0 9 3 50
+88 2 3 0 20 0 50 3 19
+89 2 3 0 20 0 47 52 24
+90 2 3 0 20 0 48 51 8
+91 2 3 0 20 0 26 43 49
+92 2 3 0 20 0 49 43 47
+93 2 3 0 20 0 47 25 49
+94 2 3 0 20 0 19 44 50
+95 2 3 0 20 0 50 44 48
+96 2 3 0 20 0 48 9 50
+97 2 3 0 20 0 52 23 24
+98 2 3 0 20 0 51 7 8
+99 2 3 0 20 0 53 51 48
+100 2 3 0 20 0 54 52 47
+101 2 3 0 20 0 53 48 44
+102 2 3 0 20 0 44 45 53
+103 2 3 0 20 0 56 53 45
+104 2 3 0 20 0 54 47 43
+105 2 3 0 20 0 43 46 54
+106 2 3 0 20 0 55 54 46
+107 2 3 0 20 0 53 56 55
+108 2 3 0 20 0 56 54 55
+109 2 3 0 20 0 46 51 55
+110 2 3 0 20 0 51 53 55
+111 2 3 0 20 0 45 52 56
+112 2 3 0 20 0 52 54 56
+113 2 3 0 20 0 57 6 23
+114 2 3 0 20 0 22 6 57
+115 2 3 0 20 0 58 2 7
+116 2 3 0 20 0 29 2 58
+117 2 3 0 20 0 57 23 52
+118 2 3 0 20 0 52 45 57
+119 2 3 0 20 0 58 7 51
+120 2 3 0 20 0 51 46 58
+121 2 3 0 20 0 22 57 21
+122 2 3 0 20 0 21 57 45
+123 2 3 0 20 0 29 58 28
+124 2 3 0 20 0 28 58 46
+$EndElements
diff --git a/benchmarks/2d/createTopo_GEO.geo b/benchmarks/2d/createTopo_GEO.geo
new file mode 100644
index 0000000000000000000000000000000000000000..27b48e73392d667167385d56d543032ae32a63f6
--- /dev/null
+++ b/benchmarks/2d/createTopo_GEO.geo
@@ -0,0 +1,11 @@
+Mesh.CharacteristicLengthFactor=4;
+
+Merge "createTopo.msh";
+CreateTopology;
+
+//Compound Line(150)={1};
+//Compound Line(160)={2,3,4};
+//Compound Line(170)={5,6,7};
+
+Compound Surface(100)={10} Boundary {{}};
+Compound Surface(200)={20} Boundary {{}};
\ No newline at end of file