diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index fa6bb9272c65a055680e8e69b006b0810ab9d639..a47ba07e0bbee6dbe4bf7145c405a6d2b862f10e 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -405,7 +405,6 @@ static MTriangle* findInvertedTri(v2t_cont &adjv, MVertex*v0, MVertex*v1, MVerte
std::map<MVertex*, SPoint3> &vCoord, double nTot)
{
- MTriangle *tri=NULL;
v2t_cont :: iterator it = adjv.find(v1);
std::vector<MElement*> vTri = it->second;
MTriangle *myTri=NULL;
@@ -774,7 +773,7 @@ void GFaceCompound::convexBoundary(double nTot) const
//start the loop
int nbInv = 0;
int nbInvTri = 0;
- for(int i = 0; i < oVert.size(); i++){
+ for(unsigned int i = 0; i < oVert.size(); i++){
MVertex *v0 = oVert[i];
MVertex *v1, *v2;
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 1781c10bb31e39dc02e68f69a663ab51f2ebbaa4..82b52c4b5b5d195effa1ffcd2d6e6cb835650c39 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -84,19 +84,19 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
vertb = Create_Vertex(gvb->tag(), gvb->x(), gvb->y(), gvb->z(),
gvb->prescribedMeshSizeAtVertex(), 1.0);
Tree_Add(gm->getGEOInternals()->Points, &vertb);
+ vertb->Typ = MSH_POINT;
+ vertb->Num = gvb->tag();
}
if (!verte){
verte = Create_Vertex(gve->tag(), gve->x(), gve->y(), gve->z(),
gve->prescribedMeshSizeAtVertex(), 1.0);
Tree_Add(gm->getGEOInternals()->Points, &verte);
+ verte->Typ = MSH_POINT;
+ verte->Num = gve->tag();
}
- List_T *temp = List_Create(2, 2, sizeof(int));
- List_Add(temp, &vertb);
- List_Add(temp, &verte);
-
if (ge->geomType() == GEntity::Line){
- c = Create_Curve(numEdge, MSH_SEGM_LINE, 1, temp, NULL, -1, -1, 0., 1.);
+ c = Create_Curve(numEdge, MSH_SEGM_LINE, 1, NULL, NULL, -1, -1, 0., 1.);
}
else if (ge->geomType() == GEntity::DiscreteCurve){
c = Create_Curve(numEdge, MSH_SEGM_DISCRETE, 1, NULL, NULL, -1, -1, 0., 1.);
@@ -104,9 +104,12 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
else if(ge->geomType() == GEntity::CompoundCurve){
c = Create_Curve(numEdge, MSH_SEGM_COMPOUND, 1, NULL, NULL, -1, -1, 0., 1.);
std::vector<GEdge*> gec = ((GEdgeCompound*)ge)->getCompounds();
- for(int i = 0; i < gec.size(); i++)
+ for(unsigned int i = 0; i < gec.size(); i++)
c->compound.push_back(gec[i]->tag());
}
+ else{
+ Msg::Info("Unknown type of curve to add to planar face ...");
+ }
c->Control_Points = List_Create(2, 1, sizeof(Vertex *));
List_Add(c->Control_Points, &vertb);
@@ -117,7 +120,6 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
Tree_Add(gm->getGEOInternals()->Curves, &c);
CreateReversedCurve(c);
- List_Delete(temp);
}
List_Add(temp, &numEdge);
}
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index b7d4b7938c9b206889ad18513d126a669ce41da8..2a69d5ab0b7f1be76d0c4b63f42119fb514e1e36 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -3503,6 +3503,11 @@ void sortEdgesInLoop(int num, List_T *edges)
while(List_Nbr(edges) < nbEdges) {
for(int i = 0; i < List_Nbr(temp); i++) {
c2 = *(Curve **)List_Pointer(temp, i);
+ //reverse loop if not ordered correctly !
+ if (c1->end == c2->end){
+ Curve *c2R = CreateReversedCurve(c2);
+ c2 = c2R;
+ }
if(c1->end == c2->beg) {
List_Add(edges, &c2->Num);
List_PSuppress(temp, i);
@@ -3518,8 +3523,8 @@ void sortEdgesInLoop(int num, List_T *edges)
}
break;
}
- }
- if(j++ > nbEdges) {
+ }
+ if(j++ > nbEdges) {
Msg::Error("Line Loop %d is wrong", num);
break;
}
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 9019575d06a6db508fb1f3dd443ff94254830ada..48226300cca770f1b92a3f39e2046c56c72c9291 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -67,43 +67,51 @@ void OCCEdge::setTrimmed(OCCFace *f)
SPoint2 OCCEdge::reparamOnFace(const GFace *face, double epar, int dir) const
{
- const TopoDS_Face *s = (TopoDS_Face*) face->getNativePtr();
- double t0, t1;
- Handle(Geom2d_Curve) c2d;
- if(dir == 1){
- c2d = BRep_Tool::CurveOnSurface(c, *s, t0, t1);
+ if (face->getNativeType() != GEntity::OpenCascadeModel){
+ const GPoint pt = point(epar);
+ SPoint3 sp(pt.x(), pt.y(), pt.z());
+ return face->parFromPoint(sp);
}
else{
- c2d = BRep_Tool::CurveOnSurface(c_rev, *s, t0, t1);
- }
+ const TopoDS_Face *s = (TopoDS_Face*) face->getNativePtr();
+ double t0, t1;
+ Handle(Geom2d_Curve) c2d;
+
+ if(dir == 1){
+ c2d = BRep_Tool::CurveOnSurface(c, *s, t0, t1);
+ }
+ else{
+ c2d = BRep_Tool::CurveOnSurface(c_rev, *s, t0, t1);
+ }
- if(c2d.IsNull()){
- Msg::Fatal("Reparam on face failed: curve %d is not on surface %d",
- tag(), face->tag());
- }
+ if(c2d.IsNull()){
+ Msg::Fatal("Reparam on face failed: curve %d is not on surface %d",
+ tag(), face->tag());
+ }
- double u, v;
- gp_Pnt2d pnt = c2d->Value(epar);
- pnt.Coord(u, v);
-
- // sometimes OCC miserably fails ...
- GPoint p1 = point(epar);
- GPoint p2 = face->point(u, v);
- const double dx = p1.x()-p2.x();
- const double dy = p1.y()-p2.y();
- const double dz = p1.z()-p2.z();
- if(sqrt(dx * dx + dy * dy + dz * dz) > 1.e-2 * CTX::instance()->lc){
- Msg::Warning("Reparam on face was inaccurate for curve %d on surface %d at point %g",
- tag(), face->tag(), epar);
- Msg::Warning("On the face %d local (%g %g) global (%g %g %g)",
- face->tag(), u, v, p2.x(), p2.y(), p2.z());
- Msg::Warning("On the edge %d local (%g) global (%g %g %g)",
- tag(), epar, p1.x(), p1.y(), p1.z());
+ double u, v;
+ gp_Pnt2d pnt = c2d->Value(epar);
+ pnt.Coord(u, v);
+
+ // sometimes OCC miserably fails ...
+ GPoint p1 = point(epar);
+ GPoint p2 = face->point(u, v);
+ const double dx = p1.x()-p2.x();
+ const double dy = p1.y()-p2.y();
+ const double dz = p1.z()-p2.z();
+ if(sqrt(dx * dx + dy * dy + dz * dz) > 1.e-2 * CTX::instance()->lc){
+ Msg::Warning("Reparam on face was inaccurate for curve %d on surface %d at point %g",
+ tag(), face->tag(), epar);
+ Msg::Warning("On the face %d local (%g %g) global (%g %g %g)",
+ face->tag(), u, v, p2.x(), p2.y(), p2.z());
+ Msg::Warning("On the edge %d local (%g) global (%g %g %g)",
+ tag(), epar, p1.x(), p1.y(), p1.z());
+ }
+ return SPoint2(u, v);
}
- return SPoint2(u, v);
-}
+}
GPoint OCCEdge::closestPoint(const SPoint3 &qp, double ¶m) const
{
gp_Pnt pnt(qp.x(), qp.y(), qp.z());
@@ -128,15 +136,17 @@ GPoint OCCEdge::closestPoint(const SPoint3 &qp, double ¶m) const
// True if the edge is a seam for the given face
bool OCCEdge::isSeam(const GFace *face) const
{
- if (face->geomType() == GEntity::CompoundSurface)return false;
+
+ if (face->geomType() == GEntity::CompoundSurface) return false;
+ if (face->getNativeType() != GEntity::OpenCascadeModel) return false;
+
const TopoDS_Face *s = (TopoDS_Face*) face->getNativePtr();
BRepAdaptor_Surface surface(*s);
// printf("asking if edge %d is a seam of face %d\n",tag(),face->tag());
// printf("periodic %d %d\n",surface.IsUPeriodic(),surface.IsVPeriodic());
// if(surface.IsUPeriodic() || surface.IsVPeriodic()){
return BRep_Tool::IsClosed(c, *s);
- // }
- // return false;
+
}
GPoint OCCEdge::point(double par) const
diff --git a/Geo/OCCVertex.cpp b/Geo/OCCVertex.cpp
index 59db6a737474d4dd27199a89154fa16735b62c4e..cac9d8e916de4e456f1eb4b284b03a00cea18bfb 100644
--- a/Geo/OCCVertex.cpp
+++ b/Geo/OCCVertex.cpp
@@ -87,14 +87,22 @@ SPoint2 OCCVertex::reparamOnFace(const GFace *gf, int dir) const
while(it != l_edges.end()){
std::list<GEdge*> l = gf->edges();
if(std::find(l.begin(), l.end(), *it) != l.end()){
- const TopoDS_Face *s = (TopoDS_Face*)gf->getNativePtr();
- const TopoDS_Edge *c = (TopoDS_Edge*)(*it)->getNativePtr();
- double s1,s0;
- Handle(Geom2d_Curve) curve2d = BRep_Tool::CurveOnSurface(*c, *s, s0, s1);
- if((*it)->getBeginVertex() == this)
- return (*it)->reparamOnFace(gf, s0, dir);
- else if((*it)->getEndVertex() == this)
- return (*it)->reparamOnFace(gf, s1, dir);
+ if (gf->getNativeType() == GEntity::OpenCascadeModel){
+ const TopoDS_Face *s = (TopoDS_Face*)gf->getNativePtr();
+ const TopoDS_Edge *c = (TopoDS_Edge*)(*it)->getNativePtr();
+ double s1,s0;
+ Handle(Geom2d_Curve) curve2d = BRep_Tool::CurveOnSurface(*c, *s, s0, s1);
+ if((*it)->getBeginVertex() == this)
+ return (*it)->reparamOnFace(gf, s0, dir);
+ else if((*it)->getEndVertex() == this)
+ return (*it)->reparamOnFace(gf, s1, dir);
+ }
+ //if not OCCFace (OK for planar faces)
+ else{
+ const GPoint pt = point();
+ SPoint3 sp(pt.x(), pt.y(), pt.z());
+ gf->parFromPoint(sp);
+ }
}
++it;
}
diff --git a/Geo/gmshFace.h b/Geo/gmshFace.h
index b70242b3ef0c5398058098fe5b41653347cff5b9..26252b0bd27679ea7e9e742ecf9c0109921101c3 100644
--- a/Geo/gmshFace.h
+++ b/Geo/gmshFace.h
@@ -31,7 +31,7 @@ class gmshFace : public GFace {
virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
virtual GEntity::GeomType geomType() const;
ModelType getNativeType() const { return GmshModel; }
- void *getNativePtr() const { return s; }
+ void *getNativePtr() const { printf("coucuo here \n"); return s; }
virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true) const;
virtual void resetMeshAttributes();
};
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 543cd12645d806695c499e6c45cb719432c0765b..3f7838de4ab7a478cd0c2ae6cc647c51cd442794 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -296,6 +296,7 @@ static void remeshUnrecoveredEdges(std::map<MVertex*, BDS_Point*> &recoverMapInv
static bool algoDelaunay2D(GFace *gf)
{
+
if(!noSeam(gf))
return false;
@@ -2014,8 +2015,7 @@ void partitionAndRemesh(GFaceCompound *gf)
GFace *gfc = gf->model()->getFaceByTag(numf + NF + i );
meshGFace mgf;
mgf(gfc);
- //gfc->lloyd(20,0);
-
+
for(unsigned int j = 0; j < gfc->triangles.size(); ++j){
MTriangle *t = gfc->triangles[j];
std::vector<MVertex *> v(3);
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index 6d31391c1a19d590afbff1ad6dd3e8f9df08e66f..11e96bc3e7065f225f6a53405b5bd87741fce33f 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -16,7 +16,7 @@ static void increaseStencil(MVertex *v, v2t_cont &adj, std::vector<MElement*> &l
std::set<MElement*> stencil;
std::set<MVertex*> vs;
stencil.insert(lt.begin(),lt.end());
- for (int i=0;i<lt.size();i++){
+ for (unsigned int i=0;i<lt.size();i++){
for (int j=0;j<lt[i]->getNumVertices();j++){
MVertex *v1 = lt[i]->getVertex(j);
if (v1 != v){
@@ -36,7 +36,7 @@ meshMetric::meshMetric(GModel *gm){
if (_dim == 2){
for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); ++fit){
- for (int i=0;i<(*fit)->getNumMeshElements();i++){
+ for (unsigned int i=0;i<(*fit)->getNumMeshElements();i++){
MElement *e = (*fit)->getMeshElement(i);
MElement *copy = e->copy(_vertexMap, newP, newD);
_elements.push_back(copy);
@@ -45,7 +45,7 @@ meshMetric::meshMetric(GModel *gm){
}
else if (_dim == 3){
for (GModel::riter rit = gm->firstRegion(); rit != gm->lastRegion(); ++rit){
- for (int i=0;i<(*rit)->getNumMeshElements();i++){
+ for (unsigned int i=0;i<(*rit)->getNumMeshElements();i++){
MElement *e = (*rit)->getMeshElement(i);
MElement *copy = e->copy(_vertexMap, newP, newD);
_elements.push_back(copy);
@@ -62,7 +62,7 @@ meshMetric::meshMetric(std::vector<MElement*> elements){
std::map<MElement*, MElement*> newP;
std::map<MElement*, MElement*> newD;
- for (int i=0;i<elements.size();i++){
+ for (unsigned int i=0;i<elements.size();i++){
MElement *e = elements[i];
MElement *copy = e->copy(_vertexMap, newP, newD);
_elements.push_back(copy);
@@ -100,7 +100,7 @@ void meshMetric::intersectMetrics(){
_nodalMetrics[ver] = setOfMetrics[0][ver];
_detMetric[ver] = setOfDetMetric[0][ver];
_nodalSizes[ver] = setOfSizes[0][ver];
- for (int i=1;i<setOfMetrics.size();i++){
+ for (unsigned int i=1;i<setOfMetrics.size();i++){
_nodalMetrics[ver] = intersection_conserve_mostaniso(_nodalMetrics[ver],setOfMetrics[i][ver]);
_nodalSizes[ver] = std::min(_nodalSizes[ver],setOfSizes[i][ver]);
}
@@ -194,7 +194,7 @@ void meshMetric::exportInfo(const char * fileendname){
meshMetric::~meshMetric(){
if (_octree) delete _octree;
- for (int i=0; i< _elements.size(); i++)
+ for (unsigned int i=0; i< _elements.size(); i++)
delete _elements[i];
}
@@ -239,7 +239,7 @@ void meshMetric::computeHessian_FE(){
MVertex *ver = itv->first;
std::vector<MElement*> vElems= itv->second;
SVector3 gradv;
- for (int i=0;i<vElems.size();i++){
+ for (unsigned int i=0;i<vElems.size();i++){
MElement *e = vElems[i];
gradv += gradElems[e];
}
@@ -289,7 +289,7 @@ void meshMetric::computeHessian_FE(){
MVertex *ver = itv->first;
std::vector<MElement*> vElems= itv->second;
STensor3 hessv(0.0);
- for (int i=0;i<vElems.size();i++){
+ for (unsigned int i=0;i<vElems.size();i++){
MElement *e = vElems[i];
hessv += hessElems[e];
}
@@ -321,7 +321,7 @@ void meshMetric::computeHessian_LS( ){
fullVector<double> ATb (DIM);
fullVector<double> result (DIM);
fullMatrix<double> f (1,1);
- for(int i = 0; i < lt.size(); i++) {
+ for(unsigned int i = 0; i < lt.size(); i++) {
MElement *e = lt[i];
int npts; IntPt *pts;
SPoint3 p;
@@ -371,8 +371,8 @@ void meshMetric::computeMetric(){
//printf("%d elements are considered in the meshMetric \n",(int)_elements.size());
computeValues();
- computeHessian_FE();
- //computeHessian_LS();
+ //computeHessian_FE();
+ computeHessian_LS();
int metricNumber = setOfMetrics.size();