From fe9615f9dbb106a1e31b79a384871779e17f1c39 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Fri, 15 Jun 2012 15:03:00 +0000
Subject: [PATCH] cleanup
---
Fltk/statisticsWindow.cpp | 6 +-
Mesh/BackgroundMesh.cpp | 272 ++++++++++++++++++--------------------
Mesh/meshGRegion.cpp | 9 +-
Solver/laplaceTerm.h | 10 +-
doc/VERSIONS.txt | 2 +-
5 files changed, 142 insertions(+), 157 deletions(-)
diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index 2faebfdb64..0c95721008 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -146,9 +146,9 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
for(int i = 0; i < 4; i++){
int ww = 3 * FL_NORMAL_SIZE;
new Fl_Box
- (FL_NO_BOX, width - 3 * ww - 2 * WB, 2 * WB + (13 + i) * BH, ww, BH, "Plot:");
+ (FL_NO_BOX, width - 3 * ww - 2 * WB, 2 * WB + (13 + i) * BH, ww, BH, "Plot");
butt[2 * i] = new Fl_Button
- (width - 2 * ww - 2 * WB, 2 * WB + (13 + i) * BH, ww, BH, "2D");
+ (width - 2 * ww - 2 * WB, 2 * WB + (13 + i) * BH, ww, BH, "X-Y");
butt[2 * i + 1] = new Fl_Button
(width - ww - 2 * WB, 2 * WB + (13 + i) * BH, ww, BH, "3D");
}
@@ -234,7 +234,7 @@ void statisticsWindow::compute(bool elementQuality)
std::map<MVertex*, int > vert2Deg;
for(unsigned int i = 0; i < entities.size(); i++){
if(entities[i]->dim() < 2 ) continue;
- // if(entities[i]->tag() < 100) continue;
+ // if(entities[i]->tag() < 100) continue;
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
MElement *e = entities[i]->getMeshElement(j);
for(unsigned int k = 0; k < e->getNumEdges(); k++){
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 6a1d7ddc56..99c175f2cf 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -35,16 +35,16 @@
// CTX::instance()->mesh.minCircPoints tells the minimum number of points per
// radius of curvature
-SMetric3 buildMetricTangentToCurve (SVector3 &t, double l_t, double l_n)
+SMetric3 buildMetricTangentToCurve(SVector3 &t, double l_t, double l_n)
{
- if (l_t == 0.0)return SMetric3(1.e-22);
+ if (l_t == 0.0) return SMetric3(1.e-22);
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);
}
@@ -57,7 +57,8 @@ SMetric3 buildMetricTangentToCurve (SVector3 &t, double l_t, double l_n)
return Metric;
}
-SMetric3 buildMetricTangentToSurface (SVector3 &t1, SVector3 &t2, double l_t1, double l_t2, double l_n)
+SMetric3 buildMetricTangentToSurface(SVector3 &t1, SVector3 &t2,
+ double l_t1, double l_t2, double l_n)
{
t1.normalize();
t2.normalize();
@@ -72,29 +73,28 @@ SMetric3 buildMetricTangentToSurface (SVector3 &t1, SVector3 &t2, double l_t1, d
return Metric;
}
-
SMetric3 max_edge_curvature_metric(const GVertex *gv)
{
SMetric3 val (1.e-12);
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;
if (gv == _myGEdge->getBeginVertex()) {
SVector3 t = _myGEdge->firstDer(range.low());
t.normalize();
- double l_t = ((2 * M_PI) /( fabs(_myGEdge->curvature(range.low()))
- * CTX::instance()->mesh.minCircPoints ));
+ double l_t = ((2 * M_PI) /( fabs(_myGEdge->curvature(range.low()))
+ * CTX::instance()->mesh.minCircPoints ));
double l_n = 1.e12;
cc = buildMetricTangentToCurve(t,l_t,l_n);
}
else {
SVector3 t = _myGEdge->firstDer(range.high());
t.normalize();
- double l_t = ((2 * M_PI) /( fabs(_myGEdge->curvature(range.high()))
- * CTX::instance()->mesh.minCircPoints ));
+ double l_t = ((2 * M_PI) /( fabs(_myGEdge->curvature(range.high()))
+ * CTX::instance()->mesh.minCircPoints ));
double l_n = 1.e12;
cc = buildMetricTangentToCurve(t,l_t,l_n);
}
@@ -107,20 +107,20 @@ SMetric3 max_edge_curvature_metric(const GEdge *ge, double u)
{
SVector3 t = ge->firstDer(u);
t.normalize();
- double l_t = ((2 * M_PI) /( fabs(ge->curvature(u))
- * CTX::instance()->mesh.minCircPoints ));
+ double l_t = ((2 * M_PI) /( fabs(ge->curvature(u))
+ * CTX::instance()->mesh.minCircPoints ));
double l_n = 1.e12;
- return buildMetricTangentToCurve(t,l_t,l_n);
+ return buildMetricTangentToCurve(t,l_t,l_n);
}
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());
@@ -146,11 +146,12 @@ static double max_surf_curvature(const GEdge *ge, double u)
return val;
}
+/*
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);
@@ -161,11 +162,12 @@ static double max_surf_curvature(const GVertex *gv)
}
return val;
}
+*/
-SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, double v,
+SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, double v,
bool surface_isotropic,
double d_normal ,
- double d_tangent_max )
+ double d_tangent_max)
{
if (gf->geomType() == GEntity::Plane)return SMetric3(1.e-12);
double cmax, cmin;
@@ -190,7 +192,7 @@ SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, double v,
lambda2 = std::min(lambda2, d_tangent_max);
SMetric3 curvMetric (1./(lambda1*lambda1),1./(lambda2*lambda2),
- 1./(lambda3*lambda3),
+ 1./(lambda3*lambda3),
dirMin, dirMax, Z );
return curvMetric;
}
@@ -198,45 +200,41 @@ SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, double v,
static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u, bool iso_surf)
{
const GEdgeCompound* ptrCompoundEdge = dynamic_cast<const GEdgeCompound*>(ge);
- if (ptrCompoundEdge)
- {
- double cmax, cmin;
- SVector3 dirMax,dirMin;
- cmax = ptrCompoundEdge->curvatures(u,&dirMax, &dirMin, &cmax,&cmin);
- if (cmin == 0)cmin =1.e-5;
- if (cmax == 0)cmax =1.e-5;
- double lambda2 = ((2 * M_PI) /( fabs(cmax) * CTX::instance()->mesh.minCircPoints ) );
- double lambda1 = ((2 * M_PI) /( fabs(cmin) * CTX::instance()->mesh.minCircPoints ) );
- SVector3 Z = crossprod(dirMax,dirMin);
-
- lambda1 = std::max(lambda1, CTX::instance()->mesh.lcMin);
- lambda2 = std::max(lambda2, CTX::instance()->mesh.lcMin);
- lambda1 = std::min(lambda1, CTX::instance()->mesh.lcMax);
- lambda2 = std::min(lambda2, CTX::instance()->mesh.lcMax);
-
- SMetric3 curvMetric (1./(lambda1*lambda1),1./(lambda2*lambda2),1.e-5,dirMin, dirMax, Z );
- return curvMetric;
+ if (ptrCompoundEdge){
+ double cmax, cmin;
+ SVector3 dirMax,dirMin;
+ cmax = ptrCompoundEdge->curvatures(u,&dirMax, &dirMin, &cmax,&cmin);
+ if (cmin == 0)cmin =1.e-5;
+ if (cmax == 0)cmax =1.e-5;
+ double lambda2 = ((2 * M_PI) /( fabs(cmax) * CTX::instance()->mesh.minCircPoints ) );
+ double lambda1 = ((2 * M_PI) /( fabs(cmin) * CTX::instance()->mesh.minCircPoints ) );
+ SVector3 Z = crossprod(dirMax,dirMin);
+
+ lambda1 = std::max(lambda1, CTX::instance()->mesh.lcMin);
+ lambda2 = std::max(lambda2, CTX::instance()->mesh.lcMin);
+ lambda1 = std::min(lambda1, CTX::instance()->mesh.lcMax);
+ lambda2 = std::min(lambda2, CTX::instance()->mesh.lcMax);
+
+ SMetric3 curvMetric (1. / (lambda1 * lambda1), 1. / (lambda2 * lambda2),
+ 1.e-5, dirMin, dirMax, Z);
+ return curvMetric;
}
-
- else
- {
+ else{
SMetric3 mesh_size(1.e-12);
std::list<GFace *> faces = ge->faces();
std::list<GFace *>::iterator it = faces.begin();
// we choose the metric eigenvectors to be the ones
// related to the edge ...
SMetric3 curvMetric = max_edge_curvature_metric(ge, u);
- while(it != faces.end())
- {
- if (((*it)->geomType() != GEntity::CompoundSurface) &&
- ((*it)->geomType() != GEntity::DiscreteSurface)){
- SPoint2 par = ge->reparamOnFace((*it), u, 1);
- SMetric3 m = metric_based_on_surface_curvature (*it, par.x(), par.y(), iso_surf);
- curvMetric = intersection_conserveM1(curvMetric,m);
- }
- ++it;
+ while(it != faces.end()){
+ if (((*it)->geomType() != GEntity::CompoundSurface) &&
+ ((*it)->geomType() != GEntity::DiscreteSurface)){
+ SPoint2 par = ge->reparamOnFace((*it), u, 1);
+ SMetric3 m = metric_based_on_surface_curvature (*it, par.x(), par.y(), iso_surf);
+ curvMetric = intersection_conserveM1(curvMetric,m);
}
-
+ ++it;
+ }
return curvMetric;
}
}
@@ -245,7 +243,7 @@ static SMetric3 metric_based_on_surface_curvature(const GVertex *gv, bool iso_su
{
SMetric3 mesh_size(1.e-15);
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);
@@ -254,18 +252,16 @@ static SMetric3 metric_based_on_surface_curvature(const GVertex *gv, bool iso_su
// This is because we want to call the function
// metric_based_on_surface_curvature(const GEdge *ge, double u) for the case when
// ge is a compound edge
- if (_myGEdge->geomType() == GEntity::CompoundCurve)
- {
- if (gv == _myGEdge->getBeginVertex())
- mesh_size = intersection
- (mesh_size,
- metric_based_on_surface_curvature(_myGEdge, bounds.low(), iso_surf));
- else
- mesh_size = intersection
- (mesh_size,
- metric_based_on_surface_curvature(_myGEdge, bounds.high(), iso_surf));
+ if (_myGEdge->geomType() == GEntity::CompoundCurve){
+ if (gv == _myGEdge->getBeginVertex())
+ mesh_size = intersection
+ (mesh_size,
+ metric_based_on_surface_curvature(_myGEdge, bounds.low(), iso_surf));
+ else
+ mesh_size = intersection
+ (mesh_size,
+ metric_based_on_surface_curvature(_myGEdge, bounds.high(), iso_surf));
}
-
}
return mesh_size;
}
@@ -282,15 +278,15 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
switch(ge->dim()){
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);
+ // Crv = std::max(max_surf_curvature((const GVertex *)ge), Crv);
+ // Crv = max_surf_curvature((const GVertex *)ge);
break;
case 1:
{
GEdge *ged = (GEdge *)ge;
Crv = ged->curvature(U);
Crv = std::max(Crv, max_surf_curvature(ged, U));
- // Crv = max_surf_curvature(ged, U);
+ // Crv = max_surf_curvature(ged, U);
}
break;
case 2:
@@ -337,15 +333,15 @@ static double LC_MVertex_PNTS(GEntity *ge, double U, double V)
GVertex *v2 = ged->getEndVertex();
if (v1 && v2){
Range<double> range = ged->parBounds(0);
- double a = (U - range.low()) / (range.high() - range.low());
+ 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;
@@ -353,22 +349,22 @@ 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)
double l1 = CTX::instance()->lc;
-
+
// 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
double l3 = MAX_LC;
if(CTX::instance()->mesh.lcFromCurvature && ge->dim() < 3)
l3 = LC_MVertex_CURV(ge, U, V);
-
+
// lc from fields
double l4 = MAX_LC;
FieldManager *fields = ge->model()->getFields();
@@ -376,12 +372,12 @@ double BGM_MeshSize(GEntity *ge, double U, double V,
Field *f = fields->get(fields->getBackgroundField());
if(f) l4 = (*f)(X, Y, Z, ge);
}
-
+
// take the minimum, then constrain by lcMin and lcMax
double lc = std::min(std::min(std::min(l1, l2), l3), l4);
lc = std::max(lc, CTX::instance()->mesh.lcMin);
lc = std::min(lc, CTX::instance()->mesh.lcMax);
-
+
if(lc <= 0.){
Msg::Error("Wrong mesh element size lc = %g (lcmin = %g, lcmax = %g)",
lc, CTX::instance()->mesh.lcMin, CTX::instance()->mesh.lcMax);
@@ -398,19 +394,19 @@ double BGM_MeshSize(GEntity *ge, double U, double V,
// anisotropic version of the background field - 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;
const double max_lc = CTX::instance()->mesh.lcMax;
- // 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)
@@ -424,11 +420,6 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
if(f){
if (!f->isotropic()){
(*f)(X, Y, Z, l4,ge);
- /*
- if (ge->tag() == 3){
- printf("X = %12.5E , l4 = %12.5E %12.5E %12.5E l2 = %12.5E\n",X,l4(0,0),l4(1,1),l4(0,1),l2);
- }
- */
}
else{
double L = (*f)(X, Y, Z, ge);
@@ -436,7 +427,7 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
}
}
}
-
+
// take the minimum, then constrain by lcMin and lcMax
double lc = std::min(l1, l2);
lc = std::max(lc, CTX::instance()->mesh.lcMin);
@@ -447,7 +438,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_conserveM1(intersection_conserveM1 (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));
@@ -460,7 +451,7 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
l3.eig(V,S,true);
printf("%g %g %g\n",S(0),S(1),S(2));
l4.eig(V,S,true);
- printf("%g %g %g\n",S(0),S(1),S(2));
+ printf("%g %g %g\n",S(0),S(1),S(2));
}
}
@@ -507,7 +498,7 @@ backgroundMesh::backgroundMesh(GFace *_gf)
MVertex *newv =0;
if (it == _3Dto2D.end()){
SPoint2 p;
- bool success = reparamMeshVertexOnFace(v, _gf, p);
+ reparamMeshVertexOnFace(v, _gf, p);
newv = new MVertex (p.x(), p.y(), 0.0);
_vertices.push_back(newv);
_3Dto2D[v] = newv;
@@ -519,7 +510,7 @@ backgroundMesh::backgroundMesh(GFace *_gf)
}
_triangles.push_back(new MTriangle(news[0],news[1],news[2]));
}
-
+
#if defined(HAVE_ANN)
//printf("creating uv kdtree %d \n", myBCNodes.size());
index = new ANNidx[2];
@@ -539,8 +530,8 @@ backgroundMesh::backgroundMesh(GFace *_gf)
#endif
// build a search structure
- _octree = new MElementOctree(_triangles);
-
+ _octree = new MElementOctree(_triangles);
+
// compute the mesh sizes at nodes
if (CTX::instance()->mesh.lcFromPoints)
propagate1dMesh(_gf);
@@ -552,10 +543,10 @@ backgroundMesh::backgroundMesh(GFace *_gf)
}
// ensure that other criteria are fullfilled
updateSizes(_gf);
-
+
// compute optimal mesh orientations
propagatecrossField(_gf);
-
+
_3Dto2D.clear();
_2Dto3D.clear();
}
@@ -585,21 +576,21 @@ 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>;
#endif
-
+
dofManager<double> myAssembler(_lsys);
-
+
// fix boundary conditions
std::map<MVertex*, double>::iterator itv = dirichlet.begin();
for ( ; itv != dirichlet.end(); ++itv){
myAssembler.fixVertex(itv->first, 0, 1, itv->second);
}
-
-
+
+
// Number vertices
std::set<MVertex*> vs;
for (unsigned int k = 0; k < _gf->triangles.size(); k++)
@@ -614,12 +605,12 @@ static void propagateValuesOnFace(GFace *_gf,
reparamMeshVertexOnFace ( *it, _gf, p);
theMap[*it] = SPoint3((*it)->x(),(*it)->y(),(*it)->z());
(*it)->setXYZ(p.x(),p.y(),0.0);
- }
+ }
}
for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it)
myAssembler.numberVertex(*it, 0, 1);
-
+
// Assemble
simpleFunction<double> ONE(1.0);
laplaceTerm l(0, 1, &ONE);
@@ -628,10 +619,10 @@ static void propagateValuesOnFace(GFace *_gf,
SElement se(t);
l.addToMatrix(myAssembler, &se);
}
-
+
// Solve
_lsys->systemSolve();
-
+
// save solution
for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
myAssembler.getDofValue(*it, 0, 1, dirichlet[*it]);
@@ -653,7 +644,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++ ){
@@ -668,15 +659,15 @@ 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));
}
}
}
}
-
+
propagateValuesOnFace(_gf, sizes);
-
+
std::map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
for ( ; itv2 != _2Dto3D.end(); ++itv2){
MVertex *v_2D = itv2->first;
@@ -688,7 +679,7 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
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;
@@ -703,20 +694,20 @@ crossField2d::crossField2d(MVertex* v, GEdge* ge)
void backgroundMesh::propagatecrossField(GFace *_gf)
{
std::map<MVertex*,double> _cosines4,_sines4;
-
+
std::list<GEdge*> e;
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;
- bool success = reparamMeshEdgeOnFace(v[0],v[1],_gf,p1,p2);
+ SPoint2 p1,p2;
+ 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() );
//double angle = atan2 ( v0->y()-v1->y() , v0->x()- v1->x() );
@@ -725,8 +716,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);
@@ -736,7 +727,7 @@ void backgroundMesh::propagatecrossField(GFace *_gf)
}
}
}
-
+
// force smooth transition
const int nbSmooth = 0;
const double threshold_angle = 2. * M_PI/180.;
@@ -767,10 +758,10 @@ void backgroundMesh::propagatecrossField(GFace *_gf)
}
}
}
-
+
propagateValuesOnFace(_gf,_cosines4,false);
propagateValuesOnFace(_gf,_sines4,false);
-
+
std::map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
for ( ; itv2 != _2Dto3D.end(); ++itv2){
MVertex *v_2D = itv2->first;
@@ -797,20 +788,20 @@ 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);
+ 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);
- }
+ }
// do not allow large variations in the size field
// (Int. J. Numer. Meth. Engng. 43, 1143-1165 (1998) MESH GRADATION
// CONTROL, BOROUCHAKI, HECHT, FREY)
std::set<MEdge,Less_Edge> edges;
- for (int i=0;i < _triangles.size();i++){
- for (int j=0;j< _triangles[i]->getNumEdges();j++){
+ for (unsigned int i = 0; i < _triangles.size(); i++){
+ for (int j = 0; j < _triangles[i]->getNumEdges(); j++){
edges.insert(_triangles[i]->getEdge(j));
}
}
@@ -838,16 +829,13 @@ double backgroundMesh::operator() (double u, double v, double w) const
if (!e) {
#if defined(HAVE_ANN)
//printf("BGM octree not found --> find in kdtree \n");
- double pt[3] = {u,v,0.0};
-
+ double pt[3] = {u, v, 0.0};
uv_kdtree->annkSearch(pt, 2, index, dist);
SPoint3 p1(nodes[index[0]][0], nodes[index[0]][1], nodes[index[0]][2]);
SPoint3 p2(nodes[index[1]][0], nodes[index[1]][1], nodes[index[1]][2]);
SPoint3 pnew; double d;
- signedDistancePointLine(p1,p2, SPoint3(u,v,0.0), d, pnew);
- double uvw[3]={pnew.x(),pnew.y(), 0.0};
- double UV[3];
- e = _octree->find(pnew.x(), pnew.y(), 0.0, 2, true);
+ signedDistancePointLine(p1, p2, SPoint3(u, v, 0.), d, pnew);
+ e = _octree->find(pnew.x(), pnew.y(), 0.0, 2, true);
#endif
if(!e){
Msg::Error("BGM octree: cannot find UVW=%g %g %g", u, v, w);
@@ -864,12 +852,12 @@ double backgroundMesh::operator() (double u, double v, double w) const
double backgroundMesh::getAngle(double u, double v, double w) const
{
- // HACK FOR LEWIS
+ // HACK FOR LEWIS
// h = 1+30(y-x^2)^2 + (1-x)^2
// double x = u;
// double y = v;
- // double dhdx = 30 * 2 * (y-x*x) * (-2*x) - 2 * (1-x);
- // double dhdy = 30 * 2 * (y-x*x);
+ // double dhdx = 30 * 2 * (y-x*x) * (-2*x) - 2 * (1-x);
+ // double dhdy = 30 * 2 * (y-x*x);
// double angles = atan2(y,x*x);
// crossField2d::normalizeAngle (angles);
// return angles;
@@ -885,10 +873,8 @@ double backgroundMesh::getAngle(double u, double v, double w) const
SPoint3 p1(nodes[index[0]][0], nodes[index[0]][1], nodes[index[0]][2]);
SPoint3 p2(nodes[index[1]][0], nodes[index[1]][1], nodes[index[1]][2]);
SPoint3 pnew; double d;
- signedDistancePointLine(p1,p2, SPoint3(u,v,0.0), d, pnew);
- double uvw[3]={pnew.x(),pnew.y(), 0.0};
- double UV[3];
- e = _octree->find(pnew.x(), pnew.y(), 0.0, 2,true);
+ signedDistancePointLine(p1, p2, SPoint3(u, v, 0.), d, pnew);
+ e = _octree->find(pnew.x(), pnew.y(), 0., 2, true);
#endif
if(!e){
Msg::Error("BGM octree angle: cannot find UVW=%g %g %g", u, v, w);
@@ -899,16 +885,16 @@ double backgroundMesh::getAngle(double u, double v, double w) const
std::map<MVertex*,double>::const_iterator itv1 = _angles.find(e->getVertex(0));
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);
-
+
return angle;
}
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 325e598694..64d68cb87b 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -71,12 +71,12 @@ void printVoronoi(GRegion *gr, std::set<SPoint3> &candidates)
//print voronoi poles
FILE *outfile;
- smooth_normals *snorm = gr->model()->normals;
+ //smooth_normals *snorm = gr->model()->normals;
outfile = fopen("nodes.pos", "w");
fprintf(outfile, "View \"Voronoi poles\" {\n");
std::map<MVertex*, std::set<MTetrahedron*> >::iterator itmap = node2Tet.begin();
for(; itmap != node2Tet.end(); itmap++){
- MVertex *v = itmap->first;
+ //MVertex *v = itmap->first;
std::set<MTetrahedron*> allTets = itmap->second;
std::set<MTetrahedron*>::const_iterator it = allTets.begin();
MTetrahedron *poleTet = *it;
@@ -98,7 +98,7 @@ void printVoronoi(GRegion *gr, std::set<SPoint3> &candidates)
double x[3] = {pc.x(), pc.y(), pc.z()};
double uvw[3];
poleTet->xyz2uvw(x, uvw);
- bool inside = poleTet->isInside(uvw[0], uvw[1], uvw[2]);
+ //bool inside = poleTet->isInside(uvw[0], uvw[1], uvw[2]);
//if (inside){
fprintf(outfile,"SP(%g,%g,%g) {%g};\n",
pc.x(), pc.y(), pc.z(), maxRadius);
@@ -520,13 +520,12 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
std::vector<MVertex*> numberedV;
char opts[128];
buildTetgenStructure(gr, in, numberedV);
- //if (Msg::GetVerbosity() == 20) sprintf(opts, "peVvS0");
if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_DEL ||
CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_HEX ||
CTX::instance()->mesh.algo3d == ALGO_3D_MMG3D ||
CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD ||
CTX::instance()->mesh.algo2d == ALGO_2D_BAMG){
- sprintf(opts, "-q1.5pY", (Msg::GetVerbosity() < 3) ? 'Q':
+ sprintf(opts, "-q1.5pY%c", (Msg::GetVerbosity() < 3) ? 'Q':
(Msg::GetVerbosity() > 6) ? 'V': '\0');
}
else {
diff --git a/Solver/laplaceTerm.h b/Solver/laplaceTerm.h
index 19db897d3e..e6e1891fdb 100644
--- a/Solver/laplaceTerm.h
+++ b/Solver/laplaceTerm.h
@@ -8,19 +8,19 @@
#include "helmholtzTerm.h"
-// \nabla \cdot k \nabla U
+// \nabla \cdot k \nabla U
class laplaceTerm : public helmholtzTerm<double> {
protected:
- std::map<MVertex*, SPoint3> *_coordView;
const int _iField;
+ std::map<MVertex*, SPoint3> *_coordView;
public:
- laplaceTerm(GModel *gm, int iField, simpleFunction<double> *k,
+ laplaceTerm(GModel *gm, int iField, simpleFunction<double> *k,
std::map<MVertex*, SPoint3> *coord=NULL)
: helmholtzTerm<double>(gm, iField, iField, k, 0), _iField(iField), _coordView(coord) {}
void elementVector(SElement *se, fullVector<double> &m) const
{
MElement *e = se->getMeshElement();
- int nbSF = e->getNumShapeFunctions();
+ int nbSF = e->getNumShapeFunctions();
fullMatrix<double> *mat;
mat = new fullMatrix<double>(nbSF, nbSF);
@@ -42,7 +42,7 @@ class laplaceTerm : public helmholtzTerm<double> {
}
};
-// a \nabla U
+// a \nabla U
class massTerm : public helmholtzTerm<double> {
public:
massTerm(GModel *gm, int iField, simpleFunction<double> *a)
diff --git a/doc/VERSIONS.txt b/doc/VERSIONS.txt
index ae30bebe87..ff6955882f 100644
--- a/doc/VERSIONS.txt
+++ b/doc/VERSIONS.txt
@@ -7,7 +7,7 @@ field export; new experimental stereo+camera visualization mode; experimental
BAMG & MMG3D support for anisotropic mesh generation; new OCC cut&merge
algorithm imported from Salome; new ability to connect extruded meshes to
tetrahedral grids using pyramids; new homology solver; Abaqus (INP) mesh export;
-various bug fixes and improvements.
+bug fixes and small improvements all over the place.
2.5.0 (Oct 15, 2010): new compound geometrical entities (for remeshing and/or
trans-patch meshing); improved mesh reclassification tool; new client/server
--
GitLab