meshMetric and adaptMesh improved

Common/Octree.cpp

@@ -5,7 +5,7 @@
 
 #include <stdlib.h>
 #include <stdio.h>
-#include <list>
+#include <vector>
 #include "Octree.h"
 
 Octree* Octree_Create(int maxElements, double origin[3], double size[3],   
@@ -69,7 +69,8 @@ void Octree_Insert(void *element, Octree *myOctree)
 void Octree_Arrange(Octree *myOctree)
 {
   if(!myOctree) return;
-  std::list<void *>::iterator iter;
+  //std::list<void *>::iterator iter;
+  std::vector<void *>::iterator iter;
   double minPt[3], maxPt[3];
   for(iter = myOctree->info->listAllElements.begin(); iter!= 
       myOctree->info->listAllElements.end(); iter++) {
@@ -86,7 +87,7 @@ void *Octree_Search(double *pt, Octree *myOctree)
                        myOctree->function_BB, myOctree->function_inElement);
 }
 
-void Octree_SearchAll(double *pt, Octree *myOctree, std::list<void*> *output)
+void Octree_SearchAll(double *pt, Octree *myOctree, std::vector<void*> *output)
 {
   if(!myOctree) return;
   searchAllElements(myOctree->root, pt, myOctree->info, myOctree->function_BB,

Common/Octree.h

@@ -6,7 +6,7 @@
 #ifndef _OCTREE_H_
 #define _OCTREE_H_
 
-#include <list>
+#include <vector>
 #include "OctreeInternals.h"
 
 Octree* Octree_Create(int maxElements, // max. num of elts allowed in an octant
@@ -20,6 +20,6 @@ void Octree_Delete(Octree *);
 void Octree_Insert(void *, Octree *);
 void Octree_Arrange(Octree *);
 void *Octree_Search(double *, Octree *);
-void Octree_SearchAll(double *, Octree *, std::list<void *> *);
+void Octree_SearchAll(double *, Octree *, std::vector<void *> *);
 
 #endif

Common/OctreeInternals.cpp

@@ -283,7 +283,7 @@ void *searchElement(octantBucket *_buckets_head, double *_pt, globalInfo *_globa
   int flag;
   octantBucket *ptrBucket;
   ELink ptr1;  
-  std::list<void*>::iterator iter;
+  std::vector<void*>::iterator iter;
   void * ptrToEle = _globalPara->ptrToPrevElement;
 
   if (ptrToEle) {
@@ -369,7 +369,8 @@ void insertOneBB(void* _region, double *_minPt, double *_maxPt, octantBucket *_b
       ptr = ptr->next;
     }
 
-    _bucket->listBB.insert(_bucket->listBB.end(),_region);      
+    //_bucket->listBB.insert(_bucket->listBB.end(),_region); 
+    _bucket->listBB.push_back(_region);      
     return;
   }
         
@@ -380,11 +381,11 @@ void insertOneBB(void* _region, double *_minPt, double *_maxPt, octantBucket *_b
 
 void *searchAllElements(octantBucket *_buckets_head, double *_pt, globalInfo *_globalPara,
                         BBFunction BBElement, InEleFunction xyzInElement, 
-                        std::list<void*> *_elements)
+                        std::vector<void*> *_elements)
 {
   int flag, flag1;
   octantBucket *ptrBucket;
-  std::list<void*>::iterator iter;
+  std::vector<void*>::iterator iter;
 
   ptrBucket = findElementBucket(_buckets_head, _pt);
   if (ptrBucket == NULL) {

Common/OctreeInternals.h

@@ -6,7 +6,7 @@
 #ifndef _OCTREE_INTERNALS_H_
 #define _OCTREE_INTERNALS_H_
 
-#include <list>
+#include <vector>
 
 // file of function prototypes and macro constants 
 typedef void (*BBFunction)(void *, double*, double*);
@@ -30,7 +30,7 @@ struct bucket {
   int numElements; // number of elements contained by bucket 
   int precision;   // the level of precision of the bucket 
   ELink lhead; // list of elements in bucket, if NULL -> no elements 
-  std::list<void *> listBB; // list of elements in bucket by Bounding Box       
+  std::vector<void *> listBB; // list of elements in bucket by Bounding Box       
   struct bucket *next; // link to ragged digit extensions to bucket array 
   struct bucket *parent; // link to the parent bucket 
 };
@@ -44,7 +44,7 @@ struct global {
   double origin[3];   // smallest x,y, z of model's bounding box  
   double size[3];    // size in x, y, z of model bounding box 
   void * ptrToPrevElement;      
-  std::list<void *> listAllElements;
+  std::vector<void *> listAllElements;
 };
 typedef struct global globalInfo;
 
@@ -76,6 +76,6 @@ int xyzInElementBB(double *xyz, void *region, BBFunction BBElement);
 void insertOneBB(void *, double *, double *, octantBucket *);
 void *searchAllElements(octantBucket *_buckets_head, double *_pt, 
                         globalInfo *_globalPara, BBFunction BBElement, 
-                        InEleFunction xyzInElement, std::list<void *> *_elements);
+                        InEleFunction xyzInElement, std::vector<void *> *_elements);
 
 #endif

Geo/GModel.cpp

@@ -523,52 +523,50 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub
 {
 #if defined(HAVE_MESH)
 
-  int ITER = 0;
-  int nbElemsOld = 0;
+  if (getNumMeshElements() == 0) mesh(getDim());
+  int nbElemsOld = getNumMeshElements();
   int nbElems;
   int niter = parameters.size() >=4 ? (int) parameters[3] : 3;
 
+  FieldManager *fields = getFields();
+  fields->reset();
+
+  int ITER = 0;
   if (meshAll){
-  
-    FieldManager *fields = getFields();
-    fields->reset();
+
     while(1){
       Msg::Info("-- adaptMesh (allDim) ITER =%d ", ITER);
 
+      fields->reset();
+      int id = fields->newId();
+      (*fields)[id] = new meshMetric(this, technique, f, parameters);;
+      fields->background_field = id;    
+
+      opt_mesh_lc_integration_precision(0, GMSH_SET, 1.e-4);
       opt_mesh_algo2d(0, GMSH_SET, 7.0); //bamg
       opt_mesh_algo3d(0, GMSH_SET, 7.0); //mmg3d
       opt_mesh_lc_from_points(0, GMSH_SET, 0.0); //do not mesh lines with lc
-      GenerateMesh(this, getDim());
-      nbElems = getNumMeshElements();
-
-      writeMSH("meshAdapt.msh");
 
-      fields->reset();
-      if (++ITER >= niter)  break;
-      if (ITER > 5 && fabs((double)(nbElems - nbElemsOld)) < 0.01 * nbElemsOld) break;
-	
-      int id = fields->newId();
-      (*fields)[id] = new meshMetric(this, technique, f, parameters);;
-      fields->background_field = id;
-   
-            
       std::for_each(firstEdge(), lastEdge(), deMeshGEdge());
       std::for_each(firstFace(), lastFace(), deMeshGFace());
       std::for_each(firstRegion(), lastRegion(), deMeshGRegion());
 
-      nbElemsOld = nbElems;
-      
+      GenerateMesh(this, getDim());
+      nbElems = getNumMeshElements();
+
+      char name[256];
+      sprintf(name, "meshAdapt-%d.msh", ITER);
+      writeMSH(name);
+      //exit(1);
+            
+      if (ITER++ >= niter)  break;
+      if (ITER > 5 && fabs((double)(nbElems - nbElemsOld)) < 0.0 * nbElemsOld) break;
+	   
+      nbElemsOld = nbElems;  
     }
-    fields->reset();
   }
-
   else{
 
-    if (getNumMeshElements() == 0) mesh(getDim());
-    //meshMetric *mm; 
-    FieldManager *fields = getFields();
-    fields->reset();
-
     while(1) {
       Msg::Info("-- adaptMesh ITER =%d ", ITER);
       std::vector<MElement*> elements;
@@ -592,13 +590,11 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub
 
       nbElems = elements.size();
       if (nbElems == 0)return -1;
-
  
+      fields->reset();
       int id = fields->newId();
-      (*fields)[id] = new meshMetric(this, technique, f, parameters);;
+      (*fields)[id] = new meshMetric(this, technique, f, parameters);
       fields->background_field = id;
-      //mm = new meshMetric(this, technique, f, parameters);
-      //mm->setAsBackgroundMesh (this);
 
       if (getDim() == 2){
 	for (fiter fit = firstFace(); fit != lastFace(); ++fit){
@@ -617,16 +613,16 @@ int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<doub
 	  _octree = 0;
 	}
       }
-      fields->reset();
-      //delete mm;
+     
       if (++ITER >= niter) break;
       if (fabs((double)(nbElems - nbElemsOld)) < 0.01 * nbElemsOld) break;
 
       nbElemsOld = nbElems;
-
     }
   }
 
+  fields->reset();
+
   return 0;
 #else
   Msg::Error("Mesh module not compiled");

Geo/GRbf.cpp

@@ -211,12 +211,12 @@ void GRbf::buildOctree(double radius)
   Octree_Arrange(oct);
 
   for (int i= 0; i< nbNodes; i++){
-    std::list<void*> l;
+    std::vector<void*> l;
     double P[3] = {centers(i,0),centers(i,1), centers(i,2)};
     Octree_SearchAll(P, oct, &l);
     nodesInSphere[i].push_back(i);
     if (l.size() == 1) printf("*** WARNING: Found only %d sphere ! \n", (int)l.size());
-    for (std::list<void*>::iterator it = l.begin(); it != l.end(); it++) {
+    for (std::vector<void*>::iterator it = l.begin(); it != l.end(); it++) {
       Sphere *sph = (Sphere*) *it;
       std::vector<int> &pts = nodesInSphere[i];
       if (sph->index != i)  nodesInSphere[i].push_back(sph->index);

Geo/MElementOctree.cpp

@@ -126,10 +126,10 @@ MElementOctree::~MElementOctree()
 std::vector<MElement *> MElementOctree::findAll(double x, double y, double z, int dim)
 {
   double P[3] = {x, y, z};
-  std::list<void*> v;
+  std::vector<void*> v;
   std::vector<MElement*> e;
   Octree_SearchAll(P, _octree,&v);
-  for (std::list<void*>::iterator it = v.begin();
+  for (std::vector<void*>::iterator it = v.begin();
        it != v.end(); ++it){
     MElement *el = (MElement*) *it;
     if (dim == -1 || el->getDim() == dim)e.push_back(el);
@@ -143,10 +143,10 @@ MElement *MElementOctree::find(double x, double y, double z, int dim, bool stric
   MElement *e = (MElement*)Octree_Search(P, _octree);
   if (e && (dim == -1 || e->getDim() == dim))
     return e;
-  std::list<void*> l;
+  std::vector<void*> l;
   if (e && e->getDim() != dim) {
     Octree_SearchAll(P, _octree, &l);
-    for (std::list<void*>::iterator it = l.begin(); it != l.end(); it++) {
+    for (std::vector<void*>::iterator it = l.begin(); it != l.end(); it++) {
       MElement *el = (MElement*) *it;
       if (el->getDim() == dim) {
         return el;

Mesh/meshGFaceLloyd.cpp

@@ -971,7 +971,7 @@ double lpcvt::total_area(){
   double total;
   SPoint2 p1,p2,p3;
   voronoi_vertex v1,v2,v3;
-  std::list<voronoi_element>::iterator it;
+  std::vector<voronoi_element>::iterator it;
   total = 0.0;
   for(it=clipped.begin();it!=clipped.end();it++){
     v1 = it->get_v1();
@@ -988,7 +988,7 @@ double lpcvt::total_area(){
 void lpcvt::print_voronoi1(){
   SPoint2 p1,p2,p3;
   voronoi_vertex v1,v2,v3;
-  std::list<voronoi_element>::iterator it;
+  std::vector<voronoi_element>::iterator it;
   std::ofstream file("voronoi1.pos");
   file << "View \"test\" {\n";
   for(it=clipped.begin();it!=clipped.end();it++){
@@ -1066,7 +1066,7 @@ void lpcvt::compute_parameters(GFace* gf,int p){
   SPoint2 p1,p2,p3;
   voronoi_vertex v1,v2,v3;
   metric m;
-  std::list<voronoi_element>::iterator it;
+  std::vector<voronoi_element>::iterator it;
 	
   k = 1.0;	
   for(it=clipped.begin();it!=clipped.end();it++){
@@ -1143,7 +1143,7 @@ void lpcvt::eval(DocRecord& triangulator,std::vector<SVector3>& gradients,double
   SVector3 grad1,grad2;
   SVector3 normal;
   voronoi_vertex v1,v2,v3;
-  std::list<voronoi_element>::iterator it;
+  std::vector<voronoi_element>::iterator it;
 	
   for(i=0;i<gradients.size();i++){
     gradients[i] = SVector3(0.0,0.0,0.0);
@@ -1207,7 +1207,7 @@ void lpcvt::eval(DocRecord& triangulator,std::vector<SVector3>& gradients,double
 
 void lpcvt::swap(){
   voronoi_vertex vertex;
-  std::list<voronoi_element>::iterator it;
+  std::vector<voronoi_element>::iterator it;
   for(it=clipped.begin();it!=clipped.end();it++){
 	if(J(it->get_v1().get_point(),it->get_v2().get_point(),it->get_v3().get_point())<0.0){
       vertex = it->get_v3();

Mesh/meshGFaceLloyd.h

@@ -34,7 +34,7 @@ class smoothing{
 
 class lpcvt{
  private :
-  std::list<voronoi_element> clipped;
+  std::vector<voronoi_element> clipped;
   std::queue<int> fifo;
   std::vector<segment_list> borders;
   std::vector<double> angles;

Mesh/meshMetric.cpp

@@ -85,7 +85,9 @@ meshMetric::meshMetric(GModel *gm, int technique, simpleFunction<double> *fct, s
 
 }
 meshMetric::~meshMetric(){
-  //if (_octree) delete _octree;
+  if (_octree) delete _octree;
+  for (int i=0; i< _elements.size(); i++)
+   delete _elements[i];
 }
 
 void meshMetric::computeValues( v2t_cont adj){
@@ -108,7 +110,7 @@ void meshMetric::computeHessian( v2t_cont adj){
     while (it != adj.end()) {
       std::vector<MElement*> lt = it->second;      
       MVertex *ver = it->first;
-      while (lt.size() < 12) increaseStencil(ver,adj,lt); //<7
+      while (lt.size() < 9) increaseStencil(ver,adj,lt); //<7
       // if ( ver->onWhat()->dim() < _dim ){
       // 	while (lt.size() < 12){
       // 	  increaseStencil(ver,adj,lt);
@@ -331,8 +333,7 @@ void meshMetric::computeMetric(){
   //smoothMetric (sol);
   //curvatureContributionToMetric();
 
-  //putOnNewView();
-
+ 
 }
 
 double meshMetric::operator() (double x, double y, double z, GEntity *ge) {
@@ -379,31 +380,25 @@ void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEnti
   }
 }
 
-void meshMetric::setAsBackgroundMesh (GModel *gm){
-  FieldManager *fields = gm->getFields();
-  int id = fields->newId();
-  (*fields)[id] = new meshMetric(*this);
-  fields->background_field = id;
-}
-
 void meshMetric::printMetric(const char* n) const{
 
   FILE *f = fopen (n,"w");
   fprintf(f,"View \"\"{\n");
 
-  //std::map<MVertex*,SMetric3 >::const_iterator it = _hessian.begin(); 
-  std::map<MVertex*,SMetric3>::const_iterator it= _nodalMetrics.begin();
-  for (; it != _nodalMetrics.end(); ++it){
-    //for (; it != _hessian.end(); ++it){
+  std::map<MVertex*,SMetric3 >::const_iterator it = _hessian.begin(); 
+  //std::map<MVertex*,SMetric3>::const_iterator it= _nodalMetrics.begin();
+  //for (; it != _nodalMetrics.end(); ++it){
+    for (; it != _hessian.end(); ++it){
       MVertex *v =  it->first;
-      // double lapl =  getLaplacian(v);
-      // fprintf(f, "SP(%g,%g,%g){%g};\n",  it->first->x(),it->first->y(),it->first->z(),lapl);
-      fprintf(f,"TP(%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
-      	    it->first->x(),it->first->y(),it->first->z(),
-      	    it->second(0,0),it->second(0,1),it->second(0,2),
-      	    it->second(1,0),it->second(1,1),it->second(1,2),
-      	    it->second(2,0),it->second(2,1),it->second(2,2)
-      	    );
+      SMetric3 h = it->second;
+      double lapl = h(0,0)+h(1,1)+h(2,2); 
+       fprintf(f, "SP(%g,%g,%g){%g};\n",  it->first->x(),it->first->y(),it->first->z(),lapl);
+      // fprintf(f,"TP(%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
+      // 	    it->first->x(),it->first->y(),it->first->z(),
+      // 	    it->second(0,0),it->second(0,1),it->second(0,2),
+      // 	    it->second(1,0),it->second(1,1),it->second(1,2),
+      // 	    it->second(2,0),it->second(2,1),it->second(2,2)
+      // 	    );
 
   } 
   fprintf(f,"};\n");

Mesh/meshMetric.h

@@ -41,8 +41,7 @@ class meshMetric: public Field {
   void computeMetric() ;
   void computeValues( v2t_cont adj);
   void computeHessian( v2t_cont adj);
-  void setAsBackgroundMesh (GModel *gm);
-
+ 
   double getLaplacian (MVertex *v);
   virtual bool isotropic () const {return false;}
   virtual const char *getName()

benchmarks/boolean/sphere.py

-options = gmshOptions()
-options:initOptions()
-options:numberSet('Mesh', 0, 'CharacteristicLengthFactor', 0.9);
+from dgpy import *
+from gmshpy import *
 
-R = 1.0;
+GmshSetOption('Mesh', 'CharacteristicLengthFactor', 0.9)
+
+R = 0.3;
 myModel = GModel();
-myModel:addSphere(0,0,0,R);
+myModel.addSphere(0.5,0.5,0.5,R);
 
-myModel:setAsCurrent();
+myModel.setAsCurrent();
 
---myModel:mesh(2);
+myModel.mesh(2);
+myModel.save("sphere.stl");