*** empty log message ***

Mesh/BDS.cpp

-// $Id: BDS.cpp,v 1.81 2007-10-11 14:34:04 remacle Exp $
+// $Id: BDS.cpp,v 1.82 2007-10-14 17:30:42 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -645,6 +645,131 @@ void BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
 }
 
 
+void BDS_Mesh::saturate_edge(BDS_Edge * e, std::vector<BDS_Point *> &mids)
+{
+  BDS_Point *op[2];
+  BDS_Point *p1 = e->p1;
+  BDS_Point *p2 = e->p2;
+
+  BDS_Point *pts1[4];
+  e->faces(0)->getNodes(pts1);
+
+  int orientation = 0;
+  for(int i = 0; i < 3; i++) {
+    if(pts1[i] == p1) {
+      if(pts1[(i + 1) % 3] == p2)
+        orientation = 1;
+      else
+        orientation = -1;
+
+      break;
+    }
+  }
+
+  //  printf("sat %d\n",mids.size());
+
+  // we should project 
+  e->oppositeof(op);
+  BDS_GeomEntity *g1 = 0, *g2 = 0, *ge = e->g;
+
+  BDS_Edge *p1_op1 = find_edge(p1, op[0], e->faces(0));
+  BDS_Edge *op1_p2 = find_edge(op[0], p2, e->faces(0));
+  BDS_Edge *p1_op2 = find_edge(p1, op[1], e->faces(1));
+  BDS_Edge *op2_p2 = find_edge(op[1], p2, e->faces(1));
+
+  if(e->faces(0)) {
+    g1 = e->faces(0)->g;
+    del_face(e->faces(0));
+  }
+  // not a bug !!!
+  if(e->faces(0)) {
+    g2 = e->faces(0)->g;
+    del_face(e->faces(0));
+  }
+
+  //  double t_l = e->target_length;
+
+  del_edge(e);
+
+  // create all the sub-edges of e 
+  std::vector<BDS_Edge*> subs;
+  for (int i=0; i<mids.size()+1; i++)
+    {
+      BDS_Point *a,*b;
+      if (i == 0)a = p1;
+      else a = mids[i-1];
+      if (i == mids.size())b = p2;
+      else b = mids[i];
+      BDS_Edge *sub = new BDS_Edge(a, b);
+      //      printf("%g %g -> %g %g\n",a->X,a->Y,b->X,b->Y);
+      edges.push_back(sub); 
+      subs.push_back(sub);
+      sub->g = ge;
+    }
+  // create edges that connect op1 and op2
+  std::vector<BDS_Edge*> conn1,conn2;
+  for (int i=0;i<mids.size();i++)
+    {
+      BDS_Edge *c1 = new BDS_Edge(mids[i], op[0]);
+      edges.push_back(c1); 
+      conn1.push_back(c1);
+      BDS_Edge *c2 = new BDS_Edge(mids[i], op[1]);
+      edges.push_back(c2); 
+      conn2.push_back(c2);
+      c1->g = g1;
+      c2->g = g2;
+      mids[i]->g = ge;
+    }
+
+  // create the triangles
+
+  for (int i=0;i<mids.size()+1;i++)
+    {
+      BDS_Edge *e1,*e2,*e3=subs[i];
+
+      //      printf("--> %d %g %g ->%d  %g %g\n",e3->p1->iD,e3->p1->X,e3->p1->Y,e3->p2->iD,e3->p2->X,e3->p2->Y);
+
+
+      if (i==0)e1 = p1_op1;
+      else e1 = conn1[i-1];
+      if (i==mids.size())e2 = op1_p2;
+      else e2 = conn1[i];
+
+//       printf("--> %d %g %g ->%d  %g %g\n",e1->p1->iD,e1->p1->X,e1->p1->Y,e1->p2->iD,e1->p2->X,e1->p2->Y);
+//       printf("--> %d %g %g ->%d  %g %g\n",e2->p1->iD,e2->p1->X,e2->p1->Y,e2->p2->iD,e2->p2->X,e2->p2->Y);
+
+      BDS_Face *t1;
+      if(orientation == 1)
+	t1 = new BDS_Face(e3, e2, e1);
+      else 
+	t1 = new BDS_Face(e3, e1, e2);
+
+
+      if (i==0)e1 = p1_op2;
+      else e1 = conn2[i-1];
+      if (i==mids.size())e2 = op2_p2;
+      else e2 = conn2[i];
+
+      BDS_Face *t2;
+      if(orientation != 1)
+	t2 = new BDS_Face(e3, e2, e1);
+      else 
+	t2 = new BDS_Face(e3, e1, e2);
+      t1->g = g1;
+      t2->g = g2;
+
+      triangles.push_back(t1);
+      triangles.push_back(t2);
+    }
+  // config has changed
+  p1->config_modified = true;
+  p2->config_modified = true;
+  op[0]->config_modified = true;
+  op[1]->config_modified = true;
+}
+
+
+
 
 // This function does actually the swap without taking into account
 // the feasability of the operation. Those conditions have to be
@@ -682,6 +807,7 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
 
   // we test if the edge is deleted
   //  return false;
+
   if(e->deleted)
     return false;
   
@@ -692,6 +818,7 @@ bool BDS_Mesh::swap_edge(BDS_Edge * e, const BDS_SwapEdgeTest &theTest)
   if(e->g && e->g->classif_degree == 1)
     return false;
 
+
   BDS_Point *op[2];
   BDS_Point *p1 = e->p1;
   BDS_Point *p2 = e->p2;

Mesh/BDS.h

@@ -457,6 +457,7 @@ public:
   bool smooth_point_centroid(BDS_Point * p, GFace *gf);
   bool move_point(BDS_Point *p , double X, double Y, double Z);
   void split_edge(BDS_Edge *, BDS_Point *);
+  void saturate_edge(BDS_Edge *, std::vector<BDS_Point *>&);
   bool edge_constraint    ( BDS_Point *p1, BDS_Point *p2 );
   bool recombine_edge    ( BDS_Edge *e );
   // Global operators

Mesh/meshGFace.cpp

-// $Id: meshGFace.cpp,v 1.96 2007-10-11 16:08:23 remacle Exp $
+// $Id: meshGFace.cpp,v 1.97 2007-10-14 17:30:42 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -255,6 +255,14 @@ double NewGetLc(BDS_Edge *e, GFace *f)
   return 2*linearLength / (l1 + l2);
 }
 
+double NewGetLc(BDS_Point *p1,BDS_Point *p2, GFace *f)
+{
+  double linearLength = computeEdgeLinearLength(p1,p2,f);
+  double l1 = NewGetLc(p1);
+  double l2 = NewGetLc(p2);
+  return 2*linearLength / (l1 + l2);
+}
+
 bool edgeSwapTest(BDS_Edge *e,GFace *gf)
 {
   BDS_Point *op[2];
@@ -318,9 +326,15 @@ bool edgeSwapTestDelaunay(BDS_Edge *e,GFace *gf)
   double p2x[3] =  {e->p2->X,e->p2->Y,e->p2->Z};
   double op1x[3] = {op[0]->X,op[0]->Y,op[0]->Z};
   double op2x[3] = {op[1]->X,op[1]->Y,op[1]->Z};
+  double p1u[2] =  {e->p1->u,e->p1->v};
+  double p2u[2] =  {e->p2->u,e->p2->v};
+  double op1u[2] = {op[0]->u,op[0]->v};
+  double op2u[2] = {op[1]->u,op[1]->v};
   double fourth[3];
   fourthPoint(p1x,p2x,op1x,fourth);
   double result = gmsh::insphere(p1x, p2x, op1x, fourth, op2x) * gmsh::orient3d(p1x, p2x, op1x, fourth);  
+    //double result = gmsh::incircle(p1u, p2u, op1u, op2u) * gmsh::orient2d(p1u, p2u, op1u);    
+  //  printf("result = a%12.5E\n",result);
   return result > 0.;
 }
 
@@ -404,11 +418,13 @@ void swapEdgePass ( GFace *gf, BDS_Mesh &m, int &nb_swap )
 	{
 	  int result = edgeSwapTestQuality(*it,5);
 	  if (result >= 0)
-	    if(edgeSwapTestDelaunay(*it,gf) || result > 0)
-	      if (m.swap_edge ( *it , BDS_SwapEdgeTestParametric()))
-		nb_swap++;
-	  ++it;
+	    {
+	      if(edgeSwapTestDelaunay(*it,gf) || result > 0)
+		if (m.swap_edge ( *it , BDS_SwapEdgeTestParametric()))
+		  nb_swap++;
+	    }
 	}
+      ++it;
     }  
 }
 
@@ -425,13 +441,18 @@ void splitEdgePass ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
 	  double lone = NewGetLc ( *it,gf);
 	  if ((*it)->numfaces() == 2 && (lone >  MAXE_))
 	    {
+
+//  	      BDS_Point *op[2];
+//  	      (*it)->oppositeof(op);
+//  	      double lone1 = NewGetLc ( op[0],mid,gf);	      
+//  	      double lone2 = NewGetLc ( op[1],mid,gf);	      
+	      
 	      const double coord = 0.5;
 	      BDS_Point *mid ;
 	      mid  = m.add_point(++m.MAXPOINTNUMBER,
 				 coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
 				 coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v,gf);
-	      double l1;
-	      //		  if (BGMExists())
+
 	      mid->lcBGM() = BGM_MeshSize(gf,
 					  (coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u)*m.scalingU,
 					  (coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v)*m.scalingV,
@@ -440,6 +461,52 @@ void splitEdgePass ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
 	      mid->lc() = 0.5 * ( (*it)->p1->lc() +  (*it)->p2->lc() );		  
 	      m.split_edge ( *it, mid );
 	      nb_split++;
+	    }
+	}
+      ++it;
+    }
+}
+
+void saturateEdgePass ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
+{
+  int NN1 = m.edges.size();
+  int NN2 = 0;
+  std::list<BDS_Edge*>::iterator it = m.edges.begin();
+  while (1)
+    {
+      if (NN2++ >= NN1)break;
+      if (!(*it)->deleted)
+	{
+	  double lone = NewGetLc ( *it,gf);
+	  if ((*it)->numfaces() == 2 && (lone >  MAXE_))
+	    {
+	      int nbSub = (int) (lone/MAXE_) ;
+	      //	      nbSub = std::min(nbSub,2);
+	      //	      printf("%d %g\n",nbSub,lone/MAXE_);
+	      std::vector<BDS_Point*> mids;
+	      for (int i=0;i<nbSub;i++)
+		{
+		  const double coord = (double)(i+1)/(nbSub+1);
+		  BDS_Point *mid ;
+		  mid  = m.add_point(++m.MAXPOINTNUMBER,
+				     (1.-coord) * (*it)->p1->u + coord * (*it)->p2->u,
+				     (1.-coord) * (*it)->p1->v + coord * (*it)->p2->v,gf);
+		  double l1;
+		  //		  if (BGMExists())
+		  mid->lcBGM() = BGM_MeshSize(gf,
+					      ((1.-coord) * (*it)->p1->u + (coord) * (*it)->p2->u)*m.scalingU,
+					      ((1.-coord) * (*it)->p1->v + (coord) * (*it)->p2->v)*m.scalingV,
+					      mid->X,mid->Y,mid->Z);
+		  //mid->lc() = 2./ ( 1./(*it)->p1->lc() +  1./(*it)->p2->lc() );		  
+		  mid->lc() = ( (1.-coord)*(*it)->p1->lc() +  coord*(*it)->p2->lc() );		 
+		  mids.push_back(mid);
+		  //		  printf("new point %g %g lc %g\n",mid->X,mid->Y,mid->lc());
+		}
+	      //	      printf("saturating an edge with %d points %d triangles\n",mids.size(),m.triangles.size());
+	      if(nbSub>0)m.saturate_edge ( *it, mids );
+	      //	      printf("-> %d triangles\n",m.triangles.size());
+	      nb_split++;
+	      //	      if (nb_split == )break;
 	    } 
 	}
       ++it;
@@ -518,7 +585,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
 
   double t_spl=0, t_sw=0,t_col=0,t_sm=0;
 
-  const double MINE_ = 0.7, MAXE_=1.4;
+  const double MINE_ = 0.67, MAXE_=1.4;
   while (1)
     {
       // we count the number of local mesh modifs.
@@ -558,6 +625,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
       double minE = MINE_;//std::min(MINE_,minL * 1.2);
       clock_t t1 = clock();
       splitEdgePass ( gf, m, maxE, nb_split);
+      //saturateEdgePass ( gf, m, maxE, nb_split);
       clock_t t2 = clock();
       swapEdgePass ( gf, m, nb_swap);
       clock_t t3 = clock();
@@ -567,12 +635,15 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
       clock_t t5 = clock();
       smoothVertexPass ( gf, m, nb_smooth);
       clock_t t6 = clock();
+      swapEdgePass ( gf, m, nb_swap);
+      clock_t t7 = clock();
       // clean up the mesh
 
 
       t_spl += (double)(t2-t1)/CLOCKS_PER_SEC;
       t_sw  += (double)(t3-t2)/CLOCKS_PER_SEC;
       t_sw  += (double)(t5-t4)/CLOCKS_PER_SEC;
+      t_sw  += (double)(t7-t6)/CLOCKS_PER_SEC;
       t_col += (double)(t4-t3)/CLOCKS_PER_SEC;
       t_sm  += (double)(t6-t5)/CLOCKS_PER_SEC;