diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index e9e41e535a505747853a074192f9bef708df3db9..9a53fc187e6f6495e875707ce83d547dff6fff67 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -271,14 +271,6 @@ SOrientedBoundingBox GFace::getOBB()
return SOrientedBoundingBox(_obb);
}
-surface_params GFace::getSurfaceParams() const
-{
- surface_params p;
- p.radius = p.radius2 = p.height = p.cx = p.cy = p.cz = 0.;
- Msg::Error("Empty surface parameters for this type of surface");
- return p;
-}
-
std::vector<MVertex*> GFace::getEmbeddedMeshVertices() const
{
std::set<MVertex*> tmp;
@@ -1106,7 +1098,7 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[
#endif
}
-bool GFace::containsParam(const SPoint2 &pt) const
+bool GFace::containsParam(const SPoint2 &pt)
{
Range<double> uu = parBounds(0);
Range<double> vv = parBounds(1);
@@ -1279,51 +1271,6 @@ bool GFace::fillVertexArray(bool force)
return true;
}
-// by default we assume that straight lines are geodesics
-SPoint2 GFace::geodesic(const SPoint2 &pt1, const SPoint2 &pt2, double t)
-{
- if(CTX::instance()->mesh.secondOrderExperimental && geomType() != GEntity::Plane ){
- // FIXME: this is buggy -- remove the CTX option once we do it in
- // a robust manner
- GPoint gp1 = point(pt1.x(), pt1.y());
- GPoint gp2 = point(pt2.x(), pt2.y());
- SPoint2 guess = pt1 + (pt2 - pt1) * t;
- GPoint gp = closestPoint(SPoint3(gp1.x() + t * (gp2.x() - gp1.x()),
- gp1.y() + t * (gp2.y() - gp1.y()),
- gp1.z() + t * (gp2.z() - gp1.z())),
- (double*)guess);
- if (gp.g())
- return SPoint2(gp.u(), gp.v());
- else
- return pt1 + (pt2 - pt1) * t;
- }
- else{
- return pt1 + (pt2 - pt1) * t;
- }
-}
-
-
-// length of a curve drawn on a surface
-// S = (X(u,v), Y(u,v), Z(u,v) );
-// u = u(t) , v = v(t)
-// C = C ( u(t), v(t) )
-// dC/dt = dC/du du/dt + dC/dv dv/dt
-double GFace::length(const SPoint2 &pt1, const SPoint2 &pt2, int nbQuadPoints)
-{
- double *t = 0, *w = 0;
- double L = 0.0;
- gmshGaussLegendre1D(nbQuadPoints, &t, &w);
- for (int i = 0; i < nbQuadPoints; i++){
- const double ti = 0.5 * (1. + t[i]);
- SPoint2 pi = geodesic(pt1, pt2, ti);
- Pair<SVector3, SVector3> der2 = firstDer(pi);
- SVector3 der = der2.left() * (pt2.x() - pt1.x()) + der2.right() * (pt2.y() - pt1.y());
- const double d = norm(der);
- L += d * w[i] ;
- }
- return L;
-}
-
int GFace::genusGeom() const
{
int nSeams = 0;
@@ -1456,7 +1403,7 @@ void GFace::lloyd(int nbiter, int infn)
void GFace::replaceEdges(std::list<GEdge*> &new_edges)
{
- replaceEdgesInternal(new_edges);
+ // replaceEdgesInternal(new_edges);
std::list<GEdge*>::iterator it = l_edges.begin();
std::list<GEdge*>::iterator it2 = new_edges.begin();
std::list<int>::iterator it3 = l_dirs.begin();
diff --git a/Geo/GFace.h b/Geo/GFace.h
index a9ac73c2f704ed8060d8f014227540a41e112133..b55a4b71382452b1c63bb7ce005c25c21d69685c 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -26,11 +26,6 @@ class MPolygon;
class ExtrudeParams;
class GFaceCompound;
-struct surface_params
-{
- double radius, radius2, height, cx, cy, cz;
-};
-
class GRegion;
// A model face.
@@ -46,9 +41,6 @@ class GFace : public GEntity {
std::list<GVertex *> embedded_vertices;
GFaceCompound *compound; // this model face belongs to a compound
- // replace edges (for gluing) for specific modelers, we have to
- // re-create internal data
- virtual void replaceEdgesInternal(std::list<GEdge*> &){}
BoundaryLayerColumns _columns;
public: // this will become protected or private
@@ -154,9 +146,6 @@ class GFace : public GEntity {
// get the oriented bounding box
virtual SOrientedBoundingBox getOBB();
- // retrieve surface params
- virtual surface_params getSurfaceParams() const;
-
// compute the genus G of the surface
virtual int genusGeom() const;
virtual bool checkTopology() const { return true; }
@@ -165,14 +154,6 @@ class GFace : public GEntity {
virtual GPoint point(double par1, double par2) const = 0;
virtual GPoint point(const SPoint2 &pt) const { return point(pt.x(), pt.y()); }
- // compute, in parametric space, the interpolation from pt1 to pt2
- // along a geodesic of the surface
- virtual SPoint2 geodesic(const SPoint2 &pt1, const SPoint2 &pt2, double t);
-
- // compute the length of a geodesic between two points in parametric
- // space
- virtual double length(const SPoint2 &pt1, const SPoint2 &pt2, int n=4);
-
// if the mapping is a conforming mapping, i.e. a mapping that
// conserves angles, this function returns the eigenvalue of the
// metric at a given point this is a special feature for
@@ -190,7 +171,7 @@ class GFace : public GEntity {
virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true) const;
// true if the parameter value is interior to the face
- virtual bool containsParam(const SPoint2 &pt) const;
+ virtual bool containsParam(const SPoint2 &pt);
// return the point on the face closest to the given point
virtual GPoint closestPoint(const SPoint3 & queryPoint,
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index cd14f1a54a6e4549da3901932a182f25cc1d0e53..144e06aac163f3927932f2b5fbfdd570aa05e1fd 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -3561,242 +3561,6 @@ void GModel::setPhysicalNumToEntitiesInBox(int EntityDimension, int PhysicalNumb
setPhysicalNumToEntitiesInBox(EntityDimension, PhysicalNumber, box);
}
-static void computeDuplicates(GModel *model,
- std::multimap<GVertex*, GVertex*> &Unique2Duplicates,
- std::map<GVertex*, GVertex*> &Duplicates2Unique,
- const double &eps)
-{
- // FIXME: currently we use a greedy algorithm in n^2 (use e.g. MVertexRTree)
-
- // FIXME: add option to remove orphaned entities after duplicate check
- std::list<GVertex*> v;
- v.insert(v.begin(), model->firstVertex(), model->lastVertex());
-
- while(!v.empty()){
- GVertex *pv = *v.begin();
- v.erase(v.begin());
- bool found = false;
- for (std::multimap<GVertex*,GVertex*>::iterator it = Unique2Duplicates.begin();
- it != Unique2Duplicates.end(); ++it){
- GVertex *unique = it->first;
- const double d = sqrt((unique->x() - pv->x()) * (unique->x() - pv->x()) +
- (unique->y() - pv->y()) * (unique->y() - pv->y()) +
- (unique->z() - pv->z()) * (unique->z() - pv->z()));
- if (d <= eps && pv->geomType() == unique->geomType()) {
- found = true;
- Unique2Duplicates.insert(std::make_pair(unique, pv));
- Duplicates2Unique[pv] = unique;
- break;
- }
- }
- if (!found) {
- Unique2Duplicates.insert(std::make_pair(pv, pv));
- Duplicates2Unique[pv] = pv;
- }
- }
-}
-
-static void glueVerticesInEdges(GModel *model,
- std::multimap<GVertex*, GVertex*> &Unique2Duplicates,
- std::map<GVertex*, GVertex*> &Duplicates2Unique)
-{
- Msg::Debug("Gluing Edges");
- for (GModel::eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
- GEdge *e = *it;
- GVertex *v1 = e->getBeginVertex();
- GVertex *v2 = e->getEndVertex();
- GVertex *replacementOfv1 = Duplicates2Unique[v1];
- GVertex *replacementOfv2 = Duplicates2Unique[v2];
- if ((v1 != replacementOfv1) || (v2 != replacementOfv2)){
- Msg::Debug("Model Edge %d is re-build", e->tag());
- e->replaceEndingPoints (replacementOfv1, replacementOfv2);
- }
- }
-}
-
-static void computeDuplicates(GModel *model,
- std::multimap<GEdge*, GEdge*> &Unique2Duplicates,
- std::map<GEdge*,GEdge*> &Duplicates2Unique,
- const double &eps)
-{
- std::list<GEdge*> e;
- e.insert(e.begin(), model->firstEdge(), model->lastEdge());
-
- while(!e.empty()){
- GEdge *pe = *e.begin();
- e.erase(e.begin());
- bool found = false;
- for (std::multimap<GEdge*,GEdge*>::iterator it = Unique2Duplicates.begin();
- it != Unique2Duplicates.end(); ++it ){
- GEdge *unique = it->first;
- // first check edges that have same endpoints
- if (((unique->getBeginVertex() == pe->getBeginVertex() &&
- unique->getEndVertex() == pe->getEndVertex()) ||
- (unique->getEndVertex() == pe->getBeginVertex() &&
- unique->getBeginVertex() == pe->getEndVertex())) &&
- unique->geomType() == pe->geomType()){
- if ((unique->geomType() == GEntity::Line && pe->geomType() == GEntity::Line) ||
- unique->geomType() == GEntity::DiscreteCurve ||
- pe->geomType() == GEntity::DiscreteCurve ||
- unique->geomType() == GEntity::BoundaryLayerCurve ||
- pe->geomType() == GEntity::BoundaryLayerCurve){
- found = true;
- Unique2Duplicates.insert(std::make_pair(unique,pe));
- Duplicates2Unique[pe] = unique;
- break;
- }
- // compute a point
- Range<double> r = pe->parBounds(0);
- GPoint gp = pe->point(0.5 * (r.low() + r.high()));
- double t;
- GPoint gp2 = pe->closestPoint(SPoint3(gp.x(),gp.y(),gp.z()),t);
- const double d = sqrt((gp.x() - gp2.x()) * (gp.x() - gp2.x()) +
- (gp.y() - gp2.y()) * (gp.y() - gp2.y()) +
- (gp.z() - gp2.z()) * (gp.z() - gp2.z()));
- if (t >= r.low() && t <= r.high() && d <= eps) {
- found = true;
- Unique2Duplicates.insert(std::make_pair(unique,pe));
- Duplicates2Unique[pe] = unique;
- break;
- }
- }
- }
- if (!found) {
- Unique2Duplicates.insert(std::make_pair(pe,pe));
- Duplicates2Unique[pe] = pe;
- }
- }
-}
-
-static void glueEdgesInFaces(GModel *model,
- std::multimap<GEdge*, GEdge*> &Unique2Duplicates,
- std::map<GEdge*, GEdge*> &Duplicates2Unique)
-{
- Msg::Debug("Gluing Model Faces");
- for (GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it){
- GFace *f = *it;
- bool aDifferenceExists = false;
- std::list<GEdge*> old = f->edges(), enew;
- for (std::list<GEdge*>::iterator eit = old.begin(); eit !=old.end(); eit++){
- GEdge *temp = Duplicates2Unique[*eit];
- enew.push_back(temp);
- if (temp != *eit){
- aDifferenceExists = true;
- }
- }
- if (aDifferenceExists){
- Msg::Debug("Model Face %d is re-build", f->tag());
- f->replaceEdges(enew);
- }
- }
-}
-
-static void computeDuplicates(GModel *model,
- std::multimap<GFace*, GFace*> &Unique2Duplicates,
- std::map<GFace*,GFace*> &Duplicates2Unique,
- const double &eps)
-{
- std::list<GFace*> f;
- f.insert(f.begin(),model->firstFace(),model->lastFace());
-
- while(!f.empty()){
- GFace *pf = *f.begin();
- Range<double> r = pf->parBounds(0);
- Range<double> s = pf->parBounds(1);
- f.erase(f.begin());
- std::list<GEdge*> pf_edges = pf->edges();
- pf_edges.sort();
- bool found = false;
- for (std::multimap<GFace*,GFace*>::iterator it = Unique2Duplicates.begin();
- it != Unique2Duplicates.end(); ++it){
- GFace *unique = it->first;
- std::list<GEdge*> unique_edges = unique->edges();
- if (pf->geomType() == unique->geomType() &&
- unique_edges.size() == pf_edges.size()){
- unique_edges.sort();
- std::list<GEdge*>::iterator it_pf = pf_edges.begin();
- std::list<GEdge*>::iterator it_unique = unique_edges.begin();
- bool all_similar = true;
- // first check faces that have same edges
- for (; it_pf != pf_edges.end() ; ++it_pf,it_unique++){
- if (*it_pf != *it_unique) all_similar = false;
- }
- if (all_similar){
- if (unique->geomType() == GEntity::Plane && pf->geomType() == GEntity::Plane){
- found = true;
- Unique2Duplicates.insert(std::make_pair(unique,pf));
- Duplicates2Unique[pf] = unique;
- break;
- }
- double t[2]={0,0};
- // FIXME: evaluate a point on the surface (use e.g. buildRepresentationCross)
- const double d = 1.0;
- if (t[0] >= r.low() && t[0] <= r.high() &&
- t[1] >= s.low() && t[1] <= s.high() && d <= eps) {
- found = true;
- Unique2Duplicates.insert(std::make_pair(unique,pf));
- Duplicates2Unique[pf] = unique;
- break;
- }
- }
- }
- }
- if (!found) {
- Unique2Duplicates.insert(std::make_pair(pf,pf));
- Duplicates2Unique[pf] = pf;
- }
- }
-}
-
-static void glueFacesInRegions(GModel *model,
- std::multimap<GFace*, GFace*> &Unique2Duplicates,
- std::map<GFace*, GFace*> &Duplicates2Unique)
-{
- Msg::Debug("Gluing Regions");
- for (GModel::riter it = model->firstRegion(); it != model->lastRegion();++it){
- GRegion *r = *it;
- bool aDifferenceExists = false;
- std::list<GFace*> old = r->faces(), fnew;
- for (std::list<GFace*>::iterator fit = old.begin(); fit != old.end(); fit++){
- std::map<GFace*, GFace*>::iterator itR = Duplicates2Unique.find(*fit);
- if (itR == Duplicates2Unique.end()){
- Msg::Error("Error in the gluing process");
- return;
- }
- GFace *temp = itR->second;;
- fnew.push_back(temp);
- if (temp != *fit){
- aDifferenceExists = true;
- }
- }
- if (aDifferenceExists){
- Msg::Debug("Model Region %d is re-build", r->tag());
- r->replaceFaces (fnew);
- }
- }
-}
-
-void GModel::glue(double eps)
-{
- {
- std::multimap<GVertex*,GVertex*> Unique2Duplicates;
- std::map<GVertex*,GVertex*> Duplicates2Unique;
- computeDuplicates(this, Unique2Duplicates, Duplicates2Unique, eps);
- glueVerticesInEdges(this, Unique2Duplicates, Duplicates2Unique);
- }
- {
- std::multimap<GEdge*,GEdge*> Unique2Duplicates;
- std::map<GEdge*,GEdge*> Duplicates2Unique;
- computeDuplicates(this, Unique2Duplicates, Duplicates2Unique, eps);
- glueEdgesInFaces(this, Unique2Duplicates, Duplicates2Unique);
- }
- {
- std::multimap<GFace*,GFace*> Unique2Duplicates;
- std::map<GFace*,GFace*> Duplicates2Unique;
- computeDuplicates(this, Unique2Duplicates, Duplicates2Unique, eps);
- glueFacesInRegions(this, Unique2Duplicates, Duplicates2Unique);
- }
-}
GEdge *getNewModelEdge(GFace *gf1, GFace *gf2,
std::map<std::pair<int, int>, GEdge*> &newEdges)
diff --git a/Geo/GModel.h b/Geo/GModel.h
index b37ed3715174d4a968d805eae56ba4c54b2bdd43..bce806399c748ae0d89dd3fe0ee14dd9980124a4 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -493,7 +493,7 @@ class GModel {
// vertices that are too close, then merge edges, faces and
// regions). Warning: the gluer changes the geometric model, so that
// some pointers could become invalid.
- void glue(double eps);
+ // void glue(double eps);
// change the entity creation factory
void setFactory(std::string name);
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 5562e4f98132ef63eb7b6bfda76657efa7ffea32..23ac5442e5e62df8baffb697967d6c8bce9a0f2f 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -3285,8 +3285,9 @@ bool OCC_Internals::_makeFaceSTL(TopoDS_Face s,
#if (OCC_VERSION_MAJOR >= 7)
BRepMesh_FastDiscret::Parameters parameters;
parameters.Deflection = 0.1;
- parameters.Angle = 0.35;
- parameters.Relative = Standard_True;
+ parameters.Angle = 0.1;
+ // parameters.InternalVerticesMode = Standard_False;
+ parameters.Relative = Standard_False;
BRepMesh_FastDiscret aMesher(aBox, parameters);
#else
BRepMesh_FastDiscret aMesher(0.1, 0.35, aBox, Standard_False, Standard_False,
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index d2674dbc0097d3f359db1043fdbf714300340af5..ea980a29670301215a6d59462fdbf8dedbe96629 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -124,20 +124,34 @@ SPoint2 OCCEdge::reparamOnFace(const GFace *face, double epar, int dir) const
if (CTX::instance()->geom.reparamOnFaceRobust){
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();
+ double dx = p1.x()-p2.x();
+ double dy = p1.y()-p2.y();
+ double dz = p1.z()-p2.z();
if(sqrt(dx * dx + dy * dy + dz * dz) > CTX::instance()->geom.tolerance){
- Msg::Warning("Reparam on face was inaccurate for curve %d on surface %d at point %g",
+ Msg::Debug("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)",
+ Msg::Debug("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)",
+ Msg::Debug("On the edge %d local (%g) global (%g %g %g)",
tag(), epar, p1.x(), p1.y(), p1.z());
double guess [2] = {u,v};
GPoint pp = face->closestPoint(SPoint3(p1.x(),p1.y(),p1.z()), guess);
u = pp.u();
v = pp.v();
+
+ GPoint p2 = face->point(u, v);
+ dx = p1.x()-p2.x();
+ dy = p1.y()-p2.y();
+ dz = p1.z()-p2.z();
+ if(sqrt(dx * dx + dy * dy + dz * dz) > CTX::instance()->geom.tolerance){
+ Msg::Error("Closest Point was inaccurate for curve %d on surface %d at point %g",
+ tag(), face->tag(), epar);
+ Msg::Error("On the face %d local (%g %g) global (%g %g %g)",
+ face->tag(), u, v, p2.x(), p2.y(), p2.z());
+ Msg::Error("On the edge %d local (%g) global (%g %g %g)",
+ tag(), epar, p1.x(), p1.y(), p1.z());
+ }
+ else Msg::Info("OK");
}
}
return SPoint2(u, v);
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index c9628090372666b55bac8fc90159c88105961e01..a33daa00033738a85b03738427b92bc3bd51462e 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -13,6 +13,7 @@
#include "OCCFace.h"
#include "Numeric.h"
#include "Context.h"
+#include "robustPredicates.h"
#if defined(HAVE_OCC)
@@ -37,9 +38,10 @@
#include <TopExp_Explorer.hxx>
#include <TopoDS.hxx>
#include <gp_Pln.hxx>
+#include <gp_Sphere.hxx>
OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num)
- : GFace(m, num), s(_s)
+: GFace(m, num), s(_s),_radius(-1)
{
setup();
if(model()->getOCCInternals())
@@ -111,7 +113,11 @@ void OCCFace::setup()
BRepAdaptor_Surface surface(s);
_periodic[0] = surface.IsUPeriodic();
_periodic[1] = surface.IsVPeriodic();
+ if (_periodic[0]) _period[0]=surface.UPeriod();
+ if (_periodic[1]) _period[1]=surface.VPeriod();
+
+
ShapeAnalysis::GetFaceUVBounds(s, umin, umax, vmin, vmax);
Msg::Debug("OCC Face %d with %d parameter bounds (%g,%g)(%g,%g)",
tag(), l_edges.size(), umin, umax, vmin, vmax);
@@ -138,6 +144,14 @@ void OCCFace::setup()
embedded_vertices.push_back(v);
}
}
+
+ if (geomType()==GEntity::Sphere){
+ BRepAdaptor_Surface surface(s);
+ gp_Sphere sphere = surface.Sphere();
+ _radius = sphere.Radius();
+ gp_Pnt loc = sphere.Location();
+ _center = SPoint3(loc.X(),loc.Y(),loc.Z());
+ }
}
SBoundingBox3d OCCFace::bounds() const
@@ -389,24 +403,9 @@ bool OCCFace::containsPoint(const SPoint3 &pt) const
return false;
}
-surface_params OCCFace::getSurfaceParams() const
-{
- surface_params p;
- switch(geomType()){
- case GEntity::Cylinder:
- p.radius = Handle(Geom_CylindricalSurface)::DownCast(occface)->Radius();
- break;
- case GEntity::Sphere:
- p.radius = Handle(Geom_SphericalSurface)::DownCast(occface)->Radius();
- break;
- default:
- break;
- }
- return p;
-}
-
bool OCCFace::buildSTLTriangulation(bool force)
{
+
if(stl_triangles.size()){
if(force){
stl_vertices.clear();
@@ -456,112 +455,14 @@ bool OCCFace::buildSTLTriangulation(bool force)
return true;
}
-void OCCFace::replaceEdgesInternal(std::list<GEdge*> &new_edges)
-{
- // we simply replace old edges by new edges in the structure
- Handle(IntTools_Context) myContext = new IntTools_Context;
-
- // make a copy of s
- TopoDS_Face copy_of_s_forward = s;
- copy_of_s_forward.Orientation(TopAbs_FORWARD);
- // make a copy of occface
- TopLoc_Location location;
- Handle(Geom_Surface) copy_of_occface = BRep_Tool::Surface(copy_of_s_forward, location);
- // check periodicity
- bool bIsUPeriodic = _periodic[0];
- // get tolerance
- double tolerance = BRep_Tool::Tolerance(copy_of_s_forward);
-
- BRep_Builder aBB;
- TopoDS_Face newFace;
- aBB.MakeFace(newFace, copy_of_occface, location, tolerance);
- // expolore the face
- TopExp_Explorer aExpW, aExpE;
- aExpW.Init(copy_of_s_forward, TopAbs_WIRE);
- for (; aExpW.More(); aExpW.Next()) {
- TopoDS_Wire newWire;
- aBB.MakeWire(newWire);
- const TopoDS_Wire& aW=TopoDS::Wire(aExpW.Current());
- aExpE.Init(aW, TopAbs_EDGE);
- for (; aExpE.More(); aExpE.Next()) {
- const TopoDS_Edge& aE=TopoDS::Edge(aExpE.Current());
- std::list<GEdge*>::iterator it = l_edges.begin();
- std::list<GEdge*>::iterator it2 = new_edges.begin();
- TopoDS_Edge aER;
- Msg::Debug("trying to replace %d by %d",(*it)->tag(),(*it2)->tag());
- for ( ; it != l_edges.end(); ++it, ++it2){
- OCCEdge *occEd = dynamic_cast<OCCEdge*>(*it);
- TopoDS_Edge olde = occEd->getTopoDS_Edge();
- if (olde.IsSame(aE)){
- aER = *((TopoDS_Edge*)(*it2)->getNativePtr());
- }
- else {
- olde = occEd->getTopoDS_EdgeOld();
- if (olde.IsSame(aE)){
- aER = *((TopoDS_Edge*)(*it2)->getNativePtr());
- }
- }
- }
- if (aER.IsNull()){
- Msg::Error("cannot find an edge for gluing a face");
- }
- aER.Orientation(TopAbs_FORWARD);
- if (!BRep_Tool::Degenerated(aER)) {
- if (bIsUPeriodic) {
- Standard_Real aT1, aT2, aTx, aUx;
- BRep_Builder aBB_;
- Handle(Geom2d_Curve) aC2D =
- BRep_Tool::CurveOnSurface(aER, copy_of_s_forward, aT1, aT2);
- if (!aC2D.IsNull()) {
- if (BRep_Tool::IsClosed(aER, copy_of_s_forward)) {
- continue;
- }
- else{
- aTx = BOPTools_AlgoTools2D::IntermediatePoint(aT1, aT2);
- gp_Pnt2d aP2D;
- aC2D->D0(aTx, aP2D);
- aUx=aP2D.X();
- if (aUx < umin || aUx > umax) {
- // need to rebuild
- Handle(Geom2d_Curve) aC2Dx;
- aBB_.UpdateEdge(aER, aC2Dx, copy_of_s_forward , BRep_Tool::Tolerance(aE));
- }
- }
- }
- }
- BOPTools_AlgoTools2D::BuildPCurveForEdgeOnFace(aER, copy_of_s_forward);
- // orient image
- Standard_Boolean bIsToReverse =
- BOPTools_AlgoTools::IsSplitToReverse(aER, aE, myContext);
- if (bIsToReverse) {
- aER.Reverse();
- }
- }
- else {
- aER.Orientation(aE.Orientation());
- }
- aBB.Add(newWire, aER);
- }
- aBB.Add(newFace, newWire);
- }
- _replaced = s;
- s = newFace;
-
- setup();
- if(model()->getOCCInternals())
- model()->getOCCInternals()->bind(_replaced, tag());
-}
bool OCCFace::isSphere (double &radius, SPoint3 ¢er) const
{
switch(geomType()){
case GEntity::Sphere:
{
- radius = Handle(Geom_SphericalSurface)::DownCast(occface)->Radius();
- gp_Ax3 pos = Handle(Geom_SphericalSurface)::DownCast(occface)->Position();
- gp_Ax1 axis = pos.Axis();
- gp_Pnt loc = axis.Location();
- center = SPoint3(loc.X(),loc.Y(),loc.Z());
+ radius = _radius;
+ center = _center;
}
return true;
default:
@@ -569,4 +470,47 @@ bool OCCFace::isSphere (double &radius, SPoint3 ¢er) const
}
}
+bool OCCFace::containsParam(const SPoint2 &pt) {
+ // return GFace::containsParam(pt);
+ if (! buildSTLTriangulation(false) ) {
+ Msg::Warning ("Inacurate computation in OCCFace::containsParam");
+ return GFace::containsParam(pt);
+ }
+ // FILE *F = fopen("HOP.pos","w");
+ // fprintf(F,"View \" \"{\n");
+ /// fprintf(F,"SP(%g,%g,%g){2,2,2};\n",pt.x(),pt.y(),1.0);
+ SPoint2 mine = pt;
+
+ // bool ok = false;
+
+ for(unsigned int i = 0; i < stl_triangles.size(); i += 3){
+ SPoint2 gp1 = stl_vertices[stl_triangles[i]];
+ SPoint2 gp2 = stl_vertices[stl_triangles[i + 1]];
+ SPoint2 gp3 = stl_vertices[stl_triangles[i + 2]];
+
+ double s1 = robustPredicates::orient2d (gp1,gp2,mine);
+ double s2 = robustPredicates::orient2d (gp2,gp3,mine);
+ double s3 = robustPredicates::orient2d (gp3,gp1,mine);
+ /*
+ fprintf(F,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1};\n",
+ gp1.x(),gp1.y(),0.0,
+ gp2.x(),gp2.y(),0.0,
+ gp3.x(),gp3.y(),0.0);
+
+ printf("%g %g %g\n",s1,s2,s3);
+ */
+ if (s1*s2 >= 0 && s1*s3 >=0){
+ //ok = true;
+ return true;
+ }
+ }
+ // printf("gasp\n");
+ // fprintf(F,"};\n");
+ // fclose(F);
+ // return ok;
+ return false;
+
+}
+
+
#endif
diff --git a/Geo/OCCFace.h b/Geo/OCCFace.h
index b8d7850e7b5d2741ff9ad0db7792a8fde888f030..734b9032b420fccc91b177747bd353efb6c845c1 100644
--- a/Geo/OCCFace.h
+++ b/Geo/OCCFace.h
@@ -22,10 +22,10 @@ class OCCFace : public GFace {
Handle(Geom_Surface) occface;
double umin, umax, vmin, vmax;
bool _periodic[2];
+ double _period[2];
bool buildSTLTriangulation(bool force=false);
- void replaceEdgesInternal (std::list<GEdge*> &);
+ // void replaceEdgesInternal (std::list<GEdge*> &);
void setup();
- bool _isSphere;
double _radius;
SPoint3 _center;
public:
@@ -47,11 +47,14 @@ class OCCFace : public GFace {
virtual double curvatureMax(const SPoint2 ¶m) const;
virtual double curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
double *curvMax, double *curvMin) const;
- surface_params getSurfaceParams() const;
TopoDS_Face getTopoDS_Face () { return s; }
TopoDS_Face getTopoDS_FaceOld () { return _replaced; }
// tells if it's a sphere, and if it is, returns parameters
virtual bool isSphere (double &radius, SPoint3 ¢er) const;
+ virtual bool periodic(int dim) const { return _periodic[dim]; }
+ virtual double period(int dim) const { return _period[dim]; }
+ // true if the parameter value is interior to the face (taking into account boundaries)
+ virtual bool containsParam(const SPoint2 &pt) ;
};
#endif
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 68ca5d558581556b99f8f9fee40b7eeecff1debd..fdbbbaeb1dd5c6c7d64fdc4e1fc6e0fc5f68dc12 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -21,7 +21,7 @@
#endif
gmshFace::gmshFace(GModel *m, Surface *face)
- : GFace(m, face->Num), isSphere(false), radius(0.)
+ : GFace(m, face->Num)
{
resetNativePtr(face);
resetMeshAttributes();
@@ -105,7 +105,6 @@ void gmshFace::resetNativePtr(Surface *face)
// the bounding vertices)
if(s->Typ == MSH_SURF_PLAN) computeMeanPlane();
- isSphere = IsRuledSurfaceASphere(s, center, radius);
}
double gmshFace::getMetricEigenvalue(const SPoint2 &pt)
diff --git a/Geo/gmshFace.h b/Geo/gmshFace.h
index 86eae0a2e937921ee724920d64b481fc33ba8499..076d079c72a5e66f5e5bfaba133806cb09087d50 100644
--- a/Geo/gmshFace.h
+++ b/Geo/gmshFace.h
@@ -13,9 +13,6 @@ class Surface;
class gmshFace : public GFace {
protected:
Surface *s;
- bool isSphere;
- SPoint3 center;
- double radius;
bool buildSTLTriangulation(bool force);
public:
gmshFace(GModel *m, Surface *face);
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index a2a0a639960d7d1cf3f4ed67d4da97ef6771fd24..82a54fa1816e5d0f39bf85cdc262f4dd2110a018 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -103,13 +103,10 @@ void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace
gf->curvatureDiv(SPoint2(pts[3]->u, pts[3]->v)));
}
else{
- fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
pts[0]->X, pts[0]->Y, pts[0]->Z,
pts[1]->X, pts[1]->Y, pts[1]->Z,
- pts[2]->X, pts[2]->Y, pts[2]->Z,
- gf->curvatureDiv(SPoint2(pts[0]->u, pts[0]->v)),
- gf->curvatureDiv(SPoint2(pts[1]->u, pts[1]->v)),
- gf->curvatureDiv(SPoint2(pts[2]->u, pts[2]->v)));
+ pts[2]->X, pts[2]->Y, pts[2]->Z,pts[0]->iD, pts[1]->iD, pts[2]->iD);
}
}
}
@@ -202,7 +199,7 @@ BDS_Point *BDS_Mesh::add_point(int num, double x, double y, double z)
BDS_Point *BDS_Mesh::add_point(int num, double u, double v, GFace *gf)
{
- GPoint gp = gf->point(u*scalingU,v*scalingV);
+ GPoint gp = gf->point(u,v);
BDS_Point *pp = new BDS_Point(num, gp.x(), gp.y(), gp.z());
pp->u = u;
pp->v = v;
@@ -1227,7 +1224,9 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
}
/* TEST */
- bool isSphere = gf->geomType() == GEntity::Sphere;
+ double radius;
+ SPoint3 center;
+ bool isSphere = gf->isSphere(radius, center);
double XX=0,YY=0,ZZ=0;
double U = 0;
@@ -1269,11 +1268,13 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
V = gp.v();
}
else
- gp = gf->point(U * scalingU, V * scalingV);
+ gp = gf->point(U , V );
if (!gp.succeeded()){
return false;
}
+ // if (!gf->containsParam(SPoint2(U,V)))return false;
+
const double oldX = p->X;
const double oldY = p->Y;
const double oldZ = p->Z;
@@ -1299,23 +1300,30 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
// printf("%22.15E %22.15E\n", snew, sold);
if(snew < .1 * sold) return false;
- if(1 || test_quality){
- p->X = gp.x();
- p->Y = gp.y();
- p->Z = gp.z();
- newWorst = std::min(newWorst, qmTriangle::gamma(*it));
- double norm1[3],norm2[3];
- normal_triangle(n[0], n[1], n[2], norm1);
- p->X = oldX;
- p->Y = oldY;
- p->Z = oldZ;
- normal_triangle(n[0], n[1], n[2], norm2);
- oldWorst = std::min(oldWorst, qmTriangle::gamma(*it));
- double ps;
+ p->X = gp.x();
+ p->Y = gp.y();
+ p->Z = gp.z();
+ newWorst = std::min(newWorst, qmTriangle::gamma(*it));
+ double norm1[3],norm2[3];
+ normal_triangle(n[0], n[1], n[2], norm1);
+ p->X = oldX;
+ p->Y = oldY;
+ p->Z = oldZ;
+ normal_triangle(n[0], n[1], n[2], norm2);
+ oldWorst = std::min(oldWorst, qmTriangle::gamma(*it));
+ double ps;
+ if (isSphere){
+ double dx = center.x() - gp.x();
+ double dy = center.y() - gp.y();
+ double dz = center.z() - gp.z();
+ ps = dx*norm1[0]+dy*norm1[1]+dz*norm1[2];
+ if (ps < 0)return false;
+ }
+ else{
prosca(norm1, norm2, &ps);
- double threshold = (isSphere ? 0.95 : 0.5);
+ double threshold = 0.5;
if(ps < threshold){
- return false;
+ return false;
}
}
++it;
@@ -1386,7 +1394,7 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point *p, GFace *gf)
++it;
}
- GPoint gp = gf->point(U * scalingU, V * scalingV);
+ GPoint gp = gf->point(U, V);
if (!gp.succeeded())return false;
p->u = U;
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 4144e84a9a0ee2ff37be36f6546cc7ef09363173..c6aa237557710a92f97600a4ebbb39f97bfff355 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -405,8 +405,7 @@ class BDS_Mesh
public:
int MAXPOINTNUMBER, SNAP_SUCCESS, SNAP_FAILURE;
double Min[3], Max[3], LC;
- double scalingU, scalingV;
- BDS_Mesh(int _MAXX = 0) : MAXPOINTNUMBER(_MAXX),scalingU(1),scalingV(1){}
+ BDS_Mesh(int _MAXX = 0) : MAXPOINTNUMBER(_MAXX){}
void load(GVertex *gv); // load in BDS all the meshes of the vertex
void load(GEdge *ge); // load in BDS all the meshes of the edge
void load(GFace *gf); // load in BDS all the meshes of the surface
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 99fe9736fa15468f466e92313c9e315b2b34643b..f1dc9d3130c3d511991ecf623aa7e35188e21abd 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -599,7 +599,7 @@ void BDS2GMSH(BDS_Mesh *m, GFace *gf,
BDS_Point *p = *itp;
if(recoverMap.find(p) == recoverMap.end()){
MVertex *v = new MFaceVertex
- (p->X, p->Y, p->Z, gf, m->scalingU * p->u, m->scalingV * p->v);
+ (p->X, p->Y, p->Z, gf, p->u, p->v);
recoverMap[p] = v;
gf->mesh_vertices.push_back(v);
}
@@ -1097,8 +1097,6 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// Buid a BDS_Mesh structure that is convenient for doing the actual
// meshing procedure
BDS_Mesh *m = new BDS_Mesh;
- m->scalingU = 1;
- m->scalingV = 1;
std::vector<BDS_Point*> points(all_vertices.size());
SBoundingBox3d bbox;
@@ -1858,8 +1856,8 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
if(pp == 0){
double U, V;
SPoint2 param = coords[i];
- U = param.x() / m->scalingU;
- V = param.y() / m->scalingV;
+ U = param.x();
+ V = param.y();
pp = m->add_point(count + countTot, U, V, gf);
if (ge->dim() == 0){
pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
@@ -1912,8 +1910,6 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// Buid a BDS_Mesh structure that is convenient for doing the actual
// meshing procedure
BDS_Mesh *m = new BDS_Mesh;
- m->scalingU = 1;
- m->scalingV = 1;
std::vector<std::vector<BDS_Point*> > edgeLoops_BDS;
SBoundingBox3d bbox;
@@ -2306,7 +2302,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
BDS_Point *p = *itp;
if(recoverMap.find(p) == recoverMap.end()){
MVertex *v = new MFaceVertex
- (p->X, p->Y, p->Z, gf, m->scalingU * p->u, m->scalingV * p->v);
+ (p->X, p->Y, p->Z, gf, p->u, p->v);
recoverMap[p] = v;
gf->mesh_vertices.push_back(v);
}
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index c403a14215a974dfbdda855b59c05365efb7955b..16ccdeb59cef132730a3fe526d462428f881c485 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -35,19 +35,10 @@ double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2)
extern double lengthInfniteNorm(const double p[2], const double q[2],
const double quadAngle);
-inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f,
- double SCALINGU, double SCALINGV)
+inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f)
{
- // FIXME SUPER HACK
- // if (CTX::instance()->mesh.recombineAll || f->meshAttributes.recombine || 1){
- // double quadAngle = backgroundMesh::current()->getAngle (0.5 * (p1->u + p2->u) * SCALINGU, 0.5 * (p1->v + p2->v) * SCALINGV,0);
- // const double a [2] = {p1->u,p1->v};
- // const double b [2] = {p2->u,p2->v};
- // return lengthInfniteNorm(a,b, quadAngle);
- // }
-
- GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
- 0.5 * (p1->v + p2->v) * SCALINGV));
+ GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u),
+ 0.5 * (p1->v + p2->v)));
if (!GP.succeeded())
return computeEdgeLinearLength(p1,p2);
@@ -63,16 +54,15 @@ inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f,
return l1 + l2;
}
-inline double computeEdgeLinearLength_new(BDS_Point *p1, BDS_Point *p2, GFace *f,
- double SCALINGU, double SCALINGV)
+inline double computeEdgeLinearLength_new(BDS_Point *p1, BDS_Point *p2, GFace *f)
{
const int nbSb = 10;
GPoint GP[nbSb];
for (int i = 1; i < nbSb; i++){
double xi = (double)i / nbSb;
- GP[i-1] = f->point(SPoint2(((1-xi) * p1->u + xi * p2->u) * SCALINGU,
- ((1-xi) * p1->v + xi * p2->v) * SCALINGV));
+ GP[i-1] = f->point(SPoint2(((1-xi) * p1->u + xi * p2->u),
+ ((1-xi) * p1->v + xi * p2->v)));
if (!GP[i-1].succeeded())
return computeEdgeLinearLength(p1,p2);
}
@@ -100,14 +90,13 @@ inline double computeEdgeLinearLength_new(BDS_Point *p1, BDS_Point *p2, GFace *
return l;
}
-inline double computeEdgeMiddleCoord(BDS_Point *p1, BDS_Point *p2, GFace *f,
- double SCALINGU, double SCALINGV)
+inline double computeEdgeMiddleCoord(BDS_Point *p1, BDS_Point *p2, GFace *f)
{
if (f->geomType() == GEntity::Plane)
return 0.5;
- GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
- 0.5 * (p1->v + p2->v) * SCALINGV));
+ GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u),
+ 0.5 * (p1->v + p2->v)));
const double dx1 = p1->X - GP.x();
const double dy1 = p1->Y - GP.y();
@@ -124,14 +113,13 @@ inline double computeEdgeMiddleCoord(BDS_Point *p1, BDS_Point *p2, GFace *f,
return 0.25 * (3 * l2 - l1) / l2;
}
-inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f,
- double SCALINGU, double SCALINGV)
+inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f)
{
// FIXME !!!
if (f->geomType() == GEntity::Plane)
return e->length();
else
- return computeEdgeLinearLength(e->p1, e->p2, f, SCALINGU, SCALINGV);
+ return computeEdgeLinearLength(e->p1, e->p2, f);
}
double NewGetLc(BDS_Point *p)
@@ -141,8 +129,7 @@ double NewGetLc(BDS_Point *p)
std::min(p->lc(), p->lcBGM()) : p->lcBGM();
}
-static double correctLC_(BDS_Point *p1,BDS_Point *p2, GFace *f,
- double SCALINGU, double SCALINGV)
+static double correctLC_(BDS_Point *p1,BDS_Point *p2, GFace *f)
{
double l1 = NewGetLc(p1);
double l2 = NewGetLc(p2);
@@ -152,14 +139,11 @@ static double correctLC_(BDS_Point *p1,BDS_Point *p2, GFace *f,
double U = coord * p1->u + (1 - coord) * p2->u;
double V = coord * p1->v + (1 - coord) * p2->v;
- GPoint gpp = f->point(SCALINGU*U, SCALINGV*V);
+ GPoint gpp = f->point(U, V);
double lmid = BGM_MeshSize(f, U, V, gpp.x(), gpp.y(), gpp.z());
l = std::min(l, lmid);
if(CTX::instance()->mesh.lcFromCurvature){
- // GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
- // 0.5 * (p1->v + p2->v) * SCALINGV));
- // double l3 = BGM_MeshSize(f,GP.u(),GP.v(),GP.x(),GP.y(),GP.z());
double l3 = l;
double lcmin = std::min(std::min(l1, l2), l3);
l1 = std::min(lcmin*1.2,l1);
@@ -170,18 +154,18 @@ static double correctLC_(BDS_Point *p1,BDS_Point *p2, GFace *f,
return l;
}
-double NewGetLc(BDS_Edge *e, GFace *f, double SCALINGU, double SCALINGV)
+double NewGetLc(BDS_Edge *e, GFace *f)
{
- double linearLength = computeEdgeLinearLength(e, f, SCALINGU, SCALINGV);
- double l = correctLC_ (e->p1,e->p2,f, SCALINGU, SCALINGV);
+ double linearLength = computeEdgeLinearLength(e, f);
+ double l = correctLC_ (e->p1,e->p2,f);
//printf("BDS correct lc =%g lreal=%g \n", l,linearLength);
return linearLength / l;
}
-double NewGetLc(BDS_Point *p1, BDS_Point *p2, GFace *f, double su, double sv)
+double NewGetLc(BDS_Point *p1, BDS_Point *p2, GFace *f)
{
- double linearLength = computeEdgeLinearLength(p1, p2, f, su, sv);
- double l = correctLC_ (p1,p2,f, su, sv);
+ double linearLength = computeEdgeLinearLength(p1, p2, f);
+ double l = correctLC_ (p1,p2,f);
return linearLength / l;
}
@@ -199,7 +183,7 @@ void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
double sqr2 = sqrt(2.);
while (it != m.edges.end()){
if (!(*it)->deleted){
- double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
+ double lone = NewGetLc(*it, gf);
if (lone > oneoversqr2 && lone < sqr2) GS++;
avg += lone >1 ? (1. / lone) - 1. : lone - 1.;
max_e = std::max(max_e, lone);
@@ -213,7 +197,7 @@ void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
// SWAP TESTS i.e. tell if swap should be done
-bool edgeSwapTestAngle(BDS_Edge *e, double min_cos)
+static bool edgeSwapTestAngle(BDS_Edge *e, double min_cos)
{
BDS_Face *f1 = e->faces(0);
BDS_Face *f2 = e->faces(1);
@@ -229,7 +213,7 @@ bool edgeSwapTestAngle(BDS_Edge *e, double min_cos)
return cosa > min_cos;
}
-bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
+static bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
{
if(e->numfaces() != 2) return false;
@@ -268,9 +252,9 @@ bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
// bool smoothCouldSwap = !(cosb < cosa * .5);
double la = computeEdgeLinearLength(p11, p12);
- double la_ = computeEdgeLinearLength(p11, p12, gf, m.scalingU, m.scalingV);
+ double la_ = computeEdgeLinearLength(p11, p12, gf);
double lb = computeEdgeLinearLength(p13, p23);
- double lb_ = computeEdgeLinearLength(p13, p23, gf, m.scalingU, m.scalingV);
+ double lb_ = computeEdgeLinearLength(p13, p23, gf);
double LA = (la_ -la) / la_;
double LB = (lb_ -lb) / lb_;
@@ -300,7 +284,7 @@ bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
return false;
}
-bool edgeSwapTestDelaunay(BDS_Edge *e, GFace *gf)
+static bool edgeSwapTestDelaunay(BDS_Edge *e, GFace *gf)
{
BDS_Point *op[2];
@@ -324,16 +308,7 @@ bool edgeSwapTestDelaunay(BDS_Edge *e, GFace *gf)
return result > 1e-12;
}
-bool edgeSwapTestHighOrder(BDS_Edge *e,GFace *gf)
-{
- // must evaluate the swap with the perspectve of
- // the generation of 2 high order elements
- // The rationale is to consider the edges as
- // exactly matching curves and surfaces
- return false;
-}
-
-bool edgeSwapTestDelaunayAniso(BDS_Edge *e, GFace *gf, std::set<swapquad> &configs)
+static bool edgeSwapTestDelaunayAniso(BDS_Edge *e, GFace *gf, std::set<swapquad> &configs)
{
BDS_Point *op[2];
@@ -362,65 +337,48 @@ bool edgeSwapTestDelaunayAniso(BDS_Edge *e, GFace *gf, std::set<swapquad> &confi
return true;
}
-bool evalSwap(BDS_Edge *e, double &qa, double &qb)
+static int edgeSwapTest(GFace *gf, BDS_Edge *e)
{
BDS_Point *op[2];
- if(e->numfaces() != 2) return false;
- e->oppositeof (op);
- double qa1 = qmTriangle::gamma(e->p1, e->p2, op[0]);
- double qa2 = qmTriangle::gamma(e->p1, e->p2, op[1]);
- double qb1 = qmTriangle::gamma(e->p1, op[0], op[1]);
- double qb2 = qmTriangle::gamma(e->p2, op[0], op[1]);
- qa = std::min(qa1, qa2);
- qb = std::min(qb1, qb2);
- return true;
-}
-
-int edgeSwapTestQuality(BDS_Edge *e, double fact=1.1, bool force=false)
-{
- BDS_Point *op[2];
-
- if (!force)
- if(!e->p1->config_modified && ! e->p2->config_modified) return 0;
-
+ if(!e->p1->config_modified && ! e->p2->config_modified) return 0;
if(e->numfaces() != 2) return 0;
-
+
e->oppositeof (op);
- if (!force)
- if (!edgeSwapTestAngle(e, cos(CTX::instance()->mesh.allowSwapEdgeAngle * M_PI / 180.)))
- return -1;
-
+ if (!edgeSwapTestAngle(e, cos(CTX::instance()->mesh.allowSwapEdgeAngle * M_PI / 180.)))
+ return -1;
+
double qa1 = qmTriangle::gamma(e->p1, e->p2, op[0]);
double qa2 = qmTriangle::gamma(e->p1, e->p2, op[1]);
double qb1 = qmTriangle::gamma(e->p1, op[0], op[1]);
double qb2 = qmTriangle::gamma(e->p2, op[0], op[1]);
double qa = std::min(qa1, qa2);
double qb = std::min(qb1, qb2);
- if(qb > fact * qa) return 1;
- if(qb < qa / fact) return -1;
+ if(qb > qa) return 1;
+ if(qb < qa) return -1;
return 0;
}
void swapEdgePass(GFace *gf, BDS_Mesh &m, int &nb_swap)
{
+ return;
int NN1 = m.edges.size();
int NN2 = 0;
std::list<BDS_Edge*>::iterator it = m.edges.begin();
while (1){
if (NN2++ >= NN1) break;
- // result = -1 => forbid swap because too badly shaped elements
- // result = 0 => whatever
- // result = 1 => oblige to swap because the quality is greatly improved
if (!(*it)->deleted){
- double qual = (CTX::instance()->mesh.algo2d == ALGO_2D_MESHADAPT_OLD) ? 1 : 5;
- int result = edgeSwapTestQuality(*it, qual);
if (CTX::instance()->mesh.algo2d == ALGO_2D_MESHADAPT_OLD){
if (m.swap_edge(*it, BDS_SwapEdgeTestQuality(true))) nb_swap++;
}
- else if (result >= 0 && edgeSwapTestDelaunay(*it, gf)){
- if (m.swap_edge(*it, BDS_SwapEdgeTestQuality(false))) nb_swap++;
+ else if (edgeSwapTestDelaunay(*it, gf)){
+ int result = edgeSwapTest(gf,*it);
+ // result = -1 => forbid swap because too badly shaped elements
+ // result = 0 => whatever
+ // result = 1 => oblige to swap because the quality is greatly improved
+ if (result >= 0)
+ if (m.swap_edge(*it, BDS_SwapEdgeTestQuality(false))) nb_swap++;
}
}
++it;
@@ -429,6 +387,7 @@ void swapEdgePass(GFace *gf, BDS_Mesh &m, int &nb_swap)
void delaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap)
{
+ // return;
nb_swap = 0;
std::set<swapquad> configs;
while(1){
@@ -452,43 +411,6 @@ void delaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap)
}
}
-/*
-static void midpointsphere(GFace *gf, double u1, double v1, double u2, double v2,
- double &u12, double &v12, SPoint3 ¢er, double r)
-{
- GPoint p1 = gf->point(u1, v1);
- GPoint p2 = gf->point(u2, v2);
-
- SVector3 DIR ((p1.x()+p2.x())/2.0 - center.x(),
- (p1.y()+p2.y())/2.0 - center.y(),
- (p1.z()+p2.z())/2.0 - center.z());
- DIR.normalize();
-
- double X,Y,Z;
-
- double guess [2] = {0.5 * (u1 + u2), 0.5 * (v1 + v2)};
- u12 = guess[0];
- v12 = guess[1];
- if ( (v1 > 0.7*M_PI/2 && v2 > 0.7*M_PI/2) ||
- (v1 < -0.7*M_PI/2 && v2 < -0.7*M_PI/2) ){
- X = center.x() + DIR.x() * r;
- Y = center.y() + DIR.y() * r;
- Z = center.z() + DIR.z() * r;
- }
- else{
- return;
- X = center.x() - DIR.x() * r;
- Y = center.y() - DIR.y() * r;
- Z = center.z() - DIR.z() * r;
- }
-
- GPoint proj = gf->closestPoint(SPoint3(X, Y, Z), guess);
- if (proj.succeeded()){
- u12 = proj.u();
- v12 = proj.v();
- }
-}
-*/
bool edges_sort(std::pair<double, BDS_Edge*> a, std::pair<double, BDS_Edge*> b)
{
@@ -510,49 +432,39 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
{
std::list<BDS_Edge*>::iterator it = m.edges.begin();
std::vector<std::pair<double, BDS_Edge*> > edges;
-
+
while (it != m.edges.end()){
if(!(*it)->deleted && (*it)->numfaces() == 2 && (*it)->g->classif_degree == 2){
- double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
- if(lone > MAXE_){
- edges.push_back(std::make_pair(-lone, *it));
- }
+ double lone = NewGetLc(*it, gf);
+ if(lone > MAXE_)edges.push_back(std::make_pair(-lone, *it));
}
++it;
}
std::sort(edges.begin(), edges.end(), edges_sort);
+ bool isPeriodic = gf->periodic(0) || gf->periodic(1) ;
+ SPoint3 center;
+ double radius;
+ bool isSphere = gf->isSphere (radius, center);
+
for (unsigned int i = 0; i < edges.size(); ++i){
BDS_Edge *e = edges[i].second;
if (!e->deleted){
- const double coord = 0.5;
BDS_Point *mid ;
- double U, V;
- U = coord * e->p1->u + (1 - coord) * e->p2->u;
- V = coord * e->p1->v + (1 - coord) * e->p2->v;
-
- GPoint gpp = gf->point(m.scalingU*U,m.scalingV*V);
- if (gpp.succeeded()){
+ double U = 0.5*(e->p1->u+e->p2->u);
+ double V = 0.5*(e->p1->v+e->p2->v);
+ GPoint gpp = gf->point(U,V);
+
+ if ((!isPeriodic || gf->containsParam(SPoint2(U,V))) && gpp.succeeded()){
mid = m.add_point(++m.MAXPOINTNUMBER, gpp.x(),gpp.y(),gpp.z());
mid->u = U;
mid->v = V;
- if (backgroundMesh::current() && 0){
- mid->lc() = mid->lcBGM() =
- backgroundMesh::current()->operator()
- ((coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
- (coord * e->p1->v + (1 - coord) * e->p2->v)*m.scalingV,
- 0.0);
- }
- else {
- mid->lcBGM() = BGM_MeshSize
- (gf,
- (coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
- (coord * e->p1->v + (1 - coord) * e->p2->v)*m.scalingV,
- mid->X,mid->Y,mid->Z);
- mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
- }
- if(!m.split_edge(e, mid)) m.del_point(mid);
+ mid->lcBGM() = BGM_MeshSize (gf,U,V, mid->X,mid->Y,mid->Z);
+ mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
+
+ // if (isSphere && !canWeSplitAnEdgeOnASphere (e, mid, center,radius))m.del_point(mid);
+ if(!m.split_edge(e, mid)) m.del_point(mid);
else nb_split++;
}
}
@@ -579,7 +491,7 @@ double getMaxLcWhenCollapsingEdge(GFace *gf, BDS_Mesh &m, BDS_Edge *e, BDS_Point
}
if(!newP1 || !newP2) break; // error
BDS_Edge collapsedEdge = BDS_Edge(newP1, newP2);
- maxLc = std::max(maxLc, NewGetLc(&collapsedEdge, gf, m.scalingU, m.scalingV));
+ maxLc = std::max(maxLc, NewGetLc(&collapsedEdge, gf));
newP1->del(&collapsedEdge);
newP2->del(&collapsedEdge);
++eit;
@@ -595,7 +507,8 @@ void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_c
while (it != m.edges.end()){
if(!(*it)->deleted && (*it)->numfaces() == 2 && (*it)->g->classif_degree == 2){
- double lone = NewGetLc(*it, gf,m.scalingU,m.scalingV);
+
+ double lone = NewGetLc(*it, gf);
if(lone < MINE_){
edges.push_back (std::make_pair(lone, *it));
}
@@ -654,7 +567,7 @@ void collapseEdgePassUnSorted(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP,
if(NN2++ >= NN1) break;
if(!(*it)->deleted){
- double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
+ double lone = NewGetLc(*it, gf);
if(!(*it)->deleted && (*it)->numfaces() == 2 && lone < MINE_){
bool res = false;
if((*it)->p1->iD > MAXNP)
@@ -672,7 +585,7 @@ void collapseEdgePassUnSorted(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP,
void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q)
{
// FIXME SUPER HACK
- // return;
+ //return;
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while(itp != m.points.end()){
if(m.smooth_point_centroid(*itp, gf,q))
@@ -756,7 +669,7 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
if (!(*it)->deleted){
(*it)->p1->config_modified = false;
(*it)->p2->config_modified = false;
- double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
+ double lone = NewGetLc(*it, gf);
maxL = std::max(maxL, lone);
minL = std::min(minL, lone);
}
@@ -902,10 +815,13 @@ int solveInvalidPeriodic(GFace *gf, BDS_Mesh &m,
coord * e->p1->u + (1 - coord) * e->p2->u,
coord * e->p1->v + (1 - coord) * e->p2->v, gf);
mid->lcBGM() = BGM_MeshSize(gf,
- (coord * e->p1->u + (1 - coord) * e->p2->u) * m.scalingU,
- (coord * e->p1->v + (1 - coord) * e->p2->v) * m.scalingV,
+ (coord * e->p1->u + (1 - coord) * e->p2->u),
+ (coord * e->p1->v + (1 - coord) * e->p2->v),
mid->X, mid->Y, mid->Z);
mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
+
+ // printf("coucou\n");
+
if(!m.split_edge(e, mid)) m.del_point(mid);
}
}
diff --git a/tutorial/t17.geo b/tutorial/t17.geo
index c5213165f13f9bfb9d7aad430659b9ee41a1aa38..8734c0526a642aae800df11f0e18b68126ec0e9e 100644
--- a/tutorial/t17.geo
+++ b/tutorial/t17.geo
@@ -32,4 +32,7 @@ Merge "t17.pos";
Background Mesh View[0];
// Use bamg
+Mesh.SmoothRatio = 3;
+Mesh.AnisoMax = 1000;
+
Mesh.Algorithm=7;
\ No newline at end of file