diff --git a/Geo/GEdgeCompound.h b/Geo/GEdgeCompound.h
index ea5476df7f9bb20a1ec765c72a9e9663cfca245f..d7db02a46be44eaa5a4fcfd152df21096a37c0b6 100644
--- a/Geo/GEdgeCompound.h
+++ b/Geo/GEdgeCompound.h
@@ -16,9 +16,9 @@ class GEdgeCompound : public GEdge {
std::vector<GEdge*> _compound;
std::vector<int> _orientation;
std::vector<double> _pars;
- void parametrize();
void orderEdges();
public:
+ void parametrize();
void getLocalParameter(const double &t, int &iEdge, double & tLoc) const;
void getCompoundParameter(GEdge *ge, const double &tLoc, double &t) const;
GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound);
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index de9f1ce6509b96ebbeb8d41ecde90650588433b4..e269e660edfe87d622f04bee83fb667ad7349155 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -431,9 +431,10 @@ void GFaceCompound::one2OneMap() const
bool GFaceCompound::parametrize() const
{
-
+
bool paramOK = true;
+ if(oct) return paramOK;
if(trivial()) return paramOK;
coordinates.clear();
@@ -716,6 +717,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
#endif
}
nbSplit = 0;
+ checked = false;
for(std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
if(!(*it)){
@@ -1218,6 +1220,9 @@ void GFaceCompound::computeNormals() const
double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
{
+
+ if(!oct) parametrize();
+
if(trivial()){
return (*(_compound.begin()))->curvatureMax(param);
}
@@ -1226,34 +1231,57 @@ double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
GFaceCompoundTriangle *lt;
getTriangle(param.x(), param.y(), <, U,V);
if(!lt){
- return 0.0;
+ return 0.0;
}
if(lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
return lt->gf->curvatureMax(pv);
}
- // emi fait qqch ici ...
+ else if (lt->gf->geomType() == GEntity::DiscreteSurface) {
+ //printf("!!!! compute here curvatureMax \n");
+ double curv= 0.;
+ curv = curvature(lt->tri,U,V);
+ return curv;
+ }
return 0.;
}
double GFaceCompound::curvature(MTriangle *t, double u, double v) const
{
+
SVector3 n1 = _normals[t->getVertex(0)];
SVector3 n2 = _normals[t->getVertex(1)];
SVector3 n3 = _normals[t->getVertex(2)];
double val[9] = {n1.x(), n2.x(), n3.x(),
n1.y(), n2.y(), n3.y(),
n1.z(), n2.z(), n3.z()};
+
return fabs(t->interpolateDiv(val, u, v, 0.0));
+
return 0.;
}
+SPoint2 GFaceCompound::parFromPoint(const SPoint3 &p) const
+{
+ if(!oct) parametrize();
+
+ std::map<SPoint3,SPoint3>::const_iterator it = _coordPoints.find(p);
+ SPoint3 sp = it->second;
+ printf("return u,v=%g %g \n", sp.x(), sp.y());
+
+ return SPoint2(sp.x(), sp.y());
+
+}
+
GPoint GFaceCompound::point(double par1, double par2) const
{
+
if(trivial()){
return (*(_compound.begin()))->point(par1,par2);
}
+ if(!oct) parametrize();
+
double U,V;
double par[2] = {par1,par2};
GFaceCompoundTriangle *lt;
@@ -1269,7 +1297,7 @@ GPoint GFaceCompound::point(double par1, double par2) const
// return lt->gf->point(pv.x(),pv.y());
// }
- const bool LINEARMESH = true; //false
+ const bool LINEARMESH = false; //false
if(LINEARMESH){
@@ -1295,18 +1323,41 @@ GPoint GFaceCompound::point(double par1, double par2) const
b300 = lt->v1;
b030 = lt->v2;
b003 = lt->v3;
- w12 = dot(lt->v2 - lt->v1, n1);
- w21 = dot(lt->v1 - lt->v2, n2);
- w23 = dot(lt->v3 - lt->v2, n2);
- w32 = dot(lt->v2 - lt->v3, n3);
- w31 = dot(lt->v1 - lt->v3, n3);
- w13 = dot(lt->v3 - lt->v1, n1);
- b210 = (2*lt->v1 + lt->v2-w12*n1)*0.333;
- b120 = (2*lt->v2 + lt->v1-w21*n2)*0.333;
- b021 = (2*lt->v2 + lt->v3-w23*n2)*0.333;
- b012 = (2*lt->v3 + lt->v2-w32*n3)*0.333;
- b102 = (2*lt->v3 + lt->v1-w31*n3)*0.333;
- b201 = (2*lt->v1 + lt->v3-w13*n1)*0.333;
+
+// w12 = dot(lt->v2 - lt->v1, n1);
+// w21 = dot(lt->v1 - lt->v2, n2);
+// w23 = dot(lt->v3 - lt->v2, n2);
+// w32 = dot(lt->v2 - lt->v3, n3);
+// w31 = dot(lt->v1 - lt->v3, n3);
+// w13 = dot(lt->v3 - lt->v1, n1);
+// b210 = (2*lt->v1 + lt->v2-w12*n1)*0.333;
+// b120 = (2*lt->v2 + lt->v1-w21*n2)*0.333;
+// b021 = (2*lt->v2 + lt->v3-w23*n2)*0.333;
+// b012 = (2*lt->v3 + lt->v2-w32*n3)*0.333;
+// b102 = (2*lt->v3 + lt->v1-w31*n3)*0.333;
+// b201 = (2*lt->v1 + lt->v3-w13*n1)*0.333;
+
+ //tagged PN trinagles (sigma=1)
+ double theta = 0.0;
+ SVector3 d1 = lt->v1+.33*(1-theta)*(lt->v2-lt->v1);
+ SVector3 d2 = lt->v2+.33*(1-theta)*(lt->v1-lt->v2);
+ SVector3 X1 = 1/norm(n1)*n1;
+ SVector3 X2 = 1/norm(n2)*n2;
+ b210 = d1 - dot(X1,d1-lt->v1)*X1;
+ b120 = d2 - dot(X2,d2-lt->v2)*X2;
+ SVector3 d3 = lt->v2+.33*(1-theta)*(lt->v3-lt->v2);
+ SVector3 d4 = lt->v3+.33*(1-theta)*(lt->v2-lt->v3);
+ SVector3 X3 = 1/norm(n2)*n2;
+ SVector3 X4 = 1/norm(n3)*n3;
+ b021 = d3 - dot(X3,d3-lt->v2)*X3;
+ b012 = d4 - dot(X4,d4-lt->v3)*X4;
+ SVector3 d5 = lt->v3+.33*(1-theta)*(lt->v1-lt->v3);
+ SVector3 d6 = lt->v1+.33*(1-theta)*(lt->v3-lt->v1);
+ SVector3 X5 = 1/norm(n3)*n3;
+ SVector3 X6 = 1/norm(n1)*n1;
+ b102 = d5 - dot(X5,d5-lt->v3)*X5;
+ b201 = d6 - dot(X6,d6-lt->v1)*X6;
+
E=(b210+b120+b021+b012+b102+b201)*0.16667;
VV=(lt->v1+lt->v2+lt->v3)*0.333;
b111=E+(E-VV)*0.5;
@@ -1324,6 +1375,8 @@ GPoint GFaceCompound::point(double par1, double par2) const
Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
{
+ if(!oct) parametrize();
+
if(trivial()){
return (*(_compound.begin()))->firstDer(param);
}
@@ -1351,7 +1404,56 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
void GFaceCompound::secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
{
- Msg::Error("Computation of the second derivatives not implemented for compound faces");
+
+ if(!oct) parametrize();
+
+ //use central differences
+
+ //EMI: TODO should take size of two or three triangles
+// double eps = 1e+2;
+
+// double u = param.x();
+// double v = param.y();
+// Pair<SVector3,SVector3> Der_u, Der_ueps, Der_v, Der_veps;
+
+// if(u - eps < 0.0) {
+// Der_u = firstDer(SPoint2(u,v));
+// Der_ueps = firstDer(SPoint2(u+eps,v));
+// }
+// else {
+// Der_u = firstDer(SPoint2(u-eps,v));
+// Der_ueps = firstDer(SPoint2(u,v));
+// }
+
+// if(v - eps < 0.0) {
+// Der_v = firstDer(SPoint2(u,v));
+// Der_veps = firstDer(SPoint2(u,v+eps));
+// }
+// else {
+// Der_v = firstDer(SPoint2(u,v-eps));
+// Der_veps = firstDer(SPoint2(u,v));
+// }
+
+// SVector3 dXdu_u = Der_u.first();
+// SVector3 dXdv_u = Der_u.second();
+// SVector3 dXdu_ueps = Der_ueps.first();
+// SVector3 dXdv_ueps = Der_ueps.second();
+// SVector3 dXdu_v = Der_v.first();
+// SVector3 dXdv_v = Der_v.second();
+// SVector3 dXdu_veps = Der_veps.first();
+// SVector3 dXdv_veps = Der_veps.second();
+
+// double inveps = 1./eps;
+// *dudu = inveps * (dXdu_u - dXdu_ueps) ;
+// *dvdv = inveps * (dXdv_v - dXdv_veps) ;
+// *dudv = inveps * (dXdu_v - dXdu_veps) ;
+
+ //printf("der second dudu = %g %g %g \n", dudu->x(), dudu->y(), dudu->z());
+ //printf("der second dvdv = %g %g %g \n", dvdv->x(), dvdv->y(), dvdv->z());
+ //printf("der second dudv = %g %g %g \n", dudv->x(), dudv->y(), dudv->z());
+
+ Msg::Error("Computation of the second derivatives is not implemented for compound faces");
+
}
static void GFaceCompoundBB(void *a, double*mmin, double*mmax)
@@ -1448,8 +1550,8 @@ void GFaceCompound::buildOct() const
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
for(int j = 0; j < 3; j++){
- std::map<MVertex*,SPoint3>::const_iterator itj =
- coordinates.find(t->getVertex(j));
+ std::map<MVertex*,SPoint3>::const_iterator itj = coordinates.find(t->getVertex(j));
+ _coordPoints.insert(std::make_pair(t->getVertex(j)->point(), itj->second));
bb += SPoint3(itj->second.x(),itj->second.y(),0.0);
}
count++;
@@ -1506,32 +1608,33 @@ bool GFaceCompound::checkTopology() const
//if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)return true;
bool correctTopo = true;
-
+
int Nb = _interior_loops.size();
int G = genusGeom() ;
if( G != 0 || Nb < 1) correctTopo = false;
-
+
if (!correctTopo){
Msg::Warning("Wrong topology: Genus=%d and N boundary=%d", G, Nb);
- nbSplit = G+2; //std::max(G+2, 2);
+ nbSplit = G+2;//std::max(G+2, 2);
Msg::Info("Split surface %d in %d parts with Mesh partitioner", tag(), nbSplit);
}
else
Msg::Info("Correct topology: Genus=%d and N boundary=%d", G, Nb);
-
+
return correctTopo;
+
}
bool GFaceCompound::checkAspectRatio() const
{
- if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)
- return true;
+ //if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)
+ // return true;
bool paramOK = true;
if(allNodes.empty()) buildAllNodes();
- double limit = 1.e15;
+ double limit = 1.e10;
double areaMax = 0.0;
int nb = 0;
std::list<GFace*>::const_iterator it = _compound.begin();
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 7efa430cb5e0f844e78b1b37f8b558eb754435d7..3e4f7ea3f709c49ac7a07b25b47af7f45af08991 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -65,6 +65,7 @@ class GFaceCompound : public GFace {
mutable v2t_cont adjv;
mutable bool mapv2Tri;
mutable std::map<MVertex*, SPoint3> coordinates;
+ mutable std::map<SPoint3,SPoint3 > _coordPoints;
mutable std::map<MVertex*, SVector3> _normals;
mutable std::list<MTriangle*> fillTris;
void buildOct() const ;
@@ -100,6 +101,7 @@ class GFaceCompound : public GFace {
Range<double> parBounds(int i) const
{ return trivial() ? (*(_compound.begin()))->parBounds(i) : Range<double>(-1, 1); }
virtual GPoint point(double par1, double par2) const;
+ SPoint2 parFromPoint(const SPoint3 &p) const;
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
virtual GEntity::GeomType geomType() const { return CompoundSurface; }
@@ -114,6 +116,7 @@ class GFaceCompound : public GFace {
void partitionFaceCM();
virtual std::list<GFace*> getCompounds() const {return _compound;};
mutable int nbSplit;
+ mutable bool checked;
private:
typeOfIsomorphism _type;
typeOfMapping _mapping;
@@ -134,12 +137,9 @@ class GFaceCompound : public GFace {
Msg::Error("Gmsh has to be compiled with solver support to use GFaceCompounds");
}
virtual ~GFaceCompound() {}
- virtual GPoint point(double par1, double par2) const { return GPoint(); }
- virtual Pair<SVector3, SVector3> firstDer(const SPoint2 ¶m) const
- {
- return Pair<SVector3, SVector3>(SVector3(0,0,0), SVector3(0,0,0));
- }
- virtual void secondDer(const SPoint2 ¶m,
+ GPoint point(double par1, double par2) const ;
+ Pair<SVector3, SVector3> firstDer(const SPoint2 ¶m) const;
+ void secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const{}
virtual SPoint2 getCoordinates(MVertex *v) const { return SPoint2(); }
void parametrize() const {}
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 4c37dedfb90fcc44d947f9949aef7c5d7b4267a7..636a0a2e47514a8ff09e928def97d419ab6f47f6 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1061,7 +1061,7 @@ void GModel::createTopologyFromMesh()
it != discRegions.end(); it++)
(*it)->setBoundFaces();
- //for each discreteFace, createTopology
+ //for each discreteFace, createTopology
std::vector<discreteFace*> discFaces;
for(fiter it = firstFace(); it != lastFace(); it++)
if((*it)->geomType() == GEntity::DiscreteSurface)
@@ -1119,7 +1119,7 @@ void GModel::createTopologyFromMesh()
// }
- //create Topo Fro Faces
+ //create Topo From Faces
createTopologyFromFaces(discFaces);
}
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index e911e3ac5cf1d219e1573e61d07f75ef986688bc..d1fe8ce0f3dcfe73b62c28c5c478a20bf8a0cdbf 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -143,10 +143,11 @@ int GModel::importGEOInternals()
ge->physicals.end())
ge->physicals.push_back(pnum);
}
+
+
// the physical is a compound i.e. we allow the meshes
// not to conform internal MEdges of the compound
- // the physical is a compound i.e. we allow the meshes
- // not to conform internal MEdges of the compound
+
if (p->Typ == MSH_PHYSICAL_LINE && p->Boundaries[0]){
GEdge *ge = getEdgeByTag(abs(p->Num));
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index cb7eb3b2ba2490937a74969e4390ef508708d5d1..8766d3ceedaa792b9d08c7b1a9251191f6373424 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -3395,6 +3395,7 @@ void setVolumeSurfaces(Volume *v, List_T *loops)
List_Add(v->SurfacesByTag, &is);
}
else{
+ printf("surface %d is not in GModel !\n", abs(is));
Msg::Error("Unknown surface %d", is);
return;
}
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 01c2d9fd16ba287b907a4e45a7cc4ac7559e813a..0b1aadbdd7dff1aa9238b8c6dd4c16681c6c4817 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -375,6 +375,24 @@ void discreteEdge::computeNormals () const
for ( ; itn != _normals.end(); ++itn){
itn->second.normalize();
}
+
+// //smooth normals
+// smooth_normals *normals = 0;
+// normals = new smooth_normals(180.);
+
+// for ( itn = _normals.begin(); itn != _normals.end(); ++itn){
+// MVertex *v = itn->first;
+// SVector3 n = itn->second;
+// normals->add(v->x(),v->y(),v->z(),n.x(),n.y(),n.z());
+// }
+// for ( itn = _normals.begin(); itn != _normals.end(); ++itn){
+// double nx,ny,nz;
+// MVertex *v = itn->first;
+// normals->get(v->x(),v->y(),v->z(),nx,ny,nz);
+// itn->second = SVector3(nx,ny,nz);
+// }
+// delete normals;
+
}
void discreteEdge::getLocalParameter(const double &t, int &iLine,
@@ -400,7 +418,7 @@ GPoint discreteEdge::point(double par) const
MVertex *vB = lines[iEdge]->getVertex(0);
MVertex *vE = lines[iEdge]->getVertex(1);
- const bool LINEARMESH = true; //false;
+ const bool LINEARMESH = false; //false; //false;
if (LINEARMESH){
//linear Lagrange mesh
@@ -423,10 +441,20 @@ GPoint discreteEdge::point(double par) const
b300 = v1;
b030 = v2;
+
w12 = dot(v2 - v1, n1);
w21 = dot(v1 - v2, n2);
b210 = (2*v1 + v2 -w12*n1)*0.333;
- b120 = (2*v2 + v1-w21*n2)*0.333;
+ b120 = (2*v2 + v1 -w21*n2)*0.333;
+
+// //tagged PN trinagles (sigma=1)
+ double theta = 0.0;
+ SVector3 d1 = v1+.33*(1-theta)*(v2-v1);
+ SVector3 d2 = v2+.33*(1-theta)*(v1-v2);
+ SVector3 X1 = 1/norm(n1)*n1;
+ SVector3 X2 = 1/norm(n2)*n2;
+ b210 = d1 - dot(X1,d1-v1)*X1;
+ b120 = d2 - dot(X2,d2-v2)*X2;
double U = tLoc;
double W = 1-U;
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 5edb000bac1af9fb17b5b73890d794747c479f4b..619415efe615bdc25f0f9a396daec22200d94fe5 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -10,6 +10,7 @@
#include "MTriangle.h"
#include "MEdge.h"
#include "Geo.h"
+#include "GFaceCompound.h"
discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
{
@@ -70,24 +71,66 @@ GPoint discreteFace::point(double par1, double par2) const
SPoint2 discreteFace::parFromPoint(const SPoint3 &p) const
{
- Msg::Error("Cannot compute parametric coordinates on discrete face");
- return SPoint2();
+
+ if (getCompound()){
+ printf("par from point in GFaceCompound %d \n", getCompound()->tag());
+ return getCompound()->parFromPoint(p);
+ }
+ else{
+ Msg::Error("Cannot compute parametric coordinates on discrete face");
+ return SPoint2();
+ }
}
SVector3 discreteFace::normal(const SPoint2 ¶m) const
{
- Msg::Error("Cannot evaluate normal on discrete face");
- return SVector3();
+
+ if (getCompound()){
+ return getCompound()->normal(param);
+ }
+ else{
+ Msg::Error("Cannot evaluate normal on discrete face");
+ return SVector3();
+ }
+
+}
+
+double discreteFace::curvatureMax(const SPoint2 ¶m) const
+{
+
+ if (getCompound()){
+ return getCompound()->curvatureMax(param);
+ }
+ else{
+ Msg::Error("Cannot evaluate curvature on discrete face");
+ return false;
+ }
+
}
Pair<SVector3, SVector3> discreteFace::firstDer(const SPoint2 ¶m) const
{
- Msg::Error("Cannot evaluate derivative on discrete face");
- return Pair<SVector3, SVector3>();
+
+ if (getCompound()){
+ return getCompound()->firstDer(param);
+ }
+ else{
+ Msg::Error("Cannot evaluate derivative on discrete face");
+ return Pair<SVector3, SVector3>();
+ }
+
}
void discreteFace::secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
{
- Msg::Error("Cannot evaluate derivative on discrete face");
+
+ if (getCompound()){
+ return getCompound()->secondDer(param, dudu, dvdv, dudv);
+ }
+ else{
+ Msg::Error("Cannot evaluate second derivative on discrete face");
+ return;
+ }
+
}
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index 87ed1db371b2153eefba5257e1fe3318e3e3baab..47095847ca208d389c1ecfdaeafecb615dc8f21b 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -15,10 +15,11 @@ class discreteFace : public GFace {
public:
discreteFace(GModel *model, int num);
virtual ~discreteFace() {}
- virtual GPoint point(double par1, double par2) const;
- virtual SPoint2 parFromPoint(const SPoint3 &p) const;
- virtual SVector3 normal(const SPoint2 ¶m) const;
- virtual GEntity::GeomType geomType() const { return DiscreteSurface; }
+ GPoint point(double par1, double par2) const;
+ SPoint2 parFromPoint(const SPoint3 &p) const;
+ SVector3 normal(const SPoint2 ¶m) const;
+ double curvatureMax(const SPoint2 ¶m) const;
+ GEntity::GeomType geomType() const { return DiscreteSurface; }
virtual Pair<SVector3, SVector3> firstDer(const SPoint2 ¶m) const;
virtual void secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 0087e967a32fb32152b239b9443804bdb518b584..65f17842934c68261ac8c34d48d2b25f1eb8c475 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -39,10 +39,13 @@ static double max_surf_curvature(const GEdge *ge, double u)
double val = 0;
std::list<GFace *> faces = ge->faces();
std::list<GFace *>::iterator it = faces.begin();
+ //printf("for gedge=%d facesnb=%d \n", ge->tag(), faces.size());
while(it != faces.end()){
if ((*it)->geomType() != GEntity::CompoundSurface){
+ printf("for gedge %d face=%d \n",ge->tag(), (*it)->tag());
SPoint2 par = ge->reparamOnFace((*it), u, 1);
double cc = (*it)->curvature(par);
+ printf("for face compute curvature =%g \n", cc);
val = std::max(cc, val);
}
++it;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 3962f0ec2e815c0000f7d5b891f218c75ffe6e46..09d1ca0abbb67b9b74469cb02f0efdc11b40ab4d 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1346,32 +1346,37 @@ void partitionAndRemesh(GFaceCompound *gf)
int nume = gf->model()->maxEdgeNum() + 1;
int numf = gf->model()->maxFaceNum() + 1;
std::vector<discreteFace*> pFaces;
- createPartitionFaces(gf->model(), gf, NF, pFaces); //WARNING CREATE NB and not NF faces
+ createPartitionFaces(gf->model(), gf, NF, pFaces);
gf->model()->createTopologyFromFaces(pFaces);
- CreateOutputFile("toto.msh", CTX::instance()->mesh.format);
- //Msg::Exit(1);
-
- //Remesh new faces (Compound Lines and Compound Surfaces)
+ Msg::Info("-----------------------------------------------------------");
+ Msg::Info("Multiscale Partition SUCCESSFULLY PERFORMED : %d parts", NF);
+ CreateOutputFile("multiscalePARTS.msh", CTX::instance()->mesh.format);
+ Msg::Info("-----------------------------------------------------------");
+
+ //Remesh new faces (Compound Lines and Compound Surfaces)
//-----------------------------------------------------
- //Remeshing discrete Edges
+ Msg::Info("-----------------------------------------------------------");
+ Msg::Info("Parametrize Compounds");
+ Msg::Info("-----------------------------------------------------------");
+
+ //Parametrize Compound Lines
int NE = gf->model()->maxEdgeNum() - nume + 1;
for (int i=0; i < NE; i++){
std::vector<GEdge*>e_compound;
GEdge *pe = gf->model()->getEdgeByTag(nume+i);//partition edge
e_compound.push_back(pe);
int num_gec = nume + NE + i ;
- Msg::Info("*** Remeshing discreteEdge %d with CompoundLine %d", pe->tag(), num_gec);
- GEdge *gec = new GEdgeCompound(gf->model(), num_gec, e_compound);
+ Msg::Info("Parametrize Compound Line (%d) = %d discrete edge", num_gec, pe->tag() );
+ GEdgeCompound *gec = new GEdgeCompound(gf->model(), num_gec, e_compound);
gf->model()->add(gec);
-
- meshGEdge mge;
- mge(gec);//meshing 1D
+
+ gec->parametrize();
}
- //Remeshing discrete Face
+ //Parametrize Compound surfaces
std::list<GEdge*> b[4];
std::set<MVertex*> allNod;
for (int i=0; i < NF; i++){
@@ -1379,10 +1384,28 @@ void partitionAndRemesh(GFaceCompound *gf)
GFace *pf = gf->model()->getFaceByTag(numf+i);//partition face
int num_gfc = numf + NF + i ;
f_compound.push_back(pf);
- Msg::Info("*** Remeshing discreteFace %d with CompoundSurface %d", pf->tag(), num_gfc);
- GFace *gfc = new GFaceCompound(gf->model(), num_gfc, f_compound, b[0], b[1], b[2], b[3]);
+ Msg::Info("Parametrize Compound Surface (%d) = %d discrete face", num_gfc, pf->tag() );
+ GFaceCompound *gfc = new GFaceCompound(gf->model(), num_gfc, f_compound, b[0], b[1], b[2], b[3]);
gf->model()->add(gfc);
+ gfc->parametrize();
+ }
+
+ Msg::Info("-----------------------------------------------------------");
+ Msg::Info("Mesh Compounds");
+ Msg::Info("-----------------------------------------------------------");
+
+ //Mesh 1D and 2D
+ for (int i=0; i < NE; i++){
+ GEdge *gec = gf->model()->getEdgeByTag(nume + NE + i);
+ meshGEdge mge;
+ mge(gec);//meshing 1D
+ }
+
+ Msg::Info("*** Starting Mesh of surface %d ...", gf->tag());
+
+ for (int i=0; i < NF; i++){
+ GFace *gfc = gf->model()->getFaceByTag(numf + NF + i );
meshGFace mgf;
mgf(gfc);//meshing 2D
@@ -1460,6 +1483,7 @@ void partitionAndRemesh(GFaceCompound *gf)
}
Msg::Info("*** Mesh of surface %d done by assembly remeshed faces", gf->tag());
+ Msg::Info("-----------------------------------------------------------");
gf->coherenceNormals();
gf->meshStatistics.status = GFace::DONE;
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index 3dfeeb33a2d920fc668a3b55c748b03cd31ca1fc..dda34a8a845659a29f94ce26cafdab678cfe909c 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -90,7 +90,7 @@ static int getGenus (std::vector<MElement *> &elements,
}
static int getAspectRatio (std::vector<MElement *> &elements,
- std::vector<std::vector<MEdge> > &boundaries)
+ std::vector<std::vector<MEdge> > &boundaries)
{
std::set<MVertex*> vs;
@@ -109,9 +109,9 @@ static int getAspectRatio (std::vector<MElement *> &elements,
bb +=pt;
}
-// obb = SOrientedBoundingBox::buildOBB(vertices);
-// double H = obbox.getMaxSize();
- double H = norm(SVector3(bb.max(), bb.min()));
+ //obbox = SOrientedBoundingBox::buildOBB(vertices);
+ //double H = obbox.getMaxSize();
+ double H = norm(SVector3(bb.max(), bb.min()));
double D = 0.0;
if (boundaries.size() == 0){
@@ -244,11 +244,6 @@ int multiscalePartition::assembleAllPartitions(){
}
}
- Msg::Info("-----------------------------------------------------------");
- Msg::Info("Multiscale Partition SUCCESSFULLY PERFORMED : %d parts", nbParts-1);
- Msg::Info("-----------------------------------------------------------");
-
-
return nbParts-1;
diff --git a/benchmarks/3d/durufle.geo b/benchmarks/3d/durufle.geo
index ddd8ec7e34f1b36b10c77f27cefce56fde0cfb5b..1e0f8bb789cce761da45580e5d159378cf7f5b16 100644
--- a/benchmarks/3d/durufle.geo
+++ b/benchmarks/3d/durufle.geo
@@ -1,5 +1,5 @@
// fails only for lc=0.55 - otherwise works great!
-lc = 0.55;
+lc = 0.95;
r = 2.0;
e = 0.8;
rb = r+e;
diff --git a/benchmarks/extrude/sphere_boundary_layer.geo b/benchmarks/extrude/sphere_boundary_layer.geo
index 1f52d6249e4ae14dd7706253d38063fc94e69818..66014e14b307ddaeb5d31f7bce7a2f8ea9982c3e 100644
--- a/benchmarks/extrude/sphere_boundary_layer.geo
+++ b/benchmarks/extrude/sphere_boundary_layer.geo
@@ -1,4 +1,5 @@
lc = .2;
+
Point(1) = {0.0,0.0,0.0,lc};
Point(2) = {1,0.0,0.0,lc};
Point(3) = {0,1,0.0,lc};
@@ -18,6 +19,7 @@ Circle(9) = {2,1,7};
Circle(10) = {7,1,4};
Circle(11) = {4,1,6};
Circle(12) = {6,1,2};
+
Line Loop(13) = {2,8,-10};
Ruled Surface(14) = {13};
Line Loop(15) = {10,3,7};
@@ -34,11 +36,13 @@ Line Loop(25) = {-7,4,9};
Ruled Surface(26) = {25};
Line Loop(27) = {-4,12,-6};
Ruled Surface(28) = {27};
-Surface Loop(29) = {28,26,16,14,20,24,22,18};
-Volume(30) = {29};
Extrude {
Surface{14:28:2}; Layers{10, 0.1}; // Recombine;
}
+Surface Loop(29) = {28,26,16,14,20,24,22,18};
+Volume(30) = {29};
+
+
Mesh.Algorithm3D = 4;
diff --git a/benchmarks/stl/PelvisARTHUR_CLASS_GEO.geo b/benchmarks/stl/PelvisARTHUR_CLASS_GEO.geo
index a60748d70510508a56eb170439e12f55e894a914..1991c97dedeaa04ce63a990035354bbfb4f4aee3 100644
--- a/benchmarks/stl/PelvisARTHUR_CLASS_GEO.geo
+++ b/benchmarks/stl/PelvisARTHUR_CLASS_GEO.geo
@@ -1,4 +1,4 @@
-Mesh.CharacteristicLengthFactor=0.2;
+Mesh.CharacteristicLengthFactor=0.133;
Merge "PelvisARTHUR_CLASS.msh";
CreateTopology;
@@ -12,4 +12,4 @@ Compound Line(20)={5};
Compound Surface(200)={2};
Compound Surface(400)={3};
-Compound Volume(2000) = {1003};
+//Compound Volume(2000) = {1003};