From 86a6ce05413d22691e7b6e711c86ae77602a8d32 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Tue, 23 Dec 2008 10:39:42 +0000
Subject: [PATCH] *** empty log message ***

---
 Geo/GFaceCompound.cpp   | 157 +++++++++++++++++++++++-----------------
 Geo/GFaceCompound.h     |  22 ++++--
 Geo/MVertex.cpp         |   9 +++
 Mesh/Refine.cpp         |  20 +++--
 benchmarks/3d/Torus.geo |   1 +
 5 files changed, 131 insertions(+), 78 deletions(-)

diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 726ab5289d..749d61b376 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -66,12 +66,22 @@ void GFaceCompound::parametrize() const
 
 void GFaceCompound::getBoundingEdges()
 {
-  if (_U0.size() && !_U1.size()){
+  if (_U0.size() && !_V1.size()){
     std::list<GEdge*> :: const_iterator it = _U0.begin();
     for ( ; it != _U0.end() ; ++it){
       l_edges.push_back(*it);
       (*it)->addFace(this);
     }
+    it = _U1.begin();
+    for ( ; it != _U1.end() ; ++it){
+      l_edges.push_back(*it);
+      (*it)->addFace(this);
+    }
+    it = _V0.begin();
+    for ( ; it != _V0.end() ; ++it){
+      l_edges.push_back(*it);
+      (*it)->addFace(this);
+    }
     return;
   }
 
@@ -88,11 +98,11 @@ void GFaceCompound::getBoundingEdges()
   it = _compound.begin();
   for ( ; it != _compound.end(); ++it){
     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);
-    }
-  }
+    std::list<GEdge*> :: iterator ite = ed.begin();
+    for ( ; ite != ed.end() ; ++ite){
+      if (!(*ite)->degenerate(0) && _touched.count(*ite) == 1)_unique.insert(*ite);
+    }    
+  }    
 
   std::set<GEdge*>::iterator itf = _unique.begin();
   for ( ; itf != _unique.end(); ++itf){
@@ -107,6 +117,11 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
   : GFace(m, tag), _compound(compound), _U0(U0), _U1(U1), _V0(V0), _V1(V1), oct(0)
 {
   getBoundingEdges();
+  if (!_U0.size()) _type = UNITCIRCLE;
+  else if (!_V0.size()) _type = UNITCIRCLE;
+  else if (!_U1.size()) _type = CYLINDER;
+  else if (!_V1.size()) _type = BIFURCATION;
+  else _type = SQUARE;
 }
 
 GFaceCompound::~GFaceCompound()
@@ -193,7 +208,7 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
   gmshGradientBasedDiffusivity diffusivity(coordinates);
   if (ITER > 0) diffusivity.setComponent(_isU);
   
-  if (!_U1.size()){
+  if (_type == UNITCIRCLE){
     // maps the boundary onto a circle
     std::vector<MVertex*> ordered;
     std::vector<double> coords;  
@@ -209,6 +224,19 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
       else myAssembler.fixVertex(v, 0, 1, sin(theta));
     }
   }
+  else if (_type == CYLINDER){
+    // maps the boundary onto an annulus
+    std::vector<MVertex*> ordered;
+    std::vector<double> coords;  
+    bool success = orderVertices (_U0, ordered, coords);
+    if (!success)throw;
+    for (unsigned int i = 0; i < ordered.size(); i++){
+      MVertex *v = ordered[i];
+      const double theta = 2 * M_PI * coords[i];
+      if (_isU) myAssembler.fixVertex(v, 0, 1, cos(theta));
+      else myAssembler.fixVertex(v, 0, 1, sin(theta));
+    }
+  }
   else{
     if (_isU){
       std::list<GEdge*> :: const_iterator it = _U0.begin();
@@ -258,28 +286,29 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
   Msg::Debug("Assembly done");
   lsys.systemSolve();
   Msg::Debug("System solved");
-  
   it = _compound.begin();
-  for ( ; it != _compound.end(); ++it){
+  std::set<MVertex *>allNodes;
+  for ( ; it != _compound.end() ; ++it){
     for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
       MTriangle *t = (*it)->triangles[i];
-      double uu[3], vv[3];
       for (int j = 0; j < 3; j++){
-	MVertex *v = t->getVertex(j);
-	double value = myAssembler.getDofValue(v,0,1);
-	std::map<MVertex*,SPoint2>::iterator itf = coordinates.find(v);
-	if (itf == coordinates.end()){
-	  SPoint2 p (_isU ? value : 0.0, _isU ? 0.0 : value);
-	  coordinates[v] = p;
-	}
-	else{
-	  if(_isU) itf->second[0]= value;
-	  else itf->second[1]= value;
-	  uu[j] = itf->second[0];
-	  vv[j] = itf->second[1];
-	}
-      }      
-    }   
+	allNodes.insert(t->getVertex(j));
+      }
+    }
+  }
+  
+  for (std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+    MVertex *v = *itv;
+    double value = myAssembler.getDofValue(v,0,1);
+    std::map<MVertex*,SPoint2>::iterator itf = coordinates.find(v);
+    if (itf == coordinates.end()){
+      SPoint2 p (_isU ? value:0.0,_isU ? 0.0:value);
+      coordinates[v] = p;
+    }
+    else{
+      if(_isU) itf->second[0]= value;
+      else itf->second[1]= value;	  
+    }      
   }
 }
 
@@ -317,7 +346,19 @@ double GFaceCompound::curvature(const SPoint2 &param) const
   if (!lt){
     return  0.0;
   }
-  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);
 }
 
 double GFaceCompound::curvature(MTriangle *t) const
@@ -328,7 +369,7 @@ double GFaceCompound::curvature(MTriangle *t) const
   double val[9] = {n1.x(),n2.x(),n3.x(),
 		   n1.y(),n2.y(),n3.y(),
 		   n1.z(),n2.z(),n3.z()};
-  return fabs(t->interpolateDiv (val,0,0,0.));
+  //  return fabs(t->interpolateDiv (val,0.,0.,0.));
 }
 
 GPoint GFaceCompound::point(double par1, double par2) const
@@ -339,13 +380,18 @@ GPoint GFaceCompound::point(double par1, double par2) const
   getTriangle (par1, par2, &lt, 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); 
-    return GPoint(p.x(), p.y(), p.z(), this);
+    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);
+  }
+  if (0 && lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
+    SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
+    return lt->gf->point(pv.x(),pv.y());
   }
-  p = lt->v1 * (1. - U - V) + lt->v2 * U + lt->v3 * V;
-  double par[2] = {par1, par2};
-  return GPoint(p.x(), p.y(), p.z(), this, par);
+  
+  p = lt->v1*(1.-U-V) + lt->v2*U + lt->v3*V;
+  double par[2] = {par1,par2};
+  return GPoint(p.x(),p.y(),p.z(),this,par);
 }
 
 Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 &param) const
@@ -353,12 +399,7 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 &param) const
   parametrize();
   double U,V,UDU,UDV,VDU,VDV;
   GFaceCompoundTriangle *lt;
-  // GFaceCompoundTriangle *ltdu;
-  // GFaceCompoundTriangle *ltdv;
-  getTriangle(param.x(), param.y(), &lt, U,V);
-  // getTriangle (param.x()+1.e-4, param.y(), &ltdu, UDU,VDU);
-  // getTriangle (param.x(), param.y()+1.e-4, &ltdv, UDV,VDV);
-
+  getTriangle (param.x(), param.y(), &lt, U,V);
   if (!lt)
     return Pair<SVector3, SVector3>(SVector3(1, 0, 0), SVector3(0, 1, 0));
 
@@ -367,29 +408,11 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 &param) const
   double inv[2][2];
   inv2x2(mat,inv);
  
-  // MVertex *v0 = lt->t->getVertex(0);
-  // MVertex *v1 = lt->t->getVertex(1);
-  // MVertex *v2 = lt->t->getVertex(2);
-
-  // SPoint3 p(0,0,0); 
-  // SPoint3 pdu(0,0,0); 
-  // SPoint3 pdv(0,0,0); 
-
-  // lt->t->pnt(U,V,0,p);
-  // ltdu->t->pnt(UDU,VDU,0,pdu);
-  // ltdv->t->pnt(UDV,VDV,0,pdv);
- 
   SVector3 dXdxi   (lt->v2 - lt->v1);
   SVector3 dXdeta  (lt->v3 - lt->v1);
 
-  // return Pair<SVector3, SVector3>(dXdxi * inv[0][0] + dXdeta * inv[1][0],
-  //				     dXdxi * inv[0][1] + dXdeta * inv[1][1]);
-
-  SVector3 dXdu (dXdxi * inv[0][0] + dXdeta * inv[1][0]);
-  SVector3 dXdv (dXdxi * inv[0][1] + dXdeta * inv[1][1]);
-  
-  // SVector3 dXduFD(SVector3(pdu - p) * 1.e4);
-  // SVector3 dXdvFD(SVector3(pdv - p) * 1.e4);
+  SVector3 dXdu (dXdxi*inv[0][0]+dXdeta*inv[1][0]);
+  SVector3 dXdv (dXdxi*inv[0][1]+dXdeta*inv[1][1]);
 
   return Pair<SVector3, SVector3>(dXdu,dXdv);
 } 
@@ -520,13 +543,15 @@ void GFaceCompound::buildOct() const
       _gfct[count].p1 = it0->second;
       _gfct[count].p2 = it1->second;
       _gfct[count].p3 = it2->second;
-      _gfct[count].t = t;      
-      _gfct[count].v1 = SPoint3
-        (t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z());      
-      _gfct[count].v2 = SPoint3
-        (t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z());      
-      _gfct[count].v3 = SPoint3
-        (t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z());      
+      if ((*it)->geomType() != GEntity::DiscreteSurface){
+	reparamMeshVertexOnFace(t->getVertex(0), *it, _gfct[count].gfp1); 
+	reparamMeshVertexOnFace(t->getVertex(1), *it, _gfct[count].gfp2); 
+	reparamMeshVertexOnFace(t->getVertex(2), *it, _gfct[count].gfp3); 
+      }
+      _gfct[count].v1 = SPoint3(t->getVertex(0)->x(),t->getVertex(0)->y(),t->getVertex(0)->z());      
+      _gfct[count].v2 = SPoint3(t->getVertex(1)->x(),t->getVertex(1)->y(),t->getVertex(1)->z());      
+      _gfct[count].v3 = SPoint3(t->getVertex(2)->x(),t->getVertex(2)->y(),t->getVertex(2)->z());      
+      _gfct[count].gf = *it;      
       fprintf(xyzu,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
 	      t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
 	      t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 0be12e1bab..6b88f6c213 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -26,11 +26,15 @@ The compound can therefore be re-meshed using any surface mesh
 generator of gmsh!
 */
 
-typedef struct {
+class  GFaceCompoundTriangle {
+public:
   SPoint2 p1, p2, p3;
+  SPoint2 gfp1, gfp2, gfp3;
   SPoint3 v1, v2, v3;
-  MTriangle *t;
-} GFaceCompoundTriangle;
+  GFace *gf;
+  GFaceCompoundTriangle () : gf(0)
+  {}
+} ;
 
 class Octree;
 
@@ -51,9 +55,13 @@ class GFaceCompound : public GFace {
                    double &_u, double &_v) const;
   virtual double curvature(MTriangle *t) const;
 public:
-  GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
-		std::list<GEdge*> &U0, std::list<GEdge*> &U1,
-		std::list<GEdge*> &V0, std::list<GEdge*> &V1);
+  typedef enum {UNITCIRCLE, CYLINDER, BIFURCATION, SQUARE} typeOfIsomorphism;
+  GFaceCompound(GModel *m, int tag, 
+		std::list<GFace*> &compound,
+		std::list<GEdge*> &U0,
+		std::list<GEdge*> &U1,
+		std::list<GEdge*> &V0,
+		std::list<GEdge*> &V1);
   virtual ~GFaceCompound();
   Range<double> parBounds(int i) const { return Range<double>(0, 1); } 
   virtual GPoint point(double par1, double par2) const; 
@@ -64,6 +72,8 @@ public:
   SPoint2 getCoordinates(MVertex *v) const { parametrize() ; return coordinates[v]; }
   virtual bool buildRepresentationCross(){ return false; }
   virtual double curvature(const SPoint2 &param) const;
+private:
+  typeOfIsomorphism _type;
 };
 
 #endif
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 9659fbcfa0..23881b11f4 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -208,6 +208,14 @@ MVertex::linearSearch(std::set<MVertex*, MVertexLessThanLexicographic> &pos)
 static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> &params)
 {
   params.clear();
+
+  if (gf->geomType() == GEntity::CompoundSurface &&
+      v->onWhat()->dim() < 2){
+    GFaceCompound *gfc = (GFaceCompound*) gf;
+    params.push_back(gfc->getCoordinates(v));
+    return;
+  }
+
   if(v->onWhat()->dim() == 0){
     GVertex *gv = (GVertex*)v->onWhat();
     std::list<GEdge*> ed = gv->edges();
@@ -250,6 +258,7 @@ static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> &params
 bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf, 
                            SPoint2 &param1, SPoint2 &param2)
 {
+
   std::vector<SPoint2> p1, p2;
   getAllParameters(v1, gf, p1);
   getAllParameters(v2, gf, p2);
diff --git a/Mesh/Refine.cpp b/Mesh/Refine.cpp
index 36ebf2dc15..040b80d931 100644
--- a/Mesh/Refine.cpp
+++ b/Mesh/Refine.cpp
@@ -86,13 +86,21 @@ static void Subdivide(GFace *gf, bool splitTrianglesIntoQuads)
       MTriangle *t = gf->triangles[i];
       if(t->getNumVertices() == 6){
 	SPoint2 pt, temp;
-	for(int k = 0; k < 6; k++){
-	  reparamMeshVertexOnFace(t->getVertex(k), gf, temp);
-	  pt[0] += temp[0] / 6.;
-	  pt[1] += temp[1] / 6.;
+	SPoint3 ptx; t->pnt(0.5,0.5,0,ptx);
+	bool reparamOK = true;
+	for(int k = 0; k<6; k++){
+	  reparamOK &= reparamMeshVertexOnFace(t->getVertex(k), gf, temp);
+	  pt[0] += temp[0]/6.;
+	  pt[1] += temp[1]/6.;
+	}
+	MVertex *newv;
+	if (reparamOK){
+	  GPoint gp = gf->point(pt);		
+	  newv = new MFaceVertex (gp.x(),gp.y(),gp.z(),gf,pt[0],pt[1]);
+	}
+	else {
+	  newv = new MVertex (ptx.x(),ptx.y(),ptx.z(),gf);
 	}
-	GPoint gp = gf->point(pt);	
-	MFaceVertex *newv = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, pt[0], pt[1]);
 	gf->mesh_vertices.push_back(newv);
 	quadrangles2.push_back
 	  (new MQuadrangle(t->getVertex(0), t->getVertex(3), newv, t->getVertex(5)));
diff --git a/benchmarks/3d/Torus.geo b/benchmarks/3d/Torus.geo
index e878d86047..5052cc330c 100644
--- a/benchmarks/3d/Torus.geo
+++ b/benchmarks/3d/Torus.geo
@@ -13,3 +13,4 @@ Line Loop(5) = {4,1,2,3};
 Plane Surface(6) = {5};         
        
 Extrude Surface{6, {0.0,1,0}, {0,0.0,0.0}, 1*3.14159/2};         
+Recombine Surface {6, 27, 23, 15, 19, 28};
-- 
GitLab