diff --git a/Mesh/Levy3D.cpp b/Mesh/Levy3D.cpp
index 1569811400f89d1474aa34af01794e52f1352359..a6fe2d97b96538483f00e0b233a3e472c64a592e 100644
--- a/Mesh/Levy3D.cpp
+++ b/Mesh/Levy3D.cpp
@@ -869,7 +869,7 @@ void LpCVT::compute_parameters(){
double LpCVT::get_size(double x,double y,double z){
//if outside domain return 1.0 (or other value > 0.0)
- return 0.25;
+ return Size_field::search(x,y,z);
}
Tensor LpCVT::get_tensor(double x,double y,double z){
@@ -1491,6 +1491,8 @@ void LpSmoother::improve_region(GRegion* gr)
alglib::real_1d_array alglib_scales;
Frame_field::init_region(gr);
+ Size_field::init_region(gr);
+ Size_field::solve(gr);
octree = new MElementOctree(gr->model());
for(i=0;i<gr->getNumMeshElements();i++){
@@ -1620,6 +1622,7 @@ void LpSmoother::improve_region(GRegion* gr)
delete octree;
free(initial_conditions);
free(scales);
+ Size_field::clear();
Frame_field::clear();
#endif
}
diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp
index ecb078db78186b8e1933db9d95569779d37eb95d..abe57b4835cf9c11c6792f8f540b7af87eeb7abd 100644
--- a/Mesh/directions3D.cpp
+++ b/Mesh/directions3D.cpp
@@ -518,10 +518,12 @@ void Size_field::solve(GRegion* gr){
#if defined(HAVE_PETSC)
linearSystem<double>* system = 0;
system = new linearSystemPETSc<double>;
- //system = new linearSystemFull<double>;
size_t i;
+ int count;
+ int count2;
double val;
+ double volume;
std::set<MVertex*> interior;
std::set<MVertex*>::iterator it;
std::map<MVertex*,double>::iterator it2;
@@ -530,9 +532,13 @@ void Size_field::solve(GRegion* gr){
dofManager<double> assembler(system);
+ count = 0;
for(it2=boundary.begin();it2!=boundary.end();it2++){
assembler.fixVertex(it2->first,0,1,it2->second);
+ count++;
}
+ //printf("\n");
+ //printf("number of boundary vertices = %d\n",count);
for(i=0;i<gr->tetrahedra.size();i++){
interior.insert(gr->tetrahedra[i]->getVertex(0));
@@ -547,13 +553,23 @@ void Size_field::solve(GRegion* gr){
assembler.numberVertex(*it,0,1);
}
}
+
+ for(i=0;i<gr->tetrahedra.size();i++){
+ gr->tetrahedra[i]->setVolumePositive();
+ }
+ count2 = 0;
+ volume = 0.0;
simpleFunction<double> ONE(1.0);
laplaceTerm term(0,1,&ONE);
for(i=0;i<gr->tetrahedra.size();i++){
SElement se(gr->tetrahedra[i]);
term.addToMatrix(assembler,&se);
+ count2++;
+ volume = volume + gr->tetrahedra[i]->getVolume();
}
+ //printf("number of tetrahedra = %d\n",count2);
+ //printf("volume = %f\n",volume);
system->systemSolve();
@@ -617,10 +633,19 @@ double Size_field::get_ratio(GFace* gf,SPoint2 point){
return val;
}
-void Size_field::print_field(){
+void Size_field::print_field(GRegion* gr){
+ size_t i;
double min,max;
double x,y,z;
+ double x1,y1,z1,h1;
+ double x2,y2,z2,h2;
+ double x3,y3,z3,h3;
+ double x4,y4,z4,h4;
std::map<MVertex*,double>::iterator it;
+ std::map<MVertex*,double>::iterator it1;
+ std::map<MVertex*,double>::iterator it2;
+ std::map<MVertex*,double>::iterator it3;
+ std::map<MVertex*,double>::iterator it4;
min = 1000000000.0;
max = -1000000000.0;
@@ -638,15 +663,58 @@ void Size_field::print_field(){
min = it->second;
}
- printf("x = %f, y = %f, z = %f, mesh size = %f\n",x,y,z,it->second);
+ //printf("x = %f, y = %f, z = %f, mesh size = %f\n",x,y,z,it->second);
}
- printf("\n");
+ //printf("\n");
printf("min mesh size = %f\n",min);
printf("max mesh size = %f\n",max);
- printf("%zu\n",boundary.size());
+ printf("total number of vertices = %zu\n",boundary.size());
+ printf("\n");
+
+ std::ofstream file("laplace.pos");
+ file << "View \"test\" {\n";
+
+ for(i=0;i<gr->tetrahedra.size();i++){
+ x1 = gr->tetrahedra[i]->getVertex(0)->x();
+ y1 = gr->tetrahedra[i]->getVertex(0)->y();
+ z1 = gr->tetrahedra[i]->getVertex(0)->z();
+ it1 = boundary.find(gr->tetrahedra[i]->getVertex(0));
+
+ x2 = gr->tetrahedra[i]->getVertex(1)->x();
+ y2 = gr->tetrahedra[i]->getVertex(1)->y();
+ z2 = gr->tetrahedra[i]->getVertex(1)->z();
+ it2 = boundary.find(gr->tetrahedra[i]->getVertex(1));
+
+ x3 = gr->tetrahedra[i]->getVertex(2)->x();
+ y3 = gr->tetrahedra[i]->getVertex(2)->y();
+ z3 = gr->tetrahedra[i]->getVertex(2)->z();
+ it3 = boundary.find(gr->tetrahedra[i]->getVertex(2));
+
+ x4 = gr->tetrahedra[i]->getVertex(3)->x();
+ y4 = gr->tetrahedra[i]->getVertex(3)->y();
+ z4 = gr->tetrahedra[i]->getVertex(3)->z();
+ it4 = boundary.find(gr->tetrahedra[i]->getVertex(3));
+
+ if(it1!=boundary.end() && it2!=boundary.end() && it3!=boundary.end() && it4!=boundary.end()){
+ h1 = it1->second;
+ h2 = it2->second;
+ h3 = it3->second;
+ h4 = it4->second;
+
+ file << "SS ("
+ << x1 << ", " << y1 << ", " << z1 << ", "
+ << x2 << ", " << y2 << ", " << z2 << ", "
+ << x3 << ", " << y3 << ", " << z3 << ", "
+ << x4 << ", " << y4 << ", " << z4 << "){"
+ << h1 << ", " << h2 << ", " << h3 << ", "
+ << h4 << "};\n";
+ }
+ }
+
+ file << "};\n";
}
GRegion* Size_field::test(){
diff --git a/Mesh/directions3D.h b/Mesh/directions3D.h
index 95600ab7515c06c28e4477d8876f34d5bbe296a7..357a9fb74e2b8ec9a54ca1ae52ad07fb4389dac2 100644
--- a/Mesh/directions3D.h
+++ b/Mesh/directions3D.h
@@ -72,7 +72,7 @@ class Size_field{
static void solve(GRegion*);
static double search(double,double,double);
static double get_ratio(GFace*,SPoint2);
- static void print_field();
+ static void print_field(GRegion* gr);
static GRegion* test();
static void clear();
};
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index 4d10a5ce58b04e4fe34ab2006071f870a8ad58d5..4b163188afffedc2f6deccfba4ef5b96d5b264f7 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -364,6 +364,49 @@ void callback(const alglib::real_1d_array& x,double& func,alglib::real_1d_array&
}
}
+void verification(alglib::real_1d_array& x,void* ptr){
+ int num;
+ int index;
+ int dimension;
+ double e;
+ double func;
+ double R,L,U,D;
+ wrapper* w;
+ DocRecord* pointer;
+
+ w = static_cast<wrapper*>(ptr);
+ dimension = w->get_dimension();
+ pointer = w->get_triangulator();
+ num = pointer->numPoints;
+ srand(time(NULL));
+ index = rand()%num;
+ e = 0.0000001;
+
+ alglib::real_1d_array grad;
+ grad.setlength(dimension);
+
+ x[index] = x[index] + e;
+ callback(x,R,grad,ptr);
+ x[index] = x[index] - e;
+
+ x[index] = x[index] - e;
+ callback(x,L,grad,ptr);
+ x[index] = x[index] + e;
+
+ x[index + dimension/2] = x[index + dimension/2] + e;
+ callback(x,U,grad,ptr);
+ x[index + dimension/2] = x[index + dimension/2] - e;
+
+ x[index + dimension/2] = x[index + dimension/2] - e;
+ callback(x,D,grad,ptr);
+ x[index + dimension/2] = x[index + dimension/2] + e;
+
+ callback(x,func,grad,ptr);
+
+ printf("%f %f\n",(R-L)/(2.0*e),(U-D)/(2.0*e));
+ printf("%f %f\n",grad[index],grad[index + dimension/2]);
+}
+
/****************class smoothing****************/
smoothing::smoothing(int param1,int param2){
@@ -483,6 +526,10 @@ void smoothing::optimize_face(GFace* gf){
w.set_triangulator(&triangulator);
w.set_octree(octree);
+ /*if(num_interior>1){
+ verification(x,&w);
+ }*/
+
if(num_interior>1){
minlbfgscreate(2*num_interior,4,x,state);
minlbfgssetscale(state,scales);
diff --git a/Mesh/simple3D.cpp b/Mesh/simple3D.cpp
index 38527c783dce169b3ef72a44ba01bceec0b30123..308719153697c8f0fd3636b7084755b321743bab 100644
--- a/Mesh/simple3D.cpp
+++ b/Mesh/simple3D.cpp
@@ -346,6 +346,8 @@ void Filler::treat_region(GRegion* gr){
RTree<Node*,double,3,double> rtree;
Frame_field::init_region(gr);
+ Size_field::init_region(gr);
+ Size_field::solve(gr);
octree = new MElementOctree(gr->model());
for(i=0;i<gr->getNumMeshElements();i++){
@@ -408,12 +410,11 @@ void Filler::treat_region(GRegion* gr){
if(inside_domain(octree,x,y,z)){
compute_parameters(spawn,gr);
- if(far_from_boundary(octree,spawn)){
+ if(far_from_boundary(octree,spawn)){
wrapper.set_too_close(0);
wrapper.set_spawn(spawn);
wrapper.set_parent(parent);
rtree.Search(spawn->min,spawn->max,rtree_callback,&wrapper);
-
if(!wrapper.get_too_close()){
fifo.push(spawn);
rtree.Insert(spawn->min,spawn->max,spawn);
@@ -425,6 +426,7 @@ void Filler::treat_region(GRegion* gr){
}
if(!ok) delete spawn;
}
+
printf("%d\n",count);
count++;
}
@@ -440,6 +442,7 @@ void Filler::treat_region(GRegion* gr){
for(i=0;i<garbage.size();i++) delete garbage[i];
for(i=0;i<new_vertices.size();i++) delete new_vertices[i];
new_vertices.clear();
+ Size_field::clear();
Frame_field::clear();
#endif
}
@@ -466,7 +469,7 @@ Metric Filler::get_metric(double x,double y,double z){
}
double Filler::get_size(double x,double y,double z){
- return 0.25;
+ return Size_field::search(x,y,z);
}
double Filler::get_size(double x,double y,double z,GEntity* ge){
@@ -606,19 +609,22 @@ double Filler::improvement(GEntity* ge,MElementOctree* octree,SPoint3 point,doub
x = point.x() + h_nearer*direction.x();
y = point.y() + h_nearer*direction.y();
z = point.z() + h_nearer*direction.z();
+
if(inside_domain(octree,x,y,z)){
h_farther = get_size(x,y,z);
}
else h_farther = h_nearer;
coeffA = 1.0;
- coeffB = 0.2;
+ coeffB = 0.16;
+
if(h_farther>h_nearer){
average = coeffA*h_nearer + (1.0-coeffA)*h_farther;
}
else{
average = coeffB*h_nearer + (1.0-coeffB)*h_farther;
}
+
return average;
}
@@ -690,6 +696,4 @@ void Filler::print_node(Node* node,std::ofstream& file){
/*********static declarations*********/
-std::vector<MVertex*> Filler::new_vertices;
-
-
+std::vector<MVertex*> Filler::new_vertices;
\ No newline at end of file
diff --git a/Mesh/yamakawa.cpp b/Mesh/yamakawa.cpp
index 277584d3509f622da8f95e28648b3eb0a9c71285..620711b00e2ed736424f96222de11df7331aa681 100755
--- a/Mesh/yamakawa.cpp
+++ b/Mesh/yamakawa.cpp
@@ -2922,7 +2922,7 @@ PostOp::PostOp(){}
PostOp::~PostOp(){}
-void PostOp::execute(){
+void PostOp::execute(bool flag){
GRegion* gr;
GModel* model = GModel::current();
GModel::riter it;
@@ -2931,12 +2931,12 @@ void PostOp::execute(){
{
gr = *it;
if(gr->getNumMeshElements()>0){
- execute(gr);
+ execute(gr,flag);
}
}
}
-void PostOp::execute(GRegion* gr){
+void PostOp::execute(GRegion* gr,bool flag){
printf("................PYRAMIDS................\n");
estimate1 = 0;
estimate2 = 0;
@@ -2947,11 +2947,13 @@ void PostOp::execute(GRegion* gr){
pyramids1(gr);
rearrange(gr);
- init_markings(gr);
- build_vertex_to_tetrahedra(gr);
- build_vertex_to_pyramids(gr);
- pyramids2(gr);
- rearrange(gr);
+ if(flag){
+ init_markings(gr);
+ build_vertex_to_tetrahedra(gr);
+ build_vertex_to_pyramids(gr);
+ pyramids2(gr);
+ rearrange(gr);
+ }
statistics(gr);
}
diff --git a/Mesh/yamakawa.h b/Mesh/yamakawa.h
index c19407497546bfd26d02d623ddd714a99164d42a..6775221e5ebf3e8f5ff18110953a7bbfa84133a9 100755
--- a/Mesh/yamakawa.h
+++ b/Mesh/yamakawa.h
@@ -242,8 +242,8 @@ class PostOp{
PostOp();
~PostOp();
- void execute();
- void execute(GRegion*);
+ void execute(bool);
+ void execute(GRegion*,bool);
void init_markings(GRegion*);
void pyramids1(GRegion*);