*** empty log message ***

Mesh/meshGRegion.cpp

@@ -162,6 +162,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, std::vector<MV
 	      //	      SPoint2 pp = gf->parFromPoint(SPoint3 (v[j]->x(),v[j]->y(),v[j]->z()));
 	      SPoint2 pp (0,0);
 	      MFaceVertex *v1b = new MFaceVertex (v[j]->x(),v[j]->y(),v[j]->z(),gf,pp.x(),pp.y() );
+	      gf->mesh_vertices.push_back(v1b);
 	      numberedV[out.trifacelist[i * 3 + j] -1] = v1b;
 	      delete v[j];
 	      v[j] = v1b;

Mesh/meshGRegionDelaunayInsertion.cpp

@@ -3,6 +3,7 @@
 #include "Numeric.h"
 #include "Message.h"
 #include <set>
+#include <map>
 #include <algorithm>
 
 int MTet4::inCircumSphere ( const double *p ) const
@@ -104,11 +105,10 @@ void recurFindCavity ( std::list<faceXtet> & shell,
     }
 }
 
-int III;
-
 bool insertVertex (MVertex *v , 
 		   MTet4   *t ,
-		   std::set<MTet4*,compareTet4Ptr> &allTets)
+		   std::set<MTet4*,compareTet4Ptr> &allTets,
+		   std::vector<double> & vSizes)
 {
   std::list<faceXtet>  shell;
   std::list<MTet4*>  cavity; 
@@ -132,7 +132,6 @@ bool insertVertex (MVertex *v ,
   double oldVolume = 0;
 
   //  char name2[245];
-  //  sprintf(name2,"testv%d.pos",III);
   //  FILE *ff2 = fopen (name2,"w");
   //  fprintf(ff2,"View\"test\"{\n");
 
@@ -198,7 +197,7 @@ bool insertVertex (MVertex *v ,
 // 		   it->v[2]->y(),
 // 		   it->v[2]->z());
       
-      MTet4 *t4 = new MTet4 ( t ); 
+      MTet4 *t4 = new MTet4 ( t , vSizes); 
       newTets[k++]=t4;
       // all new tets are pushed front in order to
       // ba able to destroy them if the cavity is not
@@ -242,45 +241,58 @@ bool insertVertex (MVertex *v ,
       return false;
     }
 }
+
+static void setLcs ( MTetrahedron *t, std::map<MVertex*,double> &vSizes)
+{
+  for (int i=0;i<4;i++)
+    {
+      for (int j=i+1;j<4;j++)
+	{
+	  MVertex *vi = t->getVertex(i);
+	  MVertex *vj = t->getVertex(j);
+	  double dx = vi->x()-vj->x();
+	  double dy = vi->y()-vj->y();
+	  double dz = vi->z()-vj->z();
+	  double l = sqrt(dx*dx+dy*dy+dz*dz);
+	  std::map<MVertex*,double>::iterator iti = vSizes.find(vi);	  
+	  std::map<MVertex*,double>::iterator itj = vSizes.find(vj);	  
+	  if (iti==vSizes.end() || iti->second > l)vSizes[vi] = l;
+	  if (itj==vSizes.end() || itj->second > l)vSizes[vj] = l;
+	}
+    }
+}
+
+
 void insertVerticesInRegion (GRegion *gr) 
 {
 
   std::set<MTet4*,compareTet4Ptr> allTets;
+  std::map<MVertex*,double> vSizesMap;
+  std::vector<double> vSizes;
 
-  for (int i=0;i<gr->tetrahedra.size();i++)
+  for (int i=0;i<gr->tetrahedra.size();i++)setLcs ( gr->tetrahedra[i] , vSizesMap);
+  
+  int NUM=0;
+  for (std::map<MVertex*,double>::iterator it = vSizesMap.begin();it!=vSizesMap.end();++it)
     {
-      allTets.insert ( new MTet4 ( gr->tetrahedra[i] ) );
+      it->first->setNum(NUM++);
+      vSizes.push_back(it->second);
     }
+
+  for (int i=0;i<gr->tetrahedra.size();i++)
+    allTets.insert ( new MTet4 ( gr->tetrahedra[i] ,vSizes ) );
+
   gr->tetrahedra.clear();
   connectTets ( allTets.begin(), allTets.end() );      
 
-//   for (std::set<MTet4*,compareTet4Ptr>::iterator it = allTets.begin();it!=allTets.end();++it)
-//     {
-//       printf("tet %d %d %d %d neigh %p %p %p %p\n",
-// 	     (*it)->tet()->getVertex(0)->getNum(),
-// 	     (*it)->tet()->getVertex(1)->getNum(),
-// 	     (*it)->tet()->getVertex(2)->getNum(),
-// 	     (*it)->tet()->getVertex(3)->getNum(),
-// 	     (*it)->getNeigh(0),(*it)->getNeigh(1),(*it)->getNeigh(2),(*it)->getNeigh(3));
-//     }
-  
-
   Msg(INFO,"All %d tets were connected",allTets.size());
 
   // here the classification should be done
 
-  III = 0;
+  int ITER = 0;
 
   while (1)
     {
-      //      connectTets ( allTets.begin(), allTets.end() );      
-      //       for (std::set<MTet4*,compareTet4Ptr>::iterator tit =  allTets.begin();tit != allTets.end();++tit)
-      // 	if(!(*tit)->assertNeigh())throw;
-      //      double vv = 0;
-      //      for (std::set<MTet4*,compareTet4Ptr>::iterator tit =  allTets.begin();tit != allTets.end();++tit)
-      //	if (!(*tit)->isDeleted())vv+=fabs((*tit)->getVolume());      
-
-
       MTet4 *worst = *allTets.begin();
 
       if (worst->isDeleted())
@@ -292,8 +304,9 @@ void insertVerticesInRegion (GRegion *gr)
 	}
       else
 	{
-	  if(III%1000 ==0)Msg(INFO,"%d points created -- Worst tet radius is %g",gr->mesh_vertices.size(),worst->getRadius());
-	  if (worst->getRadius() < .015) break;
+	  if(ITER++%5000 ==0)
+	    Msg(INFO,"%d points created -- Worst tet radius is %g",vSizes.size(),worst->getRadius());
+	  if (worst->getRadius() < 1) break;
 	  double center[3];
 	  worst->tet()->circumcenter(center);
 	  double uvw[3];
@@ -301,19 +314,25 @@ void insertVerticesInRegion (GRegion *gr)
 	  if (inside)
 	    {
 	      MVertex *v = new MVertex (center[0],center[1],center[2]);
+	      v->setNum(NUM++);
+	      double lc = 
+		(1-uvw[0]-uvw[1]-uvw[2])*vSizes[worst->tet()->getVertex(0)->getNum()] +
+		(uvw[0])*vSizes[worst->tet()->getVertex(1)->getNum()] +
+		(uvw[1])*vSizes[worst->tet()->getVertex(2)->getNum()] +
+		(uvw[2])*vSizes[worst->tet()->getVertex(3)->getNum()];
+	      vSizes.push_back(lc);
 	      // compute mesh spacing there
-	      III++;
-	      if (!insertVertex ( v  , worst, allTets))
+	      if (!insertVertex ( v  , worst, allTets,vSizes))
 		{
-		  // this one does not work
-		  //		  Msg(INFO,"Point is %g %g %g %d",uvw[0], uvw[1], uvw[2],worst->inCircumSphere ( v ));
 		  allTets.erase(allTets.begin());
 		  worst->forceRadius(0);
 		  allTets.insert(worst);		  
 		  delete v;
 		}
-	      else gr->mesh_vertices.push_back(v);
-	      //	      if(III == 15)  break;
+	      else 
+		{
+		  gr->mesh_vertices.push_back(v);
+		}
 	    }
 	  else
 	    {

Mesh/meshGRegionDelaunayInsertion.h

@@ -17,7 +17,7 @@ class MTet4
   void   forceRadius (double r){circum_radius=r;}
   double getRadius ()const {return circum_radius;}
   
-  MTet4 ( MTetrahedron * t) : deleted(false), base (t)
+  MTet4 ( MTetrahedron * t, std::vector<double> & sizes) : deleted(false), base (t)
   {
     neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
     double center[3];
@@ -26,6 +26,12 @@ class MTet4
     const double dy = base->getVertex(0)->y() - center[1];
     const double dz = base->getVertex(0)->z() - center[2];
     circum_radius = sqrt ( dx*dx + dy*dy + dz*dz);
+
+    double lc = 0.25*(sizes [base->getVertex(0)->getNum()]+
+		      sizes [base->getVertex(1)->getNum()]+
+		      sizes [base->getVertex(2)->getNum()]+
+		      sizes [base->getVertex(3)->getNum()]);
+    circum_radius /= lc;
   } 
 
   inline MTetrahedron * tet() const {return base;}

benchmarks/3d/Cube-04.geo

@@ -7,7 +7,7 @@ Extrude Point {1, {1,0.0,0} };
 Extrude Line {1, {0.0,0.0,1} };    
 aa[] = Extrude Surface {5, {0,1,0} };;
   
-Point(100) = {0.3,0.3,0.3,.02};              
+Point(100) = {0.3,0.3,0.3,.1};              
 Extrude Point {100, {.4,0.0,0} };                   
 Extrude Line {28, {0,0.4,0} };   
 Coherence;   

benchmarks/3d/Torus.geo

-lc = .2;          
+lc = .1;          
 Point(1) = {2.0,0.0,0.0,lc};          
 Point(2) = {2.0,1,0.0,lc};          
 Point(3) = {1,0,0.0,lc};