From dd8bb86a3d88c4c0da1a75277f4e876771c40ae5 Mon Sep 17 00:00:00 2001
From: Gauthier Becker <gauthierbecker@gmail.com>
Date: Sun, 29 May 2011 20:21:24 +0000
Subject: [PATCH] mpi in progress (implementation is finished remain to debug)
---
Mesh/BackgroundMesh.cpp | 105 +++++++++++++--------------
Mesh/meshGFaceLloyd.cpp | 152 ++++++++++++++++++++--------------------
Solver/dofManager.h | 4 +-
3 files changed, 131 insertions(+), 130 deletions(-)
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index a38167eb9a..b01b7d529f 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -26,6 +26,7 @@
#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
@@ -39,10 +40,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);
}
@@ -63,10 +64,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()) {
@@ -88,17 +89,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());
@@ -119,7 +120,7 @@ static double max_surf_curvature(const GEdge *ge, double u)
val = std::max(cc, val);
}
++it;
- }
+ }
return val;
}
@@ -127,7 +128,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);
@@ -162,7 +163,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;
}
@@ -182,7 +183,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));
@@ -193,7 +194,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);
@@ -203,7 +204,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;
@@ -219,7 +220,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);
@@ -229,7 +230,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:
@@ -239,7 +240,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;
}
@@ -274,16 +275,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;
@@ -291,7 +292,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)
@@ -299,7 +300,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
@@ -335,18 +336,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)
@@ -366,7 +367,7 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
}
}
}
-
+
// take the minimum, then constrain by lcMin and lcMax
double lc = std::min(l1, l2);
@@ -383,7 +384,7 @@ 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));
@@ -432,7 +433,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;
@@ -445,7 +446,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)
@@ -456,7 +457,7 @@ backgroundMesh::backgroundMesh(GFace *_gf)
_sizes[itv2->first] = MAX_LC;
}
}
- // ensure that other criteria are fullfilled
+ // ensure that other criteria are fullfilled
// updateSizes(_gf);
// compute optimal mesh orientations
@@ -473,7 +474,7 @@ backgroundMesh::~backgroundMesh()
delete _octree;
}
-static void propagateValuesOnFace (GFace *_gf,
+static void propagateValuesOnFace (GFace *_gf,
std::map<MVertex*,double> &dirichlet)
{
linearSystem<double> *_lsys = 0;
@@ -483,7 +484,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>;
@@ -513,7 +514,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
@@ -523,7 +524,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;
}
@@ -532,7 +533,7 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
std::list<GEdge*> e = _gf->edges();
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++ ){
@@ -546,9 +547,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));
- }
+ }
}
}
}
@@ -563,11 +564,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;
@@ -587,14 +588,14 @@ void backgroundMesh::propagatecrossField(GFace *_gf)
std::map<MVertex*,double> _cosines4,_sines4;
std::list<GEdge*> e = _gf->edges();
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() );
@@ -603,8 +604,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);
@@ -633,7 +634,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;
@@ -646,14 +647,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
@@ -688,11 +689,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/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index 7dde29512e..06a238f157 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;
@@ -165,16 +165,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 +188,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 +215,8 @@ void lloydAlgorithm::operator () (GFace *gf)
if(flag) count++;
}
}
- triangulator.remove_all();
-
+ triangulator.remove_all();
+
triangulator.Voronoi();
double delta;
delta = 0.01;
@@ -227,7 +227,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 +249,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 +268,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 +288,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 +310,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 +336,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 +346,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 +359,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 +376,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,10 +386,10 @@ 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){
@@ -430,7 +430,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 +670,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 +697,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 +727,7 @@ void lpcvt::step2(DocRecord& triangulator,GFace* gf){
temp[i].add_vertex(vertex);
}
}
- }
+ }
}
void lpcvt::step3(DocRecord& triangulator,GFace* gf){
@@ -844,7 +844,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 +896,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 +951,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 +1037,7 @@ void lpcvt::print_voronoi1(){
p2 = v2.get_point();
p3 = v3.get_point();
print_segment(p2,p3,file);
- }
+ }
file << "};\n";
}
@@ -1061,7 +1061,7 @@ void lpcvt::print_voronoi2(){
}
file << "};\n";
}
-
+
void lpcvt::print_delaunay(DocRecord& triangulator){
int i;
int j;
@@ -1089,7 +1089,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,9 +1101,9 @@ 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();
@@ -1112,7 +1112,7 @@ void lpcvt::compute_metrics(DocRecord& triangulator){
cosinus = cos(angle);
sinus = sin(angle);
m = metric(cosinus,-sinus,sinus,cosinus);
- metrics[i] = m;
+ metrics[i] = m;
}
}
@@ -1121,7 +1121,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 +1148,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 +1190,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 +1223,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 +1251,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 +1379,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 +1412,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 +1445,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 +1462,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 +1491,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 +1520,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 +1562,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 +1581,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 +1594,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/Solver/dofManager.h b/Solver/dofManager.h
index 895a2dc079..22a39f7a15 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -496,7 +496,7 @@ class dofManager{
}
int sizeOfR() const { return _isParallel ? _localSize : unknown.size(); }
int sizeOfF() const { return fixed.size(); }
- void systemSolve(){ _current->systemSolve(); }
+ virtual void systemSolve(){ _current->systemSolve(); }
void systemClear()
{
_current->zeroMatrix();
@@ -513,7 +513,7 @@ class dofManager{
throw;
}
}
- virtual linearSystem<dataMat> *getLinearSystem(std::string &name)
+ linearSystem<dataMat> *getLinearSystem(std::string &name)
{
typename std::map<const std::string, linearSystem<dataMat>*>::iterator it =
_linearSystems.find(name);
--
GitLab