Skip to content
Snippets Groups Projects
Commit 516edae1 authored by Tristan Carrier Baudouin's avatar Tristan Carrier Baudouin
Browse files

changes

parent 1aaad4fa
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 1
− 1
View file @ 516edae1
...@@ -85,7 +85,7 @@ set(GMSH_API ...@@ -85,7 +85,7 @@ set(GMSH_API
Geo/SOrientedBoundingBox.h Geo/CellComplex.h Geo/ChainComplex.h Geo/Cell.h Geo/SOrientedBoundingBox.h Geo/CellComplex.h Geo/ChainComplex.h Geo/Cell.h
Geo/Homology.h Geo/Homology.h
Mesh/meshGEdge.h Mesh/meshGFace.h Mesh/meshGFaceOptimize.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/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
Solver/crossConfTerm.h Solver/orthogonalTerm.h Solver/crossConfTerm.h Solver/orthogonalTerm.h
Solver/linearSystem.h Solver/linearSystemGMM.h Solver/linearSystemCSR.h Solver/linearSystem.h Solver/linearSystemGMM.h Solver/linearSystemCSR.h
......
Fltk/statisticsWindow.cpp
+ 92
− 98
View file @ 516edae1
...@@ -186,106 +186,100 @@ statisticsWindow::statisticsWindow(int deltaFontSize) ...@@ -186,106 +186,100 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
void statisticsWindow::compute(bool elementQuality) void statisticsWindow::compute(bool elementQuality)
{ {
//emi hack - MINIMUM ANGLES //emi hack - MINIMUM ANGLES
// double minAngle = 120.0; /*double minAngle = 120.0;
// double meanAngle = 0.0; double meanAngle = 0.0;
// int count = 0; int count = 0;
// std::vector<GEntity*> entities; std::vector<GEntity*> entities;
// GModel::current()->getEntities(entities); GModel::current()->getEntities(entities);
// std::map<int, std::vector<double> > d; std::map<int, std::vector<double> > d;
// for(unsigned int i = 0; i < entities.size(); i++){ for(unsigned int i = 0; i < entities.size(); i++){
// if(entities[i]->dim() < 2) continue; if(entities[i]->dim() < 2) continue;
// for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){ for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
// MElement *e = entities[i]->getMeshElement(j); MElement *e = entities[i]->getMeshElement(j);
// double angle = e->angleShapeMeasure(); double angle = e->angleShapeMeasure();
// minAngle = std::min(minAngle, angle); minAngle = std::min(minAngle, angle);
// meanAngle += angle; meanAngle += angle;
// count++; count++;
// } }
// } }
// meanAngle = meanAngle / count; meanAngle = meanAngle / count;
// printf("Angles = min=%g av=%g \n", minAngle, meanAngle); printf("Angles = min=%g av=%g \n", minAngle, meanAngle);*/
//hack emi //hack emi
//Emi hack - MESH DEGREE VERTICES //Emi hack - MESH DEGREE VERTICES
// std::vector<GEntity*> entities; std::vector<GEntity*> entities;
// std::set<MEdge, Less_Edge> edges; std::set<MEdge, Less_Edge> edges;
// GModel::current()->getEntities(entities); GModel::current()->getEntities(entities);
// std::map<MVertex*, int > vert2Deg; std::map<MVertex*, int > vert2Deg;
// for(unsigned int i = 0; i < entities.size(); i++){ for(unsigned int i = 0; i < entities.size(); i++){
// printf("entity dim =%d tag=%d \n", entities[i]->dim() , entities[i]->tag()); printf("entity dim =%d tag=%d \n", entities[i]->dim() , entities[i]->tag());
// if(entities[i]->dim() < 2 ) continue; if(entities[i]->dim() < 2 ) continue;
// if(entities[i]->tag() != 100) continue; if(entities[i]->tag() != 10) continue;
// printf("continuing \n"); printf("continuing \n");
// for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){ for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
// MElement *e = entities[i]->getMeshElement(j); MElement *e = entities[i]->getMeshElement(j);
// for(unsigned int k = 0; k < e->getNumEdges(); k++){ for(unsigned int k = 0; k < e->getNumEdges(); k++){
// edges.insert(e->getEdge(k)); edges.insert(e->getEdge(k));
// } }
// for(unsigned int k = 0; k < e->getNumVertices(); k++){ for(unsigned int k = 0; k < e->getNumVertices(); k++){
// MVertex *v = e->getVertex(k); MVertex *v = e->getVertex(k);
// if (v->onWhat()->dim() < 2) continue; if (v->onWhat()->dim() < 2) continue;
// std::map<MVertex*, int >::iterator it = vert2Deg.find(v); std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
// if (it == vert2Deg.end()) { if (it == vert2Deg.end()){
// vert2Deg.insert(std::make_pair(v,1)); vert2Deg.insert(std::make_pair(v,1));
// } }
// else{ else{
// int nbE = it->second+1; int nbE = it->second+1;
// it->second = nbE; it->second = nbE;
// } }
// } }
// } }
// } }
// int dMin = 10; int dMin = 10;
// int dMax = 0; int dMax = 0;
// int d4 = 0; int d4 = 0;
// int nbElems = vert2Deg.size(); int nbElems = vert2Deg.size();
// std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin(); std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
// for(; itmap !=vert2Deg.end(); itmap++){ for(; itmap !=vert2Deg.end(); itmap++){
// MVertex *v = itmap->first; MVertex *v = itmap->first;
// int nbE = itmap->second; int nbE = itmap->second;
// dMin = std::min(nbE, dMin); dMin = std::min(nbE, dMin);
// dMax = std::max(nbE, dMax); dMax = std::max(nbE, dMax);
// if (nbE == 4) d4 += 1; if (nbE == 4) d4 += 1;
// } }
// if (nbElems > 0) if (nbElems > 0) printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n", dMin, dMax, (double)d4/nbElems);
// 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);
// FieldManager *fields = GModel::current()->getFields(); int nbEdges = edges.size();
// Field *f = fields->get(fields->background_field); printf("nb edges =%d \n", nbEdges);
// int nbEdges = edges.size(); system("rm qualEdges.txt");
// printf("nb edges =%d \n", nbEdges); FILE *fp = fopen("qualEdges.txt", "w");
// system("rm qualEdges.txt"); std::vector<int> qualE;
// FILE *fp = fopen("qualEdges.txt", "w"); int nbS = 50;
// std::vector<int> qualE; qualE.resize(nbS);
// int nbS = 50; if(fields->background_field > 0){
// qualE.resize(nbS); printf("found field \n");
// if(fields->background_field > 0){ std::set<MEdge, Less_Edge>::iterator it = edges.begin();
// printf("found field \n"); double sum = 0;
// std::set<MEdge, Less_Edge>::iterator it = edges.begin(); for (; it !=edges.end();++it){
// double sum = 0; MVertex *v0 = it->getVertex(0);
// for (; it !=edges.end();++it){ MVertex *v1 = it->getVertex(1);
// MVertex *v0 = it->getVertex(0); 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()));
// MVertex *v1 = it->getVertex(1); double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),0.5*(v0->z()+v1->z()),v0->onWhat());
// double l = sqrt((v0->x()-v1->x())*(v0->x()-v1->x())+ double el = l/lf;
// (v0->y()-v1->y())*(v0->y()-v1->y())+ int index = (int) ceil(el*nbS*0.5);
// (v0->z()-v1->z())*(v0->z()-v1->z())); qualE[index]+= 1;
// double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),0.5*(v0->z()+v1->z()),v0->onWhat()); double e = (l>lf) ? lf/l : l/lf;
// double el = l/lf; sum += e - 1.0;
// int index = (int) ceil(el*nbS*0.5); }
// qualE[index]+= 1; double tau = exp ((1./edges.size()) * sum);
// double e = (l>lf) ? lf/l : l/lf; printf("N edges = %d tau = %g\n",(int)edges.size(),tau);
// sum += e - 1.0; double ibegin = 2./(2*nbS);
// } double inext = 2./nbS;
// double tau = exp ((1./edges.size()) * sum); for (int i= 0; i< qualE.size(); i++){
// printf("N edges = %d tau = %g\n",(int)edges.size(),tau); fprintf(fp, "0 0 0 0 %g 0 0 %g \n", ibegin+i*inext , (double)qualE[i]/nbEdges);
}
// double ibegin = 2./(2*nbS); }
// double inext = 2./nbS; fclose(fp);
// 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 //emi end hack
int num = 0; int num = 0;
......
Mesh/BackgroundMesh.cpp
+ 1
− 1
View file @ 516edae1
...@@ -727,7 +727,7 @@ double backgroundMesh::operator() (double u, double v, double w) const ...@@ -727,7 +727,7 @@ double backgroundMesh::operator() (double u, double v, double w) const
MElement *e = _octree->find(u, v, w, 2, true); MElement *e = _octree->find(u, v, w, 2, true);
if (!e) { if (!e) {
Msg::Error("cannot find %g %g", u, v); Msg::Error("cannot find %g %g", u, v);
return -1000.0; return -1000.0;//0.4;
} }
e->xyz2uvw(uv, uv2); e->xyz2uvw(uv, uv2);
std::map<MVertex*,double>::const_iterator itv1 = _sizes.find(e->getVertex(0)); std::map<MVertex*,double>::const_iterator itv1 = _sizes.find(e->getVertex(0));
......
Mesh/CMakeLists.txt
+ 1
− 0
View file @ 516edae1
...@@ -30,6 +30,7 @@ set(SRC ...@@ -30,6 +30,7 @@ set(SRC
multiscalePartition.cpp multiscalePartition.cpp
QuadTriUtils.cpp QuadTriUtils.cpp
QuadTriExtruded2D.cpp QuadTriExtruded3D.cpp QuadTriTransfinite3D.cpp QuadTriExtruded2D.cpp QuadTriExtruded3D.cpp QuadTriTransfinite3D.cpp
Levy3D.cpp
) )
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h) file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
......
Mesh/meshGFaceLloyd.cpp
+ 26
− 30
View file @ 516edae1
...@@ -1014,6 +1014,7 @@ void lpcvt::print_segment(SPoint2 p1,SPoint2 p2,std::ofstream& file){ ...@@ -1014,6 +1014,7 @@ void lpcvt::print_segment(SPoint2 p1,SPoint2 p2,std::ofstream& file){
void lpcvt::compute_parameters(GFace* gf,int p){ void lpcvt::compute_parameters(GFace* gf,int p){
double h1,h2,h3; double h1,h2,h3;
double rho1,rho2,rho3;
double k; double k;
double ratio; double ratio;
double angle; double angle;
...@@ -1038,13 +1039,16 @@ void lpcvt::compute_parameters(GFace* gf,int p){ ...@@ -1038,13 +1039,16 @@ void lpcvt::compute_parameters(GFace* gf,int p){
h1 = k*backgroundMesh::current()->operator()(p1.x(),p1.y(),0.0)*ratio; h1 = k*backgroundMesh::current()->operator()(p1.x(),p1.y(),0.0)*ratio;
h2 = k*backgroundMesh::current()->operator()(p2.x(),p2.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; 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); angle = -backgroundMesh::current()->getAngle(p1.x(),p1.y(),0.0);
cosinus = cos(angle); cosinus = cos(angle);
sinus = sin(angle); sinus = sin(angle);
m = metric(cosinus,-sinus,sinus,cosinus); m = metric(cosinus,-sinus,sinus,cosinus);
v1.set_h(h1); v1.set_rho(rho1);
v2.set_h(h2); v2.set_rho(rho2);
v3.set_h(h3); v3.set_rho(rho3);
it->set_v1(v1); it->set_v1(v1);
it->set_v2(v2); it->set_v2(v2);
it->set_v3(v3); it->set_v3(v3);
...@@ -1207,7 +1211,7 @@ double lpcvt::F(voronoi_element element,int p){ ...@@ -1207,7 +1211,7 @@ double lpcvt::F(voronoi_element element,int p){
y = Ty(u,v,generator,C1,C2); y = Ty(u,v,generator,C1,C2);
point = SPoint2(x,y); point = SPoint2(x,y);
weight = gauss_weights(i,0); 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 = energy + weight*rho*f(generator,point,m,p);
} }
energy = J(generator,C1,C2)*energy; energy = J(generator,C1,C2)*energy;
...@@ -1247,7 +1251,7 @@ SVector3 lpcvt::simple(voronoi_element element,int p){ ...@@ -1247,7 +1251,7 @@ SVector3 lpcvt::simple(voronoi_element element,int p){
y = Ty(u,v,generator,C1,C2); y = Ty(u,v,generator,C1,C2);
point = SPoint2(x,y); point = SPoint2(x,y);
weight = gauss_weights(i,0); 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_x = comp_x + weight*rho*df_dx(generator,point,m,p);
comp_y = comp_y + weight*rho*df_dy(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){ ...@@ -1292,7 +1296,7 @@ SVector3 lpcvt::dF_dC1(voronoi_element element,int p){
y = Ty(u,v,generator,C1,C2); y = Ty(u,v,generator,C1,C2);
point = SPoint2(x,y); point = SPoint2(x,y);
weight = gauss_weights(i,0); 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_dx = element.get_drho_dx();
drho_dy = element.get_drho_dy(); drho_dy = element.get_drho_dy();
distance = f(point,generator,m,p); distance = f(point,generator,m,p);
...@@ -1342,7 +1346,7 @@ SVector3 lpcvt::dF_dC2(voronoi_element element,int p){ ...@@ -1342,7 +1346,7 @@ SVector3 lpcvt::dF_dC2(voronoi_element element,int p){
y = Ty(u,v,generator,C1,C2); y = Ty(u,v,generator,C1,C2);
point = SPoint2(x,y); point = SPoint2(x,y);
weight = gauss_weights(i,0); 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_dx = element.get_drho_dx();
drho_dy = element.get_drho_dy(); drho_dy = element.get_drho_dy();
distance = f(point,generator,m,p); distance = f(point,generator,m,p);
...@@ -1586,8 +1590,8 @@ bool voronoi_vertex::get_duplicate(){ ...@@ -1586,8 +1590,8 @@ bool voronoi_vertex::get_duplicate(){
return duplicate; return duplicate;
} }
double voronoi_vertex::get_h(){ double voronoi_vertex::get_rho(){
return h; return rho;
} }
void voronoi_vertex::set_point(SPoint2 new_point){ void voronoi_vertex::set_point(SPoint2 new_point){
...@@ -1614,8 +1618,8 @@ void voronoi_vertex::set_duplicate(bool new_duplicate){ ...@@ -1614,8 +1618,8 @@ void voronoi_vertex::set_duplicate(bool new_duplicate){
duplicate = new_duplicate; duplicate = new_duplicate;
} }
void voronoi_vertex::set_h(double new_h){ void voronoi_vertex::set_rho(double new_rho){
h = new_h; rho = new_rho;
} }
...@@ -1644,18 +1648,16 @@ voronoi_vertex voronoi_element::get_v3(){ ...@@ -1644,18 +1648,16 @@ voronoi_vertex voronoi_element::get_v3(){
return v3; return v3;
} }
double voronoi_element::get_rho(double u,double v,int p){ double voronoi_element::get_rho(double u,double v){
double h1; double rho1;
double h2; double rho2;
double h3; double rho3;
double h;
double rho; double rho;
h1 = v1.get_h(); rho1 = v1.get_rho();
h2 = v2.get_h(); rho2 = v2.get_rho();
h3 = v3.get_h(); rho3 = v3.get_rho();
h = h1*(1.0-u-v) + h2*u + h3*v; rho = rho1*(1.0-u-v) + rho2*u + rho3*v;
rho = compute_rho(h,p);
return rho; return rho;
} }
...@@ -1688,9 +1690,6 @@ void voronoi_element::set_metric(metric new_m){ ...@@ -1688,9 +1690,6 @@ void voronoi_element::set_metric(metric new_m){
} }
void voronoi_element::deriv_rho(int p){ void voronoi_element::deriv_rho(int p){
double h1;
double h2;
double h3;
double rho1; double rho1;
double rho2; double rho2;
double rho3; double rho3;
...@@ -1709,12 +1708,9 @@ void voronoi_element::deriv_rho(int p){ ...@@ -1709,12 +1708,9 @@ void voronoi_element::deriv_rho(int p){
SPoint2 p2; SPoint2 p2;
SPoint2 p3; SPoint2 p3;
h1 = v1.get_h(); rho1 = v1.get_rho();
h2 = v2.get_h(); rho2 = v2.get_rho();
h3 = v3.get_h(); rho3 = v3.get_rho();
rho1 = compute_rho(h1,p);
rho2 = compute_rho(h2,p);
rho3 = compute_rho(h3,p);
p1 = v1.get_point(); p1 = v1.get_point();
p2 = v2.get_point(); p2 = v2.get_point();
p3 = v3.get_point(); p3 = v3.get_point();
......
Mesh/meshGFaceLloyd.h
+ 4
− 4
View file @ 516edae1
...@@ -126,7 +126,7 @@ class voronoi_vertex{ ...@@ -126,7 +126,7 @@ class voronoi_vertex{
int index3; int index3;
SVector3 normal; SVector3 normal;
bool duplicate; bool duplicate;
double h; double rho;
public : public :
voronoi_vertex(SPoint2); voronoi_vertex(SPoint2);
voronoi_vertex(); voronoi_vertex();
...@@ -137,14 +137,14 @@ class voronoi_vertex{ ...@@ -137,14 +137,14 @@ class voronoi_vertex{
int get_index3(); int get_index3();
SVector3 get_normal(); SVector3 get_normal();
bool get_duplicate(); bool get_duplicate();
double get_h(); double get_rho();
void set_point(SPoint2); void set_point(SPoint2);
void set_index1(int); void set_index1(int);
void set_index2(int); void set_index2(int);
void set_index3(int); void set_index3(int);
void set_normal(SVector3); void set_normal(SVector3);
void set_duplicate(bool); void set_duplicate(bool);
void set_h(double); void set_rho(double);
}; };
class voronoi_element{ class voronoi_element{
...@@ -162,7 +162,7 @@ class voronoi_element{ ...@@ -162,7 +162,7 @@ class voronoi_element{
voronoi_vertex get_v1(); voronoi_vertex get_v1();
voronoi_vertex get_v2(); voronoi_vertex get_v2();
voronoi_vertex get_v3(); voronoi_vertex get_v3();
double get_rho(double,double,int); double get_rho(double,double);
double get_drho_dx(); double get_drho_dx();
double get_drho_dy(); double get_drho_dy();
metric get_metric(); metric get_metric();
......
gmshpy/gmshMesh.i
+ 2
− 0
View file @ 516edae1
...@@ -10,6 +10,7 @@ ...@@ -10,6 +10,7 @@
#include "meshGFaceLloyd.h" #include "meshGFaceLloyd.h"
#include "meshGFaceOptimize.h" #include "meshGFaceOptimize.h"
#include "meshPartitionOptions.h" #include "meshPartitionOptions.h"
#include "Levy3D.h"
%} %}
%include "GmshConfig.h" %include "GmshConfig.h"
...@@ -18,3 +19,4 @@ ...@@ -18,3 +19,4 @@
%include "meshGFaceLloyd.h" %include "meshGFaceLloyd.h"
%include "meshGFaceOptimize.h" %include "meshGFaceOptimize.h"
%include "meshPartitionOptions.h" %include "meshPartitionOptions.h"
%include "Levy3D.h"
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment