diff --git a/CMakeLists.txt b/CMakeLists.txt
index cba64a8d9bf94813959ba1915fb7972116be5584..da4c3923be19be2a05909bbb58a33596bb1f1d79 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -494,7 +494,6 @@ if(ENABLE_BFGS)
endif(ENABLE_BFGS)
if(ENABLE_RTREE)
- add_subdirectory(contrib/rtree)
include_directories(contrib/rtree)
set_config_option(HAVE_RTREE "RTree")
endif(ENABLE_RTREE)
diff --git a/Mesh/Levy3D.cpp b/Mesh/Levy3D.cpp
index c6663545c77ad879cccf842a71e43b051304cc49..b4734f44523adbbd6bded2cebfffc3c4a4fc3995 100755
--- a/Mesh/Levy3D.cpp
+++ b/Mesh/Levy3D.cpp
@@ -10,6 +10,12 @@
#if defined(HAVE_BFGS)
+#include "ap.h"
+#include "alglibinternal.h"
+#include "alglibmisc.h"
+#include "linalg.h"
+#include "optimization.h"
+
#include "Levy3D.h"
#include "polynomialBasis.h"
#include "GModel.h"
@@ -121,7 +127,7 @@ void call_back(const alglib::real_1d_array& x,double& func,alglib::real_1d_array
std::vector<SPoint3> bank;
std::vector<int> movability;
std::vector<SVector3> gradients;
-
+
w = static_cast<Wrap*>(ptr);
p = w->get_p();
dimension = w->get_dimension();
@@ -133,14 +139,14 @@ void call_back(const alglib::real_1d_array& x,double& func,alglib::real_1d_array
octree = w->get_octree();
error1 = 0;
error2 = 0;
-
+
bank.resize(size);
movability.resize(size);
for(i=0;i<size;i++){
bank[i] = w->get_bank(i);
movability[i] = w->get_movability(i);
}
-
+
for(i=0;i<dimension/3;i++){
bank[i+offset] = SPoint3(x[i],x[i+(dimension/3)],x[i+(2*dimension/3)]);
flag = inside_domain(octree,x[i],x[i+(dimension/3)],x[i+(2*dimension/3)]);
@@ -149,12 +155,12 @@ void call_back(const alglib::real_1d_array& x,double& func,alglib::real_1d_array
printf("Vertices outside domain.\n");
}
}
-
+
if(iteration>max_iteration){
error2 = 1;
printf("Maximum number of iterations reached.\n");
}
-
+
if(!error1 && !error2){
gradients.resize(dimension/3);
obj.get_gauss();
@@ -172,7 +178,7 @@ void call_back(const alglib::real_1d_array& x,double& func,alglib::real_1d_array
grad[i] = 0.0;
}
}
-
+
if(initial_energy>0.0 && !error1 && !error2){
printf("%d %.9f\n",iteration,100.0*(initial_energy-energy)/initial_energy);
w->set_iteration(iteration+1);
@@ -190,7 +196,7 @@ VoronoiVertex::VoronoiVertex(){
index3 = -1;
index4 = -1;
normal1 = SVector3(0.0,0.0,0.0);
- normal2 = SVector3(0.0,0.0,0.0);
+ normal2 = SVector3(0.0,0.0,0.0);
}
VoronoiVertex::~VoronoiVertex(){}
@@ -278,7 +284,7 @@ Tensor::Tensor(){
t23 = 0.0;
t31 = 0.0;
t32 = 0.0;
- t33 = 1.0;
+ t33 = 1.0;
}
Tensor::~Tensor(){}
@@ -423,7 +429,7 @@ double VoronoiElement::get_h(double u,double v,double w){
double h3;
double h4;
double h;
-
+
h1 = v1.get_h();
h2 = v2.get_h();
h3 = v3.get_h();
@@ -464,7 +470,7 @@ void VoronoiElement::deriv_h(){
SPoint3 p2;
SPoint3 p3;
SPoint3 p4;
-
+
p1 = v1.get_point();
p2 = v2.get_point();
p3 = v3.get_point();
@@ -509,7 +515,7 @@ void VoronoiElement::deriv_h(){
du_dz = b31/jacobian2;
dv_dz = b32/jacobian2;
dw_dz = b33/jacobian2;
-
+
h1 = v1.get_h();
h2 = v2.get_h();
h3 = v3.get_h();
@@ -517,7 +523,7 @@ void VoronoiElement::deriv_h(){
dh_du = h2-h1;
dh_dv = h3-h1;
dh_dw = h4-h1;
-
+
dh_dx = dh_du*du_dx + dh_dv*dv_dx + dh_dw*dw_dx;
dh_dy = dh_du*du_dy + dh_dv*dv_dy + dh_dw*dw_dy;
dh_dz = dh_du*du_dz + dh_dv*dv_dz + dh_dw*dw_dz;
@@ -532,7 +538,7 @@ void VoronoiElement::compute_jacobian(){
SPoint3 p2;
SPoint3 p3;
SPoint3 p4;
-
+
p1 = v1.get_point();
p2 = v2.get_point();
p3 = v3.get_point();
@@ -576,22 +582,22 @@ double VoronoiElement::get_quality(){
double quality;
double min_l,max_l;
double l[6];
-
+
l[0] = v1.get_point().distance(v2.get_point());
l[1] = v1.get_point().distance(v3.get_point());
l[2] = v1.get_point().distance(v4.get_point());
l[3] = v2.get_point().distance(v3.get_point());
l[4] = v2.get_point().distance(v4.get_point());
l[5] = v3.get_point().distance(v4.get_point());
-
+
min_l = 1000000.0;
max_l = -1000000.0;
-
+
for(i=0;i<6;i++){
min_l = std::min(min_l,l[i]);
max_l = std::max(max_l,l[i]);
}
-
+
quality = min_l/max_l;
return quality;
}
@@ -700,36 +706,36 @@ void LpCVT::verification(std::vector<SPoint3>& bank,std::vector<int>& movability
double front,back;
double e;
std::vector<SVector3> gradients;
-
+
gradients.resize(bank.size()-offset);
e = 0.0001;
srand(time(NULL));
index = rand()%(bank.size()-offset) + offset;
-
+
bank[index] = SPoint3(bank[index].x()+e,bank[index].y(),bank[index].z());
eval(bank,movability,offset,gradients,right,p);
bank[index] = SPoint3(bank[index].x()-e,bank[index].y(),bank[index].z());
-
+
bank[index] = SPoint3(bank[index].x()-e,bank[index].y(),bank[index].z());
eval(bank,movability,offset,gradients,left,p);
bank[index] = SPoint3(bank[index].x()+e,bank[index].y(),bank[index].z());
-
+
bank[index] = SPoint3(bank[index].x(),bank[index].y()+e,bank[index].z());
eval(bank,movability,offset,gradients,up,p);
bank[index] = SPoint3(bank[index].x(),bank[index].y()-e,bank[index].z());
-
+
bank[index] = SPoint3(bank[index].x(),bank[index].y()-e,bank[index].z());
eval(bank,movability,offset,gradients,down,p);
bank[index] = SPoint3(bank[index].x(),bank[index].y()+e,bank[index].z());
-
+
bank[index] = SPoint3(bank[index].x(),bank[index].y(),bank[index].z()+e);
eval(bank,movability,offset,gradients,front,p);
bank[index] = SPoint3(bank[index].x(),bank[index].y(),bank[index].z()-e);
-
+
bank[index] = SPoint3(bank[index].x(),bank[index].y(),bank[index].z()-e);
eval(bank,movability,offset,gradients,back,p);
bank[index] = SPoint3(bank[index].x(),bank[index].y(),bank[index].z()+e);
-
+
eval(bank,movability,offset,gradients,energy,p);
printf("...%f %f %f\n",gradients[index-offset].x(),gradients[index-offset].y(),gradients[index-offset].z());
@@ -745,18 +751,18 @@ void LpCVT::eval(std::vector<SPoint3>& bank,std::vector<int>& movability,int off
double e;
SVector3 grad1,grad2,grad3;
clip approx;
-
+
for(i=0;i<gradients.size();i++){
gradients[i] = SVector3(0.0,0.0,0.0);
}
energy = 0.0;
e = 0.000001;
-
+
clipped.clear();
approx.execute(bank,clipped);
swap();
compute_parameters();
-
+
for(i=0;i<clipped.size();i++){
if(clipped[i].get_quality()<e) continue;
energy = energy + F(clipped[i],p);
@@ -823,12 +829,12 @@ void LpCVT::compute_parameters(){
double h1,h2,h3,h4;
Tensor t;
VoronoiVertex v1,v2,v3,v4;
-
+
for(i=0;i<clipped.size();i++){
v1 = clipped[i].get_v1();
v2 = clipped[i].get_v2();
v3 = clipped[i].get_v3();
- v4 = clipped[i].get_v4();
+ v4 = clipped[i].get_v4();
h1 = get_size(clipped[i].get_v1().get_point().x(),clipped[i].get_v1().get_point().y(),clipped[i].get_v1().get_point().z());
h2 = get_size(clipped[i].get_v2().get_point().x(),clipped[i].get_v2().get_point().y(),clipped[i].get_v2().get_point().z());
h3 = get_size(clipped[i].get_v3().get_point().x(),clipped[i].get_v3().get_point().y(),clipped[i].get_v3().get_point().z());
@@ -858,19 +864,19 @@ Tensor LpCVT::get_tensor(double x,double y,double z){
angle = atan2(z,x);
t = Tensor();
-
+
t.set_t11(1.0);
t.set_t12(0.0);
t.set_t13(0.0);
-
+
t.set_t21(0.0);
t.set_t22(1.0);
t.set_t23(0.0);
-
+
t.set_t31(0.0);
t.set_t32(0.0);
t.set_t33(1.0);
-
+
return t;
}
@@ -887,7 +893,7 @@ double LpCVT::get_drho_dx(double x,double y,double z,int p){
less1 = h_to_rho(get_size(x-e,y,z),p);
plus1 = h_to_rho(get_size(x+e,y,z),p);
plus2 = h_to_rho(get_size(x+2.0*e,y,z),p);
-
+
val = (less2 - 8.0*less1 + 8.0*plus1 - plus2)/(12.0*e);
return val;
}
@@ -899,13 +905,13 @@ double LpCVT::get_drho_dy(double x,double y,double z,int p){
double plus1;
double plus2;
double val;
-
+
e = 0.000001;
less2 = h_to_rho(get_size(x,y-2.0*e,z),p);
less1 = h_to_rho(get_size(x,y-e,z),p);
plus1 = h_to_rho(get_size(x,y+e,z),p);
plus2 = h_to_rho(get_size(x,y+2.0*e,z),p);
-
+
val = (less2 - 8.0*less1 + 8.0*plus1 - plus2)/(12.0*e);
return val;
}
@@ -917,15 +923,15 @@ double LpCVT::get_drho_dz(double x,double y,double z,int p){
double plus1;
double plus2;
double val;
-
+
e = 0.000001;
less2 = h_to_rho(get_size(x,y,z-2.0*e),p);
less1 = h_to_rho(get_size(x,y,z-e),p);
plus1 = h_to_rho(get_size(x,y,z+e),p);
plus2 = h_to_rho(get_size(x,y,z+2.0*e),p);
-
+
val = (less2 - 8.0*less1 + 8.0*plus1 - plus2)/(12.0*e);
- return val;
+ return val;
}
double LpCVT::h_to_rho(double h,int p){
@@ -935,9 +941,9 @@ double LpCVT::h_to_rho(double h,int p){
}
void LpCVT::swap(){
- int i;
+ int i;
for(i=0;i<clipped.size();i++){
- clipped[i].swap();
+ clipped[i].swap();
}
}
@@ -958,7 +964,7 @@ double LpCVT::F(VoronoiElement element,int p){
SPoint3 point,generator,C1,C2,C3;
VoronoiVertex v1,v2,v3,v4;
Tensor t;
-
+
v1 = element.get_v1();
v2 = element.get_v2();
v3 = element.get_v3();
@@ -969,7 +975,7 @@ double LpCVT::F(VoronoiElement element,int p){
C3 = v4.get_point();
energy = 0.0;
t = element.get_tensor();
-
+
for(i=0;i<gauss_num;i++){
u = gauss_points(i,0);
v = gauss_points(i,1);
@@ -997,7 +1003,7 @@ SVector3 LpCVT::simple(VoronoiElement element,int p){
SPoint3 point,generator,C1,C2,C3;
VoronoiVertex v1,v2,v3,v4;
Tensor t;
-
+
v1 = element.get_v1();
v2 = element.get_v2();
v3 = element.get_v3();
@@ -1011,7 +1017,7 @@ SVector3 LpCVT::simple(VoronoiElement element,int p){
comp_z = 0.0;
jacobian = element.get_jacobian();
t = element.get_tensor();
-
+
for(i=0;i<gauss_num;i++){
u = gauss_points(i,0);
v = gauss_points(i,1);
@@ -1046,7 +1052,7 @@ SVector3 LpCVT::dF_dC1(VoronoiElement element,int p){
SPoint3 point,generator,C1,C2,C3;
VoronoiVertex v1,v2,v3,v4;
Tensor t;
-
+
v1 = element.get_v1();
v2 = element.get_v2();
v3 = element.get_v3();
@@ -1063,7 +1069,7 @@ SVector3 LpCVT::dF_dC1(VoronoiElement element,int p){
gx = generator.x();
gy = generator.y();
gz = generator.z();
-
+
for(i=0;i<gauss_num;i++){
u = gauss_points(i,0);
v = gauss_points(i,1);
@@ -1087,7 +1093,7 @@ SVector3 LpCVT::dF_dC1(VoronoiElement element,int p){
comp_z = comp_z + weight*rho*df_dz(point,generator,t,p)*u*jacobian;
comp_z = comp_z + weight*rho*distance*((C2.x()-gx)*(C3.y()-gy) - (C3.x()-gx)*(C2.y()-gy));
comp_z = comp_z + weight*drho_dz*u*distance*jacobian;
- }
+ }
return SVector3(comp_x,comp_y,comp_z);
}
@@ -1105,7 +1111,7 @@ SVector3 LpCVT::dF_dC2(VoronoiElement element,int p){
SPoint3 point,generator,C1,C2,C3;
VoronoiVertex v1,v2,v3,v4;
Tensor t;
-
+
v1 = element.get_v1();
v2 = element.get_v2();
v3 = element.get_v3();
@@ -1122,7 +1128,7 @@ SVector3 LpCVT::dF_dC2(VoronoiElement element,int p){
gx = generator.x();
gy = generator.y();
gz = generator.z();
-
+
for(i=0;i<gauss_num;i++){
u = gauss_points(i,0);
v = gauss_points(i,1);
@@ -1146,7 +1152,7 @@ SVector3 LpCVT::dF_dC2(VoronoiElement element,int p){
comp_z = comp_z + weight*rho*df_dz(point,generator,t,p)*v*jacobian;
comp_z = comp_z + weight*rho*distance*((C3.x()-gx)*(C1.y()-gy) - (C1.x()-gx)*(C3.y()-gy));
comp_z = comp_z + weight*drho_dz*v*distance*jacobian;
- }
+ }
return SVector3(comp_x,comp_y,comp_z);
}
@@ -1164,7 +1170,7 @@ SVector3 LpCVT::dF_dC3(VoronoiElement element,int p){
SPoint3 point,generator,C1,C2,C3;
VoronoiVertex v1,v2,v3,v4;
Tensor t;
-
+
v1 = element.get_v1();
v2 = element.get_v2();
v3 = element.get_v3();
@@ -1181,7 +1187,7 @@ SVector3 LpCVT::dF_dC3(VoronoiElement element,int p){
gx = generator.x();
gy = generator.y();
gz = generator.z();
-
+
for(i=0;i<gauss_num;i++){
u = gauss_points(i,0);
v = gauss_points(i,1);
@@ -1205,7 +1211,7 @@ SVector3 LpCVT::dF_dC3(VoronoiElement element,int p){
comp_z = comp_z + weight*rho*df_dz(point,generator,t,p)*w*jacobian;
comp_z = comp_z + weight*rho*distance*((C1.x()-gx)*(C2.y()-gy) - (C2.x()-gx)*(C1.y()-gy));
comp_z = comp_z + weight*drho_dz*w*distance*jacobian;
- }
+ }
return SVector3(comp_x,comp_y,comp_z);
}
@@ -1217,7 +1223,7 @@ double LpCVT::f(SPoint3 p1,SPoint3 p2,Tensor t,int p){
double t31,t32,t33;
double val1,val2,val3;
double val;
-
+
x1 = p1.x();
y1 = p1.y();
z1 = p1.z();
@@ -1232,7 +1238,7 @@ double LpCVT::f(SPoint3 p1,SPoint3 p2,Tensor t,int p){
t23 = t.get_t23();
t31 = t.get_t31();
t32 = t.get_t32();
- t33 = t.get_t33();
+ t33 = t.get_t33();
val1 = t11*x1 + t12*y1 + t13*z1 - t11*x2 - t12*y2 - t13*z2;
val2 = t21*x1 + t22*y1 + t23*z1 - t21*x2 - t22*y2 - t23*z2;
val3 = t31*x1 + t32*y1 + t33*z1 - t31*x2 - t32*y2 - t33*z2;
@@ -1248,7 +1254,7 @@ double LpCVT::df_dx(SPoint3 p1,SPoint3 p2,Tensor t,int p){
double t31,t32,t33;
double val1,val2,val3;
double val;
-
+
x1 = p1.x();
y1 = p1.y();
z1 = p1.z();
@@ -1263,7 +1269,7 @@ double LpCVT::df_dx(SPoint3 p1,SPoint3 p2,Tensor t,int p){
t23 = t.get_t23();
t31 = t.get_t31();
t32 = t.get_t32();
- t33 = t.get_t33();
+ t33 = t.get_t33();
val1 = t11*x1 + t12*y1 + t13*z1 - t11*x2 - t12*y2 - t13*z2;
val2 = t21*x1 + t22*y1 + t23*z1 - t21*x2 - t22*y2 - t23*z2;
val3 = t31*x1 + t32*y1 + t33*z1 - t31*x2 - t32*y2 - t33*z2;
@@ -1279,7 +1285,7 @@ double LpCVT::df_dy(SPoint3 p1,SPoint3 p2,Tensor t,int p){
double t31,t32,t33;
double val1,val2,val3;
double val;
-
+
x1 = p1.x();
y1 = p1.y();
z1 = p1.z();
@@ -1294,7 +1300,7 @@ double LpCVT::df_dy(SPoint3 p1,SPoint3 p2,Tensor t,int p){
t23 = t.get_t23();
t31 = t.get_t31();
t32 = t.get_t32();
- t33 = t.get_t33();
+ t33 = t.get_t33();
val1 = t11*x1 + t12*y1 + t13*z1 - t11*x2 - t12*y2 - t13*z2;
val2 = t21*x1 + t22*y1 + t23*z1 - t21*x2 - t22*y2 - t23*z2;
val3 = t31*x1 + t32*y1 + t33*z1 - t31*x2 - t32*y2 - t33*z2;
@@ -1310,7 +1316,7 @@ double LpCVT::df_dz(SPoint3 p1,SPoint3 p2,Tensor t,int p){
double t31,t32,t33;
double val1,val2,val3;
double val;
-
+
x1 = p1.x();
y1 = p1.y();
z1 = p1.z();
@@ -1325,7 +1331,7 @@ double LpCVT::df_dz(SPoint3 p1,SPoint3 p2,Tensor t,int p){
t23 = t.get_t23();
t31 = t.get_t31();
t32 = t.get_t32();
- t33 = t.get_t33();
+ t33 = t.get_t33();
val1 = t11*x1 + t12*y1 + t13*z1 - t11*x2 - t12*y2 - t13*z2;
val2 = t21*x1 + t22*y1 + t23*z1 - t21*x2 - t22*y2 - t23*z2;
val3 = t31*x1 + t32*y1 + t33*z1 - t31*x2 - t32*y2 - t33*z2;
@@ -1339,13 +1345,13 @@ SVector3 LpCVT::bisectors3(SVector3 dIdC,SPoint3 C,SPoint3 x0,SPoint3 x1,SPoint3
fullMatrix<double> M(3,3);
fullMatrix<double> _dIdC(1,3);
fullMatrix<double> _val(1,3);
- A(0,0) = x1.x() - x0.x();
+ A(0,0) = x1.x() - x0.x();
A(0,1) = x1.y() - x0.y();
A(0,2) = x1.z() - x0.z();
- A(1,0) = x2.x() - x0.x();
+ A(1,0) = x2.x() - x0.x();
A(1,1) = x2.y() - x0.y();
A(1,2) = x2.z() - x0.z();
- A(2,0) = x3.x() - x0.x();
+ A(2,0) = x3.x() - x0.x();
A(2,1) = x3.y() - x0.y();
A(2,2) = x3.z() - x0.z();
A.invertInPlace();
@@ -1363,7 +1369,7 @@ SVector3 LpCVT::bisectors3(SVector3 dIdC,SPoint3 C,SPoint3 x0,SPoint3 x1,SPoint3
_dIdC(0,1) = dIdC.y();
_dIdC(0,2) = dIdC.z();
_dIdC.mult_naive(M,_val);
- return SVector3(_val(0,0),_val(0,1),_val(0,2));
+ return SVector3(_val(0,0),_val(0,1),_val(0,2));
}
SVector3 LpCVT::bisectors2(SVector3 dIdC,SPoint3 C,SPoint3 x0,SPoint3 x1,SPoint3 x2,SVector3 normal1){
@@ -1372,13 +1378,13 @@ SVector3 LpCVT::bisectors2(SVector3 dIdC,SPoint3 C,SPoint3 x0,SPoint3 x1,SPoint3
fullMatrix<double> M(3,3);
fullMatrix<double> _dIdC(1,3);
fullMatrix<double> _val(1,3);
- A(0,0) = x1.x() - x0.x();
+ A(0,0) = x1.x() - x0.x();
A(0,1) = x1.y() - x0.y();
A(0,2) = x1.z() - x0.z();
- A(1,0) = x2.x() - x0.x();
+ A(1,0) = x2.x() - x0.x();
A(1,1) = x2.y() - x0.y();
A(1,2) = x2.z() - x0.z();
- A(2,0) = normal1.x();
+ A(2,0) = normal1.x();
A(2,1) = normal1.y();
A(2,2) = normal1.z();
A.invertInPlace();
@@ -1396,7 +1402,7 @@ SVector3 LpCVT::bisectors2(SVector3 dIdC,SPoint3 C,SPoint3 x0,SPoint3 x1,SPoint3
_dIdC(0,1) = dIdC.y();
_dIdC(0,2) = dIdC.z();
_dIdC.mult_naive(M,_val);
- return SVector3(_val(0,0),_val(0,1),_val(0,2));
+ return SVector3(_val(0,0),_val(0,1),_val(0,2));
}
SVector3 LpCVT::bisectors1(SVector3 dIdC,SPoint3 C,SPoint3 x0,SPoint3 x1,SVector3 normal1,SVector3 normal2){
@@ -1405,13 +1411,13 @@ SVector3 LpCVT::bisectors1(SVector3 dIdC,SPoint3 C,SPoint3 x0,SPoint3 x1,SVector
fullMatrix<double> M(3,3);
fullMatrix<double> _dIdC(1,3);
fullMatrix<double> _val(1,3);
- A(0,0) = x1.x() - x0.x();
+ A(0,0) = x1.x() - x0.x();
A(0,1) = x1.y() - x0.y();
A(0,2) = x1.z() - x0.z();
- A(1,0) = normal1.x();
+ A(1,0) = normal1.x();
A(1,1) = normal1.y();
A(1,2) = normal1.z();
- A(2,0) = normal2.x();
+ A(2,0) = normal2.x();
A(2,1) = normal2.y();
A(2,2) = normal2.z();
A.invertInPlace();
@@ -1429,7 +1435,7 @@ SVector3 LpCVT::bisectors1(SVector3 dIdC,SPoint3 C,SPoint3 x0,SPoint3 x1,SVector
_dIdC(0,1) = dIdC.y();
_dIdC(0,2) = dIdC.z();
_dIdC.mult_naive(M,_val);
- return SVector3(_val(0,0),_val(0,1),_val(0,2));
+ return SVector3(_val(0,0),_val(0,1),_val(0,2));
}
void LpCVT::clear(){
@@ -1449,7 +1455,7 @@ void LpSmoother::improve_model(){
GRegion* gr;
GModel* model = GModel::current();
GModel::riter it;
-
+
for(it=model->firstRegion();it!=model->lastRegion();it++)
{
gr = *it;
@@ -1486,8 +1492,8 @@ void LpSmoother::improve_region(GRegion* gr){
alglib::real_1d_array x;
alglib::real_1d_array alglib_scales;
- octree = new MElementOctree(gr->model());
-
+ octree = new MElementOctree(gr->model());
+
for(i=0;i<gr->getNumMeshElements();i++){
element = gr->getMeshElement(i);
v1 = element->getVertex(0);
@@ -1500,7 +1506,7 @@ void LpSmoother::improve_region(GRegion* gr){
unmovable2.insert(v3);
unmovable2.insert(v4);
}
- }
+ }
for(i=0;i<gr->getNumMeshElements();i++){
element = gr->getMeshElement(i);
@@ -1520,7 +1526,7 @@ void LpSmoother::improve_region(GRegion* gr){
if(unmovable2.find(v4)==unmovable2.end()){
movable2.insert(v4);
}
- }
+ }
for(it=unmovable2.begin();it!=unmovable2.end();it++){
unmovable.push_back(*it);
@@ -1531,18 +1537,18 @@ void LpSmoother::improve_region(GRegion* gr){
}
offset = unmovable.size();
-
+
for(i=0;i<unmovable.size();i++){
point = SPoint3(unmovable[i]->x(),unmovable[i]->y(),unmovable[i]->z());
bank.push_back(point);
movability.push_back(0);
}
-
+
for(i=0;i<movable.size();i++){
point = SPoint3(movable[i]->x(),movable[i]->y(),movable[i]->z());
bank.push_back(point);
movability.push_back(1);
- }
+ }
w = Wrap();
w.set_p(norm);
@@ -1556,7 +1562,7 @@ void LpSmoother::improve_region(GRegion* gr){
for(i=0;i<bank.size();i++) w.set_bank(bank[i],i);
for(i=0;i<bank.size();i++) w.set_movability(movability[i],i);
- //if((bank.size()-offset)>1){
+ //if((bank.size()-offset)>1){
//LpCVT obj;
//obj.get_gauss();
//obj.verification(bank,movability,offset,6);
@@ -1566,7 +1572,7 @@ void LpSmoother::improve_region(GRegion* gr){
epsf = 0;
epsx = 0;
maxits = max_iter;
-
+
double initial_conditions[3*(bank.size()-offset)];
double scales[3*(bank.size()-offset)];
factor = 2.0;
@@ -1581,7 +1587,7 @@ void LpSmoother::improve_region(GRegion* gr){
x.setcontent(3*(bank.size()-offset),initial_conditions);
alglib_scales.setcontent(3*(bank.size()-offset),scales);
- if((bank.size()-offset)>1){
+ if((bank.size()-offset)>1){
minlbfgscreate(3*(bank.size()-offset),4,x,state);
minlbfgssetscale(state,alglib_scales);
minlbfgssetprecscale(state);
@@ -1600,15 +1606,15 @@ void LpSmoother::improve_region(GRegion* gr){
vertex = new MVertex(unmovable[i]->x(),unmovable[i]->y(),unmovable[i]->z(),gr,0);
interior_vertices.push_back(vertex);
}
- }
-
+ }
+
deleter(gr);
std::vector<GRegion*> regions;
regions.push_back(gr);
meshGRegion mesher(regions); //?
mesher(gr); //?
MeshDelaunayVolume(regions);
-
+
for(i=0;i<interior_vertices.size();i++) delete interior_vertices[i];
interior_vertices.clear();
delete octree;
@@ -1621,7 +1627,7 @@ double LpSmoother::get_size(double x,double y,double z){
int LpSmoother::get_nbr_interior_vertices(){
return interior_vertices.size();
}
-
+
MVertex* LpSmoother::get_interior_vertex(int i){
return interior_vertices[i];
}
diff --git a/Mesh/Levy3D.h b/Mesh/Levy3D.h
index 53d293607ef67b0aba28cff51d366274365a74fd..1baa6ac96c6cf18b1bf7331cde7f4699bf7eb408 100755
--- a/Mesh/Levy3D.h
+++ b/Mesh/Levy3D.h
@@ -6,16 +6,14 @@
// Contributor(s):
// Tristan Carrier
-#include "SVector3.h"
+#ifndef _LEVY3D_H_
+#define _LEVY3D_H_
+
#include <list>
+#include "SVector3.h"
#include "fullMatrix.h"
#include "GRegion.h"
#include "MElementOctree.h"
-#include "ap.h"
-#include "alglibinternal.h"
-#include "alglibmisc.h"
-#include "linalg.h"
-#include "optimization.h"
class VoronoiVertex{
private:
@@ -66,9 +64,9 @@ class Tensor{
void set_t31(double);
void set_t32(double);
void set_t33(double);
- double get_t11();
- double get_t12();
- double get_t13();
+ double get_t11();
+ double get_t12();
+ double get_t13();
double get_t21();
double get_t22();
double get_t23();
@@ -126,4 +124,6 @@ class LpSmoother{
double get_size(double,double,double);
static int get_nbr_interior_vertices();
static MVertex* get_interior_vertex(int);
-};
\ No newline at end of file
+};
+
+#endif
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 54089337372aa7f1019d1f12e174cb3a329ab14f..793f8fa0355bb720093efd443dc17f5eb28b885e 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -887,15 +887,17 @@ void bowyerWatson(GFace *gf)
vSizesBGM, vMetricsBGM, AllTris);
}
}
+#if defined(HAVE_ANN)
{
FieldManager *fields = gf->model()->getFields();
BoundaryLayerField *blf = 0;
if(fields->getBackgroundField() > 0){
Field *bl_field = fields->get(fields->getBoundaryLayerField());
blf = dynamic_cast<BoundaryLayerField*> (bl_field);
- if (blf && !blf->iRecombine)quadsToTriangles(gf,10000);
+ if (blf && !blf->iRecombine) quadsToTriangles(gf,10000);
}
}
+#endif
transferDataStructure(gf, AllTris, Us, Vs);
}
@@ -1206,6 +1208,7 @@ void bowyerWatsonFrontal(GFace *gf)
// _printTris (name, AllTris, Us, Vs,true);
transferDataStructure(gf, AllTris, Us, Vs);
// in case of boundary layer meshing
+#if defined(HAVE_ANN)
{
FieldManager *fields = gf->model()->getFields();
BoundaryLayerField *blf = 0;
@@ -1215,6 +1218,7 @@ void bowyerWatsonFrontal(GFace *gf)
if (blf && !blf->iRecombine)quadsToTriangles(gf,10000);
}
}
+#endif
}
void optimalPointFrontalQuad (GFace *gf,
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index a02a7709113e4878a54b0bdb26b5aecd30b032e2..8c70482d3376d60ec8828189fb9001c87a8ceca4 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -7,6 +7,7 @@
// Tristan Carrier
#include <set>
+#include <fstream>
#include "meshGFaceLloyd.h"
#include "DivideAndConquer.h"
#include "GFace.h"
@@ -15,8 +16,11 @@
#include "MTriangle.h"
#include "Context.h"
#include "meshGFace.h"
-#include "BackgroundMesh.h"
-#include <fstream>
+#include "BackgroundMesh.h"
+#include "GmshConfig.h"
+
+#if defined(HAVE_BFGS)
+
#include "ap.h"
#include "alglibinternal.h"
#include "alglibmisc.h"
@@ -37,13 +41,13 @@ class metric{
metric(double,double,double,double);
metric();
~metric();
- void set_a(double);
- void set_b(double);
- void set_c(double);
- void set_d(double);
- double get_a();
- double get_b();
- double get_c();
+ void set_a(double);
+ void set_b(double);
+ void set_c(double);
+ void set_d(double);
+ double get_a();
+ double get_b();
+ double get_c();
double get_d();
};
@@ -217,7 +221,7 @@ class lpcvt{
void print_voronoi2();
void print_delaunay(DocRecord&);
void print_segment(SPoint2,SPoint2,std::ofstream&);
-
+
void compute_parameters(GFace*,int);
double get_ratio(GFace*,SPoint2);
void write(DocRecord&,GFace*,int);
@@ -242,7 +246,7 @@ class lpcvt{
bool domain_search(MElementOctree* octree,double x,double y){
MElement* element;
-
+
element = (MElement*)octree->find(x,y,0.0,2,true);
if(element!=NULL) return 1;
else return 0;
@@ -273,7 +277,7 @@ void callback(const alglib::real_1d_array& x,double& func,alglib::real_1d_array&
MElementOctree* octree;
lpcvt obj;
std::vector<SVector3> gradients;
-
+
w = static_cast<wrapper*>(ptr);
p = w->get_p();
dimension = w->get_dimension();
@@ -288,7 +292,7 @@ void callback(const alglib::real_1d_array& x,double& func,alglib::real_1d_array&
error1 = 0;
error2 = 0;
error3 = 0;
-
+
index = 0;
for(i=0;i<num;i++){
if(obj.interior(*pointer,gf,i)){
@@ -324,7 +328,7 @@ void callback(const alglib::real_1d_array& x,double& func,alglib::real_1d_array&
func = energy;
}
}
-
+
if(error1 || error2 || error3){
energy = 1000000000.0;
for(i=0;i<num;i++){
@@ -338,11 +342,11 @@ void callback(const alglib::real_1d_array& x,double& func,alglib::real_1d_array&
}
if(error2){
printf("Maximum number of iterations reached.\n");
- }
+ }
if(error3){
printf("Boundary intersection.\n");
- }
-
+ }
+
if(start>0.0 && !error1 && !error2 && !error3){
printf("%d %.3f\n",iteration,100.0*(start-energy)/start);
w->set_iteration(iteration+1);
@@ -350,14 +354,14 @@ void callback(const alglib::real_1d_array& x,double& func,alglib::real_1d_array&
else if(!error1 && !error2 && !error3){
w->set_start(energy);
}
-
+
for(i=0;i<num;i++){
if(obj.interior(*pointer,gf,i)){
identificator = pointer->points[i].identificator;
grad[identificator] = gradients[i].x();
grad[identificator + dimension/2] = gradients[i].y();
}
- }
+ }
}
/****************class smoothing****************/
@@ -381,16 +385,16 @@ void smoothing::optimize_face(GFace* gf){
}
}
- backgroundMesh::set(gf);
+ backgroundMesh::set(gf);
// Create a triangulator
DocRecord triangulator(all.size());
-
+
Range<double> du = gf->parBounds(0) ;
Range<double> dv = gf->parBounds(1) ;
const double LC2D = sqrt((du.high()-du.low())*(du.high()-du.low()) +
- (dv.high()-dv.low())*(dv.high()-dv.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);
@@ -404,21 +408,21 @@ void smoothing::optimize_face(GFace* gf){
Msg::Error("A mesh vertex cannot be reparametrized");
return;
}
- double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
+ double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
(double)RAND_MAX;
- double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
+ double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
(double)RAND_MAX;
triangulator.x(i) = p.x() + XX;
triangulator.y(i) = p.y() + YY;
triangulator.data(i++) = (*it);
}
-
+
// compute the Voronoi diagram
triangulator.Voronoi();
//printf("hullSize = %d\n",triangulator.hullSize());
triangulator.makePosView("LloydInit.pos");
//triangulator.printMedialAxis("medialAxis.pos");
-
+
int exponent;
int num_interior;
int index;
@@ -438,13 +442,13 @@ void smoothing::optimize_face(GFace* gf){
alglib::real_1d_array scales;
wrapper w;
MElementOctree* octree;
-
+
exponent = NORM;
epsg = 0;
epsf = 0;
epsx = 0;
maxits = ITER_MAX;
-
+
num_interior = 0;
for(int i=0;i<triangulator.numPoints;i++){
if(obj.interior(triangulator,gf,i)){
@@ -466,12 +470,12 @@ void smoothing::optimize_face(GFace* gf){
index++;
}
}
-
+
x.setcontent(2*num_interior,initial_conditions);
scales.setcontent(2*num_interior,variables_scales);
- octree = backgroundMesh::current()->get_octree();
-
+ octree = backgroundMesh::current()->get_octree();
+
w.set_p(exponent);
w.set_dimension(2*num_interior);
w.set_face(gf);
@@ -495,8 +499,8 @@ void smoothing::optimize_face(GFace* gf){
obj2.swap();
obj2.compute_parameters(gf,exponent);
obj2.get_gauss();
- obj2.write(triangulator,gf,6);*/
-
+ obj2.write(triangulator,gf,6);*/
+
index = 0;
for(i=0;i<triangulator.numPoints;i++){
if(obj.interior(triangulator,gf,i)){
@@ -505,8 +509,8 @@ void smoothing::optimize_face(GFace* gf){
index++;
}
}
- triangulator.Voronoi();
-
+ triangulator.Voronoi();
+
// now create the vertices
std::vector<MVertex*> mesh_vertices;
for (int i=0; i<triangulator.numPoints;i++){
@@ -523,7 +527,7 @@ void smoothing::optimize_face(GFace* gf){
// destroy the mesh
deMeshGFace killer;
killer(gf);
-
+
// put all additional vertices in the list of
// vertices
gf->additionalVertices = mesh_vertices;
@@ -533,17 +537,17 @@ void smoothing::optimize_face(GFace* gf){
mesher(gf);
// assign those vertices to the face internal vertices
gf->mesh_vertices.insert(gf->mesh_vertices.begin(),
- gf->additionalVertices.begin(),
- gf->additionalVertices.end());
+ gf->additionalVertices.begin(),
+ gf->additionalVertices.end());
// clear the list of additional vertices
- gf->additionalVertices.clear();
+ gf->additionalVertices.clear();
}
void smoothing::optimize_model(){
GFace*gf;
GModel*model = GModel::current();
GModel::fiter it;
-
+
for(it=model->firstFace();it!=model->lastFace();it++)
{
gf = *it;
@@ -575,7 +579,7 @@ double lpcvt::angle(SPoint2 p1,SPoint2 p2,SPoint2 p3){
product = x1*x2 + y1*y2;
angle = product/(norm1*norm2);
angle = 180.0*myacos(angle)/M_PI;
- return angle;
+ return angle;
}
SVector3 lpcvt::normal(SPoint2 p1,SPoint2 p2){
@@ -788,7 +792,7 @@ SPoint2 lpcvt::seed(DocRecord& triangulator,GFace* gf){
int index2;
double x,y;
SPoint2 x0,x1,x2;
-
+
for(i=0;i<triangulator.numPoints;i++){
if(interior(triangulator,gf,i)){
num = triangulator._adjacencies[i].t_length;
@@ -814,14 +818,14 @@ void lpcvt::step1(DocRecord& triangulator,GFace* gf){
int index1,index2,index3;
bool ok_triangle1,ok_triangle2;
SPoint2 x0,x1,x2,x3;
-
+
borders.resize(triangulator.numPoints);
angles.resize(triangulator.numPoints);
for(i=0;i<triangulator.numPoints;i++){
angles[i] = 0.0;
}
temp.resize(triangulator.numPoints);
-
+
for(i=0;i<triangulator.numPoints;i++){
if(!interior(triangulator,gf,i) && !invisible(triangulator,gf,i)){
num = triangulator._adjacencies[i].t_length;
@@ -841,7 +845,7 @@ void lpcvt::step1(DocRecord& triangulator,GFace* gf){
else if(!ok_triangle1 && ok_triangle2){
borders[i].add_segment(i,index2,index3);
}
-
+
if(ok_triangle1){
angles[i] = angles[i] + angle(x0,x1,x2);
}
@@ -871,7 +875,7 @@ void lpcvt::step2(DocRecord& triangulator,GFace* gf){
temp[i].add_vertex(vertex);
}
}
- }
+ }
}
void lpcvt::step3(DocRecord& triangulator,GFace* gf){
@@ -988,7 +992,7 @@ void lpcvt::step3(DocRecord& triangulator,GFace* gf){
vertex2 = voronoi_vertex(val);
temp[i].add_vertex(vertex1);
temp[i].add_vertex(vertex2);
- }
+ }
}
else if(ok_triangle2){
C = circumcircle(triangulator,i,index2,index3);
@@ -1040,7 +1044,7 @@ void lpcvt::step4(DocRecord& triangulator,GFace* gf){
val = intersection(C1,C2,p1,p2,flag1);
if(flag1){
if(vertex1.get_index3()!=-1){
- opposite = vertex1.get_index3();
+ opposite = vertex1.get_index3();
}
else if(vertex1.get_index2()!=-1){
opposite = vertex1.get_index2();
@@ -1095,7 +1099,7 @@ void lpcvt::step5(DocRecord& triangulator,GFace* gf){
val = intersection(p1,p2,p3,p4,flag);
if(flag){
if(vertex1.get_index3()!=-1){
- opposite = vertex1.get_index3();
+ opposite = vertex1.get_index3();
}
else if(vertex1.get_index2()!=-1){
opposite = vertex1.get_index2();
@@ -1180,7 +1184,7 @@ void lpcvt::print_voronoi1(){
p2 = v2.get_point();
p3 = v3.get_point();
print_segment(p2,p3,file);
- }
+ }
file << "};\n";
}
@@ -1204,7 +1208,7 @@ void lpcvt::print_voronoi2(){
}
file << "};\n";
}
-
+
void lpcvt::print_delaunay(DocRecord& triangulator){
int i;
int j;
@@ -1232,7 +1236,7 @@ void lpcvt::print_delaunay(DocRecord& triangulator){
void lpcvt::print_segment(SPoint2 p1,SPoint2 p2,std::ofstream& file){
file << "SL (" << p1.x() << ", " << p1.y() << ", 0, "
<< p2.x() << ", " << p2.y() << ", 0){"
- << "10, 20};\n";
+ << "10, 20};\n";
}
void lpcvt::compute_parameters(GFace* gf,int p){
@@ -1248,8 +1252,8 @@ void lpcvt::compute_parameters(GFace* gf,int p){
voronoi_vertex v1,v2,v3;
metric m;
std::list<voronoi_element>::iterator it;
-
- k = 1.0;
+
+ k = 1.0;
for(it=clipped.begin();it!=clipped.end();it++){
v1 = it->get_v1();
v2 = it->get_v2();
@@ -1277,14 +1281,14 @@ void lpcvt::compute_parameters(GFace* gf,int p){
it->set_v3(v3);
it->deriv_rho(p);
it->set_metric(m);
- }
+ }
}
double lpcvt::get_ratio(GFace* gf,SPoint2 point){
double val;
double uv[2];
double tab[3];
-
+
uv[0] = point.x();
uv[1] = point.y();
buildMetric(gf,uv,tab);
@@ -1297,11 +1301,11 @@ void lpcvt::write(DocRecord& triangulator,GFace* gf,int p){
double energy;
SVector3 grad;
std::vector<SVector3> gradients(triangulator.numPoints);
-
+
eval(triangulator,gradients,energy,p);
-
+
std::ofstream file("gradient");
-
+
for(i=0;i<triangulator.numPoints;i++){
if(interior(triangulator,gf,i)){
grad = gradients[i];
@@ -1326,13 +1330,13 @@ void lpcvt::eval(DocRecord& triangulator,std::vector<SVector3>& gradients,double
SVector3 normal;
voronoi_vertex v1,v2,v3;
std::list<voronoi_element>::iterator it;
-
+
for(i=0;i<gradients.size();i++){
gradients[i] = SVector3(0.0,0.0,0.0);
}
energy = 0.0;
e = 0.000001;
-
+
for(it=clipped.begin();it!=clipped.end();it++){
if(it->get_quality()<e) continue;
v1 = it->get_v1();
@@ -1429,7 +1433,7 @@ double lpcvt::F(voronoi_element element,int p){
C2 = v3.get_point();
energy = 0.0;
m = element.get_metric();
-
+
for(i=0;i<gauss_num;i++){
u = gauss_points(i,0);
v = gauss_points(i,1);
@@ -1469,7 +1473,7 @@ SVector3 lpcvt::simple(voronoi_element element,int p){
comp_y = 0.0;
jacobian = J(generator,C1,C2);
m = element.get_metric();
-
+
for(i=0;i<gauss_num;i++){
u = gauss_points(i,0);
v = gauss_points(i,1);
@@ -1482,7 +1486,7 @@ SVector3 lpcvt::simple(voronoi_element element,int p){
comp_y = comp_y + weight*rho*df_dy(generator,point,m,p);
}
comp_x = jacobian*comp_x;
- comp_y = jacobian*comp_y;
+ comp_y = jacobian*comp_y;
return SVector3(comp_x,comp_y,0.0);
}
@@ -1503,7 +1507,7 @@ SVector3 lpcvt::dF_dC1(voronoi_element element,int p){
SPoint2 point,generator,C1,C2;
voronoi_vertex v1,v2,v3;
metric m;
-
+
v1 = element.get_v1();
v2 = element.get_v2();
v3 = element.get_v3();
@@ -1514,7 +1518,7 @@ SVector3 lpcvt::dF_dC1(voronoi_element element,int p){
comp_y = 0.0;
jacobian = J(generator,C1,C2);
m = element.get_metric();
-
+
for(i=0;i<gauss_num;i++){
u = gauss_points(i,0);
v = gauss_points(i,1);
@@ -1532,7 +1536,7 @@ SVector3 lpcvt::dF_dC1(voronoi_element element,int p){
comp_y = comp_y + weight*rho*df_dy(point,generator,m,p)*u*jacobian;
comp_y = comp_y + weight*rho*distance*(generator.x()-C2.x());
comp_y = comp_y + weight*drho_dy*u*distance*jacobian;
- }
+ }
return SVector3(comp_x,comp_y,0.0);
}
@@ -1553,7 +1557,7 @@ SVector3 lpcvt::dF_dC2(voronoi_element element,int p){
SPoint2 point,generator,C1,C2;
voronoi_vertex v1,v2,v3;
metric m;
-
+
v1 = element.get_v1();
v2 = element.get_v2();
v3 = element.get_v3();
@@ -1564,7 +1568,7 @@ SVector3 lpcvt::dF_dC2(voronoi_element element,int p){
comp_y = 0.0;
jacobian = J(generator,C1,C2);
m = element.get_metric();
-
+
for(i=0;i<gauss_num;i++){
u = gauss_points(i,0);
v = gauss_points(i,1);
@@ -1582,7 +1586,7 @@ SVector3 lpcvt::dF_dC2(voronoi_element element,int p){
comp_y = comp_y + weight*rho*df_dy(point,generator,m,p)*v*jacobian;
comp_y = comp_y + weight*rho*distance*(C1.x()-generator.x());
comp_y = comp_y + weight*drho_dy*v*distance*jacobian;
- }
+ }
return SVector3(comp_x,comp_y,0.0);
}
@@ -1598,7 +1602,7 @@ double lpcvt::f(SPoint2 p1,SPoint2 p2,metric m,int p){
double b;
double c;
double d;
-
+
x1 = p1.x();
y1 = p1.y();
x2 = p2.x();
@@ -1625,7 +1629,7 @@ double lpcvt::df_dx(SPoint2 p1,SPoint2 p2,metric m,int p){
double b;
double c;
double d;
-
+
x1 = p1.x();
y1 = p1.y();
x2 = p2.x();
@@ -1652,7 +1656,7 @@ double lpcvt::df_dy(SPoint2 p1,SPoint2 p2,metric m,int p){
double b;
double c;
double d;
-
+
x1 = p1.x();
y1 = p1.y();
x2 = p2.x();
@@ -1693,7 +1697,7 @@ SVector3 lpcvt::inner_dFdx0(SVector3 dFdC,SPoint2 C,SPoint2 x0,SPoint2 x1,SPoint
det = (x1.x()-x0.x())*(x2.y()-x0.y()) - (x1.y() - x0.y())*(x2.x() - x0.x());
A[0][0] = (x2.y() - x0.y())/det;
A[0][1] = -(x1.y() - x0.y())/det;
- A[1][0] = -(x2.x() - x0.x())/det;
+ A[1][0] = -(x2.x() - x0.x())/det;
A[1][1] = (x1.x() - x0.x())/det;
B[0][0] = C.x() - x0.x();
B[0][1] = C.y() - x0.y();
@@ -1712,9 +1716,9 @@ SVector3 lpcvt::boundary_dFdx0(SVector3 dFdC,SPoint2 C,SPoint2 x0,SPoint2 x1,SVe
fullMatrix<double> M(2,2);
fullMatrix<double> _dFdC(1,2);
fullMatrix<double> _val(1,2);
- A(0,0) = x1.x() - x0.x();
+ A(0,0) = x1.x() - x0.x();
A(0,1) = x1.y() - x0.y();
- A(1,0) = normal.x();
+ A(1,0) = normal.x();
A(1,1) = normal.y();
A.invertInPlace();
B(0,0) = C.x() - x0.x();
@@ -1725,7 +1729,7 @@ SVector3 lpcvt::boundary_dFdx0(SVector3 dFdC,SPoint2 C,SPoint2 x0,SPoint2 x1,SVe
_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);
+ return SVector3(_val(0,0),_val(0,1),0.0);
}
/****************class metric****************/
@@ -1873,7 +1877,7 @@ double voronoi_element::get_rho(double u,double v){
double rho2;
double rho3;
double rho;
-
+
rho1 = v1.get_rho();
rho2 = v2.get_rho();
rho3 = v3.get_rho();
@@ -1927,7 +1931,7 @@ void voronoi_element::deriv_rho(int p){
SPoint2 p1;
SPoint2 p2;
SPoint2 p3;
-
+
rho1 = v1.get_rho();
rho2 = v2.get_rho();
rho3 = v3.get_rho();
@@ -1962,21 +1966,21 @@ double voronoi_element::get_quality(){
double quality;
double l1,l2,l3;
double min_l,max_l;
-
+
x1 = v1.get_point().x();
y1 = v1.get_point().y();
x2 = v2.get_point().x();
y2 = v2.get_point().y();
x3 = v3.get_point().x();
y3 = v3.get_point().y();
-
+
l1 = sqrt(pow_int(x2-x1,2) + pow_int(y2-y1,2));
l2 = sqrt(pow_int(x3-x1,2) + pow_int(y3-y1,2));
l3 = sqrt(pow_int(x3-x2,2) + pow_int(y3-y2,2));
-
+
min_l = std::min(std::min(l1,l2),l3);
max_l = std::max(std::max(l1,l2),l3);
-
+
quality = min_l/max_l;
return quality;
}
@@ -2148,3 +2152,5 @@ MElementOctree* wrapper::get_octree(){
void wrapper::set_octree(MElementOctree* new_octree){
octree = new_octree;
}
+
+#endif
diff --git a/Mesh/meshGFaceLloyd.h b/Mesh/meshGFaceLloyd.h
index 318b592fb8798e50f7019c8e8573697fead4b411..546d7ed438147c1c5ff440f21ee8d26fb8c1662c 100644
--- a/Mesh/meshGFaceLloyd.h
+++ b/Mesh/meshGFaceLloyd.h
@@ -5,15 +5,10 @@
#ifndef _MESH_GFACE_LLOYD_H_
#define _MESH_GFACE_LLOYD_H_
+
+#include <queue>
#include "fullMatrix.h"
#include "DivideAndConquer.h"
-#include <queue>
-#include "ap.h"
-#include "alglibinternal.h"
-#include "alglibmisc.h"
-#include "linalg.h"
-#include "optimization.h"
-#include "MElementOctree.h"
class GFace;
@@ -26,4 +21,4 @@ class smoothing{
void optimize_model();
};
-#endif
\ No newline at end of file
+#endif