From 516edae1a7fc98a7769827af3aa207c67a11aedf Mon Sep 17 00:00:00 2001
From: Tristan Carrier Baudouin <tristan.carrier@uclouvain.be>
Date: Tue, 13 Sep 2011 14:30:10 +0000
Subject: [PATCH] changes
---
CMakeLists.txt | 2 +-
Fltk/statisticsWindow.cpp | 190 ++++++++++++++++++--------------------
Mesh/BackgroundMesh.cpp | 2 +-
Mesh/CMakeLists.txt | 1 +
Mesh/meshGFaceLloyd.cpp | 56 ++++++-----
Mesh/meshGFaceLloyd.h | 8 +-
gmshpy/gmshMesh.i | 2 +
7 files changed, 127 insertions(+), 134 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b854c62e62..090ee201b8 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -85,7 +85,7 @@ set(GMSH_API
Geo/SOrientedBoundingBox.h Geo/CellComplex.h Geo/ChainComplex.h Geo/Cell.h
Geo/Homology.h
Mesh/meshGEdge.h Mesh/meshGFace.h Mesh/meshGFaceOptimize.h
- Mesh/meshGFaceDelaunayInsertion.h
+ Mesh/meshGFaceDelaunayInsertion.h Mesh/levy3D.h
Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
Solver/crossConfTerm.h Solver/orthogonalTerm.h
Solver/linearSystem.h Solver/linearSystemGMM.h Solver/linearSystemCSR.h
diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index f17285986b..f24de26489 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -186,106 +186,100 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
void statisticsWindow::compute(bool elementQuality)
{
//emi hack - MINIMUM ANGLES
- // double minAngle = 120.0;
- // double meanAngle = 0.0;
- // int count = 0;
- // std::vector<GEntity*> entities;
- // GModel::current()->getEntities(entities);
- // std::map<int, std::vector<double> > d;
- // for(unsigned int i = 0; i < entities.size(); i++){
- // if(entities[i]->dim() < 2) continue;
- // for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
- // MElement *e = entities[i]->getMeshElement(j);
- // double angle = e->angleShapeMeasure();
- // minAngle = std::min(minAngle, angle);
- // meanAngle += angle;
- // count++;
- // }
- // }
- // meanAngle = meanAngle / count;
- // printf("Angles = min=%g av=%g \n", minAngle, meanAngle);
+ /*double minAngle = 120.0;
+ double meanAngle = 0.0;
+ int count = 0;
+ std::vector<GEntity*> entities;
+ GModel::current()->getEntities(entities);
+ std::map<int, std::vector<double> > d;
+ for(unsigned int i = 0; i < entities.size(); i++){
+ if(entities[i]->dim() < 2) continue;
+ for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+ MElement *e = entities[i]->getMeshElement(j);
+ double angle = e->angleShapeMeasure();
+ minAngle = std::min(minAngle, angle);
+ meanAngle += angle;
+ count++;
+ }
+ }
+ meanAngle = meanAngle / count;
+ printf("Angles = min=%g av=%g \n", minAngle, meanAngle);*/
//hack emi
//Emi hack - MESH DEGREE VERTICES
- // std::vector<GEntity*> entities;
- // std::set<MEdge, Less_Edge> edges;
- // GModel::current()->getEntities(entities);
- // std::map<MVertex*, int > vert2Deg;
- // for(unsigned int i = 0; i < entities.size(); i++){
- // printf("entity dim =%d tag=%d \n", entities[i]->dim() , entities[i]->tag());
- // if(entities[i]->dim() < 2 ) continue;
- // if(entities[i]->tag() != 100) continue;
- // printf("continuing \n");
- // for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
- // MElement *e = entities[i]->getMeshElement(j);
- // for(unsigned int k = 0; k < e->getNumEdges(); k++){
- // edges.insert(e->getEdge(k));
- // }
- // for(unsigned int k = 0; k < e->getNumVertices(); k++){
- // MVertex *v = e->getVertex(k);
- // if (v->onWhat()->dim() < 2) continue;
- // std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
- // if (it == vert2Deg.end()) {
- // vert2Deg.insert(std::make_pair(v,1));
- // }
- // else{
- // int nbE = it->second+1;
- // it->second = nbE;
- // }
- // }
- // }
- // }
- // int dMin = 10;
- // int dMax = 0;
- // int d4 = 0;
- // int nbElems = vert2Deg.size();
- // std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
- // for(; itmap !=vert2Deg.end(); itmap++){
- // MVertex *v = itmap->first;
- // int nbE = itmap->second;
- // dMin = std::min(nbE, dMin);
- // dMax = std::max(nbE, dMax);
- // if (nbE == 4) d4 += 1;
- // }
- // if (nbElems > 0)
- // printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n", dMin, dMax, (double)d4/nbElems);
-
- // FieldManager *fields = GModel::current()->getFields();
- // Field *f = fields->get(fields->background_field);
- // int nbEdges = edges.size();
- // printf("nb edges =%d \n", nbEdges);
- // system("rm qualEdges.txt");
- // FILE *fp = fopen("qualEdges.txt", "w");
- // std::vector<int> qualE;
- // int nbS = 50;
- // qualE.resize(nbS);
- // if(fields->background_field > 0){
- // printf("found field \n");
- // std::set<MEdge, Less_Edge>::iterator it = edges.begin();
- // double sum = 0;
- // for (; it !=edges.end();++it){
- // MVertex *v0 = it->getVertex(0);
- // MVertex *v1 = it->getVertex(1);
- // double l = sqrt((v0->x()-v1->x())*(v0->x()-v1->x())+
- // (v0->y()-v1->y())*(v0->y()-v1->y())+
- // (v0->z()-v1->z())*(v0->z()-v1->z()));
- // double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),0.5*(v0->z()+v1->z()),v0->onWhat());
- // double el = l/lf;
- // int index = (int) ceil(el*nbS*0.5);
- // qualE[index]+= 1;
- // double e = (l>lf) ? lf/l : l/lf;
- // sum += e - 1.0;
- // }
- // double tau = exp ((1./edges.size()) * sum);
- // printf("N edges = %d tau = %g\n",(int)edges.size(),tau);
-
- // double ibegin = 2./(2*nbS);
- // double inext = 2./nbS;
- // for (int i= 0; i< qualE.size(); i++){
- // fprintf(fp, "0 0 0 0 %g 0 0 %g \n", ibegin+i*inext , (double)qualE[i]/nbEdges);
- // }
-
- // }
- // fclose(fp);
+ std::vector<GEntity*> entities;
+ std::set<MEdge, Less_Edge> edges;
+ GModel::current()->getEntities(entities);
+ std::map<MVertex*, int > vert2Deg;
+ for(unsigned int i = 0; i < entities.size(); i++){
+ printf("entity dim =%d tag=%d \n", entities[i]->dim() , entities[i]->tag());
+ if(entities[i]->dim() < 2 ) continue;
+ if(entities[i]->tag() != 10) continue;
+ printf("continuing \n");
+ for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+ MElement *e = entities[i]->getMeshElement(j);
+ for(unsigned int k = 0; k < e->getNumEdges(); k++){
+ edges.insert(e->getEdge(k));
+ }
+ for(unsigned int k = 0; k < e->getNumVertices(); k++){
+ MVertex *v = e->getVertex(k);
+ if (v->onWhat()->dim() < 2) continue;
+ std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
+ if (it == vert2Deg.end()){
+ vert2Deg.insert(std::make_pair(v,1));
+ }
+ else{
+ int nbE = it->second+1;
+ it->second = nbE;
+ }
+ }
+ }
+ }
+ int dMin = 10;
+ int dMax = 0;
+ int d4 = 0;
+ int nbElems = vert2Deg.size();
+ std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
+ for(; itmap !=vert2Deg.end(); itmap++){
+ MVertex *v = itmap->first;
+ int nbE = itmap->second;
+ dMin = std::min(nbE, dMin);
+ dMax = std::max(nbE, dMax);
+ if (nbE == 4) d4 += 1;
+ }
+ if (nbElems > 0) printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n", dMin, dMax, (double)d4/nbElems);
+ FieldManager *fields = GModel::current()->getFields();
+ Field *f = fields->get(fields->background_field);
+ int nbEdges = edges.size();
+ printf("nb edges =%d \n", nbEdges);
+ system("rm qualEdges.txt");
+ FILE *fp = fopen("qualEdges.txt", "w");
+ std::vector<int> qualE;
+ int nbS = 50;
+ qualE.resize(nbS);
+ if(fields->background_field > 0){
+ printf("found field \n");
+ std::set<MEdge, Less_Edge>::iterator it = edges.begin();
+ double sum = 0;
+ for (; it !=edges.end();++it){
+ MVertex *v0 = it->getVertex(0);
+ MVertex *v1 = it->getVertex(1);
+ double l = sqrt((v0->x()-v1->x())*(v0->x()-v1->x())+(v0->y()-v1->y())*(v0->y()-v1->y())+(v0->z()-v1->z())*(v0->z()-v1->z()));
+ double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),0.5*(v0->z()+v1->z()),v0->onWhat());
+ double el = l/lf;
+ int index = (int) ceil(el*nbS*0.5);
+ qualE[index]+= 1;
+ double e = (l>lf) ? lf/l : l/lf;
+ sum += e - 1.0;
+ }
+ double tau = exp ((1./edges.size()) * sum);
+ printf("N edges = %d tau = %g\n",(int)edges.size(),tau);
+ double ibegin = 2./(2*nbS);
+ double inext = 2./nbS;
+ for (int i= 0; i< qualE.size(); i++){
+ fprintf(fp, "0 0 0 0 %g 0 0 %g \n", ibegin+i*inext , (double)qualE[i]/nbEdges);
+ }
+ }
+ fclose(fp);
//emi end hack
int num = 0;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index ad053eb098..1f4dadee43 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -727,7 +727,7 @@ double backgroundMesh::operator() (double u, double v, double w) const
MElement *e = _octree->find(u, v, w, 2, true);
if (!e) {
Msg::Error("cannot find %g %g", u, v);
- return -1000.0;
+ return -1000.0;//0.4;
}
e->xyz2uvw(uv, uv2);
std::map<MVertex*,double>::const_iterator itv1 = _sizes.find(e->getVertex(0));
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index 6429511992..329b610587 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -30,6 +30,7 @@ set(SRC
multiscalePartition.cpp
QuadTriUtils.cpp
QuadTriExtruded2D.cpp QuadTriExtruded3D.cpp QuadTriTransfinite3D.cpp
+ Levy3D.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index 5ac9777aa9..259efea9e3 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -1014,6 +1014,7 @@ void lpcvt::print_segment(SPoint2 p1,SPoint2 p2,std::ofstream& file){
void lpcvt::compute_parameters(GFace* gf,int p){
double h1,h2,h3;
+ double rho1,rho2,rho3;
double k;
double ratio;
double angle;
@@ -1038,13 +1039,16 @@ void lpcvt::compute_parameters(GFace* gf,int p){
h1 = k*backgroundMesh::current()->operator()(p1.x(),p1.y(),0.0)*ratio;
h2 = k*backgroundMesh::current()->operator()(p2.x(),p2.y(),0.0)*ratio;
h3 = k*backgroundMesh::current()->operator()(p3.x(),p3.y(),0.0)*ratio;
+ rho1 = it->compute_rho(h1,p);
+ rho2 = it->compute_rho(h2,p);
+ rho3 = it->compute_rho(h3,p);
angle = -backgroundMesh::current()->getAngle(p1.x(),p1.y(),0.0);
cosinus = cos(angle);
sinus = sin(angle);
m = metric(cosinus,-sinus,sinus,cosinus);
- v1.set_h(h1);
- v2.set_h(h2);
- v3.set_h(h3);
+ v1.set_rho(rho1);
+ v2.set_rho(rho2);
+ v3.set_rho(rho3);
it->set_v1(v1);
it->set_v2(v2);
it->set_v3(v3);
@@ -1207,7 +1211,7 @@ double lpcvt::F(voronoi_element element,int p){
y = Ty(u,v,generator,C1,C2);
point = SPoint2(x,y);
weight = gauss_weights(i,0);
- rho = element.get_rho(u,v,p);
+ rho = element.get_rho(u,v);
energy = energy + weight*rho*f(generator,point,m,p);
}
energy = J(generator,C1,C2)*energy;
@@ -1247,7 +1251,7 @@ SVector3 lpcvt::simple(voronoi_element element,int p){
y = Ty(u,v,generator,C1,C2);
point = SPoint2(x,y);
weight = gauss_weights(i,0);
- rho = element.get_rho(u,v,p);
+ rho = element.get_rho(u,v);
comp_x = comp_x + weight*rho*df_dx(generator,point,m,p);
comp_y = comp_y + weight*rho*df_dy(generator,point,m,p);
}
@@ -1292,7 +1296,7 @@ SVector3 lpcvt::dF_dC1(voronoi_element element,int p){
y = Ty(u,v,generator,C1,C2);
point = SPoint2(x,y);
weight = gauss_weights(i,0);
- rho = element.get_rho(u,v,p);
+ rho = element.get_rho(u,v);
drho_dx = element.get_drho_dx();
drho_dy = element.get_drho_dy();
distance = f(point,generator,m,p);
@@ -1342,7 +1346,7 @@ SVector3 lpcvt::dF_dC2(voronoi_element element,int p){
y = Ty(u,v,generator,C1,C2);
point = SPoint2(x,y);
weight = gauss_weights(i,0);
- rho = element.get_rho(u,v,p);
+ rho = element.get_rho(u,v);
drho_dx = element.get_drho_dx();
drho_dy = element.get_drho_dy();
distance = f(point,generator,m,p);
@@ -1586,8 +1590,8 @@ bool voronoi_vertex::get_duplicate(){
return duplicate;
}
-double voronoi_vertex::get_h(){
- return h;
+double voronoi_vertex::get_rho(){
+ return rho;
}
void voronoi_vertex::set_point(SPoint2 new_point){
@@ -1614,8 +1618,8 @@ void voronoi_vertex::set_duplicate(bool new_duplicate){
duplicate = new_duplicate;
}
-void voronoi_vertex::set_h(double new_h){
- h = new_h;
+void voronoi_vertex::set_rho(double new_rho){
+ rho = new_rho;
}
@@ -1644,18 +1648,16 @@ voronoi_vertex voronoi_element::get_v3(){
return v3;
}
-double voronoi_element::get_rho(double u,double v,int p){
- double h1;
- double h2;
- double h3;
- double h;
+double voronoi_element::get_rho(double u,double v){
+ double rho1;
+ double rho2;
+ double rho3;
double rho;
- h1 = v1.get_h();
- h2 = v2.get_h();
- h3 = v3.get_h();
- h = h1*(1.0-u-v) + h2*u + h3*v;
- rho = compute_rho(h,p);
+ rho1 = v1.get_rho();
+ rho2 = v2.get_rho();
+ rho3 = v3.get_rho();
+ rho = rho1*(1.0-u-v) + rho2*u + rho3*v;
return rho;
}
@@ -1688,9 +1690,6 @@ void voronoi_element::set_metric(metric new_m){
}
void voronoi_element::deriv_rho(int p){
- double h1;
- double h2;
- double h3;
double rho1;
double rho2;
double rho3;
@@ -1709,12 +1708,9 @@ void voronoi_element::deriv_rho(int p){
SPoint2 p2;
SPoint2 p3;
- h1 = v1.get_h();
- h2 = v2.get_h();
- h3 = v3.get_h();
- rho1 = compute_rho(h1,p);
- rho2 = compute_rho(h2,p);
- rho3 = compute_rho(h3,p);
+ rho1 = v1.get_rho();
+ rho2 = v2.get_rho();
+ rho3 = v3.get_rho();
p1 = v1.get_point();
p2 = v2.get_point();
p3 = v3.get_point();
diff --git a/Mesh/meshGFaceLloyd.h b/Mesh/meshGFaceLloyd.h
index 544f4c475e..59f62124b6 100644
--- a/Mesh/meshGFaceLloyd.h
+++ b/Mesh/meshGFaceLloyd.h
@@ -126,7 +126,7 @@ class voronoi_vertex{
int index3;
SVector3 normal;
bool duplicate;
- double h;
+ double rho;
public :
voronoi_vertex(SPoint2);
voronoi_vertex();
@@ -137,14 +137,14 @@ class voronoi_vertex{
int get_index3();
SVector3 get_normal();
bool get_duplicate();
- double get_h();
+ double get_rho();
void set_point(SPoint2);
void set_index1(int);
void set_index2(int);
void set_index3(int);
void set_normal(SVector3);
void set_duplicate(bool);
- void set_h(double);
+ void set_rho(double);
};
class voronoi_element{
@@ -162,7 +162,7 @@ class voronoi_element{
voronoi_vertex get_v1();
voronoi_vertex get_v2();
voronoi_vertex get_v3();
- double get_rho(double,double,int);
+ double get_rho(double,double);
double get_drho_dx();
double get_drho_dy();
metric get_metric();
diff --git a/gmshpy/gmshMesh.i b/gmshpy/gmshMesh.i
index 96c334324a..a18ea5cef5 100644
--- a/gmshpy/gmshMesh.i
+++ b/gmshpy/gmshMesh.i
@@ -10,6 +10,7 @@
#include "meshGFaceLloyd.h"
#include "meshGFaceOptimize.h"
#include "meshPartitionOptions.h"
+ #include "Levy3D.h"
%}
%include "GmshConfig.h"
@@ -18,3 +19,4 @@
%include "meshGFaceLloyd.h"
%include "meshGFaceOptimize.h"
%include "meshPartitionOptions.h"
+%include "Levy3D.h"
--
GitLab