From 6f9bca4e7bbc31127fe40e1a5bb3748dab5192d9 Mon Sep 17 00:00:00 2001
From: Tristan Carrier Baudouin <tristan.carrier@uclouvain.be>
Date: Tue, 31 May 2011 09:53:26 +0000
Subject: [PATCH] lpcvt
---
Geo/GFace.cpp | 6 +-
Geo/GFaceCompound.cpp | 3 +
Mesh/BackgroundMesh.cpp | 118 ++++++++++++-------------
Mesh/meshGFace.cpp | 2 +-
Mesh/meshGFaceLloyd.cpp | 182 ++++++++++++++++++++-------------------
Mesh/meshGFaceLloyd.h | 9 +-
benchmarks/stl/skull.geo | 2 +-
7 files changed, 160 insertions(+), 162 deletions(-)
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 1333a3ba86..ee6bd136ce 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1164,8 +1164,10 @@ bool GFace::fillPointCloud(double maxDist, std::vector<SPoint3> *points,
void GFace::lloyd(int nbiter, int infn)
{
#if defined(HAVE_MESH) && defined(HAVE_BFGS)
- lloydAlgorithm algo(nbiter, infn);
- algo(this);
+ smoothing s = smoothing(nbiter,infn);
+ s.optimize_face(this);
+ //lloydAlgorithm algo(nbiter, infn);
+ //algo(this);
#endif
}
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 8889469e39..2b03cb029f 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -2009,6 +2009,9 @@ bool GFaceCompound::checkTopology() const
//TODO: smthg to exit here for lloyd remeshing
+ //fixme tristan
+ //return true;
+
bool correctTopo = true;
if(allNodes.empty()) buildAllNodes();
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 65b2a4fdf2..48aa0a2a0b 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -26,7 +26,6 @@
#include "linearSystemCSR.h"
#include "linearSystemFull.h"
#include "linearSystemPETSc.h"
-#include "SElement.h"
#endif
// computes the characteristic length of the mesh at a vertex in order
@@ -40,10 +39,10 @@ SMetric3 buildMetricTangentToCurve (SVector3 &t, double curvature, double &lambd
SVector3 a;
if (t(0) <= t(1) && t(0) <= t(2)){
a = SVector3(1,0,0);
- }
+ }
else if (t(1) <= t(0) && t(1) <= t(2)){
a = SVector3(0,1,0);
- }
+ }
else{
a = SVector3(0,0,1);
}
@@ -64,10 +63,10 @@ SMetric3 max_edge_curvature_metric(const GVertex *gv, double &length)
SMetric3 val (1.e-12);
length = 1.e22;
std::list<GEdge*> l_edges = gv->edges();
- for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
+ for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
ite != l_edges.end(); ++ite){
GEdge *_myGEdge = *ite;
- Range<double> range = _myGEdge->parBounds(0);
+ Range<double> range = _myGEdge->parBounds(0);
SMetric3 cc;
double l;
if (gv == _myGEdge->getBeginVertex()) {
@@ -89,17 +88,17 @@ SMetric3 max_edge_curvature_metric(const GVertex *gv, double &length)
SMetric3 max_edge_curvature_metric(const GEdge *ge, double u, double &l){
SVector3 t = ge->firstDer(u);
t.normalize();
- return buildMetricTangentToCurve(t,ge->curvature(u),l);
+ return buildMetricTangentToCurve(t,ge->curvature(u),l);
}
static double max_edge_curvature(const GVertex *gv)
{
double val = 0;
std::list<GEdge*> l_edges = gv->edges();
- for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
+ for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
ite != l_edges.end(); ++ite){
GEdge *_myGEdge = *ite;
- Range<double> range = _myGEdge->parBounds(0);
+ Range<double> range = _myGEdge->parBounds(0);
double cc;
if (gv == _myGEdge->getBeginVertex()) cc = _myGEdge->curvature(range.low());
else cc = _myGEdge->curvature(range.high());
@@ -120,7 +119,7 @@ static double max_surf_curvature(const GEdge *ge, double u)
val = std::max(cc, val);
}
++it;
- }
+ }
return val;
}
@@ -128,7 +127,7 @@ static double max_surf_curvature(const GVertex *gv)
{
double val = 0;
std::list<GEdge*> l_edges = gv->edges();
- for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
+ for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
ite != l_edges.end(); ++ite){
GEdge *_myGEdge = *ite;
Range<double> bounds = _myGEdge->parBounds(0);
@@ -163,7 +162,7 @@ static SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, dou
Z.x(),Z.y(),Z.z(),
lambda1,lambda2);
*/
- SMetric3 curvMetric (1./(lambda1*lambda1),1./(lambda2*lambda2),1.e-5,
+ SMetric3 curvMetric (1./(lambda1*lambda1),1./(lambda2*lambda2),1.e-5,
dirMin, dirMax, Z );
return curvMetric;
}
@@ -183,7 +182,7 @@ static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u)
count++;
}
++it;
- }
+ }
double Crv = ge->curvature(u);
double lambda = ((2 * M_PI) /( fabs(Crv) * CTX::instance()->mesh.minCircPoints ) );
SMetric3 curvMetric (1./(lambda*lambda));
@@ -194,7 +193,7 @@ static SMetric3 metric_based_on_surface_curvature(const GVertex *gv)
{
SMetric3 mesh_size(1.e-5);
std::list<GEdge*> l_edges = gv->edges();
- for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
+ for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
ite != l_edges.end(); ++ite){
GEdge *_myGEdge = *ite;
Range<double> bounds = _myGEdge->parBounds(0);
@@ -204,7 +203,7 @@ static SMetric3 metric_based_on_surface_curvature(const GVertex *gv)
metric_based_on_surface_curvature(_myGEdge, bounds.low()));
else
mesh_size = intersection
- (mesh_size,
+ (mesh_size,
metric_based_on_surface_curvature(_myGEdge, bounds.high()));
}
return mesh_size;
@@ -220,7 +219,7 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
{
double Crv = 0;
switch(ge->dim()){
- case 0:
+ case 0:
Crv = max_edge_curvature((const GVertex *)ge);
Crv = std::max(max_surf_curvature((const GVertex *)ge), Crv);
//Crv = max_surf_curvature((const GVertex *)ge);
@@ -230,7 +229,7 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
GEdge *ged = (GEdge *)ge;
Crv = ged->curvature(U)*2;
Crv = std::max(Crv, max_surf_curvature(ged, U));
- //Crv = max_surf_curvature(ged, U);
+ //Crv = max_surf_curvature(ged, U);
}
break;
case 2:
@@ -240,7 +239,7 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
}
break;
}
-
+
double lc = Crv > 0 ? 2 * M_PI / Crv / CTX::instance()->mesh.minCircPoints : MAX_LC;
return lc;
}
@@ -275,16 +274,16 @@ static double LC_MVertex_PNTS(GEntity *ge, double U, double V)
GVertex *v1 = ged->getBeginVertex();
GVertex *v2 = ged->getEndVertex();
if (v1 && v2){
- Range<double> range = ged->parBounds(0);
- double a = (U - range.low()) / (range.high() - range.low());
+ Range<double> range = ged->parBounds(0);
+ double a = (U - range.low()) / (range.high() - range.low());
double lc = (1 - a) * v1->prescribedMeshSizeAtVertex() +
(a) * v2->prescribedMeshSizeAtVertex() ;
// FIXME we might want to remove this to make all lc treatment consistent
if(lc >= MAX_LC) return CTX::instance()->lc / 10.;
return lc;
}
- else
- return MAX_LC;
+ else
+ return MAX_LC;
}
default:
return MAX_LC;
@@ -292,7 +291,7 @@ static double LC_MVertex_PNTS(GEntity *ge, double U, double V)
}
// This is the only function that is used by the meshers
-double BGM_MeshSize(GEntity *ge, double U, double V,
+double BGM_MeshSize(GEntity *ge, double U, double V,
double X, double Y, double Z)
{
// default lc (mesh size == size of the model)
@@ -300,7 +299,7 @@ double BGM_MeshSize(GEntity *ge, double U, double V,
// lc from points
double l2 = MAX_LC;
- if(CTX::instance()->mesh.lcFromPoints && ge->dim() < 2)
+ if(CTX::instance()->mesh.lcFromPoints && ge->dim() < 2)
l2 = LC_MVertex_PNTS(ge, U, V);
// lc from curvature
@@ -336,18 +335,18 @@ double BGM_MeshSize(GEntity *ge, double U, double V,
// for now, only works with bamg in 2D, work
// in progress
-SMetric3 BGM_MeshMetric(GEntity *ge,
- double U, double V,
+SMetric3 BGM_MeshMetric(GEntity *ge,
+ double U, double V,
double X, double Y, double Z)
{
// default lc (mesh size == size of the model)
double l1 = CTX::instance()->lc;
- // lc from points
+ // lc from points
double l2 = MAX_LC;
- if(CTX::instance()->mesh.lcFromPoints && ge->dim() < 2)
+ if(CTX::instance()->mesh.lcFromPoints && ge->dim() < 2)
l2 = LC_MVertex_PNTS(ge, U, V);
-
+
// lc from curvature
SMetric3 l3(1./(MAX_LC*MAX_LC));
if(CTX::instance()->mesh.lcFromCurvature && ge->dim() < 3)
@@ -367,15 +366,9 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
}
}
}
-
+
// take the minimum, then constrain by lcMin and lcMax
double lc = std::min(l1, l2);
-
- if (backgroundMesh::current()){
- const double lcBG = backgroundMesh::current()->operator() (U,V,0);
- lc = std::min(lc,lcBG);
- }
-
lc = std::max(lc, CTX::instance()->mesh.lcMin);
lc = std::min(lc, CTX::instance()->mesh.lcMax);
@@ -384,14 +377,11 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
lc, CTX::instance()->mesh.lcMin, CTX::instance()->mesh.lcMax);
lc = l1;
}
-
+
SMetric3 LC(1./(lc*lc));
SMetric3 m = intersection(intersection (l4, LC),l3);
// printf("%g %g %g %g %g %g\n",m(0,0),m(1,1),m(2,2),m(0,1),m(0,2),m(1,2));
-
-
-
- return m;
+ return m;
// return lc * CTX::instance()->mesh.lcFactor;
}
@@ -433,7 +423,7 @@ backgroundMesh::backgroundMesh(GFace *_gf)
MVertex *newv =0;
if (it == _3Dto2D.end()){
SPoint2 p;
- bool success = reparamMeshVertexOnFace(e->getVertex(j), _gf, p);
+ bool success = reparamMeshVertexOnFace(e->getVertex(j), _gf, p);
newv = new MVertex (p.x(), p.y(), 0.0);
_vertices.push_back(newv);
_3Dto2D[e->getVertex(j)] = newv;
@@ -446,7 +436,7 @@ backgroundMesh::backgroundMesh(GFace *_gf)
}
// build a search structure
- _octree = new MElementOctree(_triangles);
+ _octree = new MElementOctree(_triangles);
// compute the mesh sizes at nodes
if (CTX::instance()->mesh.lcFromPoints)
@@ -457,8 +447,8 @@ backgroundMesh::backgroundMesh(GFace *_gf)
_sizes[itv2->first] = MAX_LC;
}
}
- // ensure that other criteria are fullfilled
- // updateSizes(_gf);
+ // ensure that other criteria are fullfilled
+ updateSizes(_gf);
// compute optimal mesh orientations
propagatecrossField(_gf);
@@ -474,7 +464,7 @@ backgroundMesh::~backgroundMesh()
delete _octree;
}
-static void propagateValuesOnFace (GFace *_gf,
+static void propagateValuesOnFace (GFace *_gf,
std::map<MVertex*,double> &dirichlet)
{
linearSystem<double> *_lsys = 0;
@@ -484,7 +474,7 @@ static void propagateValuesOnFace (GFace *_gf,
linearSystemGmm<double> *_lsysb = new linearSystemGmm<double>;
_lsysb->setGmres(1);
_lsys = _lsysb;
-#elif defined(HAVE_TAUCS)
+#elif defined(HAVE_TAUCS)
_lsys = new linearSystemCSRTaucs<double>;
#else
_lsys = new linearSystemFull<double>;
@@ -514,7 +504,7 @@ static void propagateValuesOnFace (GFace *_gf,
for (unsigned int k = 0; k < _gf->triangles.size(); k++){
MTriangle *t = _gf->triangles[k];
SElement se(t);
- l.addToMatrix(myAssembler, &se);
+ l.addToMatrix(myAssembler, &se);
}
// Solve
@@ -524,7 +514,7 @@ static void propagateValuesOnFace (GFace *_gf,
for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
myAssembler.getDofValue(*it, 0, 1, dirichlet[*it]);
}
-
+
delete _lsys;
}
@@ -553,7 +543,7 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
replaceMeshCompound(_gf, e);
std::list<GEdge*>::const_iterator it = e.begin();
std::map<MVertex*,double> sizes;
-
+
for( ; it != e.end(); ++it ){
if (!(*it)->isSeam(_gf)){
for(unsigned int i = 0; i < (*it)->lines.size(); i++ ){
@@ -568,9 +558,9 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
std::map<MVertex*, double>::iterator itv = sizes.find(v);
if (itv == sizes.end())
sizes[v] = log(d);
- else
+ else
itv->second = 0.5 * (itv->second + log(d));
- }
+ }
}
}
}
@@ -585,11 +575,11 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
}
}
-// C R O S S F I E L D S
+// C R O S S F I E L D S
crossField2d :: crossField2d (MVertex* v, GEdge* ge){
double p;
- bool success = reparamMeshVertexOnEdge(v, ge, p);
+ bool success = reparamMeshVertexOnEdge(v, ge, p);
if (!success){
Msg::Warning("cannot reparametrize a point in crossField");
_angle = 0;
@@ -613,14 +603,14 @@ void backgroundMesh::propagatecrossField(GFace *_gf)
replaceMeshCompound(_gf, e);
std::list<GEdge*>::const_iterator it = e.begin();
-
+
for( ; it != e.end(); ++it ){
if (!(*it)->isSeam(_gf)){
for(unsigned int i = 0; i < (*it)->lines.size(); i++ ){
MVertex *v[2];
v[0] = (*it)->lines[i]->getVertex(0);
v[1] = (*it)->lines[i]->getVertex(1);
- SPoint2 p1,p2;
+ SPoint2 p1,p2;
bool success = reparamMeshEdgeOnFace(v[0],v[1],_gf,p1,p2);
double angle = atan2 ( p1.y()-p2.y() , p1.x()-p2.x() );
//double angle = atan2 ( v[0]->y()-v[1]->y() , v[0]->x()- v[1]->x() );
@@ -630,8 +620,8 @@ void backgroundMesh::propagatecrossField(GFace *_gf)
std::map<MVertex*,double>::iterator itc = _cosines4.find(v[i]);
std::map<MVertex*,double>::iterator its = _sines4.find(v[i]);
if (itc != _cosines4.end()){
- itc->second = 0.5*(itc->second + cos(4*angle));
- its->second = 0.5*(its->second + sin(4*angle));
+ itc->second = 0.5*(itc->second + cos(4*angle));
+ its->second = 0.5*(its->second + sin(4*angle));
}
else {
_cosines4[v[i]] = cos(4*angle);
@@ -695,7 +685,7 @@ void backgroundMesh::propagatecrossField(GFace *_gf)
void backgroundMesh::updateSizes(GFace *_gf)
{
std::map<MVertex*,double>::iterator itv = _sizes.begin();
- for ( ; itv != _sizes.end(); ++itv){
+ for ( ; itv != _sizes.end(); ++itv){
SPoint2 p;
MVertex *v = _2Dto3D[itv->first];
double lc;
@@ -708,14 +698,14 @@ void backgroundMesh::updateSizes(GFace *_gf)
lc = BGM_MeshSize(v->onWhat(), u, 0, v->x(), v->y(), v->z());
}
else{
- bool success = reparamMeshVertexOnFace(v, _gf, p);
+ bool success = reparamMeshVertexOnFace(v, _gf, p);
lc = BGM_MeshSize(_gf, p.x(), p.y(), v->x(), v->y(), v->z());
}
// printf("2D -- %g %g 3D -- %g %g\n",p.x(),p.y(),v->x(),v->y());
itv->second = std::min(lc,itv->second);
itv->second = std::max(itv->second, CTX::instance()->mesh.lcMin);
itv->second = std::min(itv->second, CTX::instance()->mesh.lcMax);
- }
+ }
}
double backgroundMesh::operator() (double u, double v, double w) const
@@ -750,11 +740,11 @@ double backgroundMesh::getAngle (double u, double v, double w) const
std::map<MVertex*,double>::const_iterator itv2 = _angles.find(e->getVertex(1));
std::map<MVertex*,double>::const_iterator itv3 = _angles.find(e->getVertex(2));
- double cos4 = cos (4*itv1->second) * (1-uv2[0]-uv2[1]) +
- cos (4*itv2->second) * uv2[0] +
+ double cos4 = cos (4*itv1->second) * (1-uv2[0]-uv2[1]) +
+ cos (4*itv2->second) * uv2[0] +
cos (4*itv3->second) * uv2[1] ;
- double sin4 = sin (4*itv1->second) * (1-uv2[0]-uv2[1]) +
- sin (4*itv2->second) * uv2[0] +
+ double sin4 = sin (4*itv1->second) * (1-uv2[0]-uv2[1]) +
+ sin (4*itv2->second) * uv2[0] +
sin (4*itv3->second) * uv2[1] ;
double angle = atan2(sin4,cos4)/4.0;
crossField2d::normalizeAngle (angle);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index c848f82ab1..609a15a4e7 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -406,7 +406,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
// if necessary split compound and remesh parts
bool isMeshed = false;
if(gf->geomType() == GEntity::CompoundSurface){
- isMeshed = checkMeshCompound((GFaceCompound*) gf, edges);
+ isMeshed = checkMeshCompound((GFaceCompound*) gf, edges);
if (isMeshed) return true;
}
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index 06a238f157..ff9f9a9f00 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -7,7 +7,7 @@
#include "MTriangle.h"
#include "Context.h"
#include "meshGFace.h"
-#include "BackgroundMesh.h"
+#include "BackgroundMesh.h"
#include <fstream>
#include "ap.h"
#include "alglibinternal.h"
@@ -16,7 +16,7 @@
#include "optimization.h"
#include "polynomialBasis.h"
#include "MElementOctree.h"
-#include "meshGFaceOptimize.h"
+
/****************fonction callback****************/
@@ -47,8 +47,8 @@ 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);
+
+ w = static_cast<wrapper*>(ptr);
p = w->get_p();
dimension = w->get_dimension();
gf = w->get_face();
@@ -62,7 +62,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)){
@@ -80,7 +80,7 @@ void callback(const alglib::real_1d_array& x,double& func,alglib::real_1d_array&
if(iteration>=max){
error2 = 1;
}
-
+
if(!error1 && !error2){
pointer->Voronoi();
pointer->build_edges();
@@ -98,7 +98,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++){
@@ -110,15 +110,15 @@ void callback(const alglib::real_1d_array& x,double& func,alglib::real_1d_array&
if(error1){
printf("Vertices outside domain.\n");
}
-
+
if(error2){
printf("Oscillations.\n");
}
-
+
if(error3){
printf("Boundary intersection.\n");
- }
-
+ }
+
if(start>0.0){
printf("%d %.3f\n",iteration,100.0*(start-energy)/start);
}
@@ -126,20 +126,20 @@ void callback(const alglib::real_1d_array& x,double& func,alglib::real_1d_array&
w->set_start(energy);
}
w->set_iteration(iteration+1);
-
+
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();
}
- }
+ }
}
bool domain_search(MElementOctree* octree,GFace* gf,double x,double y){
GPoint gp;
MElement* element;
-
+
gp = gf->point(x,y);
element = (MElement*)octree->find(gp.x(),gp.y(),gp.z(),2,true);
if(element!=NULL) return 1;
@@ -148,10 +148,14 @@ bool domain_search(MElementOctree* octree,GFace* gf,double x,double y){
-/****************class lloydAlgorithm****************/
+/****************class smoothing****************/
+
+smoothing::smoothing(int param1,int param2){
+ ITER_MAX = param1;
+ NORM = param2;
+}
-void lloydAlgorithm::operator () (GFace *gf)
-{
+void smoothing::optimize_face(GFace* gf){
std::set<MVertex*> all;
// get all the points of the face ...
@@ -165,16 +169,16 @@ void lloydAlgorithm::operator () (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);
@@ -188,15 +192,15 @@ void lloydAlgorithm::operator () (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);
}
-
+
triangulator.Voronoi();
triangulator.initialize();
int index,number,count,max;
@@ -215,8 +219,8 @@ void lloydAlgorithm::operator () (GFace *gf)
if(flag) count++;
}
}
- triangulator.remove_all();
-
+ triangulator.remove_all();
+
triangulator.Voronoi();
double delta;
delta = 0.01;
@@ -227,7 +231,7 @@ void lloydAlgorithm::operator () (GFace *gf)
//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 index1 = -1;
int index2 = -1;
@@ -249,7 +253,7 @@ void lloydAlgorithm::operator () (GFace *gf)
else if(index6==-1) index6 = i;
else if(index7==-1) index7 = i;
else if(index8==-1) index8 = i;
-
+
}
}
/*triangulator.points[index1].where.h = 0.01;
@@ -268,13 +272,13 @@ void lloydAlgorithm::operator () (GFace *gf)
triangulator.points[index7].where.v = 0.003;
triangulator.points[index8].where.h = 0.530;
triangulator.points[index8].where.v = 0.004;*/
-
+
// 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;
double epsg;
@@ -288,13 +292,13 @@ void lloydAlgorithm::operator () (GFace *gf)
alglib::real_1d_array x;
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)){
@@ -310,25 +314,25 @@ void lloydAlgorithm::operator () (GFace *gf)
index++;
}
}
-
+
x.setcontent(2*num_interior,initial_conditions);
- octree = new MElementOctree(gf->model());
-
+ octree = new MElementOctree(gf->model());
+
w.set_p(exponent);
w.set_dimension(2*num_interior);
w.set_face(gf);
w.set_max(2*ITER_MAX);
w.set_triangulator(&triangulator);
w.set_octree(octree);
-
+
minlbfgscreate(2*num_interior,4,x,state);
minlbfgssetcond(state,epsg,epsf,epsx,maxits);
minlbfgsoptimize(state,callback,NULL,&w);
minlbfgsresults(state,x,rep);
- delete octree;
-
+ delete octree;
+
/*lpcvt obj2;
SPoint2 val = obj2.seed(triangulator,gf);
triangulator.concave(val.x(),val.y(),gf);
@@ -336,8 +340,8 @@ void lloydAlgorithm::operator () (GFace *gf)
obj2.clip_cells(triangulator,gf);
obj2.swap();
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)){
@@ -346,8 +350,8 @@ void lloydAlgorithm::operator () (GFace *gf)
index++;
}
}
- triangulator.Voronoi();
-
+ triangulator.Voronoi();
+
//lpcvt obj;
//SPoint2 val = obj.seed(triangulator,gf);
//triangulator.concave(val.x(),val.y(),gf);
@@ -359,7 +363,7 @@ void lloydAlgorithm::operator () (GFace *gf)
//obj.swap();
//obj.get_gauss();
//obj.write(triangulator,gf,6);
-
+
// now create the vertices
std::vector<MVertex*> mesh_vertices;
for (int i=0; i<triangulator.numPoints;i++){
@@ -376,7 +380,7 @@ void lloydAlgorithm::operator () (GFace *gf)
// destroy the mesh
deMeshGFace killer;
killer(gf);
-
+
// put all additional vertices in the list of
// vertices
gf->additionalVertices = mesh_vertices;
@@ -386,23 +390,23 @@ void lloydAlgorithm::operator () (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 lloydAlgorithm::optimize(int itermax,int norm){
- GFace*face;
+void smoothing::optimize_model(){
+ GFace*gf;
GModel*model = GModel::current();
- GModel::fiter iterator;
- lloydAlgorithm lloyd(itermax,norm);
- for(iterator = model->firstFace();iterator != model->lastFace();iterator++)
+ GModel::fiter it;
+
+ for(it=model->firstFace();it!=model->lastFace();it++)
{
- face = *iterator;
- if(face->getNumMeshElements() > 0){
- lloyd(face);
- recombineIntoQuads(face,1,1);
+ gf = *it;
+ if(gf->getNumMeshElements()>0){
+ optimize_face(gf);
+ recombineIntoQuads(gf,1,1);
}
}
}
@@ -430,7 +434,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){
@@ -670,14 +674,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;
@@ -697,7 +701,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);
}
@@ -727,7 +731,7 @@ void lpcvt::step2(DocRecord& triangulator,GFace* gf){
temp[i].add_vertex(vertex);
}
}
- }
+ }
}
void lpcvt::step3(DocRecord& triangulator,GFace* gf){
@@ -844,7 +848,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);
@@ -896,7 +900,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();
@@ -951,7 +955,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();
@@ -1037,7 +1041,7 @@ void lpcvt::print_voronoi1(){
p2 = v2.get_point();
p3 = v3.get_point();
print_segment(p2,p3,file);
- }
+ }
file << "};\n";
}
@@ -1061,7 +1065,7 @@ void lpcvt::print_voronoi2(){
}
file << "};\n";
}
-
+
void lpcvt::print_delaunay(DocRecord& triangulator){
int i;
int j;
@@ -1089,7 +1093,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_metrics(DocRecord& triangulator){
@@ -1101,18 +1105,18 @@ void lpcvt::compute_metrics(DocRecord& triangulator){
double sinus;
SPoint2 point;
metric m;
-
- metrics.resize(triangulator.numPoints);
-
+
+ metrics.resize(triangulator.numPoints);
+
for(i=0;i<triangulator.numPoints;i++){
point = convert(triangulator,i);
x = point.x();
y = point.y();
- angle = backgroundMesh::current()->getAngle(x,y,0.0); //-myatan2(y,x);
+ angle = -backgroundMesh::current()->getAngle(x,y,0.0); //-myatan2(y,x);
cosinus = cos(angle);
sinus = sin(angle);
m = metric(cosinus,-sinus,sinus,cosinus);
- metrics[i] = m;
+ metrics[i] = m;
}
}
@@ -1121,7 +1125,7 @@ double lpcvt::get_rho(SPoint2 point,int p){
double y;
double h;
double rho;
-
+
x = point.x();
y = point.y();
h = backgroundMesh::current()->operator()(x,y,0.0); //0.1;
@@ -1148,7 +1152,7 @@ double lpcvt::drho_dx(SPoint2 point,int p){
SPoint2 less1;
SPoint2 plus1;
SPoint2 plus2;
-
+
x = point.x();
y = point.y();
e = 0.00000001;
@@ -1190,7 +1194,7 @@ double lpcvt::drho_dy(SPoint2 point,int p){
SPoint2 less1;
SPoint2 plus1;
SPoint2 plus2;
-
+
x = point.x();
y = point.y();
e = 0.00000001;
@@ -1223,11 +1227,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];
@@ -1251,12 +1255,12 @@ 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;
-
+
for(it=clipped.begin();it!=clipped.end();it++){
v1 = it->get_v1();
v2 = it->get_v2();
@@ -1379,7 +1383,7 @@ SVector3 lpcvt::simple(SPoint2 generator,SPoint2 C1,SPoint2 C2,int p,int index){
comp_y = comp_y + weight*rho*df_dy(generator,point,p,index);
}
comp_x = jacobian*comp_x;
- comp_y = jacobian*comp_y;
+ comp_y = jacobian*comp_y;
return SVector3(comp_x,comp_y,0.0);
}
@@ -1412,7 +1416,7 @@ SVector3 lpcvt::dF_dC1(SPoint2 generator,SPoint2 C1,SPoint2 C2,int p,int index){
comp_y = comp_y + weight*rho*df_dy(point,generator,p,index)*u*jacobian;
comp_y = comp_y + weight*rho*f(point,generator,p,index)*(generator.x()-C2.x());
comp_y = comp_y + weight*drho_dy(point,p)*u*f(point,generator,p,index)*jacobian;
- }
+ }
return SVector3(comp_x,comp_y,0.0);
}
@@ -1445,7 +1449,7 @@ SVector3 lpcvt::dF_dC2(SPoint2 generator,SPoint2 C1,SPoint2 C2,int p,int index){
comp_y = comp_y + weight*rho*df_dy(point,generator,p,index)*v*jacobian;
comp_y = comp_y + weight*rho*f(point,generator,p,index)*(C1.x()-generator.x());
comp_y = comp_y + weight*drho_dy(point,p)*v*f(point,generator,p,index)*jacobian;
- }
+ }
return SVector3(comp_x,comp_y,0.0);
}
@@ -1462,7 +1466,7 @@ double lpcvt::f(SPoint2 p1,SPoint2 p2,int p,int index){
double c;
double d;
metric m;
-
+
x1 = p1.x();
y1 = p1.y();
x2 = p2.x();
@@ -1491,7 +1495,7 @@ double lpcvt::df_dx(SPoint2 p1,SPoint2 p2,int p,int index){
double c;
double d;
metric m;
-
+
x1 = p1.x();
y1 = p1.y();
x2 = p2.x();
@@ -1520,7 +1524,7 @@ double lpcvt::df_dy(SPoint2 p1,SPoint2 p2,int p,int index){
double c;
double d;
metric m;
-
+
x1 = p1.x();
y1 = p1.y();
x2 = p2.x();
@@ -1562,7 +1566,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();
@@ -1581,9 +1585,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();
@@ -1594,7 +1598,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);
}
diff --git a/Mesh/meshGFaceLloyd.h b/Mesh/meshGFaceLloyd.h
index 8499a317ea..5487c9287a 100644
--- a/Mesh/meshGFaceLloyd.h
+++ b/Mesh/meshGFaceLloyd.h
@@ -26,14 +26,13 @@ class metric;
void callback(const alglib::real_1d_array&,double&,alglib::real_1d_array&,void*);
bool domain_search(MElementOctree*,GFace*,double,double);
-class lloydAlgorithm {
+class smoothing{
int ITER_MAX;
int NORM;
public :
- lloydAlgorithm (int itermax, int norm) : ITER_MAX(itermax), NORM(norm) {}
- lloydAlgorithm() {}
- void operator () (GFace *);
- void optimize(int,int);
+ smoothing(int,int);
+ void optimize_face(GFace*);
+ void optimize_model();
};
class lpcvt{
diff --git a/benchmarks/stl/skull.geo b/benchmarks/stl/skull.geo
index 636cb45205..08a735d84c 100644
--- a/benchmarks/stl/skull.geo
+++ b/benchmarks/stl/skull.geo
@@ -1,4 +1,4 @@
-Mesh.RemeshParametrization=1; //(0) harmonic (1) conformal
+Mesh.RemeshParametrization=0; //(0) harmonic (1) conformal
Mesh.RemeshAlgorithm=1; //(0) nosplit (1) automatic (2) split metis
Mesh.CharacteristicLengthFactor=0.1;
--
GitLab