diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index ba854148177145f6a3f998824000fd5e5ba73ee9..876f3a7e301be420d82d04f6064ea5cd757978f2 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -55,7 +55,7 @@ SMetric3 buildMetricTangentToCurve (SVector3 &t, double curvature, double &lambd
c.normalize();
t.normalize();
lambda = ((2 * M_PI) /( fabs(curvature) * CTX::instance()->mesh.minCircPoints ) );
-
+
SMetric3 curvMetric (1./(lambda*lambda),1.e-10,1.e-10,t,b,c);
//SMetric3 curvMetric (1./(lambda*lambda));
return curvMetric;
@@ -118,13 +118,13 @@ static double max_surf_curvature(const GEdge *ge, double u)
std::list<GFace *>::iterator it = faces.begin();
while(it != faces.end()){
if ((*it)->geomType() != GEntity::CompoundSurface &&
- (*it)->geomType() != GEntity::DiscreteSurface){
+ (*it)->geomType() != GEntity::DiscreteSurface){
SPoint2 par = ge->reparamOnFace((*it), u, 1);
double cc = (*it)->curvature(par);
val = std::max(cc, val);
}
++it;
- }
+ }
return val;
}
@@ -162,13 +162,13 @@ static SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, dou
lambda2 = std::min(lambda2, CTX::instance()->mesh.lcMax);
/* if (gf->tag() == 36 || gf->tag() == 62)
printf("%g %g -- %g %g -- %g %g %g -- %g %g %g -- %g %g %g -- %g %g\n",u,v,cmax,cmin,
- dirMax.x(),dirMax.y(),dirMax.z(),
- dirMin.x(),dirMin.y(),dirMin.z(),
- Z.x(),Z.y(),Z.z(),
- lambda1,lambda2);
+ dirMax.x(),dirMax.y(),dirMax.z(),
+ dirMin.x(),dirMin.y(),dirMin.z(),
+ Z.x(),Z.y(),Z.z(),
+ lambda1,lambda2);
*/
SMetric3 curvMetric (1./(lambda1*lambda1),1./(lambda2*lambda2),1.e-5,
- dirMin, dirMax, Z );
+ dirMin, dirMax, Z );
return curvMetric;
}
@@ -188,7 +188,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));
@@ -223,10 +223,9 @@ static SMetric3 metric_based_on_surface_curvature(const GVertex *gv)
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);
@@ -246,7 +245,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;
}
@@ -284,7 +283,7 @@ 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);
+ Range<double> range = ged->parBounds(0);
double a = (U - range.low()) / (range.high() - range.low());
double lc = (1 - a) * v1->prescribedMeshSizeAtVertex() +
(a) * v2->prescribedMeshSizeAtVertex() ;
@@ -306,7 +305,7 @@ double BGM_MeshSize(GEntity *ge, double U, double V,
{
// 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)
@@ -324,12 +323,12 @@ double BGM_MeshSize(GEntity *ge, double U, double V,
Field *f = fields->get(fields->background_field);
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);
@@ -371,14 +370,14 @@ 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);
- }
- */
+ /*
+ 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);
+ double L = (*f)(X, Y, Z, ge);
l4 = SMetric3(1/(L*L));
}
}
@@ -438,7 +437,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;
@@ -449,10 +448,10 @@ backgroundMesh::backgroundMesh(GFace *_gf)
}
_triangles.push_back(new MTriangle(news[0],news[1],news[2]));
}
-
+
// build a search structure
_octree = new MElementOctree(_triangles);
-
+
// compute the mesh sizes at nodes
if (CTX::instance()->mesh.lcFromPoints)
propagate1dMesh(_gf);
@@ -464,10 +463,10 @@ backgroundMesh::backgroundMesh(GFace *_gf)
}
// ensure that other criteria are fullfilled
updateSizes(_gf);
-
+
// compute optimal mesh orientations
propagatecrossField(_gf);
-
+
_3Dto2D.clear();
_2Dto3D.clear();
}
@@ -479,7 +478,7 @@ backgroundMesh::~backgroundMesh()
delete _octree;
}
-static void propagateValuesOnFace(GFace *_gf,
+static void propagateValuesOnFace(GFace *_gf,
std::map<MVertex*,double> &dirichlet)
{
#if defined(HAVE_SOLVER)
@@ -495,16 +494,16 @@ static void propagateValuesOnFace(GFace *_gf,
#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++)
@@ -513,19 +512,19 @@ static void propagateValuesOnFace(GFace *_gf,
for (int j=0;j<4;j++)vs.insert(_gf->quadrangles[k]->getVertex(j));
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);
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
_lsys->systemSolve();
-
+
// save solution
for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
myAssembler.getDofValue(*it, 0, 1, dirichlet[*it]);
@@ -547,7 +546,7 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
for(unsigned int i = 0; i < (*it)->lines.size(); i++ ){
MVertex *v1 = (*it)->lines[i]->getVertex(0);
MVertex *v2 = (*it)->lines[i]->getVertex(1);
- // printf("%g %g %g\n",v1->x(),v1->y(),v1->z());
+ // printf("%g %g %g\n",v1->x(),v1->y(),v1->z());
double d = sqrt((v1->x() - v2->x()) * (v1->x() - v2->x()) +
(v1->y() - v2->y()) * (v1->y() - v2->y()) +
(v1->z() - v2->z()) * (v1->z() -v2->z()));
@@ -558,13 +557,13 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
sizes[v] = log(d);
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;
@@ -591,40 +590,40 @@ 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];
+ 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);
- 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() );
- crossField2d::normalizeAngle (angle);
- for (int i=0;i<2;i++){
- 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));
- }
- else {
- _cosines4[v[i]] = cos(4*angle);
- _sines4[v[i]] = sin(4*angle);
- }
- }
+ 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() );
+ //double angle = atan2 ( v0->y()-v1->y() , v0->x()- v1->x() );
+ crossField2d::normalizeAngle (angle);
+ for (int i=0;i<2;i++){
+ 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));
+ }
+ else {
+ _cosines4[v[i]] = cos(4*angle);
+ _sines4[v[i]] = sin(4*angle);
+ }
+ }
}
}
}
-
+
// force smooth transition
const int nbSmooth = 0;
const double threshold_angle = 2. * M_PI/180.;
@@ -632,33 +631,33 @@ void backgroundMesh::propagatecrossField(GFace *_gf)
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);
- double cos40 = _cosines4[v[0]];
- double cos41 = _cosines4[v[1]];
- double sin40 = _sines4[v[0]];
- double sin41 = _sines4[v[1]];
- double angle0 = atan2 (sin40,cos40)/4.;
- double angle1 = atan2 (sin41,cos41)/4.;
- if (fabs(angle0 - angle1) > threshold_angle ){
- double angle0_new = angle0 - (angle0-angle1) * 0.1;
- double angle1_new = angle1 + (angle0-angle1) * 0.1;
- // printf("%g %g -- %g %g\n",angle0,angle1,angle0_new,angle1_new);
- _cosines4[v[0]] = cos(4*angle0_new);
- _sines4[v[0]] = sin(4*angle0_new);
- _cosines4[v[1]] = cos(4*angle1_new);
- _sines4[v[1]] = sin(4*angle1_new);
- }
- }
+ 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);
+ double cos40 = _cosines4[v[0]];
+ double cos41 = _cosines4[v[1]];
+ double sin40 = _sines4[v[0]];
+ double sin41 = _sines4[v[1]];
+ double angle0 = atan2 (sin40,cos40)/4.;
+ double angle1 = atan2 (sin41,cos41)/4.;
+ if (fabs(angle0 - angle1) > threshold_angle ){
+ double angle0_new = angle0 - (angle0-angle1) * 0.1;
+ double angle1_new = angle1 + (angle0-angle1) * 0.1;
+ // printf("%g %g -- %g %g\n",angle0,angle1,angle0_new,angle1_new);
+ _cosines4[v[0]] = cos(4*angle0_new);
+ _sines4[v[0]] = sin(4*angle0_new);
+ _cosines4[v[1]] = cos(4*angle1_new);
+ _sines4[v[1]] = sin(4*angle1_new);
+ }
+ }
}
}
}
-
+
propagateValuesOnFace(_gf,_cosines4);
propagateValuesOnFace(_gf,_sines4);
-
+
std::map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
for ( ; itv2 != _2Dto3D.end(); ++itv2){
MVertex *v_2D = itv2->first;
@@ -667,13 +666,12 @@ void backgroundMesh::propagatecrossField(GFace *_gf)
crossField2d::normalizeAngle (angle);
_angles[v_2D] = angle;
}
-
}
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;
@@ -686,7 +684,7 @@ 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());
@@ -711,53 +709,53 @@ void backgroundMesh::updateSizes(GFace *_gf)
MVertex *v1 = it->getVertex(1);
MVertex *V0 = _2Dto3D[v0];
MVertex *V1 = _2Dto3D[v1];
- std::map<MVertex*,double>::iterator s0 = _sizes.find(V0);
- std::map<MVertex*,double>::iterator s1 = _sizes.find(V1);
+ std::map<MVertex*,double>::iterator s0 = _sizes.find(V0);
+ std::map<MVertex*,double>::iterator s1 = _sizes.find(V1);
if (s0->second < s1->second)s1->second = std::min(s1->second,_beta*s0->second);
else s0->second = std::min(s0->second,_beta*s1->second);
}
}
}
-double backgroundMesh::operator() (double x, double y, double z) const
+double backgroundMesh::operator() (double u, double v, double w) const
{
- double xyz[3] = {x, y, z};
- double uv[3];
- MElement *e = _octree->find(x, y, z, 2, true);
+ double uv[3] = {u, v, w};
+ double uv2[3];
+ MElement *e = _octree->find(u, v, w, 2, true);
if (!e) {
- Msg::Error("cannot find %g %g %g", x, y, z);
- return -1000.0;//0.4;
+ Msg::Error("cannot find %g %g %g", u, v, w);
+ return -1000.0;//0.4;
}
- e->xyz2uvw(xyz, uv);
+ e->xyz2uvw(uv, uv2);
std::map<MVertex*,double>::const_iterator itv1 = _sizes.find(e->getVertex(0));
std::map<MVertex*,double>::const_iterator itv2 = _sizes.find(e->getVertex(1));
std::map<MVertex*,double>::const_iterator itv3 = _sizes.find(e->getVertex(2));
- return itv1->second * (1-uv[0]-uv[1]) + itv2->second * uv[0] + itv3->second * uv[1];
+ return itv1->second * (1-uv2[0]-uv2[1]) + itv2->second * uv2[0] + itv3->second * uv2[1];
}
-double backgroundMesh::getAngle(double x, double y, double z) const
+double backgroundMesh::getAngle(double u, double v, double w) const
{
- double xyz[3] = {x, y, z};
- double uv[3];
- MElement *e = _octree->find(x, y, z, 2, true);
+ double uv[3] = {u, v, w};
+ double uv2[3];
+ MElement *e = _octree->find(u, v, w, 2, true);
if (!e) {
- Msg::Error("cannot find %g %g %g", x, y, z);
+ Msg::Error("cannot find %g %g %g", u, v, w);
return 0.0;
}
- e->xyz2uvw(xyz, uv);
+ e->xyz2uvw(uv, uv2);
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-uv[0]-uv[1]) +
- cos (4*itv2->second) * uv[0] +
- cos (4*itv3->second) * uv[1] ;
- double sin4 = sin (4*itv1->second) * (1-uv[0]-uv[1]) +
- sin (4*itv2->second) * uv[0] +
- sin (4*itv3->second) * uv[1] ;
+
+ 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] +
+ sin (4*itv3->second) * uv2[1] ;
double angle = atan2(sin4,cos4)/4.0;
crossField2d::normalizeAngle (angle);
-
+
return angle;
}
@@ -795,7 +793,6 @@ void backgroundMesh::print(const std::string &filename, GFace *gf,
v3->x(),v3->y(),v3->z(),
itv1->second,itv2->second,itv3->second);
}
-
}
fprintf(f,"};\n");
fclose(f);
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index 85b21450dde0efec2e8686403b177ce1f7be7bd2..82e09ed488ec453776dce17c0ee745b5b35cb2c0 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -54,8 +54,8 @@ class backgroundMesh : public simpleFunction<double>
void propagate1dMesh(GFace *);
void propagatecrossField(GFace *);
void updateSizes(GFace *);
- double operator () (double x, double y, double z) const; // returns mesh size
- double getAngle(double x, double y, double z) const ;
+ double operator () (double u, double v, double w) const; // returns mesh size
+ double getAngle(double u, double v, double w) const ;
void print(const std::string &filename, GFace *gf,
const std::map<MVertex*, double>&) const;
void print(const std::string &filename, GFace *gf, int choice = 0) const