From 4417faea4ea00a0fccace25bb76c30ec8cda2f7f Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Fri, 23 Apr 2010 13:53:50 +0000
Subject: [PATCH] dg Bloodflow added bcs
---
Geo/GFaceCompound.cpp | 9 ++++++---
Geo/GFaceCompound.h | 1 +
Mesh/Generator.cpp | 21 ++++++++++++++++++---
Mesh/meshGFaceLloyd.cpp | 26 ++++++++++++++++----------
Numeric/DivideAndConquer.cpp | 15 +++++++++------
Numeric/DivideAndConquer.h | 4 +++-
Solver/eigenSolver.cpp | 2 +-
7 files changed, 54 insertions(+), 24 deletions(-)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 46bcf30383..0cbbb498d8 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -519,8 +519,8 @@ bool GFaceCompound::parametrize() const
else if (_mapping == CONFORMAL){
Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
fillNeumannBCS();
- bool withoutFolding = parametrize_conformal_spectral() ;
- //bool withoutFolding = parametrize_conformal();
+ //bool withoutFolding = parametrize_conformal_spectral() ;
+ bool withoutFolding = parametrize_conformal();
if ( withoutFolding == false ){
//printStuff(); exit(1);
Msg::Warning("$$$ Parametrization switched to harmonic map");
@@ -1822,7 +1822,10 @@ bool GFaceCompound::checkTopology() const
// FIXME!!! I think those things are wrong with cross-patch reparametrization
//if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)return true;
+
+ //TODO: smthg to exit here for lloyd remeshing
+
bool correctTopo = true;
if(allNodes.empty()) buildAllNodes();
@@ -1913,7 +1916,7 @@ double GFaceCompound::checkAspectRatio() const
}
double AR = 2*3.14*area3D/(tot_length*tot_length);
- if (areaMin < limit && nb > 2) {
+ if (areaMin > 0 && areaMin < limit && nb > 2) {
Msg::Warning("Too small triangles in mapping (a_2D=%g)", areaMin);
}
else {
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index eeab87d002..b92b76d8e8 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -121,6 +121,7 @@ class GFaceCompound : public GFace {
void partitionFaceCM();
virtual std::list<GFace*> getCompounds() const {return _compound;};
mutable int nbSplit;
+ int getNbSplit() const {return nbSplit;}
int allowPartition() const{return _allowPartition;};
private:
typeOfIsomorphism _type;
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 6aa6d5d697..f9759546a9 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -435,12 +435,27 @@ static void Mesh2D(GModel *m)
// lloyd optimization
if (CTX::instance()->mesh.optimizeLloyd){
+ Msg::Info("------------------------------");
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
GFace *gf = *it;
+ if(gf->geomType() == GEntity::DiscreteSurface) continue;
+ if(gf->geomType() == GEntity::CompoundSurface ) {
+ GFaceCompound *gfc = (GFaceCompound*) gf;
+ if (gfc->getNbSplit() != 0) continue;
+ }
int recombine = gf->meshAttributes.recombine;
- Msg::Info("Lloyd optimization for face %d", gf->tag());
- gf->lloyd(40, recombine);
- if(recombine) recombineIntoQuads(gf);
+
+ Msg::StatusBar(2, true, "Lloyd optimization for face %d", gf->tag());
+ gf->lloyd(25,recombine);
+
+ if(recombine) recombineIntoQuads(gf);
+
+ // computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
+// gf->meshStatistics.average_element_shape,
+// gf->meshStatistics.best_element_shape,
+// gf->meshStatistics.nbTriangle,
+// gf->meshStatistics.nbGoodQuality);
+
}
}
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index d7e2546955..7f4b30d933 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -68,6 +68,7 @@ void lloydAlgorithm::operator () ( GFace * gf) {
fullMatrix<double> cgs(triangulator.numPoints,2);
// now iterate on internal vertices
double ENERGY = 0.0;
+ double criteria = 0.0;
for (int i=0; i<triangulator.numPoints;i++){
// get the ith vertex
PointRecord &pt = triangulator.points[i];
@@ -76,26 +77,23 @@ void lloydAlgorithm::operator () ( GFace * gf) {
// get the voronoi corners
std::vector<SPoint2> pts;
triangulator.voronoiCell (i,pts);
- double E;
+ double E, A;
SPoint2 p(pt.where.h,pt.where.v);
if (!infiniteNorm){
- centroidOfPolygon (p,pts, cgs(i,0),cgs(i,1),E,backgroundMesh::current());
+ centroidOfPolygon (p,pts, cgs(i,0),cgs(i,1),E, A, NULL); //backgroundMesh::current());
}
else {
- centroidOfOrientedBox (pts, 0.0, cgs(i,0),cgs(i,1),E);
+ centroidOfOrientedBox (pts, 0.0, cgs(i,0),cgs(i,1),E, A);
}
ENERGY += E;
+ double d = sqrt((p.x()-cgs(i,0))*(p.x()-cgs(i,0))+
+ (p.y()-cgs(i,1))*(p.y()-cgs(i,1)));
+ criteria = d/A;
}// if (v->onWhat() == gf)
else {
}
}// for all points
- if (ITER % 10 == 0){
- char name[234];
- sprintf(name,"LloydIter%d.pos",ITER);
- triangulator.makePosView(name);
- }
-
for(PointNumero i = 0; i < triangulator.numPoints; i++) {
MVertex *v = (MVertex*)triangulator.points[i].data;
if (v->onWhat() == gf && !triangulator.onHull(i)){
@@ -104,11 +102,18 @@ void lloydAlgorithm::operator () ( GFace * gf) {
}
}
- Msg::Debug("GFace %d Lloyd iteration %d energy %g",gf->tag(),ITER++,ENERGY);
+ Msg::Debug("GFace %d Lloyd-iter %d Inertia=%g Convergence=%g ",gf->tag(),ITER++,ENERGY, criteria);
if (ITER > ITER_MAX)break;
// compute the Voronoi diagram
triangulator.Voronoi();
+
+ if (ITER % 10 == 0){
+ char name[234];
+ sprintf(name,"LloydIter%d.pos",ITER);
+ triangulator.makePosView(name);
+ }
+
}
// now create the vertices
@@ -132,6 +137,7 @@ void lloydAlgorithm::operator () ( GFace * gf) {
// vertices
gf->_additional_vertices = mesh_vertices;
// remesh the face with all those vertices in
+ Msg::Info("Lloyd remeshing of face %d ", gf->tag());
meshGFace mesher;
mesher(gf);
// assign those vertices to the face internal vertices
diff --git a/Numeric/DivideAndConquer.cpp b/Numeric/DivideAndConquer.cpp
index ddc3f1c8b9..6bd532988e 100644
--- a/Numeric/DivideAndConquer.cpp
+++ b/Numeric/DivideAndConquer.cpp
@@ -592,7 +592,7 @@ void DocRecord::makePosView(std::string fileName)
}
void centroidOfOrientedBox(std::vector<SPoint2> &pts, const double &angle,
- double &xc, double &yc, double &inertia)
+ double &xc, double &yc, double &inertia, double &area)
{
const int N = pts.size();
@@ -614,13 +614,15 @@ void centroidOfOrientedBox(std::vector<SPoint2> &pts, const double &angle,
xc = XC*cosa - YC *sina;
yc = XC*sina + YC *cosa;
inertia = std::max(xmax-xmin,ymax-ymin);
+ area = (xmax-xmin) * (ymax-ymin);
}
void centroidOfPolygon(SPoint2 &pc, std::vector<SPoint2> &pts,
- double &xc, double &yc, double &inertia,
+ double &xc, double &yc, double &inertia, double &areaCell,
simpleFunction<double> *bgm)
{
double area_tot = 0;
+ areaCell = 0.0;
SPoint2 center(0,0);
for (unsigned int j = 0; j < pts.size(); j++){
SPoint2 &pa = pts[j];
@@ -629,7 +631,8 @@ void centroidOfPolygon(SPoint2 &pc, std::vector<SPoint2> &pts,
const double lc = bgm ? (*bgm)((pa.x()+pb.x()+pc.x())/3.0,
(pa.y()+pb.y()+pc.y())/3.0,
0.0) : 1.0;
- const double fact = 1./(lc*lc);
+ const double fact = 1./(lc*lc*lc*lc);//rho = 1/lc^4 (emi)
+ areaCell += area;
area_tot += area*fact;
center += ((pa+pb+pc) * (area*fact/3.0));
}
@@ -664,14 +667,14 @@ double DocRecord::Lloyd(int type)
PointRecord &pt = points[i];
std::vector<SPoint2> pts;
voronoiCell (i,pts);
- double E;
+ double E, A;
if (!points[i].data){
SPoint2 p (pt.where.h,pt.where.v);
if (type == 0)
- centroidOfPolygon (p,pts, cgs(i,0), cgs(i,1),E);
+ centroidOfPolygon (p,pts, cgs(i,0), cgs(i,1),E, A);
else
- centroidOfOrientedBox (pts, 0.0, cgs(i,0),cgs(i,1),E);
+ centroidOfOrientedBox (pts, 0.0, cgs(i,0),cgs(i,1),E, A);
}
inertia_tot += E;
}
diff --git a/Numeric/DivideAndConquer.h b/Numeric/DivideAndConquer.h
index e44784dcf9..b8093c2787 100644
--- a/Numeric/DivideAndConquer.h
+++ b/Numeric/DivideAndConquer.h
@@ -104,12 +104,14 @@ void centroidOfOrientedBox(std::vector<SPoint2> &pts,
const double &angle,
double &xc,
double &yc,
- double &inertia);
+ double &inertia,
+ double &area);
void centroidOfPolygon(SPoint2 &pc,
std::vector<SPoint2> &pts,
double &xc,
double &yc,
double &inertia,
+ double &areaCell,
simpleFunction<double> *bgm = 0);
diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp
index 353665900f..4e58cbad6c 100644
--- a/Solver/eigenSolver.cpp
+++ b/Solver/eigenSolver.cpp
@@ -53,7 +53,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
// set some default options
_try(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
- _try(EPSSetTolerances(eps,1.e-7,50));
+ _try(EPSSetTolerances(eps,1.e-6,20));//1.e-7 50
//_try(EPSSetType(eps, EPSKRYLOVSCHUR)); //default
_try(EPSSetType(eps,EPSARNOLDI));
//_try(EPSSetType(eps,EPSARPACK));
--
GitLab