From 43433727451fff6d7b3b19dace07bbcd4593c511 Mon Sep 17 00:00:00 2001
From: Tristan Carrier Baudouin <tristan.carrier@uclouvain.be>
Date: Fri, 4 Mar 2011 10:33:31 +0000
Subject: [PATCH] lpcvt
---
CMakeLists.txt | 5 +-
Common/LuaBindings.cpp | 3 +
Common/gmshpy.i | 2 +
Geo/GModel.h | 2 +-
Mesh/meshGFaceLloyd.cpp | 1385 +++++++++++++++++++++++++---------
Mesh/meshGFaceLloyd.h | 92 ++-
Mesh/meshGFaceOptimize.cpp | 182 ++++-
Mesh/meshGFaceOptimize.h | 18 +
Numeric/DivideAndConquer.cpp | 36 +-
Numeric/DivideAndConquer.h | 14 +-
10 files changed, 1366 insertions(+), 373 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7ed51dcdaa..78372cafba 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -43,7 +43,7 @@ option(ENABLE_OCC "Enable Open CASCADE geometrical models" ON)
option(ENABLE_ACIS "Enable ACIS geometrical models" ON)
option(ENABLE_OSMESA "Use OSMesa for offscreen rendering" OFF)
option(ENABLE_PARSER "Build the GEO file parser" ON)
-option(ENABLE_PETSC "Enable PETSc linear algebra solvers" ON)
+option(ENABLE_PETSC "Enable PETSc linear algebra solvers" OFF)
option(ENABLE_POST "Build the post-processing module" ON)
option(ENABLE_PLUGINS "Build the post-processing plugins" ON)
option(ENABLE_QT "Build QT GUI" OFF)
@@ -439,6 +439,9 @@ if(ENABLE_CHACO)
set_config_option(HAVE_CHACO "Chaco")
endif(ENABLE_CHACO)
+add_subdirectory(contrib/lbfgs)
+include_directories(contrib/lbfgs)
+
if(ENABLE_DINTEGRATION)
add_subdirectory(contrib/DiscreteIntegration)
include_directories(contrib/DiscreteIntegration)
diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 6b033c9c9b..6d6c4fe586 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -33,6 +33,8 @@
#include "polynomialBasis.h"
#include "Gauss.h"
#include "meshPartitionOptions.h"
+#include "meshGFaceOptimize.h"
+#include "meshGFaceLloyd.h"
#if defined(HAVE_OPENGL)
#include "drawContext.h"
@@ -418,6 +420,7 @@ binding::binding()
polynomialBasis::registerBindings(this);
gaussIntegration::registerBindings(this);
meshPartitionOptions::registerBindings(this);
+ Temporary::registerBindings(this);
#if defined(HAVE_SOLVER)
function::registerBindings(this);
linearSystem<double>::registerBindings(this);
diff --git a/Common/gmshpy.i b/Common/gmshpy.i
index 1136a7fc06..c225e6672e 100644
--- a/Common/gmshpy.i
+++ b/Common/gmshpy.i
@@ -31,6 +31,7 @@
#include "meshPartitionOptions.h"
#include "linearSystemCSR.h"
#include "elasticitySolver.h"
+ #include "meshGFaceLloyd.h"
#include "PView.h"
#include "PViewData.h"
#include "PViewFactory.h"
@@ -125,5 +126,6 @@ namespace std {
%include "SPoint2.h"
%include "GPoint.h"
%include "functionPython.h"
+%include "meshGFaceLloyd.h"
%include "DefaultOptions.h"
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 7297e5beca..4595ab2ff5 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -404,7 +404,7 @@ class GModel
double angle);
GEntity *extrude(GEntity *e, std::vector<double> p1, std::vector<double> p2);
GEntity *addPipe(GEntity *e, std::vector<GEdge *> edges);
-
+
void addRuledFaces(std::vector<std::vector<GEdge *> > edges);
GFace *addFace(std::vector<GEdge *> edges, std::vector< std::vector<double > > points);
GFace *addPlanarFace(std::vector<std::vector<GEdge *> > edges);
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index 4f3a57b818..9e865247ae 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -8,38 +8,89 @@
#include "Context.h"
#include "meshGFace.h"
#include "BackgroundMesh.h"
+#include <fstream>
+#include "ap.h"
+#include "alglibinternal.h"
+#include "alglibmisc.h"
+#include "linalg.h"
+#include "optimization.h"
-/*List of parameters to modify*/
-//Line 76 : the variable "number" contains the number of points to remove.
-//Line 165 : the variable "test" is printed by the script.
-//Line 190 : the variable "angle" specifies the orientation of the mesh.
-//Line 287 : the variable "num" contains the number of iteration of the bisection method
-//used to find the fixed point.
-//Line 447 : it is the function "optimize" that needs to be called by the script.
-//Line 455 : put "lloyd(face)" in commentary to avoid calling the l1 lloyds method.
-//Line 456 : put "recombineIntoQuads(face,1,1)" in commentary to avoid calling the quad
-//meshing method.
-//Line 461 : the code below line 461 refers to Lua and might need to be removed.
-
-void lloydAlgorithm::operator () (GFace *gf)
+
+
+/****************fonction callback****************/
+/*************************************************/
+
+void function1_grad(const alglib::real_1d_array& x,double& func,alglib::real_1d_array& grad,void* ptr){
+ int i;
+ int p;
+ int num;
+ int index;
+ int dimension;
+ int identificator;
+ GFace* gf;
+ DocRecord* pointer;
+
+ pointer = static_cast<DocRecord*>(ptr);
+ p = pointer->get_p();
+ dimension = pointer->get_dimension();
+ gf = pointer->get_face();
+ num = pointer->numPoints;
+
+ std::vector<SVector3> gradients(num);
+ std::vector<double> energies(num);
+ std::vector<SVector3> area_centroids(num);
+ std::vector<double> areas(num);
+
+ index = 0;
+ for(i=0;i<num;i++){
+ PointRecord &pt = pointer->points[i];
+ MVertex *v = (MVertex*)pt.data;
+ if(v->onWhat()==gf && !pointer->onHull(i)){
+ pointer->points[i].where.h = x[index];
+ pointer->points[i].where.v = x[index + dimension/2];
+ pointer->points[i].identificator = index;
+ index++;
+ }
+ }
+
+ pointer->Voronoi();
+
+ lloydAlgorithm lloyd;
+ lloyd.eval(*pointer,gf,gradients,energies,area_centroids,areas,p);
+ func = lloyd.total_energy(energies);
+ printf("%f\n",1E18*func);
+
+ for(i=0;i<num;i++){
+ PointRecord &pt = pointer->points[i];
+ MVertex *v = (MVertex*)pt.data;
+ if(v->onWhat()==gf && !pointer->onHull(i)){
+ identificator = pointer->points[i].identificator;
+ grad[identificator] = gradients[i].x();
+ grad[identificator + dimension/2] = gradients[i].y();
+ }
+ }
+}
+
+void topology(const alglib::real_1d_array &x,double func,void *ptr){
+}
+
+
+
+/****************fonctions principales****************/
+/*****************************************************/
+
+void lloydAlgorithm::operator () (GFace *gf)
{
std::set<MVertex*> all;
- // get all the points of the face ...
for (unsigned int i = 0; i < gf->getNumMeshElements(); i++){
MElement *e = gf->getMeshElement(i);
for (int j = 0;j<e->getNumVertices(); j++){
MVertex *v = e->getVertex(j);
- //if (v->onWhat()->dim() < 2){
- all.insert(v);
- //}
+ all.insert(v);
}
}
- backgroundMesh::set(gf);
- backgroundMesh::current()->print("bgm.pos", 0);
-
- // Create a triangulator
DocRecord triangulator(all.size());
Range<double> du = gf->parBounds(0) ;
@@ -48,9 +99,6 @@ void lloydAlgorithm::operator () (GFace *gf)
const double LC2D = sqrt((du.high()-du.low())*(du.high()-du.low()) +
(dv.high()-dv.low())*(dv.high()-dv.low()));
- //printf("Lloyd on face %d %d elements %d nodes LC %g\n", gf->tag(),
- // gf->getNumMeshElements(), (int)all.size(), LC2D);
-
int i = 0;
for (std::set<MVertex*>::iterator it = all.begin(); it != all.end(); ++it){
SPoint2 p;
@@ -72,10 +120,10 @@ void lloydAlgorithm::operator () (GFace *gf)
triangulator.Voronoi();
triangulator.initialize();
int index,number,count,max;
- bool recent;
- number = 0; //This variable contains the number of points to remove.
+ bool flag;
+ number = 0;
count = 0;
- max = 100000;
+ max = 1000;
for(int i=0;i<max;i++)
{
if(count>=number) break;
@@ -83,389 +131,1020 @@ void lloydAlgorithm::operator () (GFace *gf)
PointRecord& pt = triangulator.points[index];
MVertex* v = (MVertex*)pt.data;
if(v->onWhat()==gf && !triangulator.onHull(index)){
- recent = triangulator.remove_point(index);
- if(recent) count++;
+ flag = triangulator.remove_point(index);
+ if(flag) count++;
}
}
- triangulator.remove_all();
+ triangulator.remove_all();
+
+ /*triangulator.Voronoi();
+ double delta;
+ delta = 0.2;
+ for(int i=0;i<triangulator.numPoints;i++){
+ PointRecord& pt = triangulator.points[i];
+ MVertex* v = (MVertex*)pt.data;
+ if(v->onWhat()==gf && !triangulator.onHull(i)){
+ triangulator.points[i].where.h = delta + (1.0-2.0*delta)*((double)rand()/(double)RAND_MAX);
+ triangulator.points[i].where.v = delta + (1.0-2.0*delta)*((double)rand()/(double)RAND_MAX);
+ }
+ }*/
+
+ int exponent = 8;
- // compute the Voronoi diagram
triangulator.Voronoi();
- //printf("hullSize = %d\n",triangulator.hullSize());
triangulator.makePosView("LloydInit.pos");
- //triangulator.printMedialAxis("medialAxis.pos");
+
+ int num_interior;
+ double initial_conditions[2*triangulator.numPoints];
+ num_interior = 0;
+ for(int i=0;i<triangulator.numPoints;i++){
+ PointRecord& pt = triangulator.points[i];
+ MVertex* v = (MVertex*)pt.data;
+ if(v->onWhat()==gf && !triangulator.onHull(i)){
+ num_interior++;
+ }
+ }
+ index = 0;
+ for(int i=0;i<triangulator.numPoints;i++){
+ PointRecord& pt = triangulator.points[i];
+ MVertex* v = (MVertex*)pt.data;
+ if(v->onWhat()==gf && !triangulator.onHull(i)){
+ initial_conditions[index] = triangulator.points[i].where.h;
+ initial_conditions[index+num_interior] = triangulator.points[i].where.v;
+ index++;
+ }
+ }
+ alglib::real_1d_array x;
+ x.setcontent(2*num_interior,initial_conditions);
- // now do the Lloyd iterations
- int ITER = 0;
- while (1){
- // store the new coordinates there
- fullMatrix<double> cgs(triangulator.numPoints,2);
- // now iterate on internal vertices
- double ENERGY = 0.0;
- double criteria = 0.0;
- for (int i=0; i<triangulator.numPoints;i++){
- // get the ith vertex
- PointRecord &pt = triangulator.points[i];
- MVertex *v = (MVertex*)pt.data;
- if (v->onWhat() == gf && !triangulator.onHull(i)){
- // get the voronoi corners
- std::vector<SPoint2> pts;
- triangulator.voronoiCell (i,pts);
- double E, A;
- SPoint2 p(pt.where.h,pt.where.v);
- for(int k=0;k<pts.size();k++){
- rotate(p,pts[k]);
- }
- SPoint2 point = fixed_point(pts); //Fixed point for quad mesh generation
- unrotate(p,point);
- cgs(i,0) = point.x();
- cgs(i,1) = point.y();
- ENERGY += E;
- double d = sqrt((p.x()-cgs(i,0))*(p.x()-cgs(i,0))+
- (p.y()-cgs(i,1))*(p.y()-cgs(i,1)));
- criteria += d/A;
- }// if (v->onWhat() == gf)
- else {
- }
- }// for all points
-
- for(PointNumero i = 0; i < triangulator.numPoints; i++) {
- MVertex *v = (MVertex*)triangulator.points[i].data;
- if (v->onWhat() == gf && !triangulator.onHull(i)){
- triangulator.points[i].where.h = cgs(i,0);
- triangulator.points[i].where.v = cgs(i,1);
- }
- }
-
- Msg::Debug("GFace %d Lloyd-iter %d Inertia=%g Convergence=%g ", gf->tag(),
- ITER++, ENERGY, criteria);
- if (ITER > ITER_MAX)break;
-
- // compute the Voronoi diagram
- triangulator.Voronoi();
-
- if (ITER % 10 == 0){
- char name[234];
- sprintf(name,"LloydIter%d.pos",ITER);
- triangulator.makePosView(name);
- }
+ triangulator.set_p(exponent);
+ triangulator.set_dimension(2*num_interior);
+ triangulator.set_face(gf);
+
+ double epsg = 0;
+ double epsf = 0;
+ double epsx = 0;
+ alglib::ae_int_t maxits = 0;
+ alglib::minlbfgsstate state;
+ alglib::minlbfgsreport rep;
+
+ minlbfgscreate(2*num_interior,4,x,state);
+ minlbfgssetcond(state,epsg,epsf,epsx,maxits);
+ minlbfgsoptimize(state,function1_grad,NULL,&triangulator);
+ minlbfgsresults(state,x,rep);
+
+ index = 0;
+ for(i=0;i<triangulator.numPoints;i++){
+ PointRecord &pt = triangulator.points[i];
+ MVertex *v = (MVertex*)pt.data;
+ if(v->onWhat()==gf && !triangulator.onHull(i)){
+ triangulator.points[i].where.h = x[index];
+ triangulator.points[i].where.v = x[index + num_interior];
+ index++;
+ }
}
+ triangulator.Voronoi();
- // now create the vertices
test = 0;
std::vector<MVertex*> mesh_vertices;
for (int i=0; i<triangulator.numPoints;i++){
- // get the ith vertex
PointRecord &pt = triangulator.points[i];
MVertex *v = (MVertex*)pt.data;
if (v->onWhat() == gf && !triangulator.onHull(i)){
GPoint gp = gf->point (pt.where.h,pt.where.v);
MFaceVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
- mesh_vertices.push_back(v);
+ mesh_vertices.push_back(v);
test++;
}
}
- // destroy the mesh
deMeshGFace killer;
killer(gf);
- // put all additional vertices in the list of
- // vertices
gf->_additional_vertices = mesh_vertices;
- // remesh the face with all those vertices in
Msg::Info("Lloyd remeshing of face %d ", gf->tag());
meshGFace mesher;
mesher(gf);
- // assign those vertices to the face internal vertices
gf->mesh_vertices.insert(gf->mesh_vertices.begin(),
gf->_additional_vertices.begin(),
gf->_additional_vertices.end());
- // clear the list of additional vertices
- gf->_additional_vertices.clear();
-}
-
-double lloydAlgorithm::get_angle(double x,double y){
- double angle;
- angle = M_PI/4.0; //This variable specifies the orientation of the mesh.
- return angle;
-}
-
-//This function returns the a, b, c and d components of the rotation matrix.
-void lloydAlgorithm::get_rotation(double x,double y,double& a,double& b,double& c,double& d){
- double angle;
- angle = get_angle(x,y);
- a = cos(angle);
- b = -sin(angle);
- c = sin(angle);
- d = cos(angle);
-}
-
-//This function returns the a, b, c and d components of the inverse rotation matrix.
-void lloydAlgorithm::get_unrotation(double x,double y,double& _a,double& _b,double& _c,double& _d){
- double a,b,c,d;
- double det;
- get_rotation(x,y,a,b,c,d);
- det = 1.0/(a*d-b*c);
- _a = d/det;
- _b = -b/det;
- _c = -c/det;
- _d = a/det;
-}
-
-//This function rotates a vertice (point) belonging to a Voronoi cell.
-//The angle is specified by the location of the generator, because we want all vertices
-//of the Voronoi cell to be rotated by the same amount in case the orientation of the
-//mesh is not constant.
-//The generator itself is not rotated.
-void lloydAlgorithm::rotate(SPoint2 generator,SPoint2& point){
- double x,y;
- double new_x,new_y;
- double a,b,c,d;
- get_rotation(generator.x(),generator.y(),a,b,c,d);
- x = point.x();
- y = point.y();
- new_x = a*x + b*y;
- new_y = c*x + d*y;
- point.setPosition(new_x,new_y);
-}
-
-//This function rotates a vertice (point) belonging to a Voronoi cell in the reverse direction.
-//The angle is specified by the location of the generator, because we want all vertices
-//of the Voronoi cell to be rotated by the same amount in case the orientation of the
-//mesh is not constant.
-//The generator itself is not rotated.
-void lloydAlgorithm::unrotate(SPoint2 generator,SPoint2& point){
- double x,y;
- double new_x,new_y;
- double a,b,c,d;
- get_unrotation(generator.x(),generator.y(),a,b,c,d);
- x = point.x();
- y = point.y();
- new_x = a*x + b*y;
- new_y = c*x + d*y;
- point.setPosition(new_x,new_y);
-}
-
-//This function returns the fixed point of the voronoi cell specified in argument (polygon).
-//It finds the bounding box of the voronoi cell and then it calls the function range in x and y.
-SPoint2 lloydAlgorithm::fixed_point(const std::vector<SPoint2>& polygon){
- int i,num;
- double x,y;
- double x_min,y_min;
- double x_max,y_max;
- double solution_x,solution_y;
- SPoint2 point;
- x_min = 100000000.0;
- y_min = 100000000.0;
- x_max = -100000000.0;
- y_max = -100000000.0;
- num = polygon.size();
- for(i=0;i<num;i++){
- point = polygon[i];
- x = point.x();
- y = point.y();
- if(x < x_min) x_min = x;
- if(y < y_min) y_min = y;
- if(x > x_max) x_max = x;
- if(y > y_max) y_max = y;
- }
- solution_x = range(polygon,x_min,x_max,1);
- solution_y = range(polygon,y_min,y_max,0);
- return SPoint2(solution_x,solution_y);
-}
-
-//This function returns the x or y median of the Voronoi cell specified in argument (polygon).
-//If direction = 1, the function will return the x median.
-//If direction = 0, the function will return the y median.
-//The variables u1 and u2 should be respectively equal to the minimum and the maximum of x or y,
-//depending of the value of direction.
-//This function uses the bisection method.
-double lloydAlgorithm::range(const std::vector<SPoint2>& polygon,double u1,double u2,bool direction){
- int i,num;
- double mid,val;
- num = 16; //This variable specifies the number of iteration of the bisection method.
- for(i=0;i<num;i++){
- mid = 0.5*(u1 + u2);
- val = gradient(polygon,mid,direction);
- if(val >= 0) u2 = mid;
- else u1 = mid;
- }
- return 0.5*(u1 + u2);
-}
-
-//This function returns the x or y component of the l1 gradient of the voronoi cell
-//specified in argument (polygon).
-//If direction = 1, it will return the x component.
-//If direction = 0, it will return the y component.
-//The variable generator should be equal to the x or y component of the point to which the
-//voronoi cell belongs.
-//x component = A1-A2, A1 : area of the left part, A2 : area of the right part
-//y component = A3-A4, A3 : area of the lower part, A4 : area of the higher part
-double lloydAlgorithm::gradient(const std::vector<SPoint2>& polygon,double generator,bool direction){
- int i,num,index1,index2;
- bool intersection,side;
- double areaA,areaB;
- SPoint2 p1,p2,solution;
- std::vector<bool> category;
- std::vector<SPoint2> new_polygon;
- std::vector<SPoint2> polygonA;
- std::vector<SPoint2> polygonB;
- num = polygon.size();
- for(i=0;i<num;i++){
- index1 = i;
- index2 = (i+1)%num;
- p1 = polygon[index1];
- p2 = polygon[index2];
- solution = line_intersection(p1,p2,generator,direction,intersection);
- if(intersection){
- new_polygon.push_back(p1);
- new_polygon.push_back(solution);
- category.push_back(1);
- category.push_back(0);
- }
- else{
- new_polygon.push_back(p1);
- category.push_back(1);
+ gf->_additional_vertices.clear();
+}
+
+double lloydAlgorithm::optimize(int max,int flag){
+ GFace*face;
+ GModel*model = GModel::current();
+ GModel::fiter iterator;
+ lloydAlgorithm lloyd(max,flag);
+ for(iterator = model->firstFace();iterator != model->lastFace();iterator++)
+ {
+ face = *iterator;
+ if(face->getNumMeshElements() > 0){
+ lloyd(face);
+ recombineIntoQuads(face,1,1);
}
}
- num = new_polygon.size();
- for(i=0;i<num;i++){
- if(category[i]){
- side = point_side(new_polygon[i],generator,direction);
- if(side){
- polygonA.push_back(new_polygon[i]);
+ return lloyd.test;
+}
+
+
+
+/****************Calcul du gradient****************/
+/**************************************************/
+
+void lloydAlgorithm::eval(DocRecord& triangulator,GFace* gf,std::vector<SVector3>& gradients,std::vector<double>& energies,std::vector<SVector3>& area_centroids,std::vector<double>& areas,int p){
+ int i;
+
+ for(i=0;i<triangulator.numPoints;i++){
+ gradients[i] = SVector3(0.0,0.0,0.0);
+ energies[i] = 0.0;
+ area_centroids[i] = SVector3(0.0,0.0,0.0);
+ areas[i] = 0.0;
+ }
+
+ for(i=0;i<triangulator.numPoints;i++){
+ PointRecord &pt = triangulator.points[i];
+ MVertex *v = (MVertex*)pt.data;
+ if(v->onWhat()==gf && !triangulator.onHull(i)){
+ bool is_inside1;
+ bool is_inside2;
+ bool flag;
+ int num = triangulator._adjacencies[i].t_length;
+ for(int j=0;j<num;j++){
+ int index3 = triangulator._adjacencies[i].t[j];
+ int index2 = triangulator._adjacencies[i].t[(j+1)%num];
+ int index1 = triangulator._adjacencies[i].t[(j+2)%num];
+ double _x0[2] = {triangulator.points[i].where.h,triangulator.points[i].where.v};
+ double _x1[2] = {triangulator.points[index1].where.h,triangulator.points[index1].where.v};
+ double _x2[2] = {triangulator.points[index2].where.h,triangulator.points[index2].where.v};
+ double _x3[2] = {triangulator.points[index3].where.h,triangulator.points[index3].where.v};
+ SPoint2 x0(_x0[0],_x0[1]);
+ SPoint2 x1(_x1[0],_x1[1]);
+ SPoint2 x2(_x2[0],_x2[1]);
+ SPoint2 x3(_x3[0],_x3[1]);
+ double _C1[2];
+ double _C2[2];
+ circumCenterXY(_x0,_x1,_x2,_C1);
+ circumCenterXY(_x0,_x2,_x3,_C2);
+ SPoint2 C1(_C1[0],_C1[1]);
+ SPoint2 C2(_C2[0],_C2[1]);
+ is_inside1 = inside_domain(C1.x(),C1.y());
+ is_inside2 = inside_domain(C2.x(),C2.y());
+ if(is_inside1 && is_inside2){
+ energies[i] = energies[i] + F(x0,C1,C2,p);
+ gradients[i] = gradients[i] + simple(x0,C1,C2,p);
+ area_centroids[i] = area_centroids[i] + area_centroid(x0,C1,C2);
+ areas[i] = areas[i] + triangle_area(x0,C1,C2);
+ SVector3 grad1 = dF_dC1(x0,C1,C2,p);
+ SVector3 grad2 = dF_dC2(x0,C1,C2,p);
+ gradients[i] = gradients[i] + inner_dFdx0(grad1,C1,x0,x1,x2);
+ gradients[index1] = gradients[index1] + inner_dFdx0(grad1,C1,x1,x0,x2);
+ gradients[index2] = gradients[index2] + inner_dFdx0(grad1,C1,x2,x0,x1);
+ gradients[i] = gradients[i] + inner_dFdx0(grad2,C2,x0,x2,x3);
+ gradients[index2] = gradients[index2] + inner_dFdx0(grad2,C2,x2,x0,x3);
+ gradients[index3] = gradients[index3] + inner_dFdx0(grad2,C2,x3,x0,x2);
+ }
+ else if(is_inside1){
+ SPoint2 first;
+ SPoint2 second;
+ SVector3 vecA;
+ SVector3 vecB;
+ first = boundary_intersection(C1,C2,flag,vecA);
+ second = boundary_intersection(x0,C2,flag,vecB);
+ energies[i] = energies[i] + F(x0,C1,first,p);
+ energies[i] = energies[i] + F(x0,first,second,p);
+ gradients[i] = gradients[i] + simple(x0,C1,first,p);
+ gradients[i] = gradients[i] + simple(x0,first,second,p);
+ area_centroids[i] = area_centroids[i] + area_centroid(x0,C1,first);
+ area_centroids[i] = area_centroids[i] + area_centroid(x0,first,second);
+ areas[i] = areas[i] + triangle_area(x0,C1,first);
+ areas[i] = areas[i] + triangle_area(x0,first,second);
+ SVector3 grad1 = dF_dC1(x0,C1,first,p);
+ SVector3 grad2 = dF_dC2(x0,C1,first,p);
+ SVector3 grad3 = dF_dC1(x0,first,second,p);
+ gradients[i] = gradients[i] + inner_dFdx0(grad1,C1,x0,x1,x2);
+ gradients[index1] = gradients[index1] + inner_dFdx0(grad1,C1,x1,x0,x2);
+ gradients[index2] = gradients[index2] + inner_dFdx0(grad1,C1,x2,x0,x1);
+ gradients[i] = gradients[i] + boundary_dFdx0(grad2,first,x0,x2,vecA);
+ gradients[index2] = gradients[index2] + boundary_dFdx0(grad2,first,x2,x0,vecA);
+ gradients[i] = gradients[i] + boundary_dFdx0(grad3,first,x0,x2,vecA);
+ gradients[index2] = gradients[index2] + boundary_dFdx0(grad3,first,x2,x0,vecA);
+ }
+ else if(is_inside2){
+ SPoint2 first;
+ SPoint2 second;
+ SVector3 vecA;
+ SVector3 vecB;
+ first = boundary_intersection(C1,C2,flag,vecA);
+ second = boundary_intersection(x0,C1,flag,vecB);
+ energies[i] = energies[i] + F(x0,second,first,p);
+ energies[i] = energies[i] + F(x0,first,C2,p);
+ gradients[i] = gradients[i] + simple(x0,second,first,p);
+ gradients[i] = gradients[i] + simple(x0,first,C2,p);
+ area_centroids[i] = area_centroids[i] + area_centroid(x0,second,first);
+ area_centroids[i] = area_centroids[i] + area_centroid(x0,first,C2);
+ areas[i] = areas[i] + triangle_area(x0,second,first);
+ areas[i] = areas[i] + triangle_area(x0,first,C2);
+ SVector3 grad1 = dF_dC2(x0,second,first,p);
+ SVector3 grad2 = dF_dC1(x0,first,C2,p);
+ SVector3 grad3 = dF_dC2(x0,first,C2,p);
+ gradients[i] = gradients[i] + inner_dFdx0(grad3,C2,x0,x2,x3);
+ gradients[index2] = gradients[index2] + inner_dFdx0(grad3,C2,x2,x0,x3);
+ gradients[index3] = gradients[index3] + inner_dFdx0(grad3,C2,x3,x0,x2);
+ gradients[i] = gradients[i] + boundary_dFdx0(grad2,first,x0,x2,vecA);
+ gradients[index2] = gradients[index2] + boundary_dFdx0(grad2,first,x2,x0,vecA);
+ gradients[i] = gradients[i] + boundary_dFdx0(grad1,first,x0,x2,vecA);
+ gradients[index2] = gradients[index2] + boundary_dFdx0(grad1,first,x2,x0,vecA);
+ }
+ else{
+ SPoint2 first;
+ SPoint2 second;
+ SPoint2 third;
+ SPoint2 fourth;
+ SPoint2 median;
+ SPoint2 any;
+ SVector3 vecA;
+ SVector3 vecB;
+ SVector3 vec;
+ any = boundary_intersection(C1,C2,flag,vec);
+ if(flag){
+ median = mid(x0,x2);
+ first = boundary_intersection(median,C1,flag,vecA);
+ second = boundary_intersection(median,C2,flag,vecB);
+ third = boundary_intersection(x0,C1,flag,vec);
+ fourth = boundary_intersection(x0,C2,flag,vec);
+ energies[i] = energies[i] + F(x0,first,second,p);
+ energies[i] = energies[i] + F(x0,third,first,p);
+ energies[i] = energies[i] + F(x0,second,fourth,p);
+ gradients[i] = gradients[i] + simple(x0,first,second,p);
+ gradients[i] = gradients[i] + simple(x0,third,first,p);
+ gradients[i] = gradients[i] + simple(x0,second,fourth,p);
+ area_centroids[i] = area_centroids[i] + area_centroid(x0,first,second);
+ area_centroids[i] = area_centroids[i] + area_centroid(x0,third,first);
+ area_centroids[i] = area_centroids[i] + area_centroid(x0,second,fourth);
+ areas[i] = areas[i] + triangle_area(x0,first,second);
+ areas[i] = areas[i] + triangle_area(x0,third,first);
+ areas[i] = areas[i] + triangle_area(x0,second,fourth);
+ SVector3 grad1 = dF_dC1(x0,first,second,p);
+ SVector3 grad2 = dF_dC2(x0,first,second,p);
+ SVector3 grad3 = dF_dC2(x0,third,first,p);
+ SVector3 grad4 = dF_dC1(x0,second,fourth,p);
+ gradients[i] = gradients[i] + boundary_dFdx0(grad1,first,x0,x2,vecA);
+ gradients[index2] = gradients[index2] + boundary_dFdx0(grad1,first,x2,x0,vecA);
+ gradients[i] = gradients[i] + boundary_dFdx0(grad2,second,x0,x2,vecB);
+ gradients[index2] = gradients[index2] + boundary_dFdx0(grad2,second,x2,x0,vecB);
+ gradients[i] = gradients[i] + boundary_dFdx0(grad3,first,x0,x2,vecA);
+ gradients[index2] = gradients[index2] + boundary_dFdx0(grad3,first,x2,x0,vecA);
+ gradients[i] = gradients[i] + boundary_dFdx0(grad4,second,x0,x2,vecB);
+ gradients[index2] = gradients[index2] + boundary_dFdx0(grad4,second,x2,x0,vecB);
+ }
+ else{
+ first = boundary_intersection(x0,C1,flag,vec);
+ second = boundary_intersection(x0,C2,flag,vec);
+ energies[i] = energies[i] + F(x0,first,second,p);
+ gradients[i] = gradients[i] + simple(x0,first,second,p);
+ area_centroids[i] = area_centroids[i] + area_centroid(x0,first,second);
+ areas[i] = areas[i] + triangle_area(x0,first,second);
+ }
+ }
}
- else{
- polygonB.push_back(new_polygon[i]);
+ }
+ else {
+ double e = 0.00000001;
+ bool ok_triangle1;
+ bool ok_triangle2;
+ bool is_inside1;
+ bool is_inside2;
+ bool is_inside;
+ bool flag;
+ int num = triangulator._adjacencies[i].t_length;
+ for(int j=0;j<num;j++){
+ int index3 = triangulator._adjacencies[i].t[j];
+ int index2 = triangulator._adjacencies[i].t[(j+1)%num];
+ int index1 = triangulator._adjacencies[i].t[(j+2)%num];
+ if(index1!=index2 && index1!=index3 && index2!=index3 && i!=index1 && i!=index2 && i!=index3){
+ double _x0[2] = {triangulator.points[i].where.h,triangulator.points[i].where.v};
+ double _x1[2] = {triangulator.points[index1].where.h,triangulator.points[index1].where.v};
+ double _x2[2] = {triangulator.points[index2].where.h,triangulator.points[index2].where.v};
+ double _x3[2] = {triangulator.points[index3].where.h,triangulator.points[index3].where.v};
+ SPoint2 x0(_x0[0],_x0[1]);
+ SPoint2 x1(_x1[0],_x1[1]);
+ SPoint2 x2(_x2[0],_x2[1]);
+ SPoint2 x3(_x3[0],_x3[1]);
+ ok_triangle1 = adjacent(triangulator,index1,index2) && triangle_area(x0,x1,x2)>e;
+ ok_triangle2 = adjacent(triangulator,index2,index3) && triangle_area(x0,x2,x3)>e;
+ if(ok_triangle1 && ok_triangle2){
+ double _C1[2];
+ double _C2[2];
+ circumCenterXY(_x0,_x1,_x2,_C1);
+ circumCenterXY(_x0,_x2,_x3,_C2);
+ SPoint2 C1(_C1[0],_C1[1]);
+ SPoint2 C2(_C2[0],_C2[1]);
+ is_inside1 = inside_domain(C1.x(),C1.y());
+ is_inside2 = inside_domain(C2.x(),C2.y());
+ if(is_inside1 && is_inside2){
+ energies[i] = energies[i] + F(x0,C1,C2,p);
+ SVector3 grad1 = dF_dC1(x0,C1,C2,p);
+ SVector3 grad2 = dF_dC2(x0,C1,C2,p);
+ gradients[index1] = gradients[index1] + inner_dFdx0(grad1,C1,x1,x0,x2);
+ gradients[index2] = gradients[index2] + inner_dFdx0(grad1,C1,x2,x0,x1);
+ gradients[index2] = gradients[index2] + inner_dFdx0(grad2,C2,x2,x0,x3);
+ gradients[index3] = gradients[index3] + inner_dFdx0(grad2,C2,x3,x0,x2);
+ }
+ else if(is_inside1){
+ SPoint2 val;
+ SVector3 vec;
+ val = boundary_intersection(C1,C2,flag,vec);
+ energies[i] = energies[i] + F(x0,C1,val,p);
+ SVector3 grad1 = dF_dC1(x0,C1,val,p);
+ SVector3 grad2 = dF_dC2(x0,C1,val,p);
+ gradients[index1] = gradients[index1] + inner_dFdx0(grad1,C1,x1,x0,x2);
+ gradients[index2] = gradients[index2] + inner_dFdx0(grad1,C1,x2,x0,x1);
+ gradients[index2] = gradients[index2] + boundary_dFdx0(grad2,val,x2,x0,vec);
+ }
+ else if(is_inside2){
+ SPoint2 val;
+ SVector3 vec;
+ val = boundary_intersection(C1,C2,flag,vec);
+ energies[i] = energies[i] + F(x0,val,C2,p);
+ SVector3 grad1 = dF_dC1(x0,val,C2,p);
+ SVector3 grad2 = dF_dC2(x0,val,C2,p);
+ gradients[index2] = gradients[index2] + inner_dFdx0(grad2,C2,x2,x0,x3);
+ gradients[index3] = gradients[index3] + inner_dFdx0(grad2,C2,x3,x0,x2);
+ gradients[index2] = gradients[index2] + boundary_dFdx0(grad1,val,x2,x0,vec);
+ }
+ else{
+ SPoint2 first;
+ SPoint2 second;
+ SPoint2 median;
+ SPoint2 any;
+ SVector3 vecA;
+ SVector3 vecB;
+ SVector3 vec;
+ any = boundary_intersection(C1,C2,flag,vec);
+ if(flag){
+ median = mid(x0,x2);
+ first = boundary_intersection(median,C1,flag,vecA);
+ second = boundary_intersection(median,C2,flag,vecB);
+ energies[i] = energies[i] + F(x0,first,second,p);
+ SVector3 grad1 = dF_dC1(x0,first,second,p);
+ SVector3 grad2 = dF_dC2(x0,first,second,p);
+ gradients[index2] = gradients[index2] + boundary_dFdx0(grad1,first,x2,x0,vecA);
+ gradients[index2] = gradients[index2] + boundary_dFdx0(grad2,second,x2,x0,vecB);
+ }
+ }
+ }
+ else if(ok_triangle1){
+ double _C[2];
+ circumCenterXY(_x0,_x1,_x2,_C);
+ SPoint2 C(_C[0],_C[1]);
+ is_inside = inside_domain(C.x(),C.y());
+ if(is_inside){
+ SPoint2 val;
+ val = mid(x0,x2);
+ energies[i] = energies[i] + F(x0,C,val,p);
+ SVector3 grad1 = dF_dC1(x0,C,val,p);
+ gradients[index1] = gradients[index1] + inner_dFdx0(grad1,C,x1,x0,x2);
+ }
+ }
+ else if(ok_triangle2){
+ double _C[2];
+ circumCenterXY(_x0,_x2,_x3,_C);
+ SPoint2 C(_C[0],_C[1]);
+ is_inside = inside_domain(C.x(),C.y());
+ if(is_inside){
+ SPoint2 val;
+ val = mid(x0,x2);
+ energies[i] = energies[i] + F(x0,val,C,p);
+ SVector3 grad1 = dF_dC2(x0,val,C,p);
+ gradients[index3] = gradients[index3] + inner_dFdx0(grad1,C,x3,x0,x2);
+ }
+ }
+ }
}
+ }
+ }
+}
+
+double lloydAlgorithm::total_energy(std::vector<double> energies){
+ int i;
+ double total;
+ total = 0.0;
+ for(i=0;i<energies.size();i++){
+ total = total + energies[i];
+ }
+ return total;
+}
+
+SVector3 lloydAlgorithm::numerical(DocRecord& triangulator,GFace* gf,int index,int p){
+ double x_original;
+ double y_original;
+ double energy_right;
+ double energy_left;
+ double energy_up;
+ double energy_down;
+ double comp_x;
+ double comp_y;
+ double e;
+ std::vector<double> energies(triangulator.numPoints);
+ std::vector<SVector3> gradients(triangulator.numPoints);
+ std::vector<SVector3> area_centroids(triangulator.numPoints);
+ std::vector<double> areas(triangulator.numPoints);
+ e = 0.00001;
+ x_original = triangulator.points[index].where.h;
+ y_original = triangulator.points[index].where.v;
+
+ triangulator.points[index].where.h = x_original + e;
+ triangulator.points[index].where.v = y_original;
+ eval(triangulator,gf,gradients,energies,area_centroids,areas,p);
+ energy_right = total_energy(energies);
+
+ triangulator.points[index].where.h = x_original - e;
+ triangulator.points[index].where.v = y_original;
+ eval(triangulator,gf,gradients,energies,area_centroids,areas,p);
+ energy_left = total_energy(energies);
+
+ triangulator.points[index].where.h = x_original;
+ triangulator.points[index].where.v = y_original + e;
+ eval(triangulator,gf,gradients,energies,area_centroids,areas,p);
+ energy_up = total_energy(energies);
+
+ triangulator.points[index].where.h = x_original;
+ triangulator.points[index].where.v = y_original - e;
+ eval(triangulator,gf,gradients,energies,area_centroids,areas,p);
+ energy_down = total_energy(energies);
+
+ comp_x = (energy_right - energy_left)/(2*e);
+ comp_y = (energy_up - energy_down)/(2*e);
+
+ triangulator.points[index].where.h = x_original;
+ triangulator.points[index].where.v = y_original;
+
+ return SVector3(comp_x,comp_y,0.0);
+}
+
+void lloydAlgorithm::write(DocRecord& triangulator,GFace* gf,int p){
+ int i;
+ SVector3 grad1;
+ SVector3 grad2;
+ SVector3 ecart;
+ std::vector<SVector3> gradients(triangulator.numPoints);
+ std::vector<double> energies(triangulator.numPoints);
+ std::vector<SVector3> area_centroids(triangulator.numPoints);
+ std::vector<double> areas(triangulator.numPoints);
+
+ eval(triangulator,gf,gradients,energies,area_centroids,areas,p);
+
+ std::ofstream file("gradient");
+
+ for(i=0;i<triangulator.numPoints;i++){
+ PointRecord& pt = triangulator.points[i];
+ MVertex* v = (MVertex*)pt.data;
+ if(v->onWhat()==gf && !triangulator.onHull(i)){
+ grad1 = numerical(triangulator,gf,i,p);
+ grad2 = gradients[i];
+ ecart = grad1-grad2;
+ file << grad1.x() << " ";
+ file << grad2.x() << " ";
+ file << grad1.y() << " ";
+ file << grad2.y() << " ";
+ file << 100.0*ecart.norm()/grad1.norm() << " ";
+ file << "\n";
}
- else{
- polygonA.push_back(new_polygon[i]);
- polygonB.push_back(new_polygon[i]);
+ }
+}
+
+void lloydAlgorithm::write2(DocRecord& triangulator,GFace* gf){
+ int i;
+ SPoint2 generator;
+ SVector3 grad1;
+ SVector3 grad2;
+ SVector3 ecart;
+ std::vector<SVector3> gradients(triangulator.numPoints);
+ std::vector<double> energies(triangulator.numPoints);
+ std::vector<SVector3> area_centroids(triangulator.numPoints);
+ std::vector<double> areas(triangulator.numPoints);
+
+ eval(triangulator,gf,gradients,energies,area_centroids,areas,2);
+
+ std::ofstream file("gradient2");
+
+ for(i=0;i<triangulator.numPoints;i++){
+ PointRecord& pt = triangulator.points[i];
+ MVertex* v = (MVertex*)pt.data;
+ if(v->onWhat()==gf && !triangulator.onHull(i)){
+ generator = SPoint2(triangulator.points[i].where.h,triangulator.points[i].where.v);
+ grad1 = analytical(generator,area_centroids[i],areas[i]);
+ grad2 = gradients[i];
+ ecart = grad1-grad2;
+ file << grad1.x() << " ";
+ file << grad2.x() << " ";
+ file << grad1.y() << " ";
+ file << grad2.y() << " ";
+ file << 100.0*ecart.norm()/grad1.norm() << " ";
+ file << "\n";
}
}
- areaA = polygon_area(polygonA);
- areaB = polygon_area(polygonB);
- return areaA-areaB;
-}
-
-//The point p is a vertice of a Voronoi cell.
-//The variable generator contains the x or y component of the point to which
-//the Voronoi cell belongs. If direction = 1, it will contain the x component,
-//if direction = 0, it will contain the y component.
-//If direction = 1, this function will return 1 if p is at the left of the generator
-//and it will return 0 if p is at the right of the generator.
-//If direction = 0, this function will return 1 if p is below the generator
-//and it will return 0 if p is higher than the generator.
-bool lloydAlgorithm::point_side(SPoint2 p,double generator,bool direction){
- if(direction){
- if(p.x()<=generator) return 1;
- else return 0;
+}
+
+SVector3 lloydAlgorithm::analytical(SPoint2 generator,SVector3 area_centroid,double area){
+ SVector3 _generator;
+ SVector3 val;
+ _generator = SVector3(generator.x(),generator.y(),0.0);
+ val = 2.0*area*(_generator - (1.0/area)*area_centroid);
+ return val;
+}
+
+SVector3 lloydAlgorithm::area_centroid(SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double x;
+ double y;
+ double area;
+ area = triangle_area(p1,p2,p3);
+ x = (1.0/3.0)*(p1.x() + p2.x() + p3.x());
+ y = (1.0/3.0)*(p1.y() + p2.y() + p3.y());
+ return SVector3(area*x,area*y,0.0);
+}
+
+
+
+/****************Fonctions pour le calcul du gradient****************/
+/********************************************************************/
+
+double lloydAlgorithm::F(SPoint2 generator,SPoint2 C1,SPoint2 C2,int p){
+ int i;
+ int num;
+ int order;
+ double x;
+ double y;
+ double energy;
+ SPoint2 point;
+ fullMatrix<double> pts;
+ fullMatrix<double> weights;
+ order = 8;
+ gaussIntegration::getTriangle(order,pts,weights);
+ num = pts.size1();
+ energy = 0.0;
+ for(i=0;i<num;i++){
+ x = Tx(pts(i,0),pts(i,1),generator,C1,C2);
+ y = Ty(pts(i,0),pts(i,1),generator,C1,C2);
+ point = SPoint2(x,y);
+ energy = energy + weights(i,0)*f(generator,point,p);
}
- else{
- if(p.y()<=generator) return 1;
- else return 0;
+ energy = J(generator,C1,C2)*energy;
+ return energy;
+}
+
+SVector3 lloydAlgorithm::simple(SPoint2 generator,SPoint2 C1,SPoint2 C2,int p){
+ int i;
+ int num;
+ int order;
+ double x;
+ double y;
+ double comp_x;
+ double comp_y;
+ SPoint2 point;
+ fullMatrix<double> pts;
+ fullMatrix<double> weights;
+ order = 8;
+ gaussIntegration::getTriangle(order,pts,weights);
+ num = pts.size1();
+ comp_x = 0.0;
+ comp_y = 0.0;
+ for(i=0;i<num;i++){
+ x = Tx(pts(i,0),pts(i,1),generator,C1,C2);
+ y = Ty(pts(i,0),pts(i,1),generator,C1,C2);
+ point = SPoint2(x,y);
+ comp_x = comp_x + weights(i,0)*df_dx(generator,point,p);
+ comp_y = comp_y + weights(i,0)*df_dy(generator,point,p);
}
+ comp_x = J(generator,C1,C2)*comp_x;
+ comp_y = J(generator,C1,C2)*comp_y;
+ return SVector3(comp_x,comp_y,0.0);
}
-//This function uses the Surveyors Formula to calculate the area of the Voronoi cell
-//specified in argument (polygon).
-double lloydAlgorithm::polygon_area(const std::vector<SPoint2>& polygon){
+SVector3 lloydAlgorithm::dF_dC1(SPoint2 generator,SPoint2 C1,SPoint2 C2,int p){
int i;
int num;
+ int order;
+ double x;
+ double y;
+ double comp_x;
+ double comp_y;
+ SPoint2 point;
+ fullMatrix<double> pts;
+ fullMatrix<double> weights;
+ order = 8;
+ gaussIntegration::getTriangle(order,pts,weights);
+ num = pts.size1();
+ comp_x = 0.0;
+ comp_y = 0.0;
+ for(i=0;i<num;i++){
+ x = Tx(pts(i,0),pts(i,1),generator,C1,C2);
+ y = Ty(pts(i,0),pts(i,1),generator,C1,C2);
+ point = SPoint2(x,y);
+ comp_x = comp_x + weights(i,0)*df_dx(point,generator,p)*dTx_dp2x(pts(i,0),pts(i,1))*J(generator,C1,C2);
+ comp_x = comp_x + weights(i,0)*f(point,generator,p)*dJ_dp2x(generator,C1,C2);
+ comp_y = comp_y + weights(i,0)*df_dy(point,generator,p)*dTy_dp2y(pts(i,0),pts(i,1))*J(generator,C1,C2);
+ comp_y = comp_y + weights(i,0)*f(point,generator,p)*dJ_dp2y(generator,C1,C2);
+ }
+ return SVector3(comp_x,comp_y,0.0);
+}
+
+SVector3 lloydAlgorithm::dF_dC2(SPoint2 generator,SPoint2 C1,SPoint2 C2,int p){
+ int i;
+ int num;
+ int order;
+ double x;
+ double y;
+ double comp_x;
+ double comp_y;
+ SPoint2 point;
+ fullMatrix<double> pts;
+ fullMatrix<double> weights;
+ order = 8;
+ gaussIntegration::getTriangle(order,pts,weights);
+ num = pts.size1();
+ comp_x = 0.0;
+ comp_y = 0.0;
+ for(i=0;i<num;i++){
+ x = Tx(pts(i,0),pts(i,1),generator,C1,C2);
+ y = Ty(pts(i,0),pts(i,1),generator,C1,C2);
+ point = SPoint2(x,y);
+ comp_x = comp_x + weights(i,0)*df_dx(point,generator,p)*dTx_dp3x(pts(i,0),pts(i,1))*J(generator,C1,C2);
+ comp_x = comp_x + weights(i,0)*f(point,generator,p)*dJ_dp3x(generator,C1,C2);
+ comp_y = comp_y + weights(i,0)*df_dy(point,generator,p)*dTy_dp3y(pts(i,0),pts(i,1))*J(generator,C1,C2);
+ comp_y = comp_y + weights(i,0)*f(point,generator,p)*dJ_dp3y(generator,C1,C2);
+ }
+ return SVector3(comp_x,comp_y,0.0);
+}
+
+double lloydAlgorithm::f(SPoint2 p1,SPoint2 p2,int p){
+ double val;
+ val = exp(p1.x()-p2.x(),p) + exp(p1.y()-p2.y(),p);
+ return val;
+}
+
+double lloydAlgorithm::df_dx(SPoint2 p1,SPoint2 p2,int p){
+ double val;
+ val = ((double)p)*exp(p1.x()-p2.x(),p-1);
+ return val;
+}
+
+double lloydAlgorithm::df_dy(SPoint2 p1,SPoint2 p2,int p){
+ double val;
+ val = ((double)p)*exp(p1.y()-p2.y(),p-1);
+ return val;
+}
+
+double lloydAlgorithm::L1(double u,double v){
+ double val;
+ val = 1.0-u-v;
+ return val;
+}
+
+double lloydAlgorithm::dL1_du(){
+ return -1.0;
+}
+
+double lloydAlgorithm::dL1_dv(){
+ return -1.0;
+}
+
+double lloydAlgorithm::L2(double u,double v){
+ double val;
+ val = u;
+ return val;
+}
+
+double lloydAlgorithm::dL2_du(){
+ return 1.0;
+}
+
+double lloydAlgorithm::dL2_dv(){
+ return 0.0;
+}
+
+double lloydAlgorithm::L3(double u,double v){
+ double val;
+ val = v;
+ return val;
+}
+
+double lloydAlgorithm::dL3_du(){
+ return 0.0;
+}
+
+double lloydAlgorithm::dL3_dv(){
+ return 1.0;
+}
+
+double lloydAlgorithm::Tx(double u,double v,SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double val;
+ val = L1(u,v)*p1.x() + L2(u,v)*p2.x() + L3(u,v)*p3.x();
+ return val;
+}
+
+double lloydAlgorithm::dTx_dp1x(double u,double v){
+ return L1(u,v);
+}
+
+double lloydAlgorithm::dTx_dp2x(double u,double v){
+ return L2(u,v);
+}
+
+double lloydAlgorithm::dTx_dp3x(double u,double v){
+ return L3(u,v);
+}
+
+double lloydAlgorithm::dTx_du(SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double val;
+ val = dL1_du()*p1.x() + dL2_du()*p2.x() + dL3_du()*p3.x();
+ return val;
+}
+
+double lloydAlgorithm::dTx_dv(SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double val;
+ val = dL1_dv()*p1.x() + dL2_dv()*p2.x() + dL3_dv()*p3.x();
+ return val;
+}
+
+double lloydAlgorithm::Ty(double u,double v,SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double val;
+ val = L1(u,v)*p1.y() + L2(u,v)*p2.y() + L3(u,v)*p3.y();
+ return val;
+}
+
+double lloydAlgorithm::dTy_dp1y(double u,double v){
+ return L1(u,v);
+}
+
+double lloydAlgorithm::dTy_dp2y(double u,double v){
+ return L2(u,v);
+}
+
+double lloydAlgorithm::dTy_dp3y(double u,double v){
+ return L3(u,v);
+}
+
+double lloydAlgorithm::dTy_du(SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double val;
+ val = dL1_du()*p1.y() + dL2_du()*p2.y() + dL3_du()*p3.y();
+ return val;
+}
+
+double lloydAlgorithm::dTy_dv(SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double val;
+ val = dL1_dv()*p1.y() + dL2_dv()*p2.y() + dL3_dv()*p3.y();
+ return val;
+}
+
+double lloydAlgorithm::J(SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double val;
+ val = dTx_du(p1,p2,p3)*dTy_dv(p1,p2,p3) - dTx_dv(p1,p2,p3)*dTy_du(p1,p2,p3);
+ return val;
+}
+
+double lloydAlgorithm::dJ_dp1x(SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double val;
+ val = dL1_du()*dTy_dv(p1,p2,p3) - dL1_dv()*dTy_du(p1,p2,p3);
+ return val;
+}
+
+double lloydAlgorithm::dJ_dp2x(SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double val;
+ val = dL2_du()*dTy_dv(p1,p2,p3) - dL2_dv()*dTy_du(p1,p2,p3);
+ return val;
+}
+
+double lloydAlgorithm::dJ_dp3x(SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double val;
+ val = dL3_du()*dTy_dv(p1,p2,p3) - dL3_dv()*dTy_du(p1,p2,p3);
+ return val;
+}
+
+double lloydAlgorithm::dJ_dp1y(SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double val;
+ val = dTx_du(p1,p2,p3)*dL1_dv() - dTx_dv(p1,p2,p3)*dL1_du();
+ return val;
+}
+
+double lloydAlgorithm::dJ_dp2y(SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double val;
+ val = dTx_du(p1,p2,p3)*dL2_dv() - dTx_dv(p1,p2,p3)*dL2_du();
+ return val;
+}
+
+double lloydAlgorithm::dJ_dp3y(SPoint2 p1,SPoint2 p2,SPoint2 p3){
+ double val;
+ val = dTx_du(p1,p2,p3)*dL3_dv() - dTx_dv(p1,p2,p3)*dL3_du();
+ return val;
+}
+
+SVector3 lloydAlgorithm::inner_dFdx0(SVector3 dFdC,SPoint2 C,SPoint2 x0,SPoint2 x1,SPoint2 x2){
+ fullMatrix<double> A(2,2);
+ fullMatrix<double> B(2,2);
+ fullMatrix<double> M(2,2);
+ fullMatrix<double> _dFdC(1,2);
+ fullMatrix<double> _val(1,2);
+ A(0,0) = x1.x() - x0.x();
+ A(0,1) = x1.y() - x0.y();
+ A(1,0) = x2.x() - x0.x();
+ A(1,1) = x2.y() - x0.y();
+ A.invertInPlace();
+ B(0,0) = C.x() - x0.x();
+ B(0,1) = C.y() - x0.y();
+ B(1,0) = C.x() - x0.x();
+ B(1,1) = C.y() - x0.y();
+ A.mult_naive(B,M);
+ _dFdC(0,0) = dFdC.x();
+ _dFdC(0,1) = dFdC.y();
+ _dFdC.mult_naive(M,_val);
+ return SVector3(_val(0,0),_val(0,1),0.0);
+}
+
+SVector3 lloydAlgorithm::boundary_dFdx0(SVector3 dFdC,SPoint2 C,SPoint2 x0,SPoint2 x1,SVector3 normal){
+ fullMatrix<double> A(2,2);
+ fullMatrix<double> B(2,2);
+ fullMatrix<double> M(2,2);
+ fullMatrix<double> _dFdC(1,2);
+ fullMatrix<double> _val(1,2);
+ A(0,0) = x1.x() - x0.x();
+ A(0,1) = x1.y() - x0.y();
+ A(1,0) = normal.x();
+ A(1,1) = normal.y();
+ A.invertInPlace();
+ B(0,0) = C.x() - x0.x();
+ B(0,1) = C.y() - x0.y();
+ B(1,0) = 0.0;
+ B(1,1) = 0.0;
+ A.mult_naive(B,M);
+ _dFdC(0,0) = dFdC.x();
+ _dFdC(0,1) = dFdC.y();
+ _dFdC.mult_naive(M,_val);
+ return SVector3(_val(0,0),_val(0,1),0.0);
+}
+
+double lloydAlgorithm::exp(double a,int b)
+{
+ int i;
+ double val = 1.0;
+ for(i=0;i<b;i++){
+ val = val*a;
+ }
+ return val;
+}
+
+
+
+/****************Calcul des intersections****************/
+/********************************************************/
+
+SVector3 lloydAlgorithm::normal(SPoint2 p1,SPoint2 p2){
+ double x;
+ double y;
+ SVector3 val;
+ x = p2.x() - p1.x();
+ y = p2.y() - p1.y();
+ val = SVector3(-y,x,0.0);
+ return val;
+}
+
+SPoint2 lloydAlgorithm::mid(SPoint2 p1,SPoint2 p2){
+ double x;
+ double y;
+ x = 0.5*(p1.x() + p2.x());
+ y = 0.5*(p1.y() + p2.y());
+ return SPoint2(x,y);
+}
+
+void lloydAlgorithm::print_segment(SPoint2 p1,SPoint2 p2,std::ofstream& file){
+ file << "SL (" << p1.x() << ", "
+ << p1.y() << ", 0, "
+ << p2.x() << ", "
+ << p2.y() << ", 0){" << "10, 20};\n";
+}
+
+bool lloydAlgorithm::adjacent(const DocRecord& triangulator,int index1,int index2){
+ int num = triangulator._adjacencies[index1].t_length;
+ for(int i=0;i<num;i++){
+ if(triangulator._adjacencies[index1].t[i] == index2){
+ return 1;
+ }
+ }
+ return 0;
+}
+
+double lloydAlgorithm::triangle_area(SPoint2 p1,SPoint2 p2,SPoint2 p3){
double area;
- num = polygon.size();
- if(num<3) return 0.0;
+ double x1,y1;
+ double x2,y2;
+ x1 = p2.x() - p1.x();
+ y1 = p2.y() - p1.y();
+ x2 = p3.x() - p1.x();
+ y2 = p3.y() - p1.y();
+ area = 0.5*fabs(x1*y2 - x2*y1);
+ return area;
+}
+
+SPoint2 lloydAlgorithm::boundary_intersection(SPoint2 p1,SPoint2 p2,bool& flag,SVector3& vec){
+ int i;
+ int num;
+ bool flag2;
+ SPoint2 p3;
+ SPoint2 p4;
+ SPoint2 val;
+ boundary_edge edge;
+ num = get_number_boundary_edges();
+ for(i=0;i<num;i++){
+ edge = get_boundary_edge(i);
+ p3 = edge.get_p1();
+ p4 = edge.get_p2();
+ val = intersection(p1,p2,p3,p4,flag2);
+ if(flag2){
+ flag = 1;
+ vec = normal(p3,p4);
+ return val;
+ }
+ }
+ flag = 0;
+ vec = SVector3(0.0,0.0,0.0);
+ return SPoint2(0.0,0.0);
+}
+
+bool lloydAlgorithm::inside_domain(double x,double y){
+ if(x<=1.0 && x>=0.0 && y<=1.0 && y>=0.0){
+ return 1;
+ }
else{
- area = 0.0;
- for(i=0;i<num-1;i++){
- area = area + 0.5*(polygon[i].x()*polygon[i+1].y() - polygon[i+1].x()*polygon[i].y());
- }
- area = area + 0.5*(polygon[num-1].x()*polygon[0].y() - polygon[0].x()*polygon[num-1].y());
- return fabs(area);
+ return 0;
}
}
-//This function determines if there is an intersection between the segment p1-p2 and the
-//infinite line passing through the generator, which is the point to which the Voronoi cell
-//belongs.
-//If there is an intersection, then the variable intersection will contain 1, otherwise it will
-//contain 0. In that case, the point that the function returns is not imporant.
-//If direction = 1, the variable generator will contain the x component, in that case,
-//the infinite line will be vertical.
-//If direction = 0, the variable generator will contain the y component, in that case,
-//the infinite line will be horizontal.
-SPoint2 lloydAlgorithm::line_intersection(SPoint2 p1,SPoint2 p2,double generator,bool direction,bool& intersection){
- double x1,x2,solution_x;
- double y1,y2,solution_y;
- double u,e;
+int lloydAlgorithm::get_number_boundary_edges(){
+ return 4;
+}
+
+boundary_edge lloydAlgorithm::get_boundary_edge(int i){
+ if(i==0) return boundary_edge(SPoint2(0.0,0.0),SPoint2(1.0,0.0));
+ else if(i==1) return boundary_edge(SPoint2(1.0,0.0),SPoint2(1.0,1.0));
+ else if(i==2) return boundary_edge(SPoint2(1.0,1.0),SPoint2(0.0,1.0));
+ else if(i==3) return boundary_edge(SPoint2(0.0,1.0),SPoint2(0.0,0.0));
+}
+
+SPoint2 lloydAlgorithm::intersection(SPoint2 p1,SPoint2 p2,SPoint2 p3,SPoint2 p4,bool& flag){
+ double x1,y1;
+ double x2,y2;
+ double x3,y3;
+ double x4,y4;
+ double ua,ub;
+ double num_ua;
+ double num_ub;
+ double denom;
+ double e;
x1 = p1.x();
y1 = p1.y();
x2 = p2.x();
y2 = p2.y();
- e = 0.0000001;
- intersection = 0;
- if(direction){
- if(fabs(x1 - generator)<e){
- intersection = 1;
- return SPoint2(x1,y1);
- }
- else if(fabs(x2 - generator)<e){
- intersection = 1;
- return SPoint2(x2,y2);
- }
- else if((x1<generator && x2>generator)||(x1>generator && x2<generator)){
- intersection = 1;
- u = (generator-x1)/(x2-x1);
- solution_x = generator;
- solution_y = y1 + u*(y2-y1);
- return SPoint2(solution_x,solution_y);
- }
+ x3 = p3.x();
+ y3 = p3.y();
+ x4 = p4.x();
+ y4 = p4.y();
+ e = 0.000001;
+ denom = (y4-y3)*(x2-x1) - (x4-x3)*(y2-y1);
+ if(fabs(denom)<e){
+ flag = 0;
+ return SPoint2(0.0,0.0);
+ }
+ num_ua = (x4-x3)*(y1-y3) - (y4-y3)*(x1-x3);
+ num_ub = (x2-x1)*(y1-y3) - (y2-y1)*(x1-x3);
+ ua = num_ua/denom;
+ ub = num_ub/denom;
+ if(ua<=1.0 && ua>=0.0 && ub<=1.0 && ub>=0.0){
+ flag = 1;
+ return SPoint2(x1+ua*(x2-x1),y1+ua*(y2-y1));
}
else{
- if(fabs(y1 - generator)<e){
- intersection = 1;
- return SPoint2(x1,y1);
- }
- else if(fabs(y2 - generator)<e){
- intersection = 1;
- return SPoint2(x2,y2);
- }
- else if((y1<generator && y2>generator)||(y1>generator && y2<generator)){
- intersection = 1;
- u = (generator-y1)/(y2-y1);
- solution_x = x1 + u*(x2-x1);
- solution_y = generator;
- return SPoint2(solution_x,solution_y);
- }
+ flag = 0;
+ return SPoint2(0.0,0.0);
}
- return SPoint2(0.0,0.0);
}
-//This is the function called by the script.
-double lloydAlgorithm::optimize(int max,int flag){
- GFace*face;
- GModel*model = GModel::current();
- GModel::fiter iterator;
- lloydAlgorithm lloyd(max,flag);
- for(iterator = model->firstFace();iterator != model->lastFace();iterator++)
- {
- face = *iterator;
- lloyd(face);
- recombineIntoQuads(face,1,1);
- }
- return lloyd.test;
+/****************Class boundary_edge****************/
+boundary_edge::boundary_edge(SPoint2 new_p1,SPoint2 new_p2){
+ p1 = new_p1;
+ p2 = new_p2;
+}
+
+boundary_edge::boundary_edge(){}
+
+boundary_edge::~boundary_edge(){}
+
+SPoint2 boundary_edge::get_p1(){
+ return p1;
}
-#include "Bindings.h"
-void lloydAlgorithm::registerBindings(binding *b){
- classBinding*cb = b->addClass<lloydAlgorithm>("lloydAlgorithm");
- cb->setDescription("This class is used to move the points according to lloyd algorithm.");
- methodBinding*cm;
- cm = cb->setConstructor<lloydAlgorithm>();
- cm->setDescription("This is the constructor.");
- cm = cb->addMethod("optimize",&lloydAlgorithm::optimize);
- cm->setDescription("This function moves the points.");
- cm->setArgNames("max","flag",NULL);
+SPoint2 boundary_edge::get_p2(){
+ return p2;
}
+
+
diff --git a/Mesh/meshGFaceLloyd.h b/Mesh/meshGFaceLloyd.h
index b86fd54b2a..1f2462771f 100644
--- a/Mesh/meshGFaceLloyd.h
+++ b/Mesh/meshGFaceLloyd.h
@@ -7,9 +7,17 @@
#define _MESH_GFACE_LLOYD_H_
#include "fullMatrix.h"
#include "DivideAndConquer.h"
-#include "Bindings.h"
+#include "ap.h"
+#include "alglibinternal.h"
+#include "alglibmisc.h"
+#include "linalg.h"
+#include "optimization.h"
class GFace;
+class boundary_edge;
+
+void function1_grad(const alglib::real_1d_array&,double&,alglib::real_1d_array&,void*);
+void topology(const alglib::real_1d_array &,double,void*);
class lloydAlgorithm {
int ITER_MAX;
@@ -20,19 +28,77 @@ class lloydAlgorithm {
: ITER_MAX(itermax), infiniteNorm(infnorm) {}
lloydAlgorithm() {}
void operator () (GFace *);
- double get_angle(double,double);
- void get_rotation(double,double,double&,double&,double&,double&);
- void get_unrotation(double,double,double&,double&,double&,double&);
- void rotate(SPoint2,SPoint2&);
- void unrotate(SPoint2,SPoint2&);
- SPoint2 fixed_point(const std::vector<SPoint2>&);
- double range(const std::vector<SPoint2>&,double,double,bool);
- double gradient(const std::vector<SPoint2>&,double,bool);
- bool point_side(SPoint2,double,bool);
- double polygon_area(const std::vector<SPoint2>&);
- SPoint2 line_intersection(SPoint2,SPoint2,double,bool,bool&);
double optimize(int,int);
- static void registerBindings(binding *b);
+
+ void eval(DocRecord&,GFace*,std::vector<SVector3>&,std::vector<double>&,std::vector<SVector3>&,std::vector<double>&,int);
+ double total_energy(std::vector<double>);
+ SVector3 numerical(DocRecord&,GFace*,int,int);
+ void write(DocRecord&,GFace*,int);
+ void write2(DocRecord&,GFace*);
+ SVector3 analytical(SPoint2,SVector3,double);
+ SVector3 area_centroid(SPoint2,SPoint2,SPoint2);
+
+ double F(SPoint2,SPoint2,SPoint2,int);
+ SVector3 simple(SPoint2,SPoint2,SPoint2,int);
+ SVector3 dF_dC1(SPoint2,SPoint2,SPoint2,int);
+ SVector3 dF_dC2(SPoint2,SPoint2,SPoint2,int);
+ double f(SPoint2,SPoint2,int);
+ double df_dx(SPoint2,SPoint2,int);
+ double df_dy(SPoint2,SPoint2,int);
+ double L1(double,double);
+ double dL1_du();
+ double dL1_dv();
+ double L2(double,double);
+ double dL2_du();
+ double dL2_dv();
+ double L3(double,double);
+ double dL3_du();
+ double dL3_dv();
+ double Tx(double,double,SPoint2,SPoint2,SPoint2);
+ double dTx_dp1x(double,double);
+ double dTx_dp2x(double,double);
+ double dTx_dp3x(double,double);
+ double dTx_du(SPoint2,SPoint2,SPoint2);
+ double dTx_dv(SPoint2,SPoint2,SPoint2);
+ double Ty(double,double,SPoint2,SPoint2,SPoint2);
+ double dTy_dp1y(double,double);
+ double dTy_dp2y(double,double);
+ double dTy_dp3y(double,double);
+ double dTy_du(SPoint2,SPoint2,SPoint2);
+ double dTy_dv(SPoint2,SPoint2,SPoint2);
+ double J(SPoint2,SPoint2,SPoint2);
+ double dJ_dp1x(SPoint2,SPoint2,SPoint2);
+ double dJ_dp2x(SPoint2,SPoint2,SPoint2);
+ double dJ_dp3x(SPoint2,SPoint2,SPoint2);
+ double dJ_dp1y(SPoint2,SPoint2,SPoint2);
+ double dJ_dp2y(SPoint2,SPoint2,SPoint2);
+ double dJ_dp3y(SPoint2,SPoint2,SPoint2);
+ SVector3 inner_dFdx0(SVector3,SPoint2,SPoint2,SPoint2,SPoint2);
+ SVector3 boundary_dFdx0(SVector3,SPoint2,SPoint2,SPoint2,SVector3);
+ double exp(double,int);
+
+ SVector3 normal(SPoint2,SPoint2);
+ SPoint2 mid(SPoint2,SPoint2);
+ void print_segment(SPoint2,SPoint2,std::ofstream&);
+ bool adjacent(const DocRecord&,int,int);
+ double triangle_area(SPoint2,SPoint2,SPoint2);
+ SPoint2 boundary_intersection(SPoint2,SPoint2,bool&,SVector3&);
+ bool inside_domain(double,double);
+ int get_number_boundary_edges();
+ boundary_edge get_boundary_edge(int);
+ SPoint2 intersection(SPoint2,SPoint2,SPoint2,SPoint2,bool&);
+};
+
+class boundary_edge{
+ private :
+ SPoint2 p1;
+ SPoint2 p2;
+ public :
+ boundary_edge(SPoint2,SPoint2);
+ boundary_edge();
+ ~boundary_edge();
+ SPoint2 get_p1();
+ SPoint2 get_p2();
};
#endif
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 20794238d7..07956971eb 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -19,6 +19,10 @@
#include "Generator.h"
#include "Context.h"
#include "OS.h"
+#include "PView.h"
+#include "PViewData.h"
+#include "SVector3.h"
+#include "SPoint3.h"
#ifdef HAVE_MATCH
extern "C" int FAILED_NODE;
@@ -1512,6 +1516,9 @@ struct RecombineTriangle
{
MElement *t1, *t2;
double angle;
+ double cost_quality; //addition for class Temporary
+ double cost_alignment; //addition for class Temporary
+ double total_cost; //addition for class Temporary
MVertex *n1, *n2, *n3, *n4;
RecombineTriangle(const MEdge &me, MElement *_t1, MElement *_t2)
: t1(_t1), t2(_t2)
@@ -1533,11 +1540,16 @@ struct RecombineTriangle
angle = fabs(90. - a1);
angle = std::max(fabs(90. - a2),angle);
angle = std::max(fabs(90. - a3),angle);
- angle = std::max(fabs(90. - a4),angle);
+ angle = std::max(fabs(90. - a4),angle);
+ cost_quality = 1.0 - std::max(1.0 - angle/90.0,0.0); //addition for class Temporary
+ cost_alignment = Temporary::compute_alignment(me,_t1,_t2); //addition for class Temporary
+ total_cost = Temporary::compute_total_cost(cost_quality,cost_alignment); //addition for class Temporary
+ total_cost = 100.0*cost_quality; //addition for class Temporary
}
bool operator < (const RecombineTriangle &other) const
{
- return angle < other.angle;
+ //return angle < other.angle;
+ return total_cost < other.total_cost; //addition for class Temporary
}
};
@@ -1640,7 +1652,8 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
for (int i=0;i<pairs.size();++i){
elist[2*i] = t2n[pairs[i].t1];
elist[2*i+1] = t2n[pairs[i].t2];
- elen [i] = (int) pairs[i].angle;
+ //elen [i] = (int) pairs[i].angle;
+ elen [i] = (int) pairs[i].total_cost; //addition for class Temporary
double angle = atan2(pairs[i].n1->y()-pairs[i].n2->y(),
pairs[i].n1->x()-pairs[i].n2->x());
@@ -2055,3 +2068,166 @@ void recombineIntoQuadsIterative(GFace *gf)
if (COUNT == 5)break;
}
}
+
+
+
+#include "Bindings.h"
+
+double Temporary::w1,Temporary::w2,Temporary::w3;
+std::vector<SVector3> Temporary::gradients;
+
+Temporary::Temporary(){}
+
+Temporary::~Temporary(){}
+
+SVector3 Temporary::compute_gradient(MElement*element){
+ double x1,y1,z1;
+ double x2,y2,z2;
+ double x3,y3,z3;
+ double x,y,z;
+ MVertex*vertex1 = element->getVertex(0);
+ MVertex*vertex2 = element->getVertex(1);
+ MVertex*vertex3 = element->getVertex(2);
+ x1 = vertex1->x();
+ y1 = vertex1->y();
+ z1 = vertex1->z();
+ x2 = vertex2->x();
+ y2 = vertex2->y();
+ z2 = vertex2->z();
+ x3 = vertex3->x();
+ y3 = vertex3->y();
+ z3 = vertex3->z();
+ x = (x1+x2+x3)/3.0;
+ y = (y1+y2+y3)/3.0;
+ z = (z1+z2+z3)/3.0;
+ return SVector3(0.0,1.0,0.0);
+}
+
+void Temporary::select_weights(double new_w1,double new_w2,double new_w3){
+ w1 = new_w1;
+ w2 = new_w2;
+ w3 = new_w3;
+}
+
+double Temporary::compute_total_cost(double f1,double f2){
+ double cost;
+ cost = w1*f1 + w2*f2 + w3*f1*f2;
+ return cost;
+}
+
+SVector3 Temporary::compute_normal(MElement*element){
+ double x1,y1,z1;
+ double x2,y2,z2;
+ double x3,y3,z3;
+ double length;
+ SVector3 vectorA;
+ SVector3 vectorB;
+ SVector3 normal;
+ MVertex*vertex1 = element->getVertex(0);
+ MVertex*vertex2 = element->getVertex(1);
+ MVertex*vertex3 = element->getVertex(2);
+ x1 = vertex1->x();
+ y1 = vertex1->y();
+ z1 = vertex1->z();
+ x2 = vertex2->x();
+ y2 = vertex2->y();
+ z2 = vertex2->z();
+ x3 = vertex3->x();
+ y3 = vertex3->y();
+ z3 = vertex3->z();
+ vectorA = SVector3(x2-x1,y2-y1,z2-z1);
+ vectorB = SVector3(x3-x1,y3-y1,z3-z1);
+ normal = crossprod(vectorA,vectorB);
+ length = norm(normal);
+ return SVector3(normal.x()/length,normal.y()/length,normal.z()/length);
+}
+
+SVector3 Temporary::compute_other_vector(MElement*element){
+ int number;
+ double length;
+ SVector3 normal;
+ SVector3 gradient;
+ SVector3 other_vector;
+ number = element->getNum();
+ normal = Temporary::compute_normal(element);
+ gradient = Temporary::compute_gradient(element);//gradients[number];
+ other_vector = crossprod(gradient,normal);
+ length = norm(other_vector);
+ return SVector3(other_vector.x()/length,other_vector.y()/length,other_vector.z()/length);
+}
+
+double Temporary::compute_alignment(const MEdge&_edge, MElement*element1, MElement*element2){
+ int number;
+ double scalar_productA,scalar_productB;
+ double alignment;
+ SVector3 gradient;
+ SVector3 other_vector;
+ SVector3 edge;
+ MVertex*vertexA;
+ MVertex*vertexB;
+ number = element1->getNum();
+ gradient = Temporary::compute_gradient(element1);//gradients[number];
+ other_vector = Temporary::compute_other_vector(element1);
+ vertexA = _edge.getVertex(0);
+ vertexB = _edge.getVertex(1);
+ edge = SVector3(vertexB->x()-vertexA->x(),vertexB->y()-vertexA->y(),vertexB->z()-vertexA->z());
+ scalar_productA = fabs(dot(gradient,edge));
+ scalar_productB = fabs(dot(other_vector,edge));
+ alignment = std::max(scalar_productA,scalar_productB) - sqrt(2.0)/2.0;
+ alignment = alignment/(1.0-sqrt(2.0)/2.0);
+ return alignment;
+}
+
+void Temporary::read_data(std::string file_name){
+ int i,j,number;
+ double x,y,z;
+ MElement*element;
+ PView*view;
+ PViewData*data;
+ PView::readMSH(file_name,-1);
+ std::vector<PView*> list = PView::list;
+ data = list[0]->getData();
+ for(i = 0;i < data->getNumEntities(0);i++)
+ {
+ if(data->skipEntity(0,i)) continue;
+ for(j = 0;j < data->getNumElements(0,i);j++)
+ {
+ if(data->skipElement(0,i,j)) continue;
+ element = data->getElement(0,i,j);
+ number = element->getNum();
+ data->getValue(0,i,j,0,x);
+ data->getValue(0,i,j,1,y);
+ data->getValue(0,i,j,2,z);
+ gradients[number] = SVector3(x,y,z);
+ }
+ }
+}
+
+void Temporary::quadrilaterize(std::string file_name,double _w1,double _w2,double _w3){
+ GFace*face;
+ GModel*model = GModel::current();
+ GModel::fiter iterator;
+ gradients.resize(model->getNumMeshElements() + 1);
+ w1 = _w1;
+ w2 = _w2;
+ w3 = _w3;
+ Temporary::read_data(file_name);
+ for(iterator = model->firstFace();iterator != model->lastFace();iterator++)
+ {
+ face = *iterator;
+ _recombineIntoQuads(face,1,1);
+ }
+}
+
+void Temporary::registerBindings(binding *b){
+ classBinding*cb = b->addClass<Temporary>("Temporary");
+ cb->setDescription("This class is used to create quad meshes from a script.");
+ methodBinding*cm;
+ cm = cb->setConstructor<Temporary>();
+ cm->setDescription("This is the constructor.");
+ cm = cb->addMethod("quadrilaterize",&Temporary::quadrilaterize);
+ cm->setDescription("This function creates the quad mesh.");
+ cm->setArgNames("file_name","w1","w2","w3",NULL);
+}
+
+
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index a94db547f4..d4f01d47be 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -126,4 +126,22 @@ struct swapquad{
}
};
+class Temporary{
+ private :
+ static double w1,w2,w3;
+ static std::vector<SVector3> gradients;
+ void read_data(std::string);
+ static SVector3 compute_normal(MElement*);
+ static SVector3 compute_other_vector(MElement*);
+ static SVector3 compute_gradient(MElement*);
+ public :
+ Temporary();
+ ~Temporary();
+ void quadrilaterize(std::string,double,double,double);
+ static double compute_total_cost(double,double);
+ static void select_weights(double,double,double);
+ static double compute_alignment(const MEdge&,MElement*,MElement*);
+ static void registerBindings(binding *b);
+};
+
#endif
diff --git a/Numeric/DivideAndConquer.cpp b/Numeric/DivideAndConquer.cpp
index 978e4937e2..f0c29562c8 100644
--- a/Numeric/DivideAndConquer.cpp
+++ b/Numeric/DivideAndConquer.cpp
@@ -28,6 +28,30 @@
#define Pred(x) ((x)->prev)
#define Succ(x) ((x)->next)
+int DocRecord::get_p(){
+ return p;
+}
+
+void DocRecord::set_p(int new_p){
+ p = new_p;
+}
+
+int DocRecord::get_dimension(){
+ return dimension;
+}
+
+void DocRecord::set_dimension(int new_dimension){
+ dimension = new_dimension;
+}
+
+GFace* DocRecord::get_face(){
+ return gf;
+}
+
+void DocRecord::set_face(GFace* new_gf){
+ gf = new_gf;
+}
+
PointNumero DocRecord::Predecessor(PointNumero a, PointNumero b)
{
DListPeek p = points[a].adjacent;
@@ -824,7 +848,7 @@ DocRecord::DocRecord(int n)
numPoints(n), points(NULL), numTriangles(0), triangles(NULL)
{
if(numPoints)
- points = new PointRecord[numPoints+1000];
+ points = new PointRecord[numPoints+3000];
}
DocRecord::~DocRecord()
@@ -898,6 +922,7 @@ void DocRecord::remove_all(){
points2[index].where.v = points[i].where.v;
points2[index].data = points[i].data;
points2[index].flag = points[i].flag;
+ points2[index].identificator = points[i].identificator;
index++;
}
}
@@ -906,6 +931,15 @@ void DocRecord::remove_all(){
numPoints = numPoints2;
}
+void DocRecord::add_point(double x,double y,GFace*face){
+ PointRecord point;
+ point.where.h = x;
+ point.where.v = y;
+ point.data = new MVertex(x,y,0.0,(GEntity*)face,2);
+ points[numPoints] = point;
+ numPoints = numPoints+1;
+}
+
#include "Bindings.h"
void DocRecord::registerBindings(binding *b)
diff --git a/Numeric/DivideAndConquer.h b/Numeric/DivideAndConquer.h
index da5f4d7322..aa5f646cfd 100644
--- a/Numeric/DivideAndConquer.h
+++ b/Numeric/DivideAndConquer.h
@@ -30,6 +30,7 @@ struct PointRecord {
DListPeek adjacent;
void *data;
int flag; //0:to be kept, 1:to be removed
+ int identificator;
PointRecord() : adjacent(0), data (0) {}
};
@@ -59,6 +60,9 @@ typedef struct{
class DocRecord{
private:
+ int p;
+ int dimension;
+ GFace* gf;
int _hullSize;
PointNumero *_hull;
PointNumero Predecessor(PointNumero a, PointNumero b);
@@ -76,15 +80,22 @@ class DocRecord{
int Insert(PointNumero a, PointNumero b);
int DListDelete(DListPeek *dlist, PointNumero oldPoint);
int Delete(PointNumero a, PointNumero b);
- int ConvertDListToTriangles();
void ConvertDListToVoronoiData();
+ int ConvertDListToTriangles();
void RemoveAllDList();
int BuildDelaunay();
int CountPointsOnHull();
void ConvexHull();
public:
+ int get_p();
+ void set_p(int);
+ int get_dimension();
+ void set_dimension(int);
+ GFace* get_face();
+ void set_face(GFace*);
STriangle *_adjacencies;
int numPoints; // number of points
+ int size_points;
PointRecord *points; // points to triangulate
int numTriangles; // number of triangles
Triangle *triangles; // 2D results
@@ -105,6 +116,7 @@ class DocRecord{
void initialize();
bool remove_point(int);
void remove_all();
+ void add_point(double,double,GFace*);
PointNumero *ConvertDlistToArray(DListPeek *dlist, int *n);
static void registerBindings(binding *b);
};
--
GitLab