From f4d7ff4e12d13739feabcd8151c0bdc54ce5e0b7 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Thu, 22 Dec 2011 14:33:36 +0000
Subject: [PATCH] meshMetric and adaptMesh improved
---
Common/Octree.cpp | 7 ++--
Common/Octree.h | 4 +--
Common/OctreeInternals.cpp | 9 ++---
Common/OctreeInternals.h | 8 ++---
Geo/GModel.cpp | 64 +++++++++++++++++-------------------
Geo/GRbf.cpp | 4 +--
Geo/MElementOctree.cpp | 8 ++---
Mesh/meshGFaceLloyd.cpp | 10 +++---
Mesh/meshGFaceLloyd.h | 2 +-
Mesh/meshMetric.cpp | 41 ++++++++++-------------
Mesh/meshMetric.h | 3 +-
benchmarks/boolean/sphere.py | 16 +++++----
12 files changed, 85 insertions(+), 91 deletions(-)
diff --git a/Common/Octree.cpp b/Common/Octree.cpp
index 58c5c35133..a0786584ac 100644
--- a/Common/Octree.cpp
+++ b/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,
diff --git a/Common/Octree.h b/Common/Octree.h
index a831bcb77c..0a11342bbd 100644
--- a/Common/Octree.h
+++ b/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
diff --git a/Common/OctreeInternals.cpp b/Common/OctreeInternals.cpp
index bfd2d2fca8..8ac3a7c08a 100644
--- a/Common/OctreeInternals.cpp
+++ b/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) {
diff --git a/Common/OctreeInternals.h b/Common/OctreeInternals.h
index c71eda2576..9958e5d39d 100644
--- a/Common/OctreeInternals.h
+++ b/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
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 884147b1dd..db5456a403 100644
--- a/Geo/GModel.cpp
+++ b/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");
diff --git a/Geo/GRbf.cpp b/Geo/GRbf.cpp
index e7aa84ef0e..ac0e582a33 100644
--- a/Geo/GRbf.cpp
+++ b/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);
diff --git a/Geo/MElementOctree.cpp b/Geo/MElementOctree.cpp
index 5c1c5b802e..924bbc67df 100644
--- a/Geo/MElementOctree.cpp
+++ b/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;
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index 896edf9c0b..318773edc1 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/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();
diff --git a/Mesh/meshGFaceLloyd.h b/Mesh/meshGFaceLloyd.h
index 8e80fc7103..f64160fa76 100644
--- a/Mesh/meshGFaceLloyd.h
+++ b/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;
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index f56c08b54a..74d3e8d31a 100644
--- a/Mesh/meshMetric.cpp
+++ b/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");
diff --git a/Mesh/meshMetric.h b/Mesh/meshMetric.h
index 48561d834c..9460b6cb3f 100644
--- a/Mesh/meshMetric.h
+++ b/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()
diff --git a/benchmarks/boolean/sphere.py b/benchmarks/boolean/sphere.py
index 6477c57e48..d268c269c0 100644
--- a/benchmarks/boolean/sphere.py
+++ b/benchmarks/boolean/sphere.py
@@ -1,11 +1,13 @@
-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");
--
GitLab