There is also a number of more specific algorithms:
<ul>
-<li>\subpage prism_3d_algo_page "for meshing prismatic shapes"</li>
+<li>\subpage prism_3d_algo_page "for meshing prismatic 3D shapes"</li>
+<li>\subpage quad_from_ma_algo_page "for meshing faces with sinuous borders"</li>
<li>\subpage projection_algos_page "for meshing by projection of another mesh"</li>
<li>\subpage import_algos_page "for meshing by importing elements from another mesh"</li>
<li>\subpage radial_prism_algo_page "for meshing geometrical objects with cavities"</li>
--- /dev/null
+/*!
+
+\page quad_from_ma_algo_page Medial Axis Projection Quadrangle meshing algorithm
+
+Medial Axis Projection algorithm can be used for meshing faces with
+sinuous borders and having channel-like shape, for which is it
+difficult to define 1D hypotheses so that generated quadrangles to be
+of good shape.
+
+\image html quad_from_ma_mesh.png "A mesh of a river model"
+
+The algorithm assures good shape of quadrangles by constructing Medial
+Axis between sinuous borders of the face and using it to
+discretize the borders.
+
+\image html quad_from_ma_medial_axis.png "Media Axis between two blue sinuous borders"
+
+The Medial Axis is used in two ways:
+<ol>
+<li>If there is a sub-mesh on either sinuous border, then the nodes of
+ this border are mapped to the opposite border via the Medial
+ Axis.</li>
+<li> If there is no sub-meshes on the sinuous borders, then a part of
+ the Medial Axis that can be mapped to both borders is discretized
+ using a hypothesis assigned to the face or its ancestor shapes,
+ and the division points are mapped from the Medial Axis to the both
+ borders.</li>
+</ol>
+
+*/
{
};
+ /*!
+ * StdMeshers_QuadFromMedialAxis_1D2D: interface of "Quadrangle (Medial Axis Projection)" algorithm
+ */
+ interface StdMeshers_QuadFromMedialAxis_1D2D : SMESH::SMESH_2D_Algo
+ {
+ };
+
/*!
* StdMeshers_Hexa_3D: interface of "Hexahedron (i,j,k)" algorithm
*/
</python-wrap>
</algorithm>
+ <algorithm type ="QuadFromMedialAxis_1D2D"
+ label-id ="Quadrangle (Medial Axis Projection)"
+ icon-id ="mesh_algo_quad.png"
+ opt-hypos="ViscousLayers2D"
+ input ="EDGE"
+ output ="QUAD"
+ dim ="2">
+ <python-wrap>
+ <algo>QuadFromMedialAxis_1D2D=Quadrangle(algo=smeshBuilder.QUAD_MA_PROJ)</algo>
+ <hypo>ViscousLayers2D=ViscousLayers2D(SetTotalThickness(),SetNumberLayers(),SetStretchFactor(),SetIgnoreEdges())</hypo>
+ </python-wrap>
+ </algorithm>
+
<algorithm type ="Hexa_3D"
label-id ="Hexahedron (i,j,k)"
icon-id ="mesh_algo_hexa.png"
#include "SMDS_IteratorOnIterators.hxx"
#include "SMDS_VolumeTool.hxx"
#include "SMESH_Block.hxx"
+#include "SMESH_HypoFilter.hxx"
#include "SMESH_MeshAlgos.hxx"
#include "SMESH_ProxyMesh.hxx"
#include "SMESH_subMesh.hxx"
return TopAbs_SHAPE;
}
+//================================================================================
+/*!
+ * \brief Returns a shape, to which a hypothesis used to mesh a given shape is assigned
+ * \param [in] hyp - the hypothesis
+ * \param [in] shape - the shape, for meshing which the \a hyp is used
+ * \param [in] mesh - the mesh
+ * \return TopoDS_Shape - the shape the \a hyp is assigned to
+ */
+//================================================================================
+
+TopoDS_Shape SMESH_MesherHelper::GetShapeOfHypothesis( const SMESHDS_Hypothesis * hyp,
+ const TopoDS_Shape& shape,
+ SMESH_Mesh* mesh)
+{
+ const SMESH_Hypothesis* h = static_cast<const SMESH_Hypothesis*>( hyp );
+ SMESH_HypoFilter hypFilter( SMESH_HypoFilter::Is( h ));
+
+ TopoDS_Shape shapeOfHyp;
+ mesh->GetHypothesis( shape, hypFilter, /*checkAncestors=*/true, &shapeOfHyp );
+ return shapeOfHyp;
+}
+
//=======================================================================
//function : IsQuadraticMesh
//purpose : Check mesh without geometry for: if all elements on this shape are quadratic,
static TopAbs_ShapeEnum GetGroupType(const TopoDS_Shape& group,
const bool avoidCompound=false);
+ static TopoDS_Shape GetShapeOfHypothesis( const SMESHDS_Hypothesis * hyp,
+ const TopoDS_Shape& shape,
+ SMESH_Mesh* mesh);
+
public:
// ---------- PUBLIC INSTANCE METHODS ----------
// constructor
SMESH_MesherHelper(SMESH_Mesh& theMesh);
- SMESH_Mesh* GetMesh() const { return myMesh; }
+ SMESH_Gen* GetGen() const { return GetMesh()->GetGen(); }
+
+ SMESH_Mesh* GetMesh() const { return myMesh; }
SMESHDS_Mesh* GetMeshDS() const { return GetMesh()->GetMeshDS(); }
SMESH_Utils.hxx
SMESH_TryCatch.hxx
SMESH_MeshAlgos.hxx
+ SMESH_MAT2d.hxx
)
# --- sources ---
SMESH_TryCatch.cxx
SMESH_File.cxx
SMESH_MeshAlgos.cxx
+ SMESH_MAT2d.cxx
)
# --- rules ---
--- /dev/null
+// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File : SMESH_MAT2d.cxx
+// Created : Thu May 28 17:49:53 2015
+// Author : Edward AGAPOV (eap)
+
+#include "SMESH_MAT2d.hxx"
+
+#include <list>
+
+#include <BRepAdaptor_CompCurve.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepBuilderAPI_MakeVertex.hxx>
+#include <BRep_Builder.hxx>
+#include <BRep_Tool.hxx>
+#include <Bnd_B2d.hxx>
+//#include <GCPnts_AbscissaPoint.hxx>
+#include <GCPnts_TangentialDeflection.hxx>
+// #include <GCPnts_UniformAbscissa.hxx>
+// #include <GCPnts_UniformDeflection.hxx>
+#include <Geom2d_Curve.hxx>
+//#include <GeomAdaptor_Curve.hxx>
+#include <Geom2dAdaptor_Curve.hxx>
+#include <Geom_Curve.hxx>
+#include <Geom_Surface.hxx>
+#include <TopExp.hxx>
+#include <TopoDS_Vertex.hxx>
+#include <TopoDS_Wire.hxx>
+
+#ifdef _DEBUG_
+#include "SMESH_File.hxx"
+#include "SMESH_Comment.hxx"
+#endif
+
+using namespace std;
+using boost::polygon::x;
+using boost::polygon::y;
+using SMESH_MAT2d::TVD;
+using SMESH_MAT2d::TVDEdge;
+using SMESH_MAT2d::TVDCell;
+using SMESH_MAT2d::TVDVertex;
+
+namespace
+{
+ // Input data for construct_voronoi()
+ // -------------------------------------------------------------------------------------
+
+ struct InPoint
+ {
+ int _a, _b;
+ double _param;
+ InPoint(int x, int y, double param) : _a(x), _b(y), _param(param) {}
+ InPoint() : _a(0), _b(0), _param(0) {}
+
+ // working data
+ list< const TVDEdge* > _edges; // MA edges of a concave InPoint in CCW order
+
+ size_t index( const vector< InPoint >& inPoints ) const { return this - &inPoints[0]; }
+ bool operator==( const InPoint& other ) const { return _a == other._a && _b == other._b; }
+ };
+ // -------------------------------------------------------------------------------------
+
+ struct InSegment
+ {
+ InPoint * _p0;
+ InPoint * _p1;
+
+ // working data
+ size_t _geomEdgeInd; // EDGE index within the FACE
+ const TVDCell* _cell;
+ list< const TVDEdge* > _edges; // MA edges in CCW order within _cell
+
+ InSegment( InPoint * p0, InPoint * p1, size_t iE)
+ : _p0(p0), _p1(p1), _geomEdgeInd(iE) {}
+ InSegment() : _p0(0), _p1(0), _geomEdgeInd(0) {}
+
+ inline bool isConnected( const TVDEdge* edge );
+
+ inline bool isExternal( const TVDEdge* edge );
+
+ static void setGeomEdgeToCell( const TVDCell* cell, size_t eID ) { cell->color( eID ); }
+
+ static size_t getGeomEdge( const TVDCell* cell ) { return cell->color(); }
+ };
+
+ // check if a TVDEdge begins at my end or ends at my start
+ inline bool InSegment::isConnected( const TVDEdge* edge )
+ {
+ return ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&&
+ Abs( edge->vertex0()->y() - _p1->_b ) < 1. ) ||
+ (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&&
+ Abs( edge->vertex1()->y() - _p0->_b ) < 1. ));
+ }
+
+ // check if a MA TVDEdge is outside of a domain
+ inline bool InSegment::isExternal( const TVDEdge* edge )
+ {
+ double dot = // x1*x2 + y1*y2; (x1,y1) - internal normal of InSegment
+ ( _p0->_b - _p1->_b ) * ( 0.5 * ( edge->vertex0()->x() + edge->vertex1()->x() ) - _p0->_a ) +
+ ( _p1->_a - _p0->_a ) * ( 0.5 * ( edge->vertex0()->y() + edge->vertex1()->y() ) - _p0->_b );
+ return dot < 0.;
+ }
+
+ // // -------------------------------------------------------------------------------------
+ // const size_t theExternMA = 111; // to mark external MA edges
+
+ // bool isExternal( const TVDEdge* edge )
+ // {
+ // return ( SMESH_MAT2d::Branch::getBndSegment( edge ) == theExternMA );
+ // }
+
+ // // mark external MA edges
+ // void markExternalEdges( const TVDEdge* edge )
+ // {
+ // if ( isExternal( edge ))
+ // return;
+ // SMESH_MAT2d::Branch::setBndSegment( theExternMA, edge );
+ // SMESH_MAT2d::Branch::setBndSegment( theExternMA, edge->twin() );
+ // if ( edge->is_primary() && edge->vertex1() )
+ // {
+ // const TVDVertex * v = edge->vertex1();
+ // edge = v->incident_edge();
+ // do {
+ // markExternalEdges( edge );
+ // edge = edge->rot_next();
+ // } while ( edge != v->incident_edge() );
+ // }
+ // }
+
+ // -------------------------------------------------------------------------------------
+#ifdef _DEBUG_
+ // writes segments into a txt file readable by voronoi_visualizer
+ void inSegmentsToFile( vector< InSegment>& inSegments)
+ {
+ if ( inSegments.size() > 1000 )
+ return;
+ const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAdebug.txt";
+ SMESH_File file(fileName, false );
+ file.openForWriting();
+ SMESH_Comment text;
+ text << "0\n"; // nb points
+ text << inSegments.size() << "\n"; // nb segments
+ for ( size_t i = 0; i < inSegments.size(); ++i )
+ {
+ text << inSegments[i]._p0->_a << " "
+ << inSegments[i]._p0->_b << " "
+ << inSegments[i]._p1->_a << " "
+ << inSegments[i]._p1->_b << "\n";
+ }
+ text << "\n";
+ file.write( text.c_str(), text.size() );
+ cout << "Write " << fileName << endl;
+ }
+ void dumpEdge( const TVDEdge* edge )
+ {
+ cout << "*Edge_" << edge;
+ if ( !edge->vertex0() )
+ cout << " ( INF, INF";
+ else
+ cout << " ( " << edge->vertex0()->x() << ", " << edge->vertex0()->y();
+ if ( !edge->vertex1() )
+ cout << ") -> ( INF, INF";
+ else
+ cout << ") -> (" << edge->vertex1()->x() << ", " << edge->vertex1()->y();
+ cout << ")\t cell=" << edge->cell()
+ << " iBnd=" << edge->color()
+ << " twin=" << edge->twin()
+ << " twin_cell=" << edge->twin()->cell()
+ << " prev=" << edge->prev() << " next=" << edge->next()
+ << ( edge->is_primary() ? " MA " : " SCND" )
+ << ( edge->is_linear() ? " LIN " : " CURV" )
+ << endl;
+ }
+ void dumpCell( const TVDCell* cell )
+ {
+ cout << "**Cell_" << cell << " GEOM=" << cell->color() << " ";
+ cout << ( cell->contains_segment() ? " SEG " : " PNT " );
+ if ( cell-> is_degenerate() )
+ cout << " degen ";
+ else
+ {
+ cout << endl;
+ const TVDEdge* edge = cell->incident_edge();
+ size_t i = 0;
+ do {
+ edge = edge->next();
+ cout << " - " << ++i << " ";
+ dumpEdge( edge );
+ } while (edge != cell->incident_edge());
+ }
+ }
+#else
+ void inSegmentsToFile( vector< InSegment>& inSegments) {}
+ void dumpEdge( const TVDedge* edge ) {}
+ void dumpCell( const TVDCell* cell ) {}
+#endif
+}
+// -------------------------------------------------------------------------------------
+
+namespace boost {
+ namespace polygon {
+
+ template <>
+ struct geometry_concept<InPoint> {
+ typedef point_concept type;
+ };
+ template <>
+ struct point_traits<InPoint> {
+ typedef int coordinate_type;
+
+ static inline coordinate_type get(const InPoint& point, orientation_2d orient) {
+ return (orient == HORIZONTAL) ? point._a : point._b;
+ }
+ };
+
+ template <>
+ struct geometry_concept<InSegment> {
+ typedef segment_concept type;
+ };
+
+ template <>
+ struct segment_traits<InSegment> {
+ typedef int coordinate_type;
+ typedef InPoint point_type;
+
+ static inline point_type get(const InSegment& segment, direction_1d dir) {
+ return *(dir.to_int() ? segment._p1 : segment._p0);
+ }
+ };
+ } // namespace polygon
+} // namespace boost
+ // -------------------------------------------------------------------------------------
+
+namespace
+{
+ const int theNoBrachID = 0; // std::numeric_limits<int>::max();
+
+ // -------------------------------------------------------------------------------------
+ /*!
+ * \brief Intermediate DS to create InPoint's
+ */
+ struct UVU
+ {
+ gp_Pnt2d _uv;
+ double _u;
+ UVU( gp_Pnt2d uv, double u ): _uv(uv), _u(u) {}
+ InPoint getInPoint( double scale[2] )
+ {
+ return InPoint( int( _uv.X() * scale[0]), int( _uv.Y() * scale[1]), _u );
+ }
+ };
+ // -------------------------------------------------------------------------------------
+ /*!
+ * \brief A segment on EDGE, used to create BndPoints
+ */
+ struct BndSeg
+ {
+ InSegment* _inSeg;
+ const TVDEdge* _edge;
+ double _uLast;
+ int _branchID; // negative ID means reverse direction
+
+ BndSeg( InSegment* seg, const TVDEdge* edge, double u ):
+ _inSeg(seg), _edge(edge), _uLast(u), _branchID( theNoBrachID ) {}
+
+ void setIndexToEdge( size_t id )
+ {
+ SMESH_MAT2d::Branch::setBndSegment( id, _edge );
+ }
+
+ int branchID() const { return Abs( _branchID ); }
+
+ size_t geomEdge() const { return _inSeg->_geomEdgeInd; }
+
+ void setBranch( int branchID, vector< BndSeg >& bndSegs )
+ {
+ _branchID = branchID;
+
+ if ( _edge ) // pass branch to an opposite BndSeg
+ {
+ size_t oppSegIndex = SMESH_MAT2d::Branch::getBndSegment( _edge->twin() );
+ if ( oppSegIndex < bndSegs.size() /*&& bndSegs[ oppSegIndex ]._branchID == theNoBrachID*/ )
+ bndSegs[ oppSegIndex ]._branchID = -branchID;
+ }
+ }
+ bool hasOppositeEdge( const size_t noEdgeID )
+ {
+ if ( !_edge ) return false;
+ return ( _inSeg->getGeomEdge( _edge->twin()->cell() ) != noEdgeID );
+ }
+
+ // check a next segment in CW order
+ bool isSameBranch( const BndSeg& seg2 )
+ {
+ if ( !_edge || !seg2._edge )
+ return true;
+
+ const TVDCell* cell1 = this->_edge->twin()->cell();
+ const TVDCell* cell2 = seg2. _edge->twin()->cell();
+ if ( cell1 == cell2 )
+ return true;
+
+ const TVDEdge* edgeMedium1 = this->_edge->twin()->next();
+ const TVDEdge* edgeMedium2 = seg2. _edge->twin()->prev();
+
+ if ( edgeMedium1->is_secondary() && edgeMedium2->is_secondary() )
+ {
+ if ( edgeMedium1->twin() == edgeMedium2 )
+ return true;
+ // edgeMedium's are edges whose twin()->cell is built on an end point of inSegment
+ // and is located between cell1 and cell2
+ if ( edgeMedium1->twin() == edgeMedium2->twin() ) // is this possible???
+ return true;
+ if ( edgeMedium1->twin() == edgeMedium2->twin()->next() &&
+ edgeMedium1->twin()->cell()->contains_point() )
+ return true;
+ }
+ else if ( edgeMedium1->is_primary() && edgeMedium2->is_primary() )
+ {
+ if ( edgeMedium1->twin() == edgeMedium2 &&
+ SMESH_MAT2d::Branch::getBndSegment( edgeMedium1 ) ==
+ SMESH_MAT2d::Branch::getBndSegment( edgeMedium2 ))
+ // this is an ignored MA edge between inSegment's on one EDGE forming a convex corner
+ return true;
+ }
+
+ return false;
+ }
+ };
+
+ //================================================================================
+ /*!
+ * \brief Computes length of of TVDEdge
+ */
+ //================================================================================
+
+ double length( const TVDEdge* edge )
+ {
+ gp_XY d( edge->vertex0()->x() - edge->vertex1()->x(),
+ edge->vertex0()->y() - edge->vertex1()->y() );
+ return d.Modulus();
+ }
+
+ //================================================================================
+ /*!
+ * \brief Compute scale to have the same 2d proportions as in 3d
+ */
+ //================================================================================
+
+ void computeProportionScale( const TopoDS_Face& face,
+ const Bnd_B2d& uvBox,
+ double scale[2])
+ {
+ scale[0] = scale[1] = 1.;
+ if ( uvBox.IsVoid() ) return;
+
+ TopLoc_Location loc;
+ Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
+
+ const int nbDiv = 30;
+ gp_XY uvMin = uvBox.CornerMin(), uvMax = uvBox.CornerMax();
+ gp_XY uvMid = 0.5 * ( uvMin + uvMax );
+ double du = ( uvMax.X() - uvMin.X() ) / nbDiv;
+ double dv = ( uvMax.Y() - uvMin.Y() ) / nbDiv;
+
+ double uLen3d = 0, vLen3d = 0;
+ gp_Pnt uPrevP = surface->Value( uvMin.X(), uvMid.Y() );
+ gp_Pnt vPrevP = surface->Value( uvMid.X(), uvMin.Y() );
+ for (int i = 1; i <= nbDiv; i++)
+ {
+ double u = uvMin.X() + du * i;
+ double v = uvMin.Y() + dv * i;
+ gp_Pnt uP = surface->Value( u, uvMid.Y() );
+ gp_Pnt vP = surface->Value( uvMid.X(), v );
+ uLen3d += uP.Distance( uPrevP );
+ vLen3d += vP.Distance( vPrevP );
+ uPrevP = uP;
+ vPrevP = vP;
+ }
+ scale[0] = uLen3d / ( uvMax.X() - uvMin.X() );
+ scale[1] = vLen3d / ( uvMax.Y() - uvMin.Y() );
+ }
+
+ //================================================================================
+ /*!
+ * \brief Fill input data for construct_voronoi()
+ */
+ //================================================================================
+
+ bool makeInputData(const TopoDS_Face& face,
+ const std::vector< TopoDS_Edge >& edges,
+ const double minSegLen,
+ vector< InPoint >& inPoints,
+ vector< InSegment>& inSegments,
+ double scale[2])
+ {
+ const double theDiscrCoef = 0.5; // to decrease minSegLen for discretization
+ TopLoc_Location loc;
+
+ // discretize the EDGEs to get 2d points and segments
+
+ vector< vector< UVU > > uvuVec( edges.size() );
+ Bnd_B2d uvBox;
+ for ( size_t iE = 0; iE < edges.size(); ++iE )
+ {
+ vector< UVU > & points = uvuVec[ iE ];
+
+ double f,l;
+ Handle(Geom_Curve) c3d = BRep_Tool::Curve ( edges[ iE ], loc, f, l );
+ Handle(Geom2d_Curve) c2d = BRep_Tool::CurveOnSurface( edges[ iE ], face, f, l );
+ if ( c2d.IsNull() ) return false;
+
+ points.push_back( UVU( c2d->Value( f ), f ));
+ uvBox.Add( points.back()._uv );
+
+ Geom2dAdaptor_Curve c2dAdaptor (c2d, f,l );
+ double curDeflect = 0.3; //0.01; //Curvature deflection
+ double angDeflect = 0.2; // 0.09; //Angular deflection
+
+ GCPnts_TangentialDeflection discret(c2dAdaptor, angDeflect, curDeflect);
+ // if ( discret.NbPoints() > 2 )
+ // {
+ // cout << endl;
+ // do
+ // {
+ // discret.Initialize( c2dAdaptor, 100, curDeflect );
+ // cout << "C " << curDeflect << " " << discret.NbPoints() << endl;
+ // curDeflect *= 1.5;
+ // }
+ // while ( discret.NbPoints() > 5 );
+ // cout << endl;
+ // do
+ // {
+ // discret.Initialize( c2dAdaptor, angDeflect, 100 );
+ // cout << "A " << angDeflect << " " << discret.NbPoints() << endl;
+ // angDeflect *= 1.5;
+ // }
+ // while ( discret.NbPoints() > 5 );
+ // }
+ gp_Pnt p, pPrev;
+ if ( !c3d.IsNull() )
+ pPrev = c3d->Value( f );
+ for ( int i = 2; i <= discret.NbPoints(); i++ ) // skip the 1st point
+ {
+ double u = discret.Parameter(i);
+ if ( !c3d.IsNull() )
+ {
+ p = c3d->Value( u );
+ int nbDiv = int( p.Distance( pPrev ) / minSegLen / theDiscrCoef );
+ double dU = ( u - points.back()._u ) / nbDiv;
+ for ( int iD = 1; iD < nbDiv; ++iD )
+ {
+ double uD = points.back()._u + dU;
+ points.push_back( UVU( c2d->Value( uD ), uD ));
+ }
+ pPrev = p;
+ }
+ points.push_back( UVU( c2d->Value( u ), u ));
+ uvBox.Add( points.back()._uv );
+ }
+ // if ( !c3d.IsNull() )
+ // {
+ // vector<double> params;
+ // GeomAdaptor_Curve c3dAdaptor( c3d,f,l );
+ // if ( useDefl )
+ // {
+ // const double deflection = minSegLen * 0.1;
+ // GCPnts_UniformDeflection discret( c3dAdaptor, deflection, f, l, true );
+ // if ( !discret.IsDone() )
+ // return false;
+ // int nbP = discret.NbPoints();
+ // for ( int i = 2; i < nbP; i++ ) // skip 1st and last points
+ // params.push_back( discret.Parameter(i) );
+ // }
+ // else
+ // {
+ // double eLen = GCPnts_AbscissaPoint::Length( c3dAdaptor );
+ // int nbSeg = Max( 1, int( eLen / minSegLen / theDiscrCoef ));
+ // double segLen = eLen / nbSeg;
+ // GCPnts_UniformAbscissa discret( c3dAdaptor, segLen, f, l );
+ // int nbP = Min( discret.NbPoints(), nbSeg + 1 );
+ // for ( int i = 2; i < nbP; i++ ) // skip 1st and last points
+ // params.push_back( discret.Parameter(i) );
+ // }
+ // for ( size_t i = 0; i < params.size(); ++i )
+ // {
+ // points.push_back( UVU( c2d->Value( params[i] ), params[i] ));
+ // uvBox.Add( points.back()._uv );
+ // }
+ // }
+ if ( points.size() < 2 )
+ {
+ points.push_back( UVU( c2d->Value( l ), l ));
+ uvBox.Add( points.back()._uv );
+ }
+ if ( edges[ iE ].Orientation() == TopAbs_REVERSED )
+ std::reverse( points.begin(), points.end() );
+ }
+
+ // make connected EDGEs have same UV at shared VERTEX
+ TopoDS_Vertex vShared;
+ for ( size_t iE = 0; iE < edges.size(); ++iE )
+ {
+ size_t iE2 = (iE+1) % edges.size();
+ if ( !TopExp::CommonVertex( edges[iE], edges[iE2], vShared ))
+ continue;
+ if ( !vShared.IsSame( TopExp::LastVertex( edges[iE], true )))
+ return false;
+ vector< UVU > & points1 = uvuVec[ iE ];
+ vector< UVU > & points2 = uvuVec[ iE2 ];
+ gp_Pnt2d & uv1 = points1.back() ._uv;
+ gp_Pnt2d & uv2 = points2.front()._uv;
+ uv1 = uv2 = 0.5 * ( uv1.XY() + uv2.XY() );
+ }
+
+ // get scale to have the same 2d proportions as in 3d
+ computeProportionScale( face, uvBox, scale );
+
+ // make scale to have coordinates precise enough when converted to int
+
+ gp_XY uvMin = uvBox.CornerMin(), uvMax = uvBox.CornerMax();
+ uvMin.ChangeCoord(1) = uvMin.X() * scale[0];
+ uvMin.ChangeCoord(2) = uvMin.Y() * scale[1];
+ uvMax.ChangeCoord(1) = uvMax.X() * scale[0];
+ uvMax.ChangeCoord(2) = uvMax.Y() * scale[1];
+ double vMax[2] = { Max( Abs( uvMin.X() ), Abs( uvMax.X() )),
+ Max( Abs( uvMin.Y() ), Abs( uvMax.Y() )) };
+ int iMax = ( vMax[0] > vMax[1] ) ? 0 : 1;
+ const double precision = 1e-5;
+ double preciScale = Min( vMax[iMax] / precision,
+ std::numeric_limits<int>::max() / vMax[iMax] );
+ preciScale /= scale[iMax];
+ double roundedScale = 10; // to ease debug
+ while ( roundedScale * 10 < preciScale )
+ roundedScale *= 10.;
+ scale[0] *= roundedScale;
+ scale[1] *= roundedScale;
+
+ // create input points and segments
+
+ inPoints.clear();
+ inSegments.clear();
+ size_t nbPnt = 0;
+ for ( size_t iE = 0; iE < uvuVec.size(); ++iE )
+ nbPnt += uvuVec[ iE ].size();
+ inPoints.resize( nbPnt );
+ inSegments.reserve( nbPnt );
+
+ size_t iP = 0;
+ if ( face.Orientation() == TopAbs_REVERSED )
+ {
+ for ( int iE = uvuVec.size()-1; iE >= 0; --iE )
+ {
+ vector< UVU > & points = uvuVec[ iE ];
+ inPoints[ iP++ ] = points.back().getInPoint( scale );
+ for ( size_t i = points.size()-1; i >= 1; --i )
+ {
+ inPoints[ iP++ ] = points[i-1].getInPoint( scale );
+ inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
+ }
+ }
+ }
+ else
+ {
+ for ( size_t iE = 0; iE < uvuVec.size(); ++iE )
+ {
+ vector< UVU > & points = uvuVec[ iE ];
+ inPoints[ iP++ ] = points[0].getInPoint( scale );
+ for ( size_t i = 1; i < points.size(); ++i )
+ {
+ inPoints[ iP++ ] = points[i].getInPoint( scale );
+ inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
+ }
+ }
+ }
+ return true;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Create MA branches and FACE boundary data
+ * \param [in] vd - voronoi diagram of \a inSegments
+ * \param [in] inPoints - FACE boundary points
+ * \param [in,out] inSegments - FACE boundary segments
+ * \param [out] branch - MA branches to fill
+ * \param [out] branchEnd - ends of MA branches to fill
+ * \param [out] boundary - FACE boundary to fill
+ */
+ //================================================================================
+
+ void makeMA( const TVD& vd,
+ vector< InPoint >& inPoints,
+ vector< InSegment > & inSegments,
+ vector< SMESH_MAT2d::Branch >& branch,
+ vector< const SMESH_MAT2d::BranchEnd* >& branchPnt,
+ SMESH_MAT2d::Boundary& boundary )
+ {
+ const size_t noEdgeID = inSegments.size() + 1; // ID of non-existent geom EDGE
+
+ // Associate MA cells with inSegments
+ for (TVD::const_cell_iterator it = vd.cells().begin(); it != vd.cells().end(); ++it)
+ {
+ const TVDCell* cell = &(*it);
+ if ( cell->contains_segment() )
+ {
+ InSegment& seg = inSegments[ cell->source_index() ];
+ seg._cell = cell;
+ seg.setGeomEdgeToCell( cell, seg._geomEdgeInd );
+ }
+ else
+ {
+ InSegment::setGeomEdgeToCell( cell, noEdgeID );
+ }
+ }
+
+ vector< bool > inPntChecked( inPoints.size(), false );
+
+ // Find MA edges of each inSegment
+
+ for ( size_t i = 0; i < inSegments.size(); ++i )
+ {
+ InSegment& inSeg = inSegments[i];
+
+ // get edges around the cell lying on MA
+ bool hasSecondary = false;
+ const TVDEdge* edge = inSeg._cell->incident_edge();
+ do {
+ edge = edge->next(); // Returns the CCW next edge within the cell.
+ if ( edge->is_primary() && !inSeg.isExternal( edge ) )
+ inSeg._edges.push_back( edge ); // edge equidistant from two InSegments
+ else
+ hasSecondary = true;
+ } while (edge != inSeg._cell->incident_edge());
+
+ // there can be several continuous MA edges but maEdges can begin in the middle of
+ // a chain of continuous MA edges. Make the chain continuous.
+ list< const TVDEdge* >& maEdges = inSeg._edges;
+ if ( maEdges.empty() )
+ continue;
+ if ( hasSecondary )
+ while ( maEdges.back()->next() == maEdges.front() )
+ maEdges.splice( maEdges.end(), maEdges, maEdges.begin() );
+
+ // remove maEdges equidistant from two neighbor InSegments of the same geom EDGE
+ list< const TVDEdge* >::iterator e = maEdges.begin();
+ while ( e != maEdges.end() )
+ {
+ const TVDCell* cell2 = (*e)->twin()->cell(); // cell on the other side of a MA edge
+ size_t geoE2 = InSegment::getGeomEdge( cell2 );
+ bool toRemove = ( inSeg._geomEdgeInd == geoE2 && inSeg.isConnected( *e ));
+ if ( toRemove )
+ e = maEdges.erase( e );
+ else
+ ++e;
+ }
+ if ( maEdges.empty() )
+ continue;
+
+ // add MA edges corresponding to concave InPoints
+ for ( int is2nd = 0; is2nd < 2; ++is2nd ) // loop on two ends of inSeg
+ {
+ InPoint& inPnt = *( is2nd ? inSeg._p1 : inSeg._p0 );
+ size_t pInd = inPnt.index( inPoints );
+ if ( inPntChecked[ pInd ] )
+ continue;
+ if ( pInd > 0 &&
+ inPntChecked[ pInd-1 ] &&
+ inPoints[ pInd-1 ] == inPnt )
+ continue;
+ inPntChecked[ pInd ] = true;
+
+ const TVDEdge* edge = // a TVDEdge passing through an end of inSeg
+ is2nd ? maEdges.front()->prev() : maEdges.back()->next();
+ while ( true )
+ {
+ if ( edge->is_primary() ) break; // this should not happen
+ const TVDEdge* edge2 = edge->twin(); // we are in a neighbor cell, add MA edges to inPnt
+ if ( inSeg.getGeomEdge( edge2->cell() ) != noEdgeID )
+ break; // cell of an InSegment
+ bool hasInfinite = false;
+ list< const TVDEdge* > pointEdges;
+ edge = edge2;
+ do
+ {
+ edge = edge->next(); // Returns the CCW next edge within the cell.
+ if ( edge->is_infinite() )
+ hasInfinite = true;
+ else if ( edge->is_primary() && !inSeg.isExternal( edge ))
+ pointEdges.push_back( edge );
+ }
+ while ( edge != edge2 && !hasInfinite );
+
+ if ( hasInfinite || pointEdges.empty() )
+ break;
+ inPnt._edges.splice( inPnt._edges.end(), pointEdges );
+ inSeg.setGeomEdgeToCell( edge->cell(), inSeg._geomEdgeInd );
+
+ edge = is2nd ? inPnt._edges.front()->prev() : inPnt._edges.back()->next();
+ }
+ } // add MA edges corresponding to concave InPoints
+
+ } // loop on inSegments to find corresponding MA edges
+
+
+ // -------------------------------------------
+ // Create Branches and BndPoints for each EDGE
+ // -------------------------------------------
+
+ if ( inPoints.front() == inPoints.back() /*&& !inPoints[0]._edges.empty()*/ )
+ {
+ inPntChecked[0] = false; // do not use the 1st point twice
+ //InSegment::setGeomEdgeToCell( inPoints[0]._edges.back()->cell(), noEdgeID );
+ inPoints[0]._edges.clear();
+ }
+
+ // Divide InSegment's into BndSeg's
+
+ vector< BndSeg > bndSegs;
+ bndSegs.reserve( inSegments.size() * 3 );
+
+ list< const TVDEdge* >::reverse_iterator e;
+ for ( size_t i = 0; i < inSegments.size(); ++i )
+ {
+ InSegment& inSeg = inSegments[i];
+
+ // segments around 1st concave point
+ size_t ip0 = inSeg._p0->index( inPoints );
+ if ( inPntChecked[ ip0 ] )
+ for ( e = inSeg._p0->_edges.rbegin(); e != inSeg._p0->_edges.rend(); ++e )
+ bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p0->_param ));
+ inPntChecked[ ip0 ] = false;
+
+ // segments of InSegment's
+ size_t nbMaEdges = inSeg._edges.size();
+ switch ( nbMaEdges ) {
+ case 0: // "around" circle center
+ bndSegs.push_back( BndSeg( &inSeg, 0, inSeg._p1->_param )); break;
+ case 1:
+ bndSegs.push_back( BndSeg( &inSeg, inSeg._edges.back(), inSeg._p1->_param )); break;
+ default:
+ vector< double > len;
+ len.push_back(0);
+ for ( e = inSeg._edges.rbegin(); e != inSeg._edges.rend(); ++e )
+ len.push_back( len.back() + length( *e ));
+
+ e = inSeg._edges.rbegin();
+ for ( size_t l = 1; l < len.size(); ++e, ++l )
+ {
+ double dl = len[l] / len.back();
+ double u = dl * inSeg._p1->_param + ( 1. - dl ) * inSeg._p0->_param;
+ bndSegs.push_back( BndSeg( &inSeg, *e, u ));
+ }
+ }
+ // segments around 2nd concave point
+ size_t ip1 = inSeg._p1->index( inPoints );
+ if ( inPntChecked[ ip1 ] )
+ for ( e = inSeg._p1->_edges.rbegin(); e != inSeg._p1->_edges.rend(); ++e )
+ bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param ));
+ inPntChecked[ ip1 ] = false;
+ }
+
+ // make TVDEdge's know it's BndSeg to enable passing branchID to
+ // an opposite BndSeg in BndSeg::setBranch()
+ for ( size_t i = 0; i < bndSegs.size(); ++i )
+ bndSegs[i].setIndexToEdge( i );
+
+
+ // Find TVDEdge's of Branches and associate them with bndSegs
+
+ vector< vector<const TVDEdge*> > branchEdges;
+ branchEdges.reserve( boundary.nbEdges() * 4 );
+
+ map< const TVDVertex*, SMESH_MAT2d::BranchEndType > endType;
+
+ int branchID = 1; // we code orientation as branchID sign
+ branchEdges.resize( branchID + 1 );
+
+ size_t i1st = 0;
+ while ( i1st < bndSegs.size() && !bndSegs[i1st].hasOppositeEdge( noEdgeID ))
+ ++i1st;
+ bndSegs[i1st].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg
+ branchEdges[ branchID ].push_back( bndSegs[i1st]._edge );
+
+ for ( size_t i = i1st+1; i < bndSegs.size(); ++i )
+ {
+ if ( bndSegs[i].branchID() )
+ {
+ branchID = bndSegs[i]._branchID; // with sign
+
+ if ( bndSegs[i]._branchID == -bndSegs[i-1]._branchID &&
+ bndSegs[i]._edge )
+ {
+ SMESH_MAT2d::BranchEndType type =
+ ( bndSegs[i]._inSeg->isConnected( bndSegs[i]._edge ) ?
+ SMESH_MAT2d::BE_ON_VERTEX :
+ SMESH_MAT2d::BE_END );
+ endType.insert( make_pair( bndSegs[i]._edge->vertex1(), type ));
+ }
+ continue;
+ }
+ if ( !bndSegs[i-1].isSameBranch( bndSegs[i] ))
+ {
+ branchEdges.resize(( branchID = branchEdges.size()) + 1 );
+ if ( bndSegs[i]._edge )
+ endType.insert( make_pair( bndSegs[i]._edge->vertex1(),
+ SMESH_MAT2d::BE_BRANCH_POINT ));
+ }
+ bndSegs[i].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg
+ if ( bndSegs[i].hasOppositeEdge( noEdgeID ))
+ branchEdges[ bndSegs[i].branchID() ].push_back( bndSegs[i]._edge );
+ }
+ // define BranchEndType of the first TVDVertex
+ if ( bndSegs.front()._branchID == -bndSegs.back()._branchID )
+ {
+ if ( bndSegs[0]._edge )
+ {
+ SMESH_MAT2d::BranchEndType type =
+ ( bndSegs[0]._inSeg->isConnected( bndSegs[0]._edge ) ?
+ SMESH_MAT2d::BE_ON_VERTEX :
+ SMESH_MAT2d::BE_END );
+ endType.insert( make_pair( bndSegs[0]._edge->vertex1(), type ));
+ }
+ else if ( bndSegs.back()._edge )
+ {
+ SMESH_MAT2d::BranchEndType type =
+ ( bndSegs.back()._inSeg->isConnected( bndSegs.back()._edge ) ?
+ SMESH_MAT2d::BE_ON_VERTEX :
+ SMESH_MAT2d::BE_END );
+ endType.insert( make_pair( bndSegs.back()._edge->vertex0(), type ));
+ }
+ }
+ // join the 1st and the last branch edges if it is the same branch
+ if ( bndSegs.back().branchID() != bndSegs.front().branchID() &&
+ bndSegs.back().isSameBranch( bndSegs.front() ))
+ {
+ vector<const TVDEdge*> & br1 = branchEdges[ bndSegs.front().branchID() ];
+ vector<const TVDEdge*> & br2 = branchEdges[ bndSegs.back().branchID() ];
+ br1.insert( br1.begin(), br2.begin(), br2.end() );
+ br2.clear();
+ }
+
+ // associate branchIDs and the input branch vector (arg)
+ vector< const SMESH_MAT2d::Branch* > branchByID( branchEdges.size() );
+ int nbBranches = 0;
+ for ( size_t i = 0; i < branchEdges.size(); ++i )
+ {
+ nbBranches += ( !branchEdges[i].empty() );
+ }
+ branch.resize( nbBranches );
+ for ( size_t iBr = 0, brID = 0; brID < branchEdges.size(); ++brID )
+ {
+ if ( !branchEdges[ brID ].empty() )
+ branchByID[ brID ] = & branch[ iBr++ ];
+ }
+
+ // Fill in BndPoints of each EDGE of the boundary
+
+ size_t iSeg = 0;
+ int edgeInd = -1, dInd = 0;
+ while ( iSeg < bndSegs.size() )
+ {
+ const size_t geomID = bndSegs[ iSeg ].geomEdge();
+ SMESH_MAT2d::BndPoints & bndPoints = boundary.getPoints( geomID );
+
+ size_t nbSegs = 0;
+ for ( size_t i = iSeg; i < bndSegs.size() && geomID == bndSegs[ i ].geomEdge(); ++i )
+ ++nbSegs;
+ size_t iSegEnd = iSeg + nbSegs;
+
+ // make TVDEdge know an index of bndSegs within BndPoints
+ for ( size_t i = iSeg; i < iSegEnd; ++i )
+ if ( bndSegs[i]._edge )
+ SMESH_MAT2d::Branch::setBndSegment( i - iSeg, bndSegs[i]._edge );
+
+ // parameters on EDGE
+
+ bndPoints._params.reserve( nbSegs + 1 );
+ bndPoints._params.push_back( bndSegs[ iSeg ]._inSeg->_p0->_param );
+
+ for ( size_t i = iSeg; i < iSegEnd; ++i )
+ bndPoints._params.push_back( bndSegs[ i ]._uLast );
+
+ // MA edges
+
+ bndPoints._maEdges.reserve( nbSegs );
+
+ for ( size_t i = iSeg; i < iSegEnd; ++i )
+ {
+ const size_t brID = bndSegs[ i ].branchID();
+ const SMESH_MAT2d::Branch* br = branchByID[ brID ];
+
+ if ( bndSegs[ i ]._edge && !branchEdges[ brID ].empty() )
+ {
+ edgeInd += dInd;
+
+ if ( edgeInd < 0 ||
+ edgeInd >= (int) branchEdges[ brID ].size() ||
+ branchEdges[ brID ][ edgeInd ] != bndSegs[ i ]._edge )
+ {
+ if ( bndSegs[ i ]._branchID < 0 &&
+ branchEdges[ brID ].back() == bndSegs[ i ]._edge )
+ {
+ edgeInd = branchEdges[ brID ].size() - 1;
+ dInd = -1;
+ }
+ else if ( bndSegs[ i ]._branchID > 0 &&
+ branchEdges[ brID ].front() == bndSegs[ i ]._edge )
+ {
+ edgeInd = 0;
+ dInd = +1;
+ }
+ else
+ {
+ for ( edgeInd = 0; edgeInd < branchEdges[ brID ].size(); ++edgeInd )
+ if ( branchEdges[ brID ][ edgeInd ] == bndSegs[ i ]._edge )
+ break;
+ dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1;
+ }
+ }
+ }
+ else
+ {
+ // no MA edge, bndSeg corresponds to an end point of a branch
+ if ( bndPoints._maEdges.empty() )
+ {
+ // should not get here according to algo design
+ edgeInd = 0;
+ }
+ else
+ {
+ edgeInd = branchEdges[ brID ].size();
+ dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1;
+ }
+ }
+
+ bndPoints._maEdges.push_back( make_pair( br, ( 1 + edgeInd ) * dInd ));
+
+ } // loop on bndSegs of an EDGE
+
+ iSeg = iSegEnd;
+
+ } // loop on all bndSegs
+
+
+ // fill the branches with MA edges
+ for ( size_t iBr = 0, brID = 0; brID < branchEdges.size(); ++brID )
+ if ( !branchEdges[brID].empty() )
+ {
+ branch[ iBr ].init( branchEdges[brID], & boundary, endType );
+ iBr++;
+ }
+ // set branches to branch ends
+ for ( size_t i = 0; i < branch.size(); ++i )
+ branch[i].setBranchesToEnds( branch );
+
+ // fill branchPnt arg
+ map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* > v2end;
+ for ( size_t i = 0; i < branch.size(); ++i )
+ {
+ if ( branch[i].getEnd(0)->_branches.size() > 2 )
+ v2end.insert( make_pair( branch[i].getEnd(0)->_vertex, branch[i].getEnd(0) ));
+ if ( branch[i].getEnd(1)->_branches.size() > 2 )
+ v2end.insert( make_pair( branch[i].getEnd(1)->_vertex, branch[i].getEnd(1) ));
+ }
+ branchPnt.resize( v2end.size() );
+ map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* >::iterator v2e = v2end.begin();
+ for ( size_t i = 0; v2e != v2end.end(); ++v2e, ++i )
+ branchPnt[ i ] = v2e->second;
+
+ } // makeMA()
+
+} // namespace
+
+//================================================================================
+/*!
+ * \brief MedialAxis constructor
+ * \param [in] face - a face to create MA for
+ * \param [in] edges - edges of the face (possibly not all) on the order they
+ * encounter in the face boundary.
+ * \param [in] minSegLen - minimal length of a mesh segment used to discretize
+ * the edges. It is used to define precision of MA approximation
+ */
+//================================================================================
+
+SMESH_MAT2d::MedialAxis::MedialAxis(const TopoDS_Face& face,
+ const std::vector< TopoDS_Edge >& edges,
+ const double minSegLen,
+ const bool ignoreCorners):
+ _face( face ), _boundary( edges.size() )
+{
+ // input to construct_voronoi()
+ vector< InPoint > inPoints;
+ vector< InSegment> inSegments;
+ if ( !makeInputData( face, edges, minSegLen, inPoints, inSegments, _scale ))
+ return;
+
+ //inSegmentsToFile( inSegments );
+
+ // build voronoi diagram
+ construct_voronoi( inSegments.begin(), inSegments.end(), &_vd );
+
+ // make MA data
+ makeMA( _vd, inPoints, inSegments, _branch, _branchPnt, _boundary );
+}
+
+//================================================================================
+/*!
+ * \brief Return UVs of ends of MA edges of a branch
+ */
+//================================================================================
+
+void SMESH_MAT2d::MedialAxis::getPoints( const Branch& branch,
+ std::vector< gp_XY >& points) const
+{
+ branch.getPoints( points, _scale );
+}
+
+//================================================================================
+/*!
+ * \brief Returns a BranchPoint corresponding to a given point on a geom EDGE
+ * \param [in] iGeomEdge - index of geom EDGE within a vector passed at construction
+ * \param [in] u - parameter of the point on EDGE curve
+ * \param [out] p - the found BranchPoint
+ * \return bool - is OK
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
+ double u,
+ BranchPoint& p ) const
+{
+ if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() )
+ return false;
+
+ const BndPoints& points = _pointsPerEdge[ iEdge ];
+ const bool edgeReverse = ( points._params[0] > points._params.back() );
+
+ if ( u < ( edgeReverse ? points._params.back() : points._params[0] ))
+ u = edgeReverse ? points._params.back() : points._params[0];
+ else if ( u > ( edgeReverse ? points._params[0] : points._params.back()) )
+ u = edgeReverse ? points._params[0] : points._params.back();
+
+ double r = ( u - points._params[0] ) / ( points._params.back() - points._params[0] );
+ int i = int( r * double( points._maEdges.size()-1 ));
+ if ( edgeReverse )
+ {
+ while ( points._params[i ] < u ) --i;
+ while ( points._params[i+1] > u ) ++i;
+ }
+ else
+ {
+ while ( points._params[i ] > u ) --i;
+ while ( points._params[i+1] < u ) ++i;
+ }
+ double edgeParam = ( u - points._params[i] ) / ( points._params[i+1] - points._params[i] );
+
+ const std::pair< const Branch*, int >& maE = points._maEdges[ i ];
+ bool maReverse = ( maE.second < 0 );
+
+ p._branch = maE.first;
+ p._iEdge = maE.second - 1; // countered from 1 to store sign
+ p._edgeParam = maReverse ? ( 1. - edgeParam ) : edgeParam;
+
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Check if a given boundary segment is a null-length segment on a concave
+ * boundary corner.
+ * \param [in] iEdge - index of a geom EDGE
+ * \param [in] iSeg - index of a boundary segment
+ * \return bool - true if the segment is on concave corner
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Boundary::IsConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const
+{
+ if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() )
+ return false;
+
+ const BndPoints& points = _pointsPerEdge[ iEdge ];
+ if ( points._params.size() >= iSeg+1 )
+ return false;
+
+ return Abs( points._params[ iEdge ] - points._params[ iEdge+1 ]) < 1e-20;
+}
+
+//================================================================================
+/*!
+ * \brief Creates a 3d curve corresponding to a Branch
+ * \param [in] branch - the Branch
+ * \return Adaptor3d_Curve* - the new curve the caller is to delete
+ */
+//================================================================================
+
+Adaptor3d_Curve* SMESH_MAT2d::MedialAxis::make3DCurve(const Branch& branch) const
+{
+ Handle(Geom_Surface) surface = BRep_Tool::Surface( _face );
+ if ( surface.IsNull() )
+ return 0;
+
+ vector< gp_XY > uv;
+ branch.getPoints( uv, _scale );
+ if ( uv.size() < 2 )
+ return 0;
+
+ vector< TopoDS_Vertex > vertex( uv.size() );
+ for ( size_t i = 0; i < uv.size(); ++i )
+ vertex[i] = BRepBuilderAPI_MakeVertex( surface->Value( uv[i].X(), uv[i].Y() ));
+
+ TopoDS_Wire aWire;
+ BRep_Builder aBuilder;
+ aBuilder.MakeWire(aWire);
+ for ( size_t i = 1; i < vertex.size(); ++i )
+ {
+ TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( vertex[i-1], vertex[i] );
+ aBuilder.Add( aWire, edge );
+ }
+
+ // if ( myEdge.size() == 2 && FirstVertex().IsSame( LastVertex() ))
+ // aWire.Closed(true); // issue 0021141
+
+ return new BRepAdaptor_CompCurve( aWire );
+}
+
+//================================================================================
+/*!
+ * \brief Copy points of an EDGE
+ */
+//================================================================================
+
+void SMESH_MAT2d::Branch::init( vector<const TVDEdge*>& maEdges,
+ const Boundary* boundary,
+ map< const TVDVertex*, BranchEndType > endType )
+{
+ if ( maEdges.empty() ) return;
+
+ _boundary = boundary;
+ _maEdges.swap( maEdges );
+
+
+ _params.reserve( _maEdges.size() + 1 );
+ _params.push_back( 0. );
+ for ( size_t i = 0; i < _maEdges.size(); ++i )
+ _params.push_back( _params.back() + length( _maEdges[i] ));
+
+ for ( size_t i = 1; i < _params.size(); ++i )
+ _params[i] /= _params.back();
+
+
+ _endPoint1._vertex = _maEdges.front()->vertex1();
+ _endPoint2._vertex = _maEdges.back ()->vertex0();
+
+ if ( endType.count( _endPoint1._vertex ))
+ _endPoint1._type = endType[ _endPoint1._vertex ];
+ if ( endType.count( _endPoint2._vertex ))
+ _endPoint2._type = endType[ _endPoint2._vertex ];
+}
+
+//================================================================================
+/*!
+ * \brief fill BranchEnd::_branches of its ends
+ */
+//================================================================================
+
+void SMESH_MAT2d::Branch::setBranchesToEnds( const vector< Branch >& branches )
+{
+ for ( size_t i = 0; i < branches.size(); ++i )
+ {
+ if ( this->_endPoint1._vertex == branches[i]._endPoint1._vertex ||
+ this->_endPoint1._vertex == branches[i]._endPoint2._vertex )
+ this->_endPoint1._branches.push_back( &branches[i] );
+
+ if ( this->_endPoint2._vertex == branches[i]._endPoint1._vertex ||
+ this->_endPoint2._vertex == branches[i]._endPoint2._vertex )
+ this->_endPoint2._branches.push_back( &branches[i] );
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
+ * \param [in] param - [0;1] normalized param on the Branch
+ * \param [out] bp1 - BoundaryPoint on EDGE with a lower index
+ * \param [out] bp2 - BoundaryPoint on EDGE with a higher index
+ * \return bool - true if the BoundaryPoint's found
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Branch::getBoundaryPoints(double param,
+ BoundaryPoint& bp1,
+ BoundaryPoint& bp2 ) const
+{
+ if ( param < _params[0] || param > _params.back() )
+ return false;
+
+ // look for an index of a MA edge by param
+ double ip = param * _params.size();
+ size_t i = size_t( Min( int( _maEdges.size()-1), int( ip )));
+
+ while ( param < _params[i ] ) --i;
+ while ( param > _params[i+1] ) ++i;
+
+ double r = ( param - _params[i] ) / ( _params[i+1] - _params[i] );
+
+ return getBoundaryPoints( i, r, bp1, bp2 );
+}
+
+//================================================================================
+/*!
+ * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
+ * \param [in] iMAEdge - index of a MA edge within this Branch
+ * \param [in] maEdgeParam - [0;1] normalized param on the \a iMAEdge
+ * \param [out] bp1 - BoundaryPoint on EDGE with a lower index
+ * \param [out] bp2 - BoundaryPoint on EDGE with a higher index
+ * \return bool - true if the BoundaryPoint's found
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Branch::getBoundaryPoints(std::size_t iMAEdge,
+ double maEdgeParam,
+ BoundaryPoint& bp1,
+ BoundaryPoint& bp2 ) const
+{
+ if ( iMAEdge > _maEdges.size() )
+ return false;
+
+ size_t iGeom1 = getGeomEdge( _maEdges[ iMAEdge ] );
+ size_t iGeom2 = getGeomEdge( _maEdges[ iMAEdge ]->twin() );
+ size_t iSeg1 = getBndSegment( _maEdges[ iMAEdge ] );
+ size_t iSeg2 = getBndSegment( _maEdges[ iMAEdge ]->twin() );
+
+ return ( _boundary->getPoint( iGeom1, iSeg1, maEdgeParam, bp1 ) &&
+ _boundary->getPoint( iGeom2, iSeg2, maEdgeParam, bp2 ));
+}
+
+//================================================================================
+/*!
+ * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Branch::getBoundaryPoints(const BranchPoint& p,
+ BoundaryPoint& bp1,
+ BoundaryPoint& bp2 ) const
+{
+ return ( p._branch ? p._branch : this )->getBoundaryPoints( p._iEdge, p._edgeParam, bp1, bp2 );
+}
+
+//================================================================================
+/*!
+ * \brief Return a parameter of a BranchPoint normalized within this Branch
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Branch::getParameter(const BranchPoint & p, double & u ) const
+{
+ if ( p._iEdge > _params.size()-1 )
+ return false;
+
+ u = ( _params[ p._iEdge ] * ( 1 - p._edgeParam ) +
+ _params[ p._iEdge+1 ] * p._edgeParam );
+
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Check type of both ends
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Branch::hasEndOfType(BranchEndType type) const
+{
+ return ( _endPoint1._type == type || _endPoint2._type == type );
+}
+
+//================================================================================
+/*!
+ * \brief Returns MA points
+ * \param [out] points - the 2d points
+ * \param [in] scale - the scale that was used to scale the 2d space of MA
+ */
+//================================================================================
+
+void SMESH_MAT2d::Branch::getPoints( std::vector< gp_XY >& points,
+ const double scale[2]) const
+{
+ points.resize( _maEdges.size() + 1 );
+
+ points[0].SetCoord( _maEdges[0]->vertex1()->x() / scale[0], // CCW order! -> vertex1 not vertex0
+ _maEdges[0]->vertex1()->y() / scale[1] );
+
+ for ( size_t i = 0; i < _maEdges.size(); ++i )
+ points[i+1].SetCoord( _maEdges[i]->vertex0()->x() / scale[0],
+ _maEdges[i]->vertex0()->y() / scale[1] );
+}
+
+//================================================================================
+/*!
+ * \brief Return indices of EDGEs equidistant from this branch
+ */
+//================================================================================
+
+void SMESH_MAT2d::Branch::getGeomEdges( std::vector< std::size_t >& edgeIDs1,
+ std::vector< std::size_t >& edgeIDs2 ) const
+{
+ edgeIDs1.push_back( getGeomEdge( _maEdges[0] ));
+ edgeIDs2.push_back( getGeomEdge( _maEdges[0]->twin() ));
+
+ for ( size_t i = 1; i < _maEdges.size(); ++i )
+ {
+ size_t ie1 = getGeomEdge( _maEdges[i] );
+ size_t ie2 = getGeomEdge( _maEdges[i]->twin() );
+
+ if ( edgeIDs1.back() != ie1 ) edgeIDs1.push_back( ie1 );
+ if ( edgeIDs2.back() != ie2 ) edgeIDs2.push_back( ie2 );
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Looks for a BranchPoint position around a concave VERTEX
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >& edgeIDs1,
+ std::vector< std::size_t >& edgeIDs2,
+ std::vector< BranchPoint >& divPoints,
+ const vector<const TVDEdge*>& maEdges,
+ const vector<const TVDEdge*>& maEdgesTwin,
+ size_t & i) const
+{
+ // if there is a concave vertex between EDGEs
+ // then position of a dividing BranchPoint is undefined, it is somewhere
+ // on an arc-shaped part of the Branch around the concave vertex.
+ // Chose this position by a VERTEX of the opposite EDGE, or put it in the middle
+ // of the arc if there is no opposite VERTEX.
+ // All null-length segments around a VERTEX belong to one of EDGEs.
+
+ BranchPoint divisionPnt;
+ divisionPnt._branch = this;
+
+ size_t ie1 = getGeomEdge( maEdges [i] );
+ size_t ie2 = getGeomEdge( maEdgesTwin[i] );
+
+ size_t iSeg1 = getBndSegment( maEdges[ i-1 ] );
+ size_t iSeg2 = getBndSegment( maEdges[ i ] );
+ bool isConcaPrev = _boundary->IsConcaveSegment( edgeIDs1.back(), iSeg1 );
+ bool isConcaNext = _boundary->IsConcaveSegment( ie1, iSeg2 );
+ if ( !isConcaNext && !isConcaPrev )
+ return false;
+
+ bool isConcaveV = false;
+
+ int iPrev = i-1, iNext = i;
+ if ( isConcaNext ) // all null-length segments follow
+ {
+ // look for a VERTEX of the opposite EDGE
+ ++iNext; // end of null-length segments
+ while ( iNext < maEdges.size() )
+ {
+ iSeg2 = getBndSegment( maEdges[ iNext ] );
+ if ( _boundary->IsConcaveSegment( ie1, iSeg2 ))
+ ++iNext;
+ else
+ break;
+ }
+ bool vertexFound = false;
+ for ( size_t iE = i+1; iE < iNext; ++iE )
+ {
+ ie2 = getGeomEdge( maEdgesTwin[iE] );
+ if ( ie2 != edgeIDs2.back() )
+ {
+ // opposite VERTEX found
+ divisionPnt._iEdge = iE;
+ divisionPnt._edgeParam = 0;
+ divPoints.push_back( divisionPnt );
+ edgeIDs1.push_back( ie1 );
+ edgeIDs2.push_back( ie2 );
+ vertexFound = true;
+ }
+ }
+ if ( vertexFound )
+ {
+ i = --iNext;
+ isConcaveV = true;
+ }
+ }
+ else if ( isConcaPrev )
+ {
+ // all null-length segments passed, find their beginning
+ while ( iPrev-1 >= 0 )
+ {
+ iSeg1 = getBndSegment( maEdges[ iPrev-1 ] );
+ if ( _boundary->IsConcaveSegment( edgeIDs1.back(), iSeg1 ))
+ --iPrev;
+ else
+ break;
+ }
+ }
+
+ if ( iPrev < i-1 || iNext > i )
+ {
+ // no VERTEX on the opposite EDGE, put the Branch Point in the middle
+ double par1 = _params[ iPrev ], par2 = _params[ iNext ];
+ double midPar = 0.5 * ( par1 + par2 );
+ divisionPnt._iEdge = iPrev;
+ while ( _params[ divisionPnt._iEdge + 1 ] < midPar )
+ ++divisionPnt._iEdge;
+ divisionPnt._edgeParam =
+ ( _params[ divisionPnt._iEdge + 1 ] - midPar ) /
+ ( _params[ divisionPnt._iEdge + 1 ] - _params[ divisionPnt._iEdge ] );
+ divPoints.push_back( divisionPnt );
+ isConcaveV = true;
+ }
+
+ return isConcaveV;
+}
+
+//================================================================================
+/*!
+ * \brief Return indices of opposite parts of EDGEs equidistant from this branch
+ * \param [out] edgeIDs1 - EDGE index opposite to the edgeIDs2[i]-th EDGE
+ * \param [out] edgeIDs2 - EDGE index opposite to the edgeIDs1[i]-th EDGE
+ * \param [out] divPoints - BranchPoint's located between two successive unique
+ * pairs of EDGE indices. A \a divPoints[i] can separate e.g. two following pairs
+ * of EDGE indices < 0, 2 > and < 0, 1 >. Number of \a divPoints is one less
+ * than number of \a edgeIDs
+ */
+//================================================================================
+
+void SMESH_MAT2d::Branch::getOppositeGeomEdges( std::vector< std::size_t >& edgeIDs1,
+ std::vector< std::size_t >& edgeIDs2,
+ std::vector< BranchPoint >& divPoints) const
+{
+ edgeIDs1.clear();
+ edgeIDs2.clear();
+ divPoints.clear();
+
+ edgeIDs1.push_back( getGeomEdge( _maEdges[0] ));
+ edgeIDs2.push_back( getGeomEdge( _maEdges[0]->twin() ));
+
+ std::vector<const TVDEdge*> twins( _maEdges.size() );
+ for ( size_t i = 0; i < _maEdges.size(); ++i )
+ twins[i] = _maEdges[i]->twin();
+
+ // size_t lastConcaE1 = _boundary.nbEdges();
+ // size_t lastConcaE2 = _boundary.nbEdges();
+
+ BranchPoint divisionPnt;
+ divisionPnt._branch = this;
+
+ for ( size_t i = 0; i < _maEdges.size(); ++i )
+ {
+ size_t ie1 = getGeomEdge( _maEdges[i] );
+ size_t ie2 = getGeomEdge( _maEdges[i]->twin() );
+
+ if ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 )
+ {
+ bool isConcaveV = false;
+ if ( edgeIDs1.back() != ie1 && edgeIDs2.back() == ie2 )
+ {
+ isConcaveV = addDivPntForConcaVertex( edgeIDs1, edgeIDs2, divPoints, _maEdges, twins, i );
+ }
+ if ( edgeIDs1.back() == ie1 && edgeIDs2.back() != ie2 )
+ {
+ isConcaveV = addDivPntForConcaVertex( edgeIDs2, edgeIDs1, divPoints, twins, _maEdges, i );
+ }
+
+ if ( isConcaveV )
+ {
+ ie1 = getGeomEdge( _maEdges[i] );
+ ie2 = getGeomEdge( _maEdges[i]->twin() );
+ }
+ if (( !isConcaveV ) ||
+ ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 ))
+ {
+ edgeIDs1.push_back( ie1 );
+ edgeIDs2.push_back( ie2 );
+ }
+ if ( divPoints.size() < edgeIDs1.size() - 1 )
+ {
+ divisionPnt._iEdge = i;
+ divisionPnt._edgeParam = 0;
+ divPoints.push_back( divisionPnt );
+ }
+
+ } // if ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 )
+ } // loop on _maEdges
+}
+
+//================================================================================
+/*!
+ * \brief Store data of boundary segments in TVDEdge
+ */
+//================================================================================
+
+void SMESH_MAT2d::Branch::setGeomEdge( std::size_t geomIndex, const TVDEdge* maEdge )
+{
+ if ( maEdge ) maEdge->cell()->color( geomIndex );
+}
+std::size_t SMESH_MAT2d::Branch::getGeomEdge( const TVDEdge* maEdge )
+{
+ return maEdge ? maEdge->cell()->color() : std::string::npos;
+}
+void SMESH_MAT2d::Branch::setBndSegment( std::size_t segIndex, const TVDEdge* maEdge )
+{
+ if ( maEdge ) maEdge->color( segIndex );
+}
+std::size_t SMESH_MAT2d::Branch::getBndSegment( const TVDEdge* maEdge )
+{
+ return maEdge ? maEdge->color() : std::string::npos;
+}
+
+//================================================================================
+/*!
+ * \brief Returns a boundary point on a given EDGE
+ * \param [in] iEdge - index of the EDGE within MedialAxis
+ * \param [in] iSeg - index of a boundary segment within this Branch
+ * \param [in] u - [0;1] normalized param within \a iSeg-th segment
+ * \param [out] bp - the found BoundaryPoint
+ * \return bool - true if the BoundaryPoint is found
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Boundary::getPoint( std::size_t iEdge,
+ std::size_t iSeg,
+ double u,
+ BoundaryPoint& bp ) const
+{
+ if ( iEdge >= _pointsPerEdge.size() )
+ return false;
+ if ( iSeg+1 >= _pointsPerEdge[ iEdge ]._params.size() )
+ return false;
+
+ // This method is called by Branch that can have an opposite orientation,
+ // hence u is inverted depending on orientation coded as a sign of _maEdge index
+ bool isReverse = ( _pointsPerEdge[ iEdge ]._maEdges[ iSeg ].second < 0 );
+ if ( isReverse )
+ u = 1. - u;
+
+ double p0 = _pointsPerEdge[ iEdge ]._params[ iSeg ];
+ double p1 = _pointsPerEdge[ iEdge ]._params[ iSeg+1 ];
+
+ bp._param = p0 * ( 1. - u ) + p1 * u;
+ bp._edgeIndex = iEdge;
+
+ return true;
+}
+
--- /dev/null
+// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File : SMESH_MAT2d.hxx
+// Created : Thu May 28 17:49:53 2015
+// Author : Edward AGAPOV (eap)
+
+#ifndef __SMESH_MAT2d_HXX__
+#define __SMESH_MAT2d_HXX__
+
+#include "SMESH_Utils.hxx"
+
+#include <TopoDS_Face.hxx>
+#include <TopoDS_Edge.hxx>
+
+#include <vector>
+#include <map>
+
+#include <boost/polygon/polygon.hpp>
+#include <boost/polygon/voronoi.hpp>
+
+class Adaptor3d_Curve;
+
+// Medial Axis Transform 2D
+namespace SMESH_MAT2d
+{
+ class MedialAxis; // MedialAxis is the entry point
+ class Branch;
+ class BranchEnd;
+ class Boundary;
+ struct BoundaryPoint;
+
+ typedef boost::polygon::voronoi_diagram<double> TVD;
+ typedef TVD::cell_type TVDCell;
+ typedef TVD::edge_type TVDEdge;
+ typedef TVD::vertex_type TVDVertex;
+
+ //-------------------------------------------------------------------------------------
+ // type of Branch end point
+ enum BranchEndType { BE_UNDEF,
+ BE_ON_VERTEX, // branch ends at a convex VRTEX
+ BE_BRANCH_POINT, // branch meats 2 or more other branches
+ BE_END // branch end equidistant from several adjacent segments
+ };
+ //-------------------------------------------------------------------------------------
+ /*!
+ * \brief End point of MA Branch
+ */
+ struct SMESHUtils_EXPORT BranchEnd
+ {
+ const TVDVertex* _vertex;
+ BranchEndType _type;
+ std::vector< const Branch* > _branches;
+
+ BranchEnd(): _vertex(0), _type( BE_UNDEF ) {}
+ };
+ //-------------------------------------------------------------------------------------
+ /*!
+ * \brief Point on MA Branch
+ */
+ struct SMESHUtils_EXPORT BranchPoint
+ {
+ const Branch* _branch;
+ std::size_t _iEdge; // MA edge index within the branch
+ double _edgeParam; // normalized param within the MA edge
+
+ BranchPoint(): _branch(0), _iEdge(0), _edgeParam(-1) {}
+ };
+ //-------------------------------------------------------------------------------------
+ /*!
+ * \brief Branch is a set of MA edges enclosed between branch points and/or MA ends.
+ * It's main feature is to return two BoundaryPoint's per a point on it.
+ */
+ class SMESHUtils_EXPORT Branch
+ {
+ public:
+ bool getBoundaryPoints(double param, BoundaryPoint& bp1, BoundaryPoint& bp2 ) const;
+ bool getBoundaryPoints(std::size_t iMAEdge, double maEdgeParam,
+ BoundaryPoint& bp1, BoundaryPoint& bp2 ) const;
+ bool getBoundaryPoints(const BranchPoint& p,
+ BoundaryPoint& bp1, BoundaryPoint& bp2 ) const;
+ bool getParameter(const BranchPoint& p, double & u ) const;
+
+ std::size_t nbEdges() const { return _maEdges.size(); }
+
+ const BranchEnd* getEnd(bool the2nd) const { return & ( the2nd ? _endPoint2 : _endPoint1 ); }
+
+ bool hasEndOfType(BranchEndType type) const;
+
+ void getPoints( std::vector< gp_XY >& points, const double scale[2]) const;
+
+ void getGeomEdges( std::vector< std::size_t >& edgeIDs1,
+ std::vector< std::size_t >& edgeIDs2 ) const;
+
+ void getOppositeGeomEdges( std::vector< std::size_t >& edgeIDs1,
+ std::vector< std::size_t >& edgeIDs2,
+ std::vector< BranchPoint >& divPoints) const;
+
+ // construction
+ void init( std::vector<const TVDEdge*>& maEdges,
+ const Boundary* boundary,
+ std::map< const TVDVertex*, BranchEndType > endType);
+ void setBranchesToEnds( const std::vector< Branch >& branches);
+
+ static void setGeomEdge( std::size_t geomIndex, const TVDEdge* maEdge );
+ static std::size_t getGeomEdge( const TVDEdge* maEdge );
+ static void setBndSegment( std::size_t segIndex, const TVDEdge* maEdge );
+ static std::size_t getBndSegment( const TVDEdge* maEdge );
+
+ private:
+
+ bool addDivPntForConcaVertex( std::vector< std::size_t >& edgeIDs1,
+ std::vector< std::size_t >& edgeIDs2,
+ std::vector< BranchPoint >& divPoints,
+ const std::vector<const TVDEdge*>& maEdges,
+ const std::vector<const TVDEdge*>& maEdgesTwin,
+ size_t & i) const;
+
+ // association of _maEdges with boundary segments is stored in this way:
+ // index of an EDGE: TVDEdge->cell()->color()
+ // index of a segment on EDGE: TVDEdge->color()
+ std::vector<const TVDEdge*> _maEdges; // MA edges ending at points located at _params
+ std::vector<double> _params; // params of points on MA, normalized [0;1] within this branch
+ const Boundary* _boundary; // face boundary
+ BranchEnd _endPoint1;
+ BranchEnd _endPoint2;
+ };
+
+ //-------------------------------------------------------------------------------------
+ /*!
+ * \brief Data of a discretized EDGE allowing to get a point on MA by a parameter on EDGE
+ */
+ struct BndPoints
+ {
+ std::vector< double > _params; // params of discretization points on an EDGE
+ std::vector< std::pair< const Branch*, int > > _maEdges; /* index of TVDEdge in branch;
+ index sign means orientation;
+ index == Branch->nbEdges() means
+ end point of a Branch */
+ };
+ //-------------------------------------------------------------------------------------
+ /*!
+ * \brief Face boundary is discretized so that each its segment to correspond to
+ * an edge of MA
+ */
+ class Boundary
+ {
+ public:
+
+ Boundary( std::size_t nbEdges ): _pointsPerEdge( nbEdges ) {}
+ BndPoints& getPoints( std::size_t iEdge ) { return _pointsPerEdge[ iEdge ]; }
+ std::size_t nbEdges() const { return _pointsPerEdge.size(); }
+
+ bool getPoint( std::size_t iEdge, std::size_t iSeg, double u, BoundaryPoint& bp ) const;
+
+ bool getBranchPoint( const std::size_t iEdge, double u, BranchPoint& p ) const;
+
+ bool IsConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const;
+
+ private:
+ std::vector< BndPoints > _pointsPerEdge;
+ };
+
+ //-------------------------------------------------------------------------------------
+ /*!
+ * \brief Point on FACE boundary
+ */
+ struct SMESHUtils_EXPORT BoundaryPoint
+ {
+ std::size_t _edgeIndex; // index of an EDGE in a sequence passed to MedialAxis()
+ double _param; // parameter of this EDGE
+ };
+ //-------------------------------------------------------------------------------------
+ /*!
+ * \brief Medial axis (MA) is defined as the loci of centres of locally
+ * maximal balls inside 2D representation of a face. This class
+ * implements a piecewise approximation of MA.
+ */
+ class SMESHUtils_EXPORT MedialAxis
+ {
+ public:
+ MedialAxis(const TopoDS_Face& face,
+ const std::vector< TopoDS_Edge >& edges,
+ const double minSegLen,
+ const bool ignoreCorners = false );
+ const Boundary& getBoundary() const { return _boundary; }
+ const std::vector< Branch >& getBranches() const { return _branch; }
+ const std::vector< const BranchEnd* >& getBranchPoints() const { return _branchPnt; }
+
+ void getPoints( const Branch& branch, std::vector< gp_XY >& points) const;
+ Adaptor3d_Curve* make3DCurve(const Branch& branch) const;
+
+ private:
+
+ private:
+ TopoDS_Face _face;
+ TVD _vd;
+ std::vector< Branch > _branch;
+ std::vector< const BranchEnd* > _branchPnt;
+ Boundary _boundary;
+ double _scale[2];
+ };
+
+}
+
+#endif
QUADRANGLE = "Quadrangle_2D"
## Algorithm type: Radial Quadrangle 1D-2D algorithm, see StdMeshersBuilder_RadialQuadrangle1D2D
RADIAL_QUAD = "RadialQuadrangle_1D2D"
+## Algorithm type: Quadrangle (Medial Axis Projection) 1D-2D algorithm, see StdMeshersBuilder_QuadMA_1D2D
+QUAD_MA_PROJ = "QuadFromMedialAxis_1D2D"
# import items of enums
for e in StdMeshers.QuadType._items: exec('%s = StdMeshers.%s'%(e,e))
algoType = PYTHON
## doc string of the method
# @internal
- docHelper = "Creates tetrahedron 3D algorithm for solids"
- ## doc string of the method
- # @internal
docHelper = "Creates segment 1D algorithm for edges"
## Private constructor.
algoType = "Projection_1D2D"
## doc string of the method
# @internal
- docHelper = "Creates projection 1D-2D algorithm for edges and faces"
+ docHelper = "Creates projection 1D-2D algorithm for faces"
## Private constructor.
# @param mesh parent mesh object algorithm is assigned to
algoType = "RadialPrism_3D"
## doc string of the method
# @internal
- docHelper = "Creates prism 3D algorithm for volumes"
+ docHelper = "Creates Raial Prism 3D algorithm for volumes"
## Private constructor.
# @param mesh parent mesh object algorithm is assigned to
algoType = RADIAL_QUAD
## doc string of the method
# @internal
- docHelper = "Creates quadrangle 1D-2D algorithm for triangular faces"
+ docHelper = "Creates quadrangle 1D-2D algorithm for faces having a shape of disk or a disk segment"
## Private constructor.
# @param mesh parent mesh object algorithm is assigned to
pass # end of StdMeshersBuilder_RadialQuadrangle1D2D class
+## Defines a Quadrangle (Medial Axis Projection) 1D-2D algorithm
+#
+# It is created by calling smeshBuilder.Mesh.Quadrangle(smeshBuilder.QUAD_MA_PROJ,geom=0)
+#
+# @ingroup l2_algos_quad_ma
+class StdMeshersBuilder_QuadMA_1D2D(Mesh_Algorithm):
+
+ ## name of the dynamic method in smeshBuilder.Mesh class
+ # @internal
+ meshMethod = "Quadrangle"
+ ## type of algorithm used with helper function in smeshBuilder.Mesh class
+ # @internal
+ algoType = QUAD_MA_PROJ
+ ## doc string of the method
+ # @internal
+ docHelper = "Creates quadrangle 1D-2D algorithm for faces"
+
+ ## Private constructor.
+ # @param mesh parent mesh object algorithm is assigned to
+ # @param geom geometry (shape/sub-shape) algorithm is assigned to;
+ # if it is @c 0 (default), the algorithm is assigned to the main shape
+ def __init__(self, mesh, geom=0):
+ Mesh_Algorithm.__init__(self)
+ self.Create(mesh, geom, self.algoType)
+ pass
+
+ pass
+
+
## Defines a Use Existing Elements 1D algorithm
#
# It is created by calling smeshBuilder.Mesh.UseExisting1DElements(geom=0)
isDefault = True
## doc string of the method
# @internal
- docHelper = "Creates 1D-2D algorithm for edges/faces with reusing of existing mesh elements"
+ docHelper = "Creates 1D-2D algorithm for faces with reusing of existing mesh elements"
## Private constructor.
# @param mesh parent mesh object algorithm is assigned to
isDefault = True
## doc string of the method
# @internal
- docHelper = "Creates body fitting 3D algorithm for volumes"
+ docHelper = "Creates Body Fitting 3D algorithm for volumes"
## Private constructor.
# @param mesh parent mesh object algorithm is assigned to
algoType = "UseExisting_1D"
## doc string of the method
# @internal
- docHelper = "Creates 1D algorithm for edges with reusing of existing mesh elements"
+ docHelper = "Creates 1D algorithm allowing batch meshing of edges"
## Private constructor.
# @param mesh parent mesh object algorithm is assigned to
algoType = "UseExisting_2D"
## doc string of the method
# @internal
- docHelper = "Creates 2D algorithm for faces with reusing of existing mesh elements"
+ docHelper = "Creates 2D algorithm allowing batch meshing of faces"
## Private constructor.
# @param mesh parent mesh object algorithm is assigned to
StdMeshers_Projection_1D2D.hxx
StdMeshers_CartesianParameters3D.hxx
StdMeshers_Cartesian_3D.hxx
+ StdMeshers_QuadFromMedialAxis_1D2D.hxx
)
IF(SALOME_SMESH_ENABLE_MEFISTO)
StdMeshers_CartesianParameters3D.cxx
StdMeshers_Cartesian_3D.cxx
StdMeshers_Adaptive1D.cxx
+ StdMeshers_QuadFromMedialAxis_1D2D.cxx
)
IF(SALOME_SMESH_ENABLE_MEFISTO)
double xmax = -1.e300;
double ymin = 1.e300;
double ymax = -1.e300;
- int nbp = 23;
+ const int nbp = 23;
scalex = 1;
scaley = 1;
ymin = p.Y();
if (p.Y() > ymax)
ymax = p.Y();
- // MESSAGE(" "<< f<<" "<<l<<" "<<param<<" "<<xmin<<" "<<xmax<<" "<<ymin<<" "<<ymax);
}
}
- // SCRUTE(xmin);
- // SCRUTE(xmax);
- // SCRUTE(ymin);
- // SCRUTE(ymax);
double xmoy = (xmax + xmin) / 2.;
double ymoy = (ymax + ymin) / 2.;
double xsize = xmax - xmin;
}
scalex = length_x / xsize;
scaley = length_y / ysize;
-// SCRUTE(xsize);
-// SCRUTE(ysize);
double xyratio = xsize*scalex/(ysize*scaley);
const double maxratio = 1.e2;
- //SCRUTE(xyratio);
if (xyratio > maxratio) {
- SCRUTE( scaley );
scaley *= xyratio / maxratio;
- SCRUTE( scaley );
}
else if (xyratio < 1./maxratio) {
- SCRUTE( scalex );
scalex *= 1 / xyratio / maxratio;
- SCRUTE( scalex );
}
- ASSERT(scalex);
- ASSERT(scaley);
}
// namespace
--- /dev/null
+// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File : StdMeshers_QuadFromMedialAxis_1D2D.cxx
+// Created : Wed Jun 3 17:33:45 2015
+// Author : Edward AGAPOV (eap)
+
+#include "StdMeshers_QuadFromMedialAxis_1D2D.hxx"
+
+#include "SMESH_Block.hxx"
+#include "SMESH_Gen.hxx"
+#include "SMESH_MAT2d.hxx"
+#include "SMESH_Mesh.hxx"
+#include "SMESH_MesherHelper.hxx"
+#include "SMESH_ProxyMesh.hxx"
+#include "SMESH_subMesh.hxx"
+#include "StdMeshers_FaceSide.hxx"
+#include "StdMeshers_Regular_1D.hxx"
+#include "StdMeshers_ViscousLayers2D.hxx"
+
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepTools.hxx>
+#include <BRep_Tool.hxx>
+#include <GeomAPI_Interpolate.hxx>
+#include <Geom_Surface.hxx>
+#include <Precision.hxx>
+#include <TColgp_HArray1OfPnt.hxx>
+#include <TopExp.hxx>
+#include <TopLoc_Location.hxx>
+#include <TopTools_MapOfShape.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopoDS_Vertex.hxx>
+#include <gp_Pnt.hxx>
+
+#include <list>
+#include <vector>
+
+//================================================================================
+/*!
+ * \brief 1D algo
+ */
+class StdMeshers_QuadFromMedialAxis_1D2D::Algo1D : public StdMeshers_Regular_1D
+{
+public:
+ Algo1D(int studyId, SMESH_Gen* gen):
+ StdMeshers_Regular_1D( gen->GetANewId(), studyId, gen )
+ {
+ }
+ void SetSegmentLength( double len )
+ {
+ _value[ BEG_LENGTH_IND ] = len;
+ _value[ PRECISION_IND ] = 1e-7;
+ _hypType = LOCAL_LENGTH;
+ }
+};
+
+//================================================================================
+/*!
+ * \brief Constructor sets algo features
+ */
+//================================================================================
+
+StdMeshers_QuadFromMedialAxis_1D2D::StdMeshers_QuadFromMedialAxis_1D2D(int hypId,
+ int studyId,
+ SMESH_Gen* gen)
+ : StdMeshers_Quadrangle_2D(hypId, studyId, gen),
+ _regular1D( 0 )
+{
+ _name = "QuadFromMedialAxis_1D2D";
+ _shapeType = (1 << TopAbs_FACE);
+ _onlyUnaryInput = true; // FACE by FACE so far
+ _requireDiscreteBoundary = false; // make 1D by myself
+ _supportSubmeshes = true; // make 1D by myself
+ _neededLowerHyps[ 1 ] = true; // suppress warning on hiding a global 1D algo
+ _neededLowerHyps[ 2 ] = true; // suppress warning on hiding a global 2D algo
+ _compatibleHypothesis.clear();
+ //_compatibleHypothesis.push_back("ViscousLayers2D");
+}
+
+//================================================================================
+/*!
+ * \brief Destructor
+ */
+//================================================================================
+
+StdMeshers_QuadFromMedialAxis_1D2D::~StdMeshers_QuadFromMedialAxis_1D2D()
+{
+ delete _regular1D;
+ _regular1D = 0;
+}
+
+//================================================================================
+/*!
+ * \brief Check if needed hypotheses are present
+ */
+//================================================================================
+
+bool StdMeshers_QuadFromMedialAxis_1D2D::CheckHypothesis(SMESH_Mesh& aMesh,
+ const TopoDS_Shape& aShape,
+ Hypothesis_Status& aStatus)
+{
+ return true; // does not require hypothesis
+}
+
+namespace
+{
+ //================================================================================
+ /*!
+ * \brief Temporary mesh
+ */
+ struct TmpMesh : public SMESH_Mesh
+ {
+ TmpMesh()
+ {
+ _myMeshDS = new SMESHDS_Mesh(/*id=*/0, /*isEmbeddedMode=*/true);
+ }
+ };
+
+ //================================================================================
+ /*!
+ * \brief Select two EDGEs from a map, either mapped to least values or to max values
+ */
+ //================================================================================
+
+ // template< class TVal2EdgesMap >
+ // void getTwo( bool least,
+ // TVal2EdgesMap& map,
+ // vector<TopoDS_Edge>& twoEdges,
+ // vector<TopoDS_Edge>& otherEdges)
+ // {
+ // twoEdges.clear();
+ // otherEdges.clear();
+ // if ( least )
+ // {
+ // TVal2EdgesMap::iterator i = map.begin();
+ // twoEdges.push_back( i->second );
+ // twoEdges.push_back( ++i->second );
+ // for ( ; i != map.end(); ++i )
+ // otherEdges.push_back( i->second );
+ // }
+ // else
+ // {
+ // TVal2EdgesMap::reverse_iterator i = map.rbegin();
+ // twoEdges.push_back( i->second );
+ // twoEdges.push_back( ++i->second );
+ // for ( ; i != map.rend(); ++i )
+ // otherEdges.push_back( i->second );
+ // }
+ // TopoDS_Vertex v;
+ // if ( TopExp::CommonVertex( twoEdges[0], twoEdges[1], v ))
+ // {
+ // twoEdges.clear(); // two EDGEs must not be connected
+ // otherEdges.clear();
+ // }
+ // }
+
+ //================================================================================
+ /*!
+ * \brief Finds out a minimal segment length given EDGEs will be divided into.
+ * This length is further used to discretize the Medial Axis
+ */
+ //================================================================================
+
+ double getMinSegLen(SMESH_MesherHelper& theHelper,
+ const vector<TopoDS_Edge>& theEdges)
+ {
+ TmpMesh tmpMesh;
+ SMESH_Mesh* mesh = theHelper.GetMesh();
+
+ vector< SMESH_Algo* > algos( theEdges.size() );
+ for ( size_t i = 0; i < theEdges.size(); ++i )
+ {
+ SMESH_subMesh* sm = mesh->GetSubMesh( theEdges[i] );
+ algos[i] = sm->GetAlgo();
+ }
+
+ const int nbSegDflt = mesh->GetGen()->GetDefaultNbSegments();
+ double minSegLen = Precision::Infinite();
+
+ for ( size_t i = 0; i < theEdges.size(); ++i )
+ {
+ SMESH_subMesh* sm = mesh->GetSubMesh( theEdges[i] );
+ if ( SMESH_Algo::IsStraight( theEdges[i], /*degenResult=*/true ))
+ continue;
+ // get algo
+ size_t iOpp = ( theEdges.size() == 4 ? (i+2)%4 : i );
+ SMESH_Algo* algo = sm->GetAlgo();
+ if ( !algo ) algo = algos[ iOpp ];
+ // get hypo
+ SMESH_Hypothesis::Hypothesis_Status status = SMESH_Hypothesis::HYP_MISSING;
+ if ( algo )
+ {
+ if ( !algo->CheckHypothesis( *mesh, theEdges[i], status ))
+ algo->CheckHypothesis( *mesh, theEdges[iOpp], status );
+ }
+ // compute
+ if ( status != SMESH_Hypothesis::HYP_OK )
+ {
+ minSegLen = Min( minSegLen, SMESH_Algo::EdgeLength( theEdges[i] ) / nbSegDflt );
+ }
+ else
+ {
+ tmpMesh.Clear();
+ tmpMesh.ShapeToMesh( TopoDS_Shape());
+ tmpMesh.ShapeToMesh( theEdges[i] );
+ try {
+ mesh->GetGen()->Compute( tmpMesh, theEdges[i], true, true ); // make nodes on VERTEXes
+ if ( !algo->Compute( tmpMesh, theEdges[i] ))
+ continue;
+ }
+ catch (...) {
+ continue;
+ }
+ SMDS_EdgeIteratorPtr segIt = tmpMesh.GetMeshDS()->edgesIterator();
+ while ( segIt->more() )
+ {
+ const SMDS_MeshElement* seg = segIt->next();
+ double len = SMESH_TNodeXYZ( seg->GetNode(0) ).Distance( seg->GetNode(1) );
+ minSegLen = Min( minSegLen, len );
+ }
+ }
+ }
+ if ( Precision::IsInfinite( minSegLen ))
+ minSegLen = mesh->GetShapeDiagonalSize() / nbSegDflt;
+
+ return minSegLen;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Returns EDGEs located between two VERTEXes at which given MA branches end
+ * \param [in] br1 - one MA branch
+ * \param [in] br2 - one more MA branch
+ * \param [in] allEdges - all EDGEs of a FACE
+ * \param [out] shortEdges - the found EDGEs
+ * \return bool - is OK or not
+ */
+ //================================================================================
+
+ bool getConnectedEdges( const SMESH_MAT2d::Branch* br1,
+ const SMESH_MAT2d::Branch* br2,
+ const vector<TopoDS_Edge>& allEdges,
+ vector<TopoDS_Edge>& shortEdges)
+ {
+ vector< size_t > edgeIDs[4];
+ br1->getGeomEdges( edgeIDs[0], edgeIDs[1] );
+ br2->getGeomEdges( edgeIDs[2], edgeIDs[3] );
+
+ // EDGEs returned by a Branch form a connected chain with a VERTEX where
+ // the Branch ends at the chain middle. One of end EDGEs of the chain is common
+ // with either end EDGE of the chain of the other Branch, or the chains are connected
+ // at a common VERTEX;
+
+ // Get indices of end EDGEs of the branches
+ bool vAtStart1 = ( br1->getEnd(0)->_type == SMESH_MAT2d::BE_ON_VERTEX );
+ bool vAtStart2 = ( br2->getEnd(0)->_type == SMESH_MAT2d::BE_ON_VERTEX );
+ size_t iEnd[4] = {
+ vAtStart1 ? edgeIDs[0].back() : edgeIDs[0][0],
+ vAtStart1 ? edgeIDs[1].back() : edgeIDs[1][0],
+ vAtStart2 ? edgeIDs[2].back() : edgeIDs[2][0],
+ vAtStart2 ? edgeIDs[3].back() : edgeIDs[3][0]
+ };
+
+ set< size_t > connectedIDs;
+ TopoDS_Vertex vCommon;
+ // look for the same EDGEs
+ for ( int i = 0; i < 2; ++i )
+ for ( int j = 2; j < 4; ++j )
+ if ( iEnd[i] == iEnd[j] )
+ {
+ connectedIDs.insert( edgeIDs[i].begin(), edgeIDs[i].end() );
+ connectedIDs.insert( edgeIDs[j].begin(), edgeIDs[j].end() );
+ i = j = 4;
+ }
+ if ( connectedIDs.empty() )
+ // look for connected EDGEs
+ for ( int i = 0; i < 2; ++i )
+ for ( int j = 2; j < 4; ++j )
+ if ( TopExp::CommonVertex( allEdges[ iEnd[i]], allEdges[ iEnd[j]], vCommon ))
+ {
+ connectedIDs.insert( edgeIDs[i].begin(), edgeIDs[i].end() );
+ connectedIDs.insert( edgeIDs[j].begin(), edgeIDs[j].end() );
+ i = j = 4;
+ }
+ if ( connectedIDs.empty() || // nothing
+ allEdges.size() - connectedIDs.size() < 2 ) // too many
+ return false;
+
+ // set shortEdges in the order as in allEdges
+ if ( connectedIDs.count( 0 ) &&
+ connectedIDs.count( allEdges.size()-1 ))
+ {
+ size_t iE = allEdges.size()-1;
+ while ( connectedIDs.count( iE-1 ))
+ --iE;
+ for ( size_t i = 0; i < connectedIDs.size(); ++i )
+ {
+ shortEdges.push_back( allEdges[ iE ]);
+ iE = ( iE + 1 ) % allEdges.size();
+ }
+ }
+ else
+ {
+ set< size_t >::iterator i = connectedIDs.begin();
+ for ( ; i != connectedIDs.end(); ++i )
+ shortEdges.push_back( allEdges[ *i ]);
+ }
+ return true;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Find EDGEs to discretize using projection from MA
+ * \param [in] theFace - the FACE to be meshed
+ * \param [in] theWire - ordered EDGEs of the FACE
+ * \param [out] theSinuEdges - the EDGEs to discretize using projection from MA
+ * \param [out] theShortEdges - other EDGEs
+ * \return bool - OK or not
+ *
+ * Is separate all EDGEs into four sides of a quadrangle connected in the order:
+ * theSinuEdges[0], theShortEdges[0], theSinuEdges[1], theShortEdges[1]
+ */
+ //================================================================================
+
+ bool getSinuousEdges( SMESH_MesherHelper& theHelper,
+ const TopoDS_Face& theFace,
+ list<TopoDS_Edge>& theWire,
+ vector<TopoDS_Edge> theSinuEdges[2],
+ vector<TopoDS_Edge> theShortEdges[2])
+ {
+ theSinuEdges[0].clear();
+ theSinuEdges[1].clear();
+ theShortEdges[0].clear();
+ theShortEdges[1].clear();
+
+ vector<TopoDS_Edge> allEdges( theWire.begin(), theWire.end() );
+ const size_t nbEdges = allEdges.size();
+ if ( nbEdges < 4 )
+ return false;
+
+ // create MedialAxis to find short edges by analyzing MA branches
+ double minSegLen = getMinSegLen( theHelper, allEdges );
+ SMESH_MAT2d::MedialAxis ma( theFace, allEdges, minSegLen );
+
+ // in an initial request case, theFace represents a part of a river with almost parallel banks
+ // so there should be two branch points
+ using SMESH_MAT2d::BranchEnd;
+ using SMESH_MAT2d::Branch;
+ const vector< const BranchEnd* >& braPoints = ma.getBranchPoints();
+ if ( braPoints.size() < 2 )
+ return false;
+ TopTools_MapOfShape shortMap;
+ size_t nbBranchPoints = 0;
+ for ( size_t i = 0; i < braPoints.size(); ++i )
+ {
+ vector< const Branch* > vertBranches; // branches with an end on VERTEX
+ for ( size_t ib = 0; ib < braPoints[i]->_branches.size(); ++ib )
+ {
+ const Branch* branch = braPoints[i]->_branches[ ib ];
+ if ( branch->hasEndOfType( SMESH_MAT2d::BE_ON_VERTEX ))
+ vertBranches.push_back( branch );
+ }
+ if ( vertBranches.size() != 2 || braPoints[i]->_branches.size() != 3)
+ continue;
+
+ // get common EDGEs of two branches
+ if ( !getConnectedEdges( vertBranches[0], vertBranches[1],
+ allEdges, theShortEdges[ nbBranchPoints > 0 ] ))
+ return false;
+
+ for ( size_t iS = 0; iS < theShortEdges[ nbBranchPoints ].size(); ++iS )
+ shortMap.Add( theShortEdges[ nbBranchPoints ][ iS ]);
+
+ ++nbBranchPoints;
+ }
+
+ if ( nbBranchPoints != 2 )
+ return false;
+
+ // add to theSinuEdges all edges that are not theShortEdges
+ vector< vector<TopoDS_Edge> > sinuEdges(1);
+ TopoDS_Vertex vCommon;
+ for ( size_t i = 0; i < allEdges.size(); ++i )
+ {
+ if ( !shortMap.Contains( allEdges[i] ))
+ {
+ if ( !sinuEdges.back().empty() )
+ if ( !TopExp::CommonVertex( sinuEdges.back().back(), allEdges[ i ], vCommon ))
+ sinuEdges.resize( sinuEdges.size() + 1 );
+
+ sinuEdges.back().push_back( allEdges[i] );
+ }
+ }
+ if ( sinuEdges.size() == 3 )
+ {
+ if ( !TopExp::CommonVertex( sinuEdges.back().back(), sinuEdges[0][0], vCommon ))
+ return false;
+ vector<TopoDS_Edge>& last = sinuEdges.back();
+ last.insert( last.end(), sinuEdges[0].begin(), sinuEdges[0].end() );
+ sinuEdges[0].swap( last );
+ sinuEdges.resize( 2 );
+ }
+ if ( sinuEdges.size() != 2 )
+ return false;
+
+ theSinuEdges[0].swap( sinuEdges[0] );
+ theSinuEdges[1].swap( sinuEdges[1] );
+
+ if ( !TopExp::CommonVertex( theSinuEdges[0].back(), theShortEdges[0][0], vCommon ) ||
+ !vCommon.IsSame( theHelper.IthVertex( 1, theSinuEdges[0].back() )))
+ theShortEdges[0].swap( theShortEdges[1] );
+
+ return ( theShortEdges[0].size() > 0 && theShortEdges[1].size() > 0 &&
+ theSinuEdges [0].size() > 0 && theSinuEdges [1].size() > 0 );
+
+ // the sinuous EDGEs can be composite and C0 continuous,
+ // therefor we use a complex criterion to find TWO short non-sinuous EDGEs
+ // and the rest EDGEs will be treated as sinuous.
+ // A short edge should have the following features:
+ // a) straight
+ // b) short
+ // c) with convex corners at ends
+ // d) far from the other short EDGE
+
+ // vector< double > isStraightEdge( nbEdges, 0 ); // criterion value
+
+ // // a0) evaluate continuity
+ // const double contiWgt = 0.5; // weight of continuity in the criterion
+ // multimap< int, TopoDS_Edge > continuity;
+ // for ( size_t i = 0; i < nbEdges; ++I )
+ // {
+ // BRepAdaptor_Curve curve( allEdges[i] );
+ // GeomAbs_Shape C = GeomAbs_CN;
+ // try:
+ // C = curve.Continuity(); // C0, G1, C1, G2, C2, C3, CN
+ // catch ( Standard_Failure ) {}
+ // continuity.insert( make_pair( C, allEdges[i] ));
+ // isStraight[i] += double( C ) / double( CN ) * contiWgt;
+ // }
+
+ // // try to choose by continuity
+ // int mostStraight = (int) continuity.rbegin()->first;
+ // int lessStraight = (int) continuity.begin()->first;
+ // if ( mostStraight != lessStraight )
+ // {
+ // int nbStraight = continuity.count( mostStraight );
+ // if ( nbStraight == 2 )
+ // {
+ // getTwo( /*least=*/false, continuity, theShortEdges, theSinuEdges );
+ // }
+ // else if ( nbStraight == 3 && nbEdges == 4 )
+ // {
+ // theSinuEdges.push_back( continuity.begin()->second );
+ // vector<TopoDS_Edge>::iterator it =
+ // std::find( allEdges.begin(), allEdges.end(), theSinuEdges[0] );
+ // int i = std::distance( allEdges.begin(), it );
+ // theSinuEdges .push_back( allEdges[( i+2 )%4 ]);
+ // theShortEdges.push_back( allEdges[( i+1 )%4 ]);
+ // theShortEdges.push_back( allEdges[( i+3 )%4 ]);
+ // }
+ // if ( theShortEdges.size() == 2 )
+ // return true;
+ // }
+
+ // // a) curvature; evaluate aspect ratio
+ // {
+ // const double curvWgt = 0.5;
+ // for ( size_t i = 0; i < nbEdges; ++I )
+ // {
+ // BRepAdaptor_Curve curve( allEdges[i] );
+ // double curvature = 1;
+ // if ( !curve.IsClosed() )
+ // {
+ // const double f = curve.FirstParameter(), l = curve.LastParameter();
+ // gp_Pnt pf = curve.Value( f ), pl = curve.Value( l );
+ // gp_Lin line( pf, pl.XYZ() - pf.XYZ() );
+ // double distMax = 0;
+ // for ( double u = f; u < l; u += (l-f)/30. )
+ // distMax = Max( distMax, line.SquareDistance( curve.Value( u )));
+ // curvature = Sqrt( distMax ) / ( pf.Distance( pl ));
+ // }
+ // isStraight[i] += curvWgt / ( curvature + 1e-20 );
+ // }
+ // }
+ // // b) length
+ // {
+ // const double lenWgt = 0.5;
+ // for ( size_t i = 0; i < nbEdges; ++I )
+ // {
+ // double length = SMESH_Algo::Length( allEdges[i] );
+ // if ( length > 0 )
+ // isStraight[i] += lenWgt / length;
+ // }
+ // }
+ // // c) with convex corners at ends
+ // {
+ // const double cornerWgt = 0.25;
+ // for ( size_t i = 0; i < nbEdges; ++I )
+ // {
+ // double convex = 0;
+ // int iPrev = SMESH_MesherHelper::WrapIndex( int(i)-1, nbEdges );
+ // int iNext = SMESH_MesherHelper::WrapIndex( int(i)+1, nbEdges );
+ // TopoDS_Vertex v = helper.IthVertex( 0, allEdges[i] );
+ // double angle = SMESH_MesherHelper::GetAngle( allEdges[iPrev], allEdges[i], theFace, v );
+ // if ( angle < M_PI ) // [-PI; PI]
+ // convex += ( angle + M_PI ) / M_PI / M_PI;
+ // v = helper.IthVertex( 1, allEdges[i] );
+ // angle = SMESH_MesherHelper::GetAngle( allEdges[iNext], allEdges[i], theFace, v );
+ // if ( angle < M_PI ) // [-PI; PI]
+ // convex += ( angle + M_PI ) / M_PI / M_PI;
+ // isStraight[i] += cornerWgt * convex;
+ // }
+ // }
+ }
+
+ //================================================================================
+ /*!
+ * \brief Creates an EDGE from a sole branch of MA
+ */
+ //================================================================================
+
+ TopoDS_Edge makeEdgeFromMA( SMESH_MesherHelper& theHelper,
+ const SMESH_MAT2d::MedialAxis& theMA )
+ {
+ if ( theMA.getBranches().size() != 1 )
+ return TopoDS_Edge();
+
+ vector< gp_XY > uv;
+ theMA.getPoints( theMA.getBranches()[0], uv );
+ if ( uv.size() < 2 )
+ return TopoDS_Edge();
+
+ TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() );
+ Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
+
+ // cout << "from salome.geom import geomBuilder" << endl;
+ // cout << "geompy = geomBuilder.New(salome.myStudy)" << endl;
+ Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, uv.size());
+ for ( size_t i = 0; i < uv.size(); ++i )
+ {
+ gp_Pnt p = surface->Value( uv[i].X(), uv[i].Y() );
+ points->SetValue( i+1, p );
+ //cout << "geompy.MakeVertex( "<< p.X()<<", " << p.Y()<<", " << p.Z()<<" )" << endl;
+ }
+
+ GeomAPI_Interpolate interpol( points, /*isClosed=*/false, gp::Resolution());
+ interpol.Perform();
+ if ( !interpol.IsDone())
+ return TopoDS_Edge();
+
+ TopoDS_Edge branchEdge = BRepBuilderAPI_MakeEdge(interpol.Curve());
+ return branchEdge;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Returns a type of shape, to which a hypothesis used to mesh a given edge is assigned
+ */
+ //================================================================================
+
+ TopAbs_ShapeEnum getHypShape( SMESH_Mesh* mesh, const TopoDS_Shape& edge )
+ {
+ TopAbs_ShapeEnum shapeType = TopAbs_SHAPE;
+
+ SMESH_subMesh* sm = mesh->GetSubMesh( edge );
+ SMESH_Algo* algo = sm->GetAlgo();
+ if ( !algo ) return shapeType;
+
+ const list <const SMESHDS_Hypothesis *> & hyps =
+ algo->GetUsedHypothesis( *mesh, edge, /*ignoreAuxiliary=*/true );
+ if ( hyps.empty() ) return shapeType;
+
+ TopoDS_Shape shapeOfHyp =
+ SMESH_MesherHelper::GetShapeOfHypothesis( hyps.front(), edge, mesh);
+
+ return SMESH_MesherHelper::GetGroupType( shapeOfHyp, /*woCompound=*/true);
+ }
+
+ //================================================================================
+ /*!
+ * \brief Discretize a sole branch of MA an returns parameters of divisions on MA
+ */
+ //================================================================================
+
+ bool divideMA( SMESH_MesherHelper& theHelper,
+ const SMESH_MAT2d::MedialAxis& theMA,
+ const vector<TopoDS_Edge>& theSinuEdges,
+ const size_t theSinuSide0Size,
+ SMESH_Algo* the1dAlgo,
+ vector<double>& theMAParams )
+ {
+ // check if all EDGEs of one size are meshed, then MA discretization is not needed
+ SMESH_Mesh* mesh = theHelper.GetMesh();
+ size_t nbComputedEdges[2] = { 0, 0 };
+ for ( size_t i = 1; i < theSinuEdges.size(); ++i )
+ {
+ bool isComputed = ( ! mesh->GetSubMesh( theSinuEdges[i] )->IsEmpty() );
+ nbComputedEdges[ i < theSinuSide0Size ] += isComputed;
+ }
+ if ( nbComputedEdges[0] == theSinuSide0Size ||
+ nbComputedEdges[1] == theSinuEdges.size() - theSinuSide0Size )
+ return true; // discretization is not needed
+
+
+ TopoDS_Edge branchEdge = makeEdgeFromMA( theHelper, theMA );
+ if ( branchEdge.IsNull() )
+ return false;
+
+ // const char* file = "/misc/dn25/salome/eap/salome/misc/tmp/MAedge.brep";
+ // BRepTools::Write( branchEdge, file);
+ // cout << "Write " << file << endl;
+
+ // look for a most local hyps assigned to theSinuEdges
+ TopoDS_Edge edge = theSinuEdges[0];
+ int mostSimpleShape = (int) getHypShape( mesh, edge );
+ for ( size_t i = 1; i < theSinuEdges.size(); ++i )
+ {
+ int shapeType = (int) getHypShape( mesh, theSinuEdges[i] );
+ if ( shapeType > mostSimpleShape )
+ edge = theSinuEdges[i];
+ }
+
+ SMESH_Algo* algo = the1dAlgo;
+ if ( mostSimpleShape != TopAbs_SHAPE )
+ {
+ algo = mesh->GetSubMesh( edge )->GetAlgo();
+ SMESH_Hypothesis::Hypothesis_Status status;
+ if ( !algo->CheckHypothesis( *mesh, edge, status ))
+ algo = the1dAlgo;
+ }
+
+ TmpMesh tmpMesh;
+ tmpMesh.ShapeToMesh( branchEdge );
+ try {
+ mesh->GetGen()->Compute( tmpMesh, branchEdge, true, true ); // make nodes on VERTEXes
+ if ( !algo->Compute( tmpMesh, branchEdge ))
+ return false;
+ }
+ catch (...) {
+ return false;
+ }
+ return SMESH_Algo::GetNodeParamOnEdge( tmpMesh.GetMeshDS(), branchEdge, theMAParams );
+ }
+
+ //================================================================================
+ /*!
+ * \brief Modifies division parameters on MA to make them coincide with projections
+ * of VERTEXes to MA for a given pair of opposite EDGEs
+ * \param [in] theEdgePairInd - index of the EDGE pair
+ * \param [in] theDivPoints - the BranchPoint's dividing MA into parts each
+ * corresponding to a unique pair of opposite EDGEs
+ * \param [in,out] theMAParams - the MA division parameters to modify
+ * \param [in,out] theParBeg - index of the 1st division point for the given EDGE pair
+ * \param [in,out] theParEnd - index of the last division point for the given EDGE pair
+ * \return bool - is OK
+ */
+ //================================================================================
+
+ bool getParamsForEdgePair( const size_t theEdgePairInd,
+ const vector< SMESH_MAT2d::BranchPoint >& theDivPoints,
+ const vector<double>& theMAParams,
+ vector<double>& theSelectedMAParams)
+ {
+ if ( theDivPoints.empty() )
+ {
+ theSelectedMAParams = theMAParams;
+ return true;
+ }
+ if ( theEdgePairInd > theDivPoints.size() )
+ return false;
+
+ // TODO
+ return false;
+ }
+
+ //--------------------------------------------------------------------------------
+ // node or node parameter on EDGE
+ struct NodePoint
+ {
+ const SMDS_MeshNode* _node;
+ double _u;
+ int _edgeInd; // index in theSinuEdges vector
+
+ NodePoint(const SMDS_MeshNode* n=0 ): _node(n), _u(0), _edgeInd(-1) {}
+ NodePoint(double u, size_t iEdge) : _node(0), _u(u), _edgeInd(iEdge) {}
+ NodePoint(const SMESH_MAT2d::BoundaryPoint& p) : _node(0), _u(p._param), _edgeInd(p._edgeIndex) {}
+ };
+
+ //================================================================================
+ /*!
+ * \brief Finds a VERTEX corresponding to a point on EDGE, which is also filled
+ * with a node on the VERTEX, present or created
+ * \param [in,out] theNodePnt - the node position on the EDGE
+ * \param [in] theSinuEdges - the sinuous EDGEs
+ * \param [in] theMeshDS - the mesh
+ * \return bool - true if the \a theBndPnt is on VERTEX
+ */
+ //================================================================================
+
+ bool findVertex( NodePoint& theNodePnt,
+ const vector<TopoDS_Edge>& theSinuEdges,
+ SMESHDS_Mesh* theMeshDS)
+ {
+ if ( theNodePnt._edgeInd >= theSinuEdges.size() )
+ return false;
+
+ double f,l;
+ BRep_Tool::Range( theSinuEdges[ theNodePnt._edgeInd ], f,l );
+
+ TopoDS_Vertex V;
+ if ( Abs( f - theNodePnt._u ))
+ V = SMESH_MesherHelper::IthVertex( 0, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false);
+ else if ( Abs( l - theNodePnt._u ))
+ V = SMESH_MesherHelper::IthVertex( 1, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false);
+
+ if ( !V.IsNull() )
+ {
+ theNodePnt._node = SMESH_Algo::VertexNode( V, theMeshDS );
+ if ( !theNodePnt._node )
+ {
+ gp_Pnt p = BRep_Tool::Pnt( V );
+ theNodePnt._node = theMeshDS->AddNode( p.X(), p.Y(), p.Z() );
+ theMeshDS->SetNodeOnVertex( theNodePnt._node, V );
+ }
+ return true;
+ }
+ return false;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Add to the map of NodePoint's those on VERTEXes
+ * \param [in,out] theHelper - the helper
+ * \param [in] theMA - Medial Axis
+ * \param [in] theDivPoints - projections of VERTEXes to MA
+ * \param [in] theSinuEdges - the sinuous EDGEs
+ * \param [in] theSideEdgeIDs - indices of sinuous EDGEs per side
+ * \param [in] theIsEdgeComputed - is sinuous EGDE is meshed
+ * \param [in,out] thePointsOnE - the map to fill
+ */
+ //================================================================================
+
+ bool projectVertices( SMESH_MesherHelper& theHelper,
+ const SMESH_MAT2d::MedialAxis& theMA,
+ const vector< SMESH_MAT2d::BranchPoint >& theDivPoints,
+ const vector<TopoDS_Edge>& theSinuEdges,
+ //const vector< int > theSideEdgeIDs[2],
+ const vector< bool >& theIsEdgeComputed,
+ map< double, pair< NodePoint, NodePoint > > & thePointsOnE)
+ {
+ if ( theDivPoints.empty() )
+ return true;
+
+ SMESHDS_Mesh* meshDS = theHelper.GetMeshDS();
+
+ double uMA;
+ SMESH_MAT2d::BoundaryPoint bp[2];
+ const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0];
+
+ for ( size_t i = 0; i < theDivPoints.size(); ++i )
+ {
+ if ( !branch.getParameter( theDivPoints[i], uMA ))
+ return false;
+ if ( !branch.getBoundaryPoints( theDivPoints[i], bp[0], bp[1] ))
+ return false;
+
+ NodePoint np[2] = { NodePoint( bp[0] ),
+ NodePoint( bp[1] ) };
+ bool isVertex[2] = { findVertex( np[0], theSinuEdges, meshDS ),
+ findVertex( np[1], theSinuEdges, meshDS )};
+
+ map< double, pair< NodePoint, NodePoint > >::iterator u2NP =
+ thePointsOnE.insert( make_pair( uMA, make_pair( np[0], np[1]))).first;
+
+ if ( isVertex[0] && isVertex[1] )
+ continue;
+
+ bool isOppComputed = theIsEdgeComputed[ np[ isVertex[0] ]._edgeInd ];
+
+ if ( !isOppComputed )
+ continue;
+
+ // a VERTEX is projected on a meshed EDGE; there are two options:
+ // - a projected point is joined with a closet node if a strip between this and neighbor
+ // projection is wide enough; joining is done by setting the same node to the BoundaryPoint
+ // - a neighbor projection is merged this this one if it too close; a node of deleted
+ // projection is set to the BoundaryPoint of this projection
+
+
+ }
+ return true;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Divide the sinuous EDGEs by projecting the division point of Medial
+ * Axis to the EGDEs
+ * \param [in] theHelper - the helper
+ * \param [in] theMA - the Medial Axis
+ * \param [in] theMAParams - parameters of division points of \a theMA
+ * \param [in] theSinuEdges - the EDGEs to make nodes on
+ * \param [in] theSinuSide0Size - the number of EDGEs in the 1st sinuous side
+ * \return bool - is OK or not
+ */
+ //================================================================================
+
+ bool computeSinuEdges( SMESH_MesherHelper& theHelper,
+ SMESH_MAT2d::MedialAxis& theMA,
+ vector<double>& theMAParams,
+ const vector<TopoDS_Edge>& theSinuEdges,
+ const size_t theSinuSide0Size)
+ {
+ if ( theMA.getBranches().size() != 1 )
+ return false;
+
+ // normalize theMAParams
+ for ( size_t i = 0; i < theMAParams.size(); ++i )
+ theMAParams[i] /= theMAParams.back();
+
+
+ SMESH_Mesh* mesh = theHelper.GetMesh();
+ SMESHDS_Mesh* meshDS = theHelper.GetMeshDS();
+ double f,l;
+
+ vector< Handle(Geom_Curve) > curves ( theSinuEdges.size() );
+ vector< int > edgeIDs( theSinuEdges.size() );
+ vector< bool > isComputed( theSinuEdges.size() );
+ //bool hasComputed = false;
+ for ( size_t i = 0; i < theSinuEdges.size(); ++i )
+ {
+ curves[i] = BRep_Tool::Curve( theSinuEdges[i], f,l );
+ if ( !curves[i] )
+ return false;
+ SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] );
+ edgeIDs [i] = sm->GetId();
+ isComputed[i] = ( !sm->IsEmpty() );
+ // if ( isComputed[i] )
+ // hasComputed = true;
+ }
+
+ const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0];
+ SMESH_MAT2d::BoundaryPoint bp[2];
+
+ vector< std::size_t > edgeIDs1, edgeIDs2;
+ vector< SMESH_MAT2d::BranchPoint > divPoints;
+ branch.getOppositeGeomEdges( edgeIDs1, edgeIDs2, divPoints );
+ for ( size_t i = 0; i < edgeIDs1.size(); ++i )
+ if ( isComputed[ edgeIDs1[i]] &&
+ isComputed[ edgeIDs2[i]])
+ return false;
+
+ // map param on MA to parameters of nodes on a pair of theSinuEdges
+ typedef map< double, pair< NodePoint, NodePoint > > TMAPar2NPoints;
+ TMAPar2NPoints pointsOnE;
+ vector<double> maParams;
+
+ // compute params of nodes on EDGEs by projecting division points from MA
+ //const double tol = 1e-5 * theMAParams.back();
+ size_t iEdgePair = 0;
+ while ( iEdgePair < edgeIDs1.size() )
+ {
+ if ( isComputed[ edgeIDs1[ iEdgePair ]] ||
+ isComputed[ edgeIDs2[ iEdgePair ]])
+ {
+ // "projection" from one side to the other
+
+ size_t iEdgeComputed = edgeIDs1[iEdgePair], iSideComputed = 0;
+ if ( !isComputed[ iEdgeComputed ])
+ ++iSideComputed, iEdgeComputed = edgeIDs2[iEdgePair];
+
+ map< double, const SMDS_MeshNode* > nodeParams; // params of existing nodes
+ if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iEdgeComputed ], /*skipMedium=*/true, nodeParams ))
+ return false;
+
+ SMESH_MAT2d::BoundaryPoint& bndPnt = bp[ 1-iSideComputed ];
+ SMESH_MAT2d::BranchPoint brp;
+ NodePoint npN, npB;
+ NodePoint& np0 = iSideComputed ? npB : npN;
+ NodePoint& np1 = iSideComputed ? npN : npB;
+
+ double maParam1st, maParamLast, maParam;
+ if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.begin()->first, brp ))
+ return false;
+ branch.getParameter( brp, maParam1st );
+ if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.rbegin()->first, brp ))
+ return false;
+ branch.getParameter( brp, maParamLast );
+
+ map< double, const SMDS_MeshNode* >::iterator u2n = nodeParams.begin(), u2nEnd = --nodeParams.end();
+ TMAPar2NPoints::iterator pos, end = pointsOnE.end();
+ TMAPar2NPoints::iterator & hint = (maParamLast > maParam1st) ? end : pos;
+ for ( ++u2n; u2n != u2nEnd; ++u2n )
+ {
+ if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, u2n->first, brp ))
+ return false;
+ if ( !branch.getBoundaryPoints( brp, bp[0], bp[1] ))
+ return false;
+ if ( !branch.getParameter( brp, maParam ))
+ return false;
+
+ npN = NodePoint( u2n->second );
+ npB = NodePoint( bndPnt );
+ pos = pointsOnE.insert( hint, make_pair( maParam, make_pair( np0, np1 )));
+ }
+
+ // move iEdgePair forward
+ while ( iEdgePair < edgeIDs1.size() )
+ if ( edgeIDs1[ iEdgePair ] == bp[0]._edgeIndex &&
+ edgeIDs2[ iEdgePair ] == bp[1]._edgeIndex )
+ break;
+ else
+ ++iEdgePair;
+ }
+ else
+ {
+ // projection from MA
+ maParams.clear();
+ if ( !getParamsForEdgePair( iEdgePair, divPoints, theMAParams, maParams ))
+ return false;
+
+ for ( size_t i = 1; i < maParams.size()-1; ++i )
+ {
+ if ( !branch.getBoundaryPoints( maParams[i], bp[0], bp[1] ))
+ return false;
+
+ pointsOnE.insert( pointsOnE.end(), make_pair( maParams[i], make_pair( NodePoint(bp[0]),
+ NodePoint(bp[1]))));
+ }
+ }
+ ++iEdgePair;
+ }
+
+ if ( !projectVertices( theHelper, theMA, divPoints, theSinuEdges, isComputed, pointsOnE ))
+ return false;
+
+ // create nodes
+ TMAPar2NPoints::iterator u2np = pointsOnE.begin();
+ for ( ; u2np != pointsOnE.end(); ++u2np )
+ {
+ NodePoint* np[2] = { & u2np->second.first, & u2np->second.second };
+ for ( int iSide = 0; iSide < 2; ++iSide )
+ {
+ if ( np[ iSide ]->_node ) continue;
+ size_t iEdge = np[ iSide ]->_edgeInd;
+ double u = np[ iSide ]->_u;
+ gp_Pnt p = curves[ iEdge ]->Value( u );
+ np[ iSide ]->_node = meshDS->AddNode( p.X(), p.Y(), p.Z() );
+ meshDS->SetNodeOnEdge( np[ iSide ]->_node, edgeIDs[ iEdge ], u );
+ }
+ }
+
+ // create mesh segments on EDGEs
+ theHelper.SetElementsOnShape( false );
+ TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() );
+ for ( size_t i = 0; i < theSinuEdges.size(); ++i )
+ {
+ SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] );
+ if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 )
+ continue;
+
+ StdMeshers_FaceSide side( face, theSinuEdges[i], mesh,
+ /*isFwd=*/true, /*skipMediumNodes=*/true );
+ vector<const SMDS_MeshNode*> nodes = side.GetOrderedNodes();
+ for ( size_t in = 1; in < nodes.size(); ++in )
+ {
+ const SMDS_MeshElement* seg = theHelper.AddEdge( nodes[in-1], nodes[in], 0, false );
+ meshDS->SetMeshElementOnShape( seg, edgeIDs[ i ] );
+ }
+ }
+
+ // update sub-meshes on VERTEXes
+ for ( size_t i = 0; i < theSinuEdges.size(); ++i )
+ {
+ mesh->GetSubMesh( theHelper.IthVertex( 0, theSinuEdges[i] ))
+ ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
+ mesh->GetSubMesh( theHelper.IthVertex( 1, theSinuEdges[i] ))
+ ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
+ }
+
+ return true;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Mesh short EDGEs
+ */
+ //================================================================================
+
+ bool computeShortEdges( SMESH_MesherHelper& theHelper,
+ const vector<TopoDS_Edge>& theShortEdges,
+ SMESH_Algo* the1dAlgo )
+ {
+ for ( size_t i = 0; i < theShortEdges.size(); ++i )
+ {
+ theHelper.GetGen()->Compute( *theHelper.GetMesh(), theShortEdges[i], true, true );
+
+ SMESH_subMesh* sm = theHelper.GetMesh()->GetSubMesh(theShortEdges[i] );
+ if ( sm->IsEmpty() )
+ {
+ try {
+ if ( !the1dAlgo->Compute( *theHelper.GetMesh(), theShortEdges[i] ))
+ return false;
+ }
+ catch (...) {
+ return false;
+ }
+ sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
+ if ( sm->IsEmpty() )
+ return false;
+ }
+ }
+ return true;
+ }
+
+ inline double area( const UVPtStruct& p1, const UVPtStruct& p2, const UVPtStruct& p3 )
+ {
+ gp_XY v1 = p2.UV() - p1.UV();
+ gp_XY v2 = p3.UV() - p1.UV();
+ return v2 ^ v1;
+ }
+
+ bool ellipticSmooth( FaceQuadStruct::Ptr quad, int nbLoops )
+ {
+ //nbLoops = 10;
+ if ( quad->uv_grid.empty() )
+ return true;
+
+ int nbhoriz = quad->iSize;
+ int nbvertic = quad->jSize;
+
+ const double dksi = 0.5, deta = 0.5;
+ const double dksi2 = dksi*dksi, deta2 = deta*deta;
+ double err = 0., g11, g22, g12;
+ int nbErr = 0;
+
+ FaceQuadStruct& q = *quad;
+ UVPtStruct pNew;
+
+ double refArea = area( q.UVPt(0,0), q.UVPt(1,0), q.UVPt(1,1) );
+
+ for ( int iLoop = 0; iLoop < nbLoops; ++iLoop )
+ {
+ err = 0;
+ for ( int i = 1; i < nbhoriz - 1; i++ )
+ for ( int j = 1; j < nbvertic - 1; j++ )
+ {
+ g11 = ( (q.U(i,j+1) - q.U(i,j-1))*(q.U(i,j+1) - q.U(i,j-1))/dksi2 +
+ (q.V(i,j+1) - q.V(i,j-1))*(q.V(i,j+1) - q.V(i,j-1))/deta2 )/4;
+
+ g22 = ( (q.U(i+1,j) - q.U(i-1,j))*(q.U(i+1,j) - q.U(i-1,j))/dksi2 +
+ (q.V(i+1,j) - q.V(i-1,j))*(q.V(i+1,j) - q.V(i-1,j))/deta2 )/4;
+
+ g12 = ( (q.U(i+1,j) - q.U(i-1,j))*(q.U(i,j+1) - q.U(i,j-1))/dksi2 +
+ (q.V(i+1,j) - q.V(i-1,j))*(q.V(i,j+1) - q.V(i,j-1))/deta2 )/(4*dksi*deta);
+
+ pNew.u = dksi2/(2*(g11+g22)) * (g11*(q.U(i+1,j) + q.U(i-1,j))/dksi2 +
+ g22*(q.U(i,j+1) + q.U(i,j-1))/dksi2
+ - 0.5*g12*q.U(i+1,j+1) + 0.5*g12*q.U(i-1,j+1) +
+ - 0.5*g12*q.U(i-1,j-1) + 0.5*g12*q.U(i+1,j-1));
+
+ pNew.v = deta2/(2*(g11+g22)) * (g11*(q.V(i+1,j) + q.V(i-1,j))/deta2 +
+ g22*(q.V(i,j+1) + q.V(i,j-1))/deta2
+ - 0.5*g12*q.V(i+1,j+1) + 0.5*g12*q.V(i-1,j+1) +
+ - 0.5*g12*q.V(i-1,j-1) + 0.5*g12*q.V(i+1,j-1));
+
+ // if (( refArea * area( q.UVPt(i-1,j-1), q.UVPt(i,j-1), pNew ) > 0 ) &&
+ // ( refArea * area( q.UVPt(i+1,j-1), q.UVPt(i+1,j), pNew ) > 0 ) &&
+ // ( refArea * area( q.UVPt(i+1,j+1), q.UVPt(i,j+1), pNew ) > 0 ) &&
+ // ( refArea * area( q.UVPt(i-1,j), q.UVPt(i-1,j-1), pNew ) > 0 ))
+ {
+ err += sqrt(( q.U(i,j) - pNew.u ) * ( q.U(i,j) - pNew.u ) +
+ ( q.V(i,j) - pNew.v ) * ( q.V(i,j) - pNew.v ));
+ q.U(i,j) = pNew.u;
+ q.V(i,j) = pNew.v;
+ }
+ // else if ( ++nbErr < 10 )
+ // {
+ // cout << i << ", " << j << endl;
+ // cout << "x = ["
+ // << "[ " << q.U(i-1,j-1) << ", " <<q.U(i,j-1) << ", " << q.U(i+1,j-1) << " ],"
+ // << "[ " << q.U(i-1,j-0) << ", " <<q.U(i,j-0) << ", " << q.U(i+1,j-0) << " ],"
+ // << "[ " << q.U(i-1,j+1) << ", " <<q.U(i,j+1) << ", " << q.U(i+1,j+1) << " ]]" << endl;
+ // cout << "y = ["
+ // << "[ " << q.V(i-1,j-1) << ", " <<q.V(i,j-1) << ", " << q.V(i+1,j-1) << " ],"
+ // << "[ " << q.V(i-1,j-0) << ", " <<q.V(i,j-0) << ", " << q.V(i+1,j-0) << " ],"
+ // << "[ " << q.V(i-1,j+1) << ", " <<q.V(i,j+1) << ", " << q.V(i+1,j+1) << " ]]" << endl<<endl;
+ // }
+ }
+
+ if ( err / ( nbhoriz - 2 ) / ( nbvertic - 2 ) < 1e-6 )
+ break;
+ }
+ //cout << " ERR " << err / ( nbhoriz - 2 ) / ( nbvertic - 2 ) << endl;
+
+ return true;
+ }
+
+
+} // namespace
+
+//================================================================================
+/*!
+ * \brief Create quadrangle elements
+ * \param [in] theHelper - the helper
+ * \param [in] theFace - the face to mesh
+ * \param [in] theSinuEdges - the sinuous EDGEs
+ * \param [in] theShortEdges - the short EDGEs
+ * \return bool - is OK or not
+ */
+//================================================================================
+
+bool StdMeshers_QuadFromMedialAxis_1D2D::computeQuads( SMESH_MesherHelper& theHelper,
+ const TopoDS_Face& theFace,
+ const vector<TopoDS_Edge> theSinuEdges[2],
+ const vector<TopoDS_Edge> theShortEdges[2])
+{
+ SMESH_Mesh* mesh = theHelper.GetMesh();
+ SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( *mesh, theFace );
+ if ( !proxyMesh )
+ return false;
+
+ StdMeshers_Quadrangle_2D::myProxyMesh = proxyMesh;
+ StdMeshers_Quadrangle_2D::myHelper = &theHelper;
+ StdMeshers_Quadrangle_2D::myNeedSmooth = false;
+ StdMeshers_Quadrangle_2D::myCheckOri = false;
+ StdMeshers_Quadrangle_2D::myQuadList.clear();
+
+ // fill FaceQuadStruct
+
+ list< TopoDS_Edge > side[4];
+ side[0].insert( side[0].end(), theShortEdges[0].begin(), theShortEdges[0].end() );
+ side[1].insert( side[1].end(), theSinuEdges[1].begin(), theSinuEdges[1].end() );
+ side[2].insert( side[2].end(), theShortEdges[1].begin(), theShortEdges[1].end() );
+ side[3].insert( side[3].end(), theSinuEdges[0].begin(), theSinuEdges[0].end() );
+
+ FaceQuadStruct::Ptr quad( new FaceQuadStruct );
+ quad->side.resize( 4 );
+ quad->face = theFace;
+ for ( int i = 0; i < 4; ++i )
+ {
+ quad->side[i] = StdMeshers_FaceSide::New( theFace, side[i], mesh, i < QUAD_TOP_SIDE,
+ /*skipMediumNodes=*/true, proxyMesh );
+ }
+ int nbNodesShort0 = quad->side[0].NbPoints();
+ int nbNodesShort1 = quad->side[2].NbPoints();
+
+ // compute UV of internal points
+ myQuadList.push_back( quad );
+ if ( !StdMeshers_Quadrangle_2D::setNormalizedGrid( quad ))
+ return false;
+
+ // elliptic smooth of internal points to get boundary cell normal to the boundary
+ ellipticSmooth( quad, 1 );
+
+ // create quadrangles
+ bool ok;
+ if ( nbNodesShort0 == nbNodesShort1 )
+ ok = StdMeshers_Quadrangle_2D::computeQuadDominant( *mesh, theFace, quad );
+ else
+ ok = StdMeshers_Quadrangle_2D::computeTriangles( *mesh, theFace, quad );
+
+ StdMeshers_Quadrangle_2D::myProxyMesh.reset();
+ StdMeshers_Quadrangle_2D::myHelper = 0;
+
+ return ok;
+}
+
+//================================================================================
+/*!
+ * \brief Generate quadrangle mesh
+ */
+//================================================================================
+
+bool StdMeshers_QuadFromMedialAxis_1D2D::Compute(SMESH_Mesh& theMesh,
+ const TopoDS_Shape& theShape)
+{
+ SMESH_MesherHelper helper( theMesh );
+ helper.SetSubShape( theShape );
+
+ TopoDS_Face F = TopoDS::Face( theShape );
+ if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD );
+
+ list< TopoDS_Edge > edges;
+ list< int > nbEdgesInWire;
+ int nbWire = SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire);
+
+ vector< TopoDS_Edge > sinuSide[2], shortSide[2];
+ if ( nbWire == 1 && getSinuousEdges( helper, F, edges, sinuSide, shortSide ))
+ {
+ vector< TopoDS_Edge > sinuEdges = sinuSide[0];
+ sinuEdges.insert( sinuEdges.end(), sinuSide[1].begin(), sinuSide[1].end() );
+ if ( sinuEdges.size() > 2 )
+ return error(COMPERR_BAD_SHAPE, "Not yet supported case" );
+
+ double minSegLen = getMinSegLen( helper, sinuEdges );
+ SMESH_MAT2d::MedialAxis ma( F, sinuEdges, minSegLen, /*ignoreCorners=*/true );
+
+ if ( !_regular1D )
+ _regular1D = new Algo1D( _studyId, _gen );
+ _regular1D->SetSegmentLength( minSegLen );
+
+ vector<double> maParams;
+ if ( ! divideMA( helper, ma, sinuEdges, sinuSide[0].size(), _regular1D, maParams ))
+ return error(COMPERR_BAD_SHAPE);
+
+ if ( !computeShortEdges( helper, shortSide[0], _regular1D ) ||
+ !computeShortEdges( helper, shortSide[1], _regular1D ))
+ return error("Failed to mesh short edges");
+
+ if ( !computeSinuEdges( helper, ma, maParams, sinuEdges, sinuSide[0].size() ))
+ return error("Failed to mesh sinuous edges");
+
+ return computeQuads( helper, F, sinuSide, shortSide );
+ }
+
+ return error(COMPERR_BAD_SHAPE, "Not implemented so far");
+}
+
+//================================================================================
+/*!
+ * \brief Predict nb of elements
+ */
+//================================================================================
+
+bool StdMeshers_QuadFromMedialAxis_1D2D::Evaluate(SMESH_Mesh & theMesh,
+ const TopoDS_Shape & theShape,
+ MapShapeNbElems& theResMap)
+{
+ return StdMeshers_Quadrangle_2D::Evaluate(theMesh,theShape,theResMap);
+}
+
--- /dev/null
+// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File : StdMeshers_QuadFromMedialAxis_1D2D.hxx
+// Created : Wed Jun 3 17:22:35 2015
+// Author : Edward AGAPOV (eap)
+
+
+#ifndef __StdMeshers_QuadFromMedialAxis_1D2D_HXX__
+#define __StdMeshers_QuadFromMedialAxis_1D2D_HXX__
+
+#include "StdMeshers_Quadrangle_2D.hxx"
+
+#include <vector>
+
+/*!
+ * \brief Quadrangle mesher using Medial Axis
+ */
+class STDMESHERS_EXPORT StdMeshers_QuadFromMedialAxis_1D2D: public StdMeshers_Quadrangle_2D
+{
+ public:
+ StdMeshers_QuadFromMedialAxis_1D2D(int hypId, int studyId, SMESH_Gen* gen);
+ virtual ~StdMeshers_QuadFromMedialAxis_1D2D();
+
+ virtual bool CheckHypothesis(SMESH_Mesh& aMesh,
+ const TopoDS_Shape& aShape,
+ Hypothesis_Status& aStatus);
+
+ virtual bool Compute(SMESH_Mesh& aMesh,
+ const TopoDS_Shape& aShape);
+
+ virtual bool Evaluate(SMESH_Mesh & aMesh,
+ const TopoDS_Shape & aShape,
+ MapShapeNbElems& aResMap);
+
+ private:
+
+ bool computeQuads( SMESH_MesherHelper& theHelper,
+ const TopoDS_Face& theFace,
+ const std::vector<TopoDS_Edge> theSinuEdges[2],
+ const std::vector<TopoDS_Edge> theShortEdges[2]);
+
+ class Algo1D;
+ Algo1D* _regular1D;
+};
+
+#endif
//================================================================================
/*!
- * \brief Rotate sides of a quad by given nb of quartes
+ * \brief Rotate sides of a quad CCW by given nb of quartes
* \param nb - number of rotation quartes
* \param ori - to keep orientation of sides as in an unit quad or not
* \param keepGrid - if \c true Side::grid is not changed, Side::from and Side::to
{
if ( nb == 0 ) return;
+ nb = nb % NB_QUAD_SIDES;
+
vector< Side > newSides( side.size() );
vector< Side* > sidePtrs( side.size() );
for (int i = QUAD_BOTTOM_SIDE; i < NB_QUAD_SIDES; ++i)
}
newSides.swap( side );
- uv_grid.clear();
+ if ( keepGrid && !uv_grid.empty() )
+ {
+ if ( nb == 2 ) // "PI"
+ {
+ std::reverse( uv_grid.begin(), uv_grid.end() );
+ }
+ else
+ {
+ FaceQuadStruct newQuad;
+ newQuad.uv_grid.resize( uv_grid.size() );
+ newQuad.iSize = jSize;
+ newQuad.jSize = iSize;
+ int i, j, iRev, jRev;
+ int *iNew = ( nb == 1 ) ? &jRev : &j;
+ int *jNew = ( nb == 1 ) ? &i : &iRev;
+ for ( i = 0, iRev = iSize-1; i < iSize; ++i, --iRev )
+ for ( j = 0, jRev = jSize-1; j < jSize; ++j, --jRev )
+ newQuad.UVPt( *iNew, *jNew ) = UVPt( i, j );
+
+ std::swap( iSize, jSize );
+ std::swap( uv_grid, newQuad.uv_grid );
+ }
+ }
+ else
+ {
+ uv_grid.clear();
+ }
}
//=======================================================================
FaceQuadStruct ( const TopoDS_Face& F = TopoDS_Face(), const std::string& nm="main" );
UVPtStruct& UVPt( int i, int j ) { return uv_grid[ i + j * iSize ]; }
+ double& U( int i, int j ) { return UVPt( i, j ).u; }
+ double& V( int i, int j ) { return UVPt( i, j ).v; }
void shift ( size_t nb, bool keepUnitOri, bool keepGrid=false );
int & nbNodeOut( int iSide ) { return side[ iSide ].nbNodeOut; }
bool findCell ( const gp_XY& uv, int & i, int & j );
<source>ICON_SMESH_TREE_ALGO_RadialQuadrangle_1D2D</source>
<translation>mesh_tree_algo_radial_quadrangle_1D2D.png</translation>
</message>
+ <message>
+ <source>ICON_SMESH_TREE_ALGO_QuadFromMedialAxis_1D2D</source>
+ <translation>mesh_tree_algo_quad.png</translation>
+ </message>
<message>
<source>ICON_SMESH_TREE_ALGO_Prism_3D</source>
<translation>mesh_tree_algo_prism.png</translation>
#include "Utils_CorbaException.hxx"
#include "utilities.h"
+#include "StdMeshers_QuadFromMedialAxis_1D2D.hxx"
+
using namespace std;
//=============================================================================
return ::StdMeshers_Quadrangle_2D::IsApplicable( S, toCheckAll );
}
+//=============================================================================
+/*!
+ * StdMeshers_QuadFromMedialAxis_1D2D_i::StdMeshers_QuadFromMedialAxis_1D2D_i
+ *
+ * Constructor
+ */
+//=============================================================================
+
+StdMeshers_QuadFromMedialAxis_1D2D_i::
+StdMeshers_QuadFromMedialAxis_1D2D_i( PortableServer::POA_ptr thePOA,
+ int theStudyId,
+ ::SMESH_Gen* theGenImpl )
+ : SALOME::GenericObj_i( thePOA ),
+ SMESH_Hypothesis_i( thePOA ),
+ SMESH_Algo_i( thePOA ),
+ SMESH_2D_Algo_i( thePOA )
+{
+ MESSAGE( "StdMeshers_QuadFromMedialAxis_1D2D_i::StdMeshers_QuadFromMedialAxis_1D2D_i" );
+ myBaseImpl = new ::StdMeshers_QuadFromMedialAxis_1D2D( theGenImpl->GetANewId(),
+ theStudyId,
+ theGenImpl );
+}
+
+//=============================================================================
+/*!
+ * StdMeshers_QuadFromMedialAxis_1D2D_i::~StdMeshers_QuadFromMedialAxis_1D2D_i
+ *
+ * Destructor
+ *
+ */
+//=============================================================================
+
+StdMeshers_QuadFromMedialAxis_1D2D_i::~StdMeshers_QuadFromMedialAxis_1D2D_i()
+{
+ MESSAGE( "StdMeshers_QuadFromMedialAxis_1D2D_i::~StdMeshers_QuadFromMedialAxis_1D2D_i" );
+}
static CORBA::Boolean IsApplicable(const TopoDS_Shape &S, CORBA::Boolean toCheckAll);
};
+// ======================================================
+// Quadrangle (Medial Axis Projection) 2d algorithm
+// ======================================================
+class STDMESHERS_I_EXPORT StdMeshers_QuadFromMedialAxis_1D2D_i:
+ public virtual POA_StdMeshers::StdMeshers_QuadFromMedialAxis_1D2D,
+ public virtual SMESH_2D_Algo_i
+{
+ public:
+ // Constructor
+ StdMeshers_QuadFromMedialAxis_1D2D_i( PortableServer::POA_ptr thePOA,
+ int theStudyId,
+ ::SMESH_Gen* theGenImpl );
+
+ // Destructor
+ virtual ~StdMeshers_QuadFromMedialAxis_1D2D_i();
+
+ // Get implementation
+ //::StdMeshers_Quadrangle_2D* GetImpl();
+
+ // Return true if the algorithm is applicable to a shape
+ //static CORBA::Boolean IsApplicable(const TopoDS_Shape &S, CORBA::Boolean toCheckAll);
+};
+
#endif
#endif
else if (strcmp(aHypName, "Quadrangle_2D") == 0)
aCreator = new StdHypothesisCreator_i<StdMeshers_Quadrangle_2D_i, StdMeshers_Quadrangle_2D_i>;
+ else if (strcmp(aHypName, "QuadFromMedialAxis_1D2D") == 0)
+ aCreator = new StdHypothesisCreator_i<StdMeshers_QuadFromMedialAxis_1D2D_i, StdMeshers_Quadrangle_2D_i>;
else if (strcmp(aHypName, "Hexa_3D") == 0)
aCreator = new StdHypothesisCreator_i<StdMeshers_Hexa_3D_i, StdMeshers_Hexa_3D_i>;
else if (strcmp(aHypName, "Projection_1D") == 0)
