diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index a2fbfc862e261dd431ff9ce48f184e9cec58fff1..808ae12cd5300f690b401213782cc99f2eac49be 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -954,8 +954,8 @@ bool GFaceCompound::parametrize() const
if(!success) {Msg::Error("Could not order vertices on boundary");exit(1);}
}
- fillNeumannBCS_Plane();
- //fillNeumannBCS();
+ //fillNeumannBCS_Plane();
+ fillNeumannBCS();
// Convex parametrization
if (_mapping == CONVEX){
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 48226300cca770f1b92a3f39e2046c56c72c9291..0fbc440d6e838ce92ac8866b4fa468ecd1d95c20 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -158,7 +158,7 @@ GPoint OCCEdge::point(double par) const
}
else if(!curve.IsNull()){
gp_Pnt pnt = curve->Value(par);
- return GPoint(pnt.X(), pnt.Y(), pnt.Z());
+ return GPoint(pnt.X(), pnt.Y(), pnt.Z(),this,par);
}
else{
Msg::Warning("OCC Curve %d is neither a 3D curve not a trimmed curve", tag());
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 02e05831927324c14e18439a994e1a1e574f9825..e11a62f94b9c72095d6f4c9245154f49f6fff91f 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -510,7 +510,8 @@ backgroundMesh::~backgroundMesh()
}
static void propagateValuesOnFace(GFace *_gf,
- std::map<MVertex*,double> &dirichlet)
+ std::map<MVertex*,double> &dirichlet,
+ bool in_parametric_plane = false)
{
#if defined(HAVE_SOLVER)
linearSystem<double> *_lsys = 0;
@@ -541,6 +542,17 @@ static void propagateValuesOnFace(GFace *_gf,
for (int j=0;j<3;j++)vs.insert(_gf->triangles[k]->getVertex(j));
for (unsigned int k = 0; k < _gf->quadrangles.size(); k++)
for (int j=0;j<4;j++)vs.insert(_gf->quadrangles[k]->getVertex(j));
+
+ std::map<MVertex*,SPoint3> theMap;
+ if ( in_parametric_plane) {
+ for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
+ SPoint2 p;
+ 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);
@@ -560,7 +572,13 @@ static void propagateValuesOnFace(GFace *_gf,
for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
myAssembler.getDofValue(*it, 0, 1, dirichlet[*it]);
}
-
+
+ if ( in_parametric_plane) {
+ for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
+ SPoint3 p = theMap[(*it)];
+ (*it)->setXYZ(p.x(),p.y(),p.z());
+ }
+ }
delete _lsys;
#endif
}
@@ -686,8 +704,8 @@ void backgroundMesh::propagatecrossField(GFace *_gf)
}
}
- propagateValuesOnFace(_gf,_cosines4);
- propagateValuesOnFace(_gf,_sines4);
+ propagateValuesOnFace(_gf,_cosines4,false);
+ propagateValuesOnFace(_gf,_sines4,false);
std::map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
for ( ; itv2 != _2Dto3D.end(); ++itv2){
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 9df5f0d22cb5880ab2648d7e79939d7e7a503903..7d1a2483a028f4c79b6f7d0322b3b1d384ded60d 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -317,7 +317,7 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
else {
int count = u0<u1? j + 1 : nPts + 1 - (j + 1);
GPoint pc = ge->point(US[count]);
- v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge,pc.u());
+ v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge,US[count]);
// printf("Edge %d(%g) %d(%g) new vertex %g %g at %g\n",v0->getNum(),u0,v1->getNum(),u1,v->x(),v->y(), US[count]);
}
temp.push_back(v);
@@ -1325,17 +1325,15 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
// printJacobians(m, "smoothness.pos");
// m->writeMSH("SMOOTHED.msh");
// FIXME !!
- if (!linear){
+ if (!linear && CTX::instance()->mesh.smoothInternalEdges){
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
std::vector<MElement*> v;
v.insert(v.begin(), (*it)->triangles.begin(), (*it)->triangles.end());
v.insert(v.end(), (*it)->quadrangles.begin(), (*it)->quadrangles.end());
- //hot.applySmoothingTo(v, (*it));
- //hot.applySmoothingTo(v, .1,0);
+ hot.applySmoothingTo(v, (*it));
}
// hot.ensureMinimumDistorsion(0.1);
- checkHighOrderTriangles("Final mesh", m, bad, worst);
- checkHighOrderTetrahedron("Final mesh", m, bad, worst);
+ checkHighOrderTriangles("Final surface mesh", m, bad, worst);
}
// m->writeMSH("CORRECTED.msh");
diff --git a/Mesh/highOrderTools.cpp b/Mesh/highOrderTools.cpp
index dd872d712d5fc5255154931afad627c3d22242ea..7ed6dce1a15525db34ab8c5bd339500f12f89d0c 100644
--- a/Mesh/highOrderTools.cpp
+++ b/Mesh/highOrderTools.cpp
@@ -50,6 +50,7 @@ void highOrderTools::moveToStraightSidedLocation(MElement *e) const
}
}
+
void highOrderTools::ensureMinimumDistorsion(MElement *e, double threshold)
{
double dist = e->distoShapeMeasure();
@@ -154,8 +155,6 @@ double highOrderTools::applySmoothingTo(GFace *gf, double tres, bool mixed)
std::vector<MElement*> v;
v.insert(v.begin(), gf->triangles.begin(),gf->triangles.end());
v.insert(v.begin(), gf->quadrangles.begin(),gf->quadrangles.end());
- Msg::Info("Smoothing high order mesh : model face %d (%d elements)", gf->tag(),
- v.size());
// applySmoothingTo(v,gf);
return applySmoothingTo(v,tres, mixed);
}
@@ -201,9 +200,10 @@ void highOrderTools::computeMetricInfo(GFace *gf,
fullVector<double> &D)
{
int nbNodes = e->getNumVertices();
+ // printf("ELEMENT --\n");
for (int j = 0; j < nbNodes; j++){
SPoint2 param;
- reparamMeshVertexOnFace(e->getVertex(j), gf, param);
+ bool success = reparamMeshVertexOnFace(e->getVertex(j), gf, param);
// printf("%g %g vs %g %g %g\n",param.x(),param.y(),
// e->getVertex(j)->x(),e->getVertex(j)->y(),e->getVertex(j)->z());
@@ -228,18 +228,18 @@ void highOrderTools::computeMetricInfo(GFace *gf,
JT(VJ,YJ) = der.second().y();
JT(VJ,ZJ) = der.second().z();
- GPoint gp = gf->point(param);
SVector3 ss = getSSL(e->getVertex(j));
- D(XJ) = gp.x() - ss.x();
- D(YJ) = gp.y() - ss.y();
- D(ZJ) = gp.z() - ss.z();
- if (e->getVertex(j)->onWhat()->dim() == 1)
- printf("point %d on %d %d dx = %g %g %g --> %g %g %g vs. %g %g %g --> %g %g %g \n",e->getVertex(j)->getNum(),
- e->getVertex(j)->onWhat()->dim(),e->getVertex(j)->onWhat()->tag(),
- D(XJ),D(YJ),D(ZJ),
- gp.x(),gp.y(),gp.z(),
- ss.x(),ss.y(),ss.z(),
- e->getVertex(j)->x(),e->getVertex(j)->y(),e->getVertex(j)->z());
+ GPoint gp = gf->point(param);
+ D(XJ) = (gp.x() - ss.x());
+ D(YJ) = (gp.y() - ss.y());
+ D(ZJ) = (gp.z() - ss.z());
+ /*
+ printf("point %d on %d %d dx = %g %g %g --> %g %g %g --> %g %g %g \n",e->getVertex(j)->getIndex(),
+ e->getVertex(j)->onWhat()->dim(),e->getVertex(j)->onWhat()->tag(),
+ D(XJ),D(YJ),D(ZJ),
+ ss.x(),ss.y(),ss.z(),
+ e->getVertex(j)->x(),e->getVertex(j)->y(),e->getVertex(j)->z());
+ */
}
}
@@ -251,25 +251,26 @@ void highOrderTools::applySmoothingTo(std::vector<MElement*> & all, GFace *gf)
#else
linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
#endif
-
+ // compute the straight sided positions of high order nodes that are
+ // on the edges of the face in the UV plane
dofManager<double> myAssembler(lsys);
elasticityTerm El(0, 1.0, .33, _tag);
std::vector<MElement*> layer, v;
double minD;
getDistordedElements(all, 0.5, v, minD);
- const int nbLayers = 322;
+ int numBad = v.size();
+ const int nbLayers = 3;
for (int i = 0; i < nbLayers; i++){
addOneLayer(all, v, layer);
v.insert(v.end(), layer.begin(), layer.end());
}
- Msg::Debug("%d elements after adding %d layers", (int)v.size(), nbLayers);
if (!v.size()) return;
+ Msg::Info("Smoothing high order mesh : model face %d (%d elements considered in the elastic analogy, worst mapping %12.5E, %3d bad elements)", gf->tag(),
+ v.size(),minD,numBad);
+
addOneLayer(all, v, layer);
std::set<MVertex*>::iterator it;
- std::map<MVertex*,SVector3>::iterator its;
- std::map<MVertex*,SVector3>::iterator itpresent;
- std::map<MVertex*,SVector3>::iterator ittarget;
std::set<MVertex*> verticesToMove;
// on the last layer, fix displacement to 0
@@ -283,12 +284,13 @@ void highOrderTools::applySmoothingTo(std::vector<MElement*> & all, GFace *gf)
// fix all vertices that cannot move
for (unsigned int i = 0; i < v.size(); i++){
- // moveToStraightSidedLocation(v[i]);
+ moveToStraightSidedLocation(v[i]);
for (int j = 0; j < v[i]->getNumVertices(); j++){
MVertex *vert = v[i]->getVertex(j);
if (vert->onWhat()->dim() < 2){
- myAssembler.fixVertex(vert, 0, _tag, 0);
- myAssembler.fixVertex(vert, 1, _tag, 0);
+ double du = 0, dv = 0;
+ myAssembler.fixVertex(vert, 0, _tag, du);
+ myAssembler.fixVertex(vert, 1, _tag, dv);
}
}
}
@@ -307,7 +309,7 @@ void highOrderTools::applySmoothingTo(std::vector<MElement*> & all, GFace *gf)
double dx = dx0;
Msg::Debug(" dx0 = %12.5E", dx0);
int iter = 0;
- while(0){
+ while(1){
double dx2 = smooth_metric_(v, gf, myAssembler, verticesToMove, El);
Msg::Debug(" dx2 = %12.5E", dx2);
if (fabs(dx2 - dx) < 1.e-4 * dx0)break;
@@ -317,21 +319,20 @@ void highOrderTools::applySmoothingTo(std::vector<MElement*> & all, GFace *gf)
for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
SPoint2 param;
- reparamMeshVertexOnFace(*it, gf, param);
- GPoint gp = gf->point(param);
if ((*it)->onWhat()->dim() == 2){
+ reparamMeshVertexOnFace(*it, gf, param);
+ GPoint gp = gf->point(param);
(*it)->x() = gp.x();
(*it)->y() = gp.y();
(*it)->z() = gp.z();
_targetLocation[*it] = SVector3(gp.x(), gp.y(), gp.z());
}
else{
- SVector3 p = _targetLocation[(*it)];
- // SVector3 p = getSSL(*it);
- printf("%12.5E %12.5E %12.5E %d %d\n",p.x(),p.y(),p.z(),(*it)->onWhat()->dim(),(*it)->onWhat()->tag());
- // (*it)->x() = p.x();
- // (*it)->y() = p.y();
- // (*it)->z() = p.z();
+ SVector3 p = getTL(*it);
+ // printf("%12.5E %12.5E %12.5E %d %d\n",p.x(),p.y(),p.z(),(*it)->onWhat()->dim(),(*it)->onWhat()->tag());
+ (*it)->x() = p.x();
+ (*it)->y() = p.y();
+ (*it)->z() = p.z();
}
}
delete lsys;
@@ -373,6 +374,7 @@ double highOrderTools::smooth_metric_(std::vector<MElement*> & v,
// K33.print("K33");
// K22.print("K22");
J23K33.mult(D3, R2);
+ // R2.print("hopla");
for (int j = 0; j < n2; j++){
Dof RDOF = El.getLocalDofR(&se, j);
myAssembler.assemble(RDOF, -R2(j));
@@ -388,16 +390,16 @@ double highOrderTools::smooth_metric_(std::vector<MElement*> & v,
for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
if ((*it)->onWhat()->dim() == 2){
- SPoint2 param;
- reparamMeshVertexOnFace((*it), gf, param);
- SPoint2 dparam;
- myAssembler.getDofValue((*it), 0, _tag, dparam[0]);
- myAssembler.getDofValue((*it), 1, _tag, dparam[1]);
- printf("%g %g -- %g %g\n",dparam[0],dparam[1],param[0],param[1]);
+ SPoint2 param;
+ reparamMeshVertexOnFace((*it), gf, param);
+ SPoint2 dparam;
+ myAssembler.getDofValue((*it), 0, _tag, dparam[0]);
+ myAssembler.getDofValue((*it), 1, _tag, dparam[1]);
+ // printf("%g %g -- %g %g\n",dparam[0],dparam[1],param[0],param[1]);
SPoint2 newp = param+dparam;
- dx += newp.x() * newp.x() + newp.y() * newp.y();
- (*it)->setParameter(0, newp.x());
- (*it)->setParameter(1, newp.y());
+ dx += newp.x() * newp.x() + newp.y() * newp.y();
+ (*it)->setParameter(0, newp.x());
+ (*it)->setParameter(1, newp.y());
}
}
myAssembler.systemClear();
@@ -420,7 +422,6 @@ highOrderTools::highOrderTools (GModel *gm, GModel *mesh, int order) : _gm(gm),
computeStraightSidedPositions();
}
-
void highOrderTools::computeStraightSidedPositions ()
{
_clean();
@@ -450,21 +451,13 @@ void highOrderTools::computeStraightSidedPositions ()
(*it)->lines[i]->getVertex(1)->y(),
(*it)->lines[i]->getVertex(1)->z());
- for (int k=0;k<2;k++){
- if (_straightSidedLocation.find((*it)->lines[i]->getVertex(k)) == _straightSidedLocation.end()){
- _straightSidedLocation [(*it)->lines[i]->getVertex(k)] = k ? p1 : p0;
- _targetLocation[(*it)->lines[i]->getVertex(k)] = k ? p1 : p0;
- }
- }
-
-
for (int j=1;j<=N;j++){
const double xi = (double)(j)/(N+1);
// printf("xi = %g\n",xi);
- const SVector3 straightSidedPoint = p0 *(1.-xi) + p1*xi;
+ const SVector3 straightSidedPoint = p0 *(1.-xi) + p1*xi;
MVertex *v = (*it)->lines[i]->getVertex(j+1);
if (_straightSidedLocation.find(v) == _straightSidedLocation.end()){
- _straightSidedLocation [v] = straightSidedPoint;
+ _straightSidedLocation [v] = straightSidedPoint;
_targetLocation[v] = SVector3(v->x(),v->y(),v->z());
}
}
@@ -562,7 +555,6 @@ void highOrderTools::computeStraightSidedPositions ()
Msg::Info("highOrderTools has been set up : %d nodes are considered",_straightSidedLocation.size());
}
-
// apply a displacement that does not create elements that are
// distorted over a value "thres"
double highOrderTools::apply_incremental_displacement (double max_incr,
diff --git a/Mesh/highOrderTools.h b/Mesh/highOrderTools.h
index 71b5cf4ab79fca0644223133b3135e01d4b2dac9..222fc7741c57ca5df855cfb206f2e7e5506b5aad 100644
--- a/Mesh/highOrderTools.h
+++ b/Mesh/highOrderTools.h
@@ -73,6 +73,12 @@ class highOrderTools
if (it != _straightSidedLocation.end())return it->second;
else return SVector3(v->x(),v->y(),v->z());
}
+ inline SVector3 getTL(MVertex *v) const
+ {
+ std::map<MVertex*,SVector3>::const_iterator it = _targetLocation.find(v);
+ if (it != _targetLocation.end())return it->second;
+ else return SVector3(v->x(),v->y(),v->z());
+ }
void makePosViewWithJacobians(const char *nm);
};
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index c8e8974b52c32e9d527cc30a692cf82adb3f76bf..f4e13e4be8eaf7973565a640f616311fd6e99d29 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -763,7 +763,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
else {
ptin = search4Triangle (worst, center, Us, Vs, AllTris,uv);
// if (!ptin){
- // printf("strange : %g %g seems to be out of the domain for face %d\n",center[0],center[1],gf->tag());
+ // printf("strange : %g %g seems to be out of the domain for face %d\n",center[0],center[1],gf->tag());
// }
if (ptin)inside = true;
}
@@ -801,7 +801,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
MTriangle *base = worst->tri();
- // Msg::Info("Point %g %g cannot be inserted because %d", center[0], center[1], p.succeeded() );
+ Msg::Info("Point %g %g cannot be inserted because %d", center[0], center[1], p.succeeded() );
AllTris.erase(it);
worst->forceRadius(-1);
@@ -817,7 +817,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
}
else {
MTriangle *base = worst->tri();
- /*
+ /*
Msg::Info("Point %g %g is outside (%g %g , %g %g , %g %g)",
center[0], center[1],
Us[base->getVertex(0)->getIndex()],
@@ -1312,9 +1312,10 @@ void optimalPointFrontalQuadB (GFace *gf,
double metric[3]){
- // optimalPointFrontalQuad (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
-
+ optimalPointFrontalQuad (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
+ return;
MTriangle *base = worst->tri();
+
int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
int ip2 = active_edge;
int ip3 = (active_edge+1)%3;
@@ -1328,10 +1329,11 @@ void optimalPointFrontalQuadB (GFace *gf,
double midpoint[2] = {0.5 * (P[0] + Q[0]), 0.5 * (P[1] + Q[1])};
// compute background mesh data
- double quadAngle = backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0);
+ double quadAngle = backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0);
+ // double quadAngle = 0;
double center[2];
circumCenterInfinite (base, quadAngle,Us,Vs,center);
-
+
// rotate the points with respect to the angle
double XP1 = 0.5*(Q[0] - P[0]);
double YP1 = 0.5*(Q[1] - P[1]);
@@ -1339,23 +1341,44 @@ void optimalPointFrontalQuadB (GFace *gf,
double yp = -XP1 * sin(quadAngle) + YP1 * cos(quadAngle);
double t = 0.0;
+ double factor;
+
+ /*
+ consider a straight edge going from -x,-y to point x,y
+ the distance in infinity norm between the two points is 2x
+ if |x| > |y|
+
+ the point Q (-x,2x-y) is in infinity norm at equal distance
+ 2x from the 2 points
+
+ The orthogonal projection of Q onto the line can be found as
+ P = t (x,y) , t \in R
+ (Q - P) . (x,y) = 0 --> (-x-tx,2x-y-ty) . (x,y) = 0 --> -x^2 -tx^2 + 2xy - y^2 -t y^2 = 0
+ --> t (x^2 + y^2) = -x^2 + 2xy - y^2 --> t = - (x-y)^2 / (x^2 + y^2) = -1 + 2 xy /(x^2+y^2)
+
+ The L2 distance between Q and the line is
+
+ h^2 = (-x - tx)^2 + (2x-y-ty)^2 = x^2 (1+t)^2 +
+
+
+ */
+
if (fabs(xp)>=fabs(yp)){
t = -1. + 2*fabs(xp*yp)/(xp*xp + yp*yp);
+ factor = fabs(xp) / sqrt(xp*xp+yp*yp);
}
else {
t = 1. - 2*fabs(xp*yp)/(xp*xp + yp*yp);
+ factor = fabs(yp) / sqrt(xp*xp+yp*yp);
}
- const double usq = 1./sqrt(2.0);
- double factor = usq + fabs(t) * (1.-usq);
-
- // printf("t = %lf %lf\n",t,factor);
SVector3 P3(base->getVertex(ip1)->x(),base->getVertex(ip1)->y(),base->getVertex(ip1)->z());
SVector3 Q3(base->getVertex(ip2)->x(),base->getVertex(ip2)->y(),base->getVertex(ip2)->z());
SVector3 O3(base->getVertex(ip3)->x(),base->getVertex(ip3)->y(),base->getVertex(ip3)->z());
SVector3 PMID = P3*(1.-t)*.5 + Q3*(t+1.)*.5;
SVector3 N = crossprod (O3-P3,Q3-P3);
SVector3 T = Q3-P3;
+ double sizeEdge = T.norm();
N.normalize();
T.normalize();
SVector3 N2 = crossprod (T,N);
@@ -1369,7 +1392,8 @@ void optimalPointFrontalQuadB (GFace *gf,
(vSizesBGM[base->getVertex(ip1)->getIndex()] +
vSizesBGM[base->getVertex(ip2)->getIndex()] ) ;// * RATIO;
const double a = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
- double d = a*factor;
+
+ double d = (a+sizeEdge*factor)*factor*.5;
if (gf->geomType() == GEntity::CompoundSurface){
GFaceCompound *gfc = dynamic_cast<GFaceCompound*> (gf);
if (gfc){
@@ -1377,8 +1401,22 @@ void optimalPointFrontalQuadB (GFace *gf,
if (gp.succeeded()){
newPoint[0] = gp.u();
newPoint[1] = gp.v();
+ buildMetric(gf, newPoint, metric);
return ;
}
+ else {
+ gp = gfc->intersectionWithCircle(N2*(-1.0),N,PMID,d,newPoint);
+ if (gp.succeeded()){
+ Msg::Warning("--- HEY HEY -----------");
+ newPoint[0] = gp.u();
+ newPoint[1] = gp.v();
+ buildMetric(gf, newPoint, metric);
+ return ;
+ }
+ else {
+ Msg::Warning("--- NEVER GO THERE -----------");
+ }
+ }
}
}
double uvt[3] = {newPoint[0],newPoint[1],0.0};
@@ -1395,7 +1433,7 @@ void optimalPointFrontalQuadB (GFace *gf,
newPoint[1] = 100000;
Msg::Warning("--- Non optimal point found -----------");
}
- buildMetric(gf, newPoint, metric);
+ // buildMetric(gf, newPoint, metric);
}
@@ -1524,7 +1562,7 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
// compute the middle point of the edge
double newPoint[2],metric[3]={1,0,1};
- if (quad)optimalPointFrontalQuad (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
+ if (quad)optimalPointFrontalQuadB (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
else optimalPointFrontalB (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
diff --git a/Solver/elasticityTerm.cpp b/Solver/elasticityTerm.cpp
index d47611f2e32a1f971c0912268b004345934564df..512144a6ce9a36e01341df8ee934d4f6175b0090 100644
--- a/Solver/elasticityTerm.cpp
+++ b/Solver/elasticityTerm.cpp
@@ -37,7 +37,7 @@ void elasticityTerm::createData(MElement *e) const
d.w.push_back(w);
d.weight.push_back(weight);
}
- printf("coucou creation of a data for %d with %d points\n",e->getTypeForMSH(),npts);
+ // printf("coucou creation of a data for %d with %d points\n",e->getTypeForMSH(),npts);
_data[e->getTypeForMSH()] = d;
}