Commit 56c18b6d by Christophe Geuzaine

Merge branch 'MeshAdaptNew' into 'master'

Mesh adapt new

See merge request !27
parents a6fdc494 f5ff9ef9
Pipeline #363 passed with stage
in 11 minutes 43 seconds
......@@ -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;
......@@ -1109,7 +1101,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);
......@@ -1282,51 +1274,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;
......@@ -1459,7 +1406,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();
......
......@@ -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,
......
......@@ -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)
......
......@@ -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);
......
......@@ -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,
......
......@@ -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);
......
......@@ -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,6 +113,8 @@ 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)",
......@@ -138,6 +142,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,22 +401,6 @@ 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()){
......@@ -456,112 +452,13 @@ 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 &center) 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 +466,46 @@ bool OCCFace::isSphere (double &radius, SPoint3 &center) 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
......@@ -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 &param) const;
virtual double curvatures(const SPoint2 &param, 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 &center) 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
......
......@@ -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)
......
......@@ -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);
......
......@@ -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",