src/BLSURFPlugin/BLSURFPlugin_BLSURF.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
   2 //
   3 //  This library is free software; you can redistribute it and/or
   4 //  modify it under the terms of the GNU Lesser General Public
   5 //  License as published by the Free Software Foundation; either
   6 //  version 2.1 of the License.
   7 //
   8 //  This library is distributed in the hope that it will be useful,
   9 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 //  Lesser General Public License for more details.
  12 //
  13 //  You should have received a copy of the GNU Lesser General Public
  14 //  License along with this library; if not, write to the Free Software
  15 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // ---
  20 // File    : BLSURFPlugin_BLSURF.cxx
  21 // Authors : Francis KLOSS (OCC) & Patrick LAUG (INRIA) & Lioka RAZAFINDRAZAKA (CEA)
  22 //           & Aurelien ALLEAUME (DISTENE)
  23 //           Size maps developement: Nicolas GEIMER (OCC) & Gilles DAVID (EURIWARE)
  24 // ---
  25 //
  26 #include "BLSURFPlugin_BLSURF.hxx"
  27 #include "BLSURFPlugin_Hypothesis.hxx"
  28
  29 #include <structmember.h>
  30
  31
  32 #include <SMESH_Gen.hxx>
  33 #include <SMESH_Mesh.hxx>
  34 #include <SMESH_ControlsDef.hxx>
  35
  36 #include <utilities.h>
  37
  38 #include <limits>
  39 #include <list>
  40 #include <vector>
  41 #include <set>
  42 #include <cstdlib>
  43
  44 #include <BRep_Tool.hxx>
  45 #include <TopExp.hxx>
  46 #include <TopExp_Explorer.hxx>
  47 #include <TopoDS.hxx>
  48 #include <NCollection_Map.hxx>
  49 #include <Standard_ErrorHandler.hxx>
  50
  51 extern "C"{
  52 #include "distene/blsurf.h"
  53 #include <distene/api.h>
  54 }
  55
  56 #include <Geom_Surface.hxx>
  57 #include <Handle_Geom_Surface.hxx>
  58 #include <Geom2d_Curve.hxx>
  59 #include <Handle_Geom2d_Curve.hxx>
  60 #include <Geom_Curve.hxx>
  61 #include <Handle_Geom_Curve.hxx>
  62 #include <TopoDS_Vertex.hxx>
  63 #include <TopoDS_Edge.hxx>
  64 #include <TopoDS_Wire.hxx>
  65 #include <TopoDS_Face.hxx>
  66 #include <TopoDS_Shape.hxx>
  67 #include <gp_Pnt2d.hxx>
  68 #include <TopTools_IndexedMapOfShape.hxx>
  69 #include <BRepTools.hxx>
  70 #include <TopTools_DataMapOfShapeInteger.hxx>
  71 #include <GProp_GProps.hxx>
  72 #include <BRepGProp.hxx>
  73
  74 #ifndef WNT
  75 #include <fenv.h>
  76 #endif
  77
  78 #include <GeomAPI_ProjectPointOnCurve.hxx>
  79 #include <GeomAPI_ProjectPointOnSurf.hxx>
  80 #include <gp_XY.hxx>
  81 #include <gp_XYZ.hxx>
  82 // #include <BRepClass_FaceClassifier.hxx>
  83 #include <TopTools_MapOfShape.hxx>
  84
  85 /* ==================================
  86  * ===========  PYTHON ==============
  87  * ==================================*/
  88
  89 typedef struct {
  90   PyObject_HEAD
  91   int softspace;
  92   std::string *out;
  93   } PyStdOut;
  94
  95 static void
  96 PyStdOut_dealloc(PyStdOut *self)
  97 {
  98   PyObject_Del(self);
  99 }
 100
 101 static PyObject *
 102 PyStdOut_write(PyStdOut *self, PyObject *args)
 103 {
 104   char *c;
 105   int l;
 106   if (!PyArg_ParseTuple(args, "t#:write",&c, &l))
 107     return NULL;
 108
 109   //std::cerr << c ;
 110   *(self->out)=*(self->out)+c;
 111
 112   Py_INCREF(Py_None);
 113   return Py_None;
 114 }
 115
 116 static PyMethodDef PyStdOut_methods[] = {
 117   {"write",  (PyCFunction)PyStdOut_write,  METH_VARARGS,
 118     PyDoc_STR("write(string) -> None")},
 119   {NULL,    NULL}   /* sentinel */
 120 };
 121
 122 static PyMemberDef PyStdOut_memberlist[] = {
 123   {"softspace", T_INT,  offsetof(PyStdOut, softspace), 0,
 124    "flag indicating that a space needs to be printed; used by print"},
 125   {NULL} /* Sentinel */
 126 };
 127
 128 static PyTypeObject PyStdOut_Type = {
 129   /* The ob_type field must be initialized in the module init function
 130    * to be portable to Windows without using C++. */
 131   PyObject_HEAD_INIT(NULL)
 132   0,                            /*ob_size*/
 133   "PyOut",                      /*tp_name*/
 134   sizeof(PyStdOut),             /*tp_basicsize*/
 135   0,                            /*tp_itemsize*/
 136   /* methods */
 137   (destructor)PyStdOut_dealloc, /*tp_dealloc*/
 138   0,                            /*tp_print*/
 139   0,                            /*tp_getattr*/
 140   0,                            /*tp_setattr*/
 141   0,                            /*tp_compare*/
 142   0,                            /*tp_repr*/
 143   0,                            /*tp_as_number*/
 144   0,                            /*tp_as_sequence*/
 145   0,                            /*tp_as_mapping*/
 146   0,                            /*tp_hash*/
 147   0,                            /*tp_call*/
 148   0,                            /*tp_str*/
 149   PyObject_GenericGetAttr,      /*tp_getattro*/
 150   /* softspace is writable:  we must supply tp_setattro */
 151   PyObject_GenericSetAttr,      /* tp_setattro */
 152   0,                            /*tp_as_buffer*/
 153   Py_TPFLAGS_DEFAULT,           /*tp_flags*/
 154   0,                            /*tp_doc*/
 155   0,                            /*tp_traverse*/
 156   0,                            /*tp_clear*/
 157   0,                            /*tp_richcompare*/
 158   0,                            /*tp_weaklistoffset*/
 159   0,                            /*tp_iter*/
 160   0,                            /*tp_iternext*/
 161   PyStdOut_methods,             /*tp_methods*/
 162   PyStdOut_memberlist,          /*tp_members*/
 163   0,                            /*tp_getset*/
 164   0,                            /*tp_base*/
 165   0,                            /*tp_dict*/
 166   0,                            /*tp_descr_get*/
 167   0,                            /*tp_descr_set*/
 168   0,                            /*tp_dictoffset*/
 169   0,                            /*tp_init*/
 170   0,                            /*tp_alloc*/
 171   0,                            /*tp_new*/
 172   0,                            /*tp_free*/
 173   0,                            /*tp_is_gc*/
 174 };
 175
 176 PyObject * newPyStdOut( std::string& out )
 177 {
 178   PyStdOut *self;
 179   self = PyObject_New(PyStdOut, &PyStdOut_Type);
 180   if (self == NULL)
 181     return NULL;
 182   self->softspace = 0;
 183   self->out=&out;
 184   return (PyObject*)self;
 185 }
 186
 187
 188 ////////////////////////END PYTHON///////////////////////////
 189
 190 //////////////////MY MAPS////////////////////////////////////////
 191 TopTools_IndexedMapOfShape FacesWithSizeMap;
 192 std::map<int,string> FaceId2SizeMap;
 193 TopTools_IndexedMapOfShape EdgesWithSizeMap;
 194 std::map<int,string> EdgeId2SizeMap;
 195 TopTools_IndexedMapOfShape VertecesWithSizeMap;
 196 std::map<int,string> VertexId2SizeMap;
 197
 198 std::map<int,PyObject*> FaceId2PythonSmp;
 199 std::map<int,PyObject*> EdgeId2PythonSmp;
 200 std::map<int,PyObject*> VertexId2PythonSmp;
 201
 202 std::map<int,std::vector<double> > FaceId2AttractorCoords;
 203
 204 TopTools_IndexedMapOfShape FacesWithEnforcedVerteces;
 205 std::map< int, std::set< std::vector<double> > > FaceId2EnforcedVertexCoords;
 206
 207 bool HasSizeMapOnFace=false;
 208 bool HasSizeMapOnEdge=false;
 209 bool HasSizeMapOnVertex=false;
 210
 211
 212 //=============================================================================
 213 /*!
 214  *
 215  */
 216 //=============================================================================
 217
 218 BLSURFPlugin_BLSURF::BLSURFPlugin_BLSURF(int hypId, int studyId,
 219                                                SMESH_Gen* gen)
 220   : SMESH_2D_Algo(hypId, studyId, gen)
 221 {
 222   MESSAGE("BLSURFPlugin_BLSURF::BLSURFPlugin_BLSURF");
 223
 224   _name = "BLSURF";
 225   _shapeType = (1 << TopAbs_FACE); // 1 bit /shape type
 226   _compatibleHypothesis.push_back("BLSURF_Parameters");
 227   _requireDescretBoundary = false;
 228   _onlyUnaryInput = false;
 229   _hypothesis = NULL;
 230
 231   smeshGen_i = SMESH_Gen_i::GetSMESHGen();
 232   CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager");
 233   SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject);
 234
 235   MESSAGE("studyid = " << _studyId);
 236
 237   myStudy = NULL;
 238   myStudy = aStudyMgr->GetStudyByID(_studyId);
 239   MESSAGE("myStudy->StudyId() = " << myStudy->StudyId());
 240
 241   /* Initialize the Python interpreter */
 242   assert(Py_IsInitialized());
 243   PyGILState_STATE gstate;
 244   gstate = PyGILState_Ensure();
 245
 246   main_mod = NULL;
 247   main_mod = PyImport_AddModule("__main__");
 248
 249   main_dict = NULL;
 250   main_dict = PyModule_GetDict(main_mod);
 251
 252   PyRun_SimpleString("from math import *");
 253   PyGILState_Release(gstate);
 254
 255   FacesWithSizeMap.Clear();
 256   FaceId2SizeMap.clear();
 257   EdgesWithSizeMap.Clear();
 258   EdgeId2SizeMap.clear();
 259   VertecesWithSizeMap.Clear();
 260   VertexId2SizeMap.clear();
 261   FaceId2PythonSmp.clear();
 262   EdgeId2PythonSmp.clear();
 263   VertexId2PythonSmp.clear();
 264   FaceId2AttractorCoords.clear();
 265   FacesWithEnforcedVerteces.Clear();
 266   FaceId2EnforcedVertexCoords.clear();
 267
 268 }
 269
 270 //=============================================================================
 271 /*!
 272  *
 273  */
 274 //=============================================================================
 275
 276 BLSURFPlugin_BLSURF::~BLSURFPlugin_BLSURF()
 277 {
 278   MESSAGE("BLSURFPlugin_BLSURF::~BLSURFPlugin_BLSURF");
 279 }
 280
 281
 282 //=============================================================================
 283 /*!
 284  *
 285  */
 286 //=============================================================================
 287
 288 bool BLSURFPlugin_BLSURF::CheckHypothesis
 289                          (SMESH_Mesh&                          aMesh,
 290                           const TopoDS_Shape&                  aShape,
 291                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 292 {
 293   _hypothesis = NULL;
 294
 295   list<const SMESHDS_Hypothesis*>::const_iterator itl;
 296   const SMESHDS_Hypothesis* theHyp;
 297
 298   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
 299   int nbHyp = hyps.size();
 300   if (!nbHyp)
 301   {
 302     aStatus = SMESH_Hypothesis::HYP_OK;
 303     return true;  // can work with no hypothesis
 304   }
 305
 306   itl = hyps.begin();
 307   theHyp = (*itl); // use only the first hypothesis
 308
 309   string hypName = theHyp->GetName();
 310
 311   if (hypName == "BLSURF_Parameters")
 312   {
 313     _hypothesis = static_cast<const BLSURFPlugin_Hypothesis*> (theHyp);
 314     ASSERT(_hypothesis);
 315     if ( _hypothesis->GetPhysicalMesh() == BLSURFPlugin_Hypothesis::DefaultSize &&
 316          _hypothesis->GetGeometricMesh() == BLSURFPlugin_Hypothesis::DefaultGeom )
 317       //  hphy_flag = 0 and hgeo_flag = 0 is not allowed (spec)
 318       aStatus = SMESH_Hypothesis::HYP_BAD_PARAMETER;
 319     else
 320       aStatus = SMESH_Hypothesis::HYP_OK;
 321   }
 322   else
 323     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 324
 325   return aStatus == SMESH_Hypothesis::HYP_OK;
 326 }
 327
 328 //=============================================================================
 329 /*!
 330  * Pass parameters to BLSURF
 331  */
 332 //=============================================================================
 333
 334 inline std::string to_string(double d)
 335 {
 336    std::ostringstream o;
 337    o << d;
 338    return o.str();
 339 }
 340
 341 inline std::string to_string(int i)
 342 {
 343    std::ostringstream o;
 344    o << i;
 345    return o.str();
 346 }
 347
 348 double _smp_phy_size;
 349 status_t size_on_surface(integer face_id, real *uv, real *size, void *user_data);
 350 status_t size_on_edge(integer edge_id, real t, real *size, void *user_data);
 351 status_t size_on_vertex(integer vertex_id, real *size, void *user_data);
 352
 353 double my_u_min=1e6,my_v_min=1e6,my_u_max=-1e6,my_v_max=-1e6;
 354
 355 typedef struct {
 356         gp_XY uv;
 357         gp_XYZ xyz;
 358 } projectionPoint;
 359 /////////////////////////////////////////////////////////
 360 projectionPoint getProjectionPoint(const TopoDS_Face& face, const gp_XYZ& point)
 361 {
 362   projectionPoint myPoint;
 363   Handle(Geom_Surface) surface = BRep_Tool::Surface(face);
 364   GeomAPI_ProjectPointOnSurf projector( point, surface );
 365   if ( !projector.IsDone() || projector.NbPoints()==0 )
 366     throw "Can't project";
 367
 368   Quantity_Parameter u,v;
 369   projector.LowerDistanceParameters(u,v);
 370   myPoint.uv = gp_XY(u,v);
 371   gp_Pnt aPnt = projector.NearestPoint();
 372   myPoint.xyz = gp_XYZ(aPnt.X(),aPnt.Y(),aPnt.Z());
 373   //return gp_XY(u,v);
 374   return myPoint;
 375 }
 376 /////////////////////////////////////////////////////////
 377
 378 /////////////////////////////////////////////////////////
 379 double getT(const TopoDS_Edge& edge, const gp_XYZ& point)
 380 {
 381   Standard_Real f,l;
 382   Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, f,l);
 383   GeomAPI_ProjectPointOnCurve projector( point, curve);
 384   if ( projector.NbPoints() == 0 )
 385     throw;
 386   return projector.LowerDistanceParameter();
 387 }
 388
 389 /////////////////////////////////////////////////////////
 390 TopoDS_Shape BLSURFPlugin_BLSURF::entryToShape(std::string entry)
 391 {
 392   MESSAGE("BLSURFPlugin_BLSURF::entryToShape"<<entry );
 393   GEOM::GEOM_Object_var aGeomObj;
 394   TopoDS_Shape S = TopoDS_Shape();
 395   SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() );
 396   SALOMEDS::GenericAttribute_var anAttr;
 397
 398   if (!aSObj->_is_nil() && aSObj->FindAttribute(anAttr, "AttributeIOR")) {
 399     SALOMEDS::AttributeIOR_var anIOR = SALOMEDS::AttributeIOR::_narrow(anAttr);
 400     CORBA::String_var aVal = anIOR->Value();
 401     CORBA::Object_var obj = myStudy->ConvertIORToObject(aVal);
 402     aGeomObj = GEOM::GEOM_Object::_narrow(obj);
 403   }
 404   if ( !aGeomObj->_is_nil() )
 405     S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
 406   return S;
 407 }
 408
 409 /////////////////////////////////////////////////////////
 410 void createEnforcedVertexOnFace(TopoDS_Shape GeomShape, BLSURFPlugin_Hypothesis::TEnforcedVertexList enforcedVertexList)
 411 {
 412   double xe, ye, ze;
 413   std::vector<double> coords;
 414   BLSURFPlugin_Hypothesis::TEnforcedVertex enforcedVertex;
 415   // enforcedVertexList = set < vector<double> >
 416   BLSURFPlugin_Hypothesis::TEnforcedVertexList::const_iterator evlIt = enforcedVertexList.begin();
 417
 418   for( ; evlIt != enforcedVertexList.end() ; ++evlIt ) {
 419     coords.clear();
 420     enforcedVertex = *evlIt;
 421     xe = enforcedVertex[0];
 422     ye = enforcedVertex[1];
 423     ze = enforcedVertex[2];
 424     MESSAGE("Enforced Vertex: " << xe << ", " << ye << ", " << ze);
 425     // Get the (u,v) values of the enforced vertex on the face
 426     projectionPoint myPoint = getProjectionPoint(TopoDS::Face(GeomShape),gp_XYZ(xe,ye,ze));
 427     gp_XY uvPoint = myPoint.uv;
 428     gp_XYZ xyzPoint = myPoint.xyz;
 429     Standard_Real u0 = uvPoint.X();
 430     Standard_Real v0 = uvPoint.Y();
 431     Standard_Real x0 = xyzPoint.X();
 432     Standard_Real y0 = xyzPoint.Y();
 433     Standard_Real z0 = xyzPoint.Z();
 434     MESSAGE("Projected Vertex: " << x0 << ", " << y0 << ", " << z0);
 435     coords.push_back(u0);
 436     coords.push_back(v0);
 437     coords.push_back(x0);
 438     coords.push_back(y0);
 439     coords.push_back(z0);
 440
 441     int key = 0;
 442     if (! FacesWithEnforcedVerteces.Contains(TopoDS::Face(GeomShape))) {
 443       key = FacesWithEnforcedVerteces.Add(TopoDS::Face(GeomShape));
 444     }
 445     else {
 446       key = FacesWithEnforcedVerteces.FindIndex(TopoDS::Face(GeomShape));
 447     }
 448
 449     if (FaceId2EnforcedVertexCoords.find(key) != FaceId2EnforcedVertexCoords.end()) {
 450       MESSAGE("Map of enf. vertex has key " << key)
 451       MESSAGE("Enf. vertex list size is: " << FaceId2EnforcedVertexCoords[key].size())
 452       FaceId2EnforcedVertexCoords[key].insert(coords);
 453       MESSAGE("New Enf. vertex list size is: " << FaceId2EnforcedVertexCoords[key].size())
 454     }
 455     else {
 456       MESSAGE("Map of enf. vertex has not key " << key << ": creating it")
 457       std::set< std::vector<double> > ens;
 458       ens.insert(coords);
 459       FaceId2EnforcedVertexCoords[key] = ens;
 460     }
 461   }
 462 }
 463
 464 /////////////////////////////////////////////////////////
 465 void createAttractorOnFace(TopoDS_Shape GeomShape, std::string AttractorFunction)
 466 {
 467   MESSAGE("Attractor function: "<< AttractorFunction);
 468   double xa, ya, za; // Coordinates of attractor point
 469   double a, b;       // Attractor parameter
 470   bool createNode=false; // To create a node on attractor projection
 471   int pos1, pos2;
 472   // atIt->second has the following pattern:
 473   // ATTRACTOR(xa;ya;za;a;b)
 474   // where:
 475   // xa;ya;za : coordinates of  attractor
 476   // a        : desired size on attractor
 477   // b        : distance of influence of attractor
 478   //
 479   // We search the parameters in the string
 480   // xa
 481   pos1 = AttractorFunction.find(";");
 482   if (pos1!=string::npos)
 483   xa = atof(AttractorFunction.substr(10, pos1-10).c_str());
 484   // ya
 485   pos2 = AttractorFunction.find(";", pos1+1);
 486   if (pos2!=string::npos) {
 487   ya = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
 488   pos1 = pos2;
 489     }
 490   // za
 491   pos2 = AttractorFunction.find(";", pos1+1);
 492   if (pos2!=string::npos) {
 493   za = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
 494   pos1 = pos2;
 495   }
 496   // a
 497   pos2 = AttractorFunction.find(";", pos1+1);
 498   if (pos2!=string::npos) {
 499   a = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
 500   pos1 = pos2;
 501   }
 502   // b
 503   pos2 = AttractorFunction.find(";", pos1+1);
 504   if (pos2!=string::npos) {
 505   b = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
 506     pos1 = pos2;
 507   }
 508   // createNode
 509   pos2 = AttractorFunction.find(")");
 510   if (pos2!=string::npos) {
 511     string createNodeStr = AttractorFunction.substr(pos1+1, pos2-pos1-1);
 512     MESSAGE("createNode: " << createNodeStr);
 513     createNode = (AttractorFunction.substr(pos1+1, pos2-pos1-1) == "True");
 514   }
 515
 516   // Get the (u,v) values of the attractor on the face
 517   projectionPoint myPoint = getProjectionPoint(TopoDS::Face(GeomShape),gp_XYZ(xa,ya,za));
 518   gp_XY uvPoint = myPoint.uv;
 519   gp_XYZ xyzPoint = myPoint.xyz;
 520   Standard_Real u0 = uvPoint.X();
 521   Standard_Real v0 = uvPoint.Y();
 522   Standard_Real x0 = xyzPoint.X();
 523   Standard_Real y0 = xyzPoint.Y();
 524   Standard_Real z0 = xyzPoint.Z();
 525   std::vector<double> coords;
 526   coords.push_back(u0);
 527   coords.push_back(v0);
 528   coords.push_back(x0);
 529   coords.push_back(y0);
 530   coords.push_back(z0);
 531   // We construct the python function
 532   ostringstream attractorFunctionStream;
 533   attractorFunctionStream << "def f(u,v): return ";
 534   attractorFunctionStream << _smp_phy_size << "-(" << _smp_phy_size <<"-" << a << ")";
 535   attractorFunctionStream << "*exp(-((u-("<<u0<<"))*(u-("<<u0<<"))+(v-("<<v0<<"))*(v-("<<v0<<")))/(" << b << "*" << b <<"))";
 536
 537   MESSAGE("Python function for attractor:" << std::endl << attractorFunctionStream.str());
 538
 539   int key;
 540   if (! FacesWithSizeMap.Contains(TopoDS::Face(GeomShape))) {
 541     key = FacesWithSizeMap.Add(TopoDS::Face(GeomShape));
 542   }
 543   else {
 544     key = FacesWithSizeMap.FindIndex(TopoDS::Face(GeomShape));
 545   }
 546   FaceId2SizeMap[key] =attractorFunctionStream.str();
 547   if (createNode) {
 548     MESSAGE("Creating node on ("<<x0<<","<<y0<<","<<z0<<")");
 549     FaceId2AttractorCoords[key] = coords;
 550   }
 551 }
 552
 553 /////////////////////////////////////////////////////////
 554
 555 void BLSURFPlugin_BLSURF::SetParameters(const BLSURFPlugin_Hypothesis* hyp, blsurf_session_t *bls)
 556 {
 557   int    _topology      = BLSURFPlugin_Hypothesis::GetDefaultTopology();
 558   int    _physicalMesh  = BLSURFPlugin_Hypothesis::GetDefaultPhysicalMesh();
 559   double _phySize       = BLSURFPlugin_Hypothesis::GetDefaultPhySize();
 560   int    _geometricMesh = BLSURFPlugin_Hypothesis::GetDefaultGeometricMesh();
 561   double _angleMeshS    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshS();
 562   double _angleMeshC    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshC();
 563   double _gradation     = BLSURFPlugin_Hypothesis::GetDefaultGradation();
 564   bool   _quadAllowed   = BLSURFPlugin_Hypothesis::GetDefaultQuadAllowed();
 565   bool   _decimesh      = BLSURFPlugin_Hypothesis::GetDefaultDecimesh();
 566   int    _verb          = BLSURFPlugin_Hypothesis::GetDefaultVerbosity();
 567
 568   if (hyp) {
 569     MESSAGE("BLSURFPlugin_BLSURF::SetParameters");
 570     _topology      = (int) hyp->GetTopology();
 571     _physicalMesh  = (int) hyp->GetPhysicalMesh();
 572     _phySize       = hyp->GetPhySize();
 573     _geometricMesh = (int) hyp->GetGeometricMesh();
 574     _angleMeshS    = hyp->GetAngleMeshS();
 575     _angleMeshC    = hyp->GetAngleMeshC();
 576     _gradation     = hyp->GetGradation();
 577     _quadAllowed   = hyp->GetQuadAllowed();
 578     _decimesh      = hyp->GetDecimesh();
 579     _verb          = hyp->GetVerbosity();
 580
 581     if ( hyp->GetPhyMin() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 582       blsurf_set_param(bls, "hphymin", to_string(hyp->GetPhyMin()).c_str());
 583     if ( hyp->GetPhyMax() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 584       blsurf_set_param(bls, "hphymax", to_string(hyp->GetPhyMax()).c_str());
 585     if ( hyp->GetGeoMin() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 586       blsurf_set_param(bls, "hgeomin", to_string(hyp->GetGeoMin()).c_str());
 587     if ( hyp->GetGeoMax() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 588       blsurf_set_param(bls, "hgeomax", to_string(hyp->GetGeoMax()).c_str());
 589
 590     const BLSURFPlugin_Hypothesis::TOptionValues & opts = hyp->GetOptionValues();
 591     BLSURFPlugin_Hypothesis::TOptionValues::const_iterator opIt;
 592     for ( opIt = opts.begin(); opIt != opts.end(); ++opIt )
 593       if ( !opIt->second.empty() ) {
 594         MESSAGE("blsurf_set_param(): " << opIt->first << " = " << opIt->second);
 595         blsurf_set_param(bls, opIt->first.c_str(), opIt->second.c_str());
 596       }
 597
 598   } else {
 599     MESSAGE("BLSURFPlugin_BLSURF::SetParameters using defaults");
 600   }
 601   _smp_phy_size = _phySize;
 602   blsurf_set_param(bls, "topo_points",       _topology > 0 ? "1" : "0");
 603   blsurf_set_param(bls, "topo_curves",       _topology > 0 ? "1" : "0");
 604   blsurf_set_param(bls, "topo_project",      _topology > 0 ? "1" : "0");
 605   blsurf_set_param(bls, "clean_boundary",    _topology > 1 ? "1" : "0");
 606   blsurf_set_param(bls, "close_boundary",    _topology > 1 ? "1" : "0");
 607   blsurf_set_param(bls, "hphy_flag",         to_string(_physicalMesh).c_str());
 608 //  blsurf_set_param(bls, "hphy_flag",         "2");
 609   if ((to_string(_physicalMesh))=="2"){
 610
 611     TopoDS_Shape GeomShape;
 612     TopAbs_ShapeEnum GeomType;
 613     //
 614     // Standard Size Maps
 615     //
 616     MESSAGE("Setting a Size Map");
 617     const BLSURFPlugin_Hypothesis::TSizeMap sizeMaps = BLSURFPlugin_Hypothesis::GetSizeMapEntries(hyp);
 618     BLSURFPlugin_Hypothesis::TSizeMap::const_iterator smIt = sizeMaps.begin();
 619     for ( ; smIt != sizeMaps.end(); ++smIt ) {
 620       if ( !smIt->second.empty() ) {
 621         MESSAGE("blsurf_set_sizeMap(): " << smIt->first << " = " << smIt->second);
 622         GeomShape = entryToShape(smIt->first);
 623         GeomType  = GeomShape.ShapeType();
 624         MESSAGE("Geomtype is " << GeomType);
 625         int key = -1;
 626         // Group Management
 627         if (GeomType == TopAbs_COMPOUND){
 628           for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
 629             // Group of faces
 630             if (it.Value().ShapeType() == TopAbs_FACE){
 631               HasSizeMapOnFace = true;
 632               if (! FacesWithSizeMap.Contains(TopoDS::Face(it.Value()))) {
 633                 key = FacesWithSizeMap.Add(TopoDS::Face(it.Value()));
 634               }
 635               else {
 636                 key = FacesWithSizeMap.FindIndex(TopoDS::Face(it.Value()));
 637 //                 MESSAGE("Face with key " << key << " already in map");
 638               }
 639               FaceId2SizeMap[key] = smIt->second;
 640             }
 641             // Group of edges
 642             if (it.Value().ShapeType() == TopAbs_EDGE){
 643               HasSizeMapOnEdge = true;
 644               HasSizeMapOnFace = true;
 645               if (! EdgesWithSizeMap.Contains(TopoDS::Edge(it.Value()))) {
 646                 key = EdgesWithSizeMap.Add(TopoDS::Edge(it.Value()));
 647               }
 648               else {
 649                 key = EdgesWithSizeMap.FindIndex(TopoDS::Edge(it.Value()));
 650 //                 MESSAGE("Edge with key " << key << " already in map");
 651               }
 652               EdgeId2SizeMap[key] = smIt->second;
 653             }
 654             // Group of verteces
 655             if (it.Value().ShapeType() == TopAbs_VERTEX){
 656               HasSizeMapOnVertex = true;
 657               HasSizeMapOnEdge = true;
 658               HasSizeMapOnFace = true;
 659               if (! VertecesWithSizeMap.Contains(TopoDS::Vertex(it.Value()))) {
 660                 key = VertecesWithSizeMap.Add(TopoDS::Vertex(it.Value()));
 661               }
 662               else {
 663                 key = VertecesWithSizeMap.FindIndex(TopoDS::Vertex(it.Value()));
 664 //                 MESSAGE("Vertex with key " << key << " already in map");
 665               }
 666               VertexId2SizeMap[key] = smIt->second;
 667             }
 668           }
 669         }
 670         // Single face
 671         if (GeomType == TopAbs_FACE){
 672           HasSizeMapOnFace = true;
 673           if (! FacesWithSizeMap.Contains(TopoDS::Face(GeomShape))) {
 674             key = FacesWithSizeMap.Add(TopoDS::Face(GeomShape));
 675           }
 676           else {
 677             key = FacesWithSizeMap.FindIndex(TopoDS::Face(GeomShape));
 678 //             MESSAGE("Face with key " << key << " already in map");
 679           }
 680           FaceId2SizeMap[key] = smIt->second;
 681         }
 682         // Single edge
 683         if (GeomType == TopAbs_EDGE){
 684           HasSizeMapOnEdge = true;
 685           HasSizeMapOnFace = true;
 686           if (! EdgesWithSizeMap.Contains(TopoDS::Edge(GeomShape))) {
 687             key = EdgesWithSizeMap.Add(TopoDS::Edge(GeomShape));
 688           }
 689           else {
 690             key = EdgesWithSizeMap.FindIndex(TopoDS::Edge(GeomShape));
 691 //             MESSAGE("Edge with key " << key << " already in map");
 692           }
 693           EdgeId2SizeMap[key] = smIt->second;
 694         }
 695         // Single vertex
 696         if (GeomType == TopAbs_VERTEX){
 697           HasSizeMapOnVertex = true;
 698           HasSizeMapOnEdge   = true;
 699           HasSizeMapOnFace   = true;
 700           if (! VertecesWithSizeMap.Contains(TopoDS::Vertex(GeomShape))) {
 701             key = VertecesWithSizeMap.Add(TopoDS::Vertex(GeomShape));
 702           }
 703           else {
 704             key = VertecesWithSizeMap.FindIndex(TopoDS::Vertex(GeomShape));
 705 //             MESSAGE("Vertex with key " << key << " already in map");
 706           }
 707           VertexId2SizeMap[key] = smIt->second;
 708         }
 709       }
 710     }
 711
 712     //
 713     // Attractors
 714     //
 715     MESSAGE("Setting Attractors");
 716     const BLSURFPlugin_Hypothesis::TSizeMap attractors = BLSURFPlugin_Hypothesis::GetAttractorEntries(hyp);
 717     BLSURFPlugin_Hypothesis::TSizeMap::const_iterator atIt = attractors.begin();
 718     for ( ; atIt != attractors.end(); ++atIt ) {
 719       if ( !atIt->second.empty() ) {
 720         MESSAGE("blsurf_set_attractor(): " << atIt->first << " = " << atIt->second);
 721         GeomShape = entryToShape(atIt->first);
 722         GeomType  = GeomShape.ShapeType();
 723         // Group Management
 724         if (GeomType == TopAbs_COMPOUND){
 725           for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
 726             if (it.Value().ShapeType() == TopAbs_FACE){
 727               HasSizeMapOnFace = true;
 728               createAttractorOnFace(it.Value(), atIt->second);
 729             }
 730           }
 731         }
 732
 733         if (GeomType == TopAbs_FACE){
 734           HasSizeMapOnFace = true;
 735           createAttractorOnFace(GeomShape, atIt->second);
 736         }
 737 /*
 738         if (GeomType == TopAbs_EDGE){
 739           HasSizeMapOnEdge = true;
 740           HasSizeMapOnFace = true;
 741         EdgeId2SizeMap[TopoDS::Edge(GeomShape).HashCode(IntegerLast())] = atIt->second;
 742         }
 743         if (GeomType == TopAbs_VERTEX){
 744           HasSizeMapOnVertex = true;
 745           HasSizeMapOnEdge   = true;
 746           HasSizeMapOnFace   = true;
 747         VertexId2SizeMap[TopoDS::Vertex(GeomShape).HashCode(IntegerLast())] = atIt->second;
 748         }
 749 */
 750       }
 751     }
 752
 753
 754     //
 755     // Enforced Verteces
 756     //
 757     MESSAGE("Setting Enforced Verteces");
 758     const BLSURFPlugin_Hypothesis::TEnforcedVertexMap enforcedVertexMap = BLSURFPlugin_Hypothesis::GetAllEnforcedVerteces(hyp);
 759     BLSURFPlugin_Hypothesis::TEnforcedVertexMap::const_iterator enfIt = enforcedVertexMap.begin();
 760     for ( ; enfIt != enforcedVertexMap.end(); ++enfIt ) {
 761       if ( !enfIt->second.empty() ) {
 762         GeomShape = entryToShape(enfIt->first);
 763         GeomType  = GeomShape.ShapeType();
 764         // Group Management
 765         if (GeomType == TopAbs_COMPOUND){
 766           for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
 767             if (it.Value().ShapeType() == TopAbs_FACE){
 768               HasSizeMapOnFace = true;
 769               createEnforcedVertexOnFace(it.Value(), enfIt->second);
 770             }
 771           }
 772         }
 773
 774         if (GeomType == TopAbs_FACE){
 775           HasSizeMapOnFace = true;
 776           createEnforcedVertexOnFace(GeomShape, enfIt->second);
 777         }
 778       }
 779     }
 780
 781 //    if (HasSizeMapOnFace){
 782     // In all size map cases (hphy_flag = 2), at least map on face must be defined
 783     MESSAGE("Setting Size Map on FACES ");
 784     blsurf_data_set_sizemap_iso_cad_face(bls, size_on_surface, &_smp_phy_size);
 785 //    }
 786
 787     if (HasSizeMapOnEdge){
 788       MESSAGE("Setting Size Map on EDGES ");
 789       blsurf_data_set_sizemap_iso_cad_edge(bls, size_on_edge, &_smp_phy_size);
 790     }
 791     if (HasSizeMapOnVertex){
 792       MESSAGE("Setting Size Map on VERTICES ");
 793       blsurf_data_set_sizemap_iso_cad_point(bls, size_on_vertex, &_smp_phy_size);
 794     }
 795   }
 796   blsurf_set_param(bls, "hphydef",           to_string(_phySize).c_str());
 797   blsurf_set_param(bls, "hgeo_flag",         to_string(_geometricMesh).c_str());
 798   blsurf_set_param(bls, "relax_size",        _decimesh ? "0": to_string(_geometricMesh).c_str());
 799   blsurf_set_param(bls, "angle_meshs",       to_string(_angleMeshS).c_str());
 800   blsurf_set_param(bls, "angle_meshc",       to_string(_angleMeshC).c_str());
 801   blsurf_set_param(bls, "gradation",         to_string(_gradation).c_str());
 802   blsurf_set_param(bls, "patch_independent", _decimesh ? "1" : "0");
 803   blsurf_set_param(bls, "element",           _quadAllowed ? "q1.0" : "p1");
 804   blsurf_set_param(bls, "verb",              to_string(_verb).c_str());
 805 }
 806
 807
 808
 809
 810
 811
 812 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data);
 813 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
 814                   real *duu, real *duv, real *dvv, void *user_data);
 815 status_t message_callback(message_t *msg, void *user_data);
 816
 817 //=============================================================================
 818 /*!
 819  *
 820  */
 821 //=============================================================================
 822
 823 bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) {
 824
 825   MESSAGE("BLSURFPlugin_BLSURF::Compute");
 826
 827   if (aShape.ShapeType() == TopAbs_COMPOUND) {
 828     MESSAGE("  the shape is a COMPOUND");
 829   }
 830   else {
 831     MESSAGE("  the shape is UNKNOWN");
 832   };
 833
 834   context_t *ctx =  context_new();
 835   context_set_message_callback(ctx, message_callback, &_comment);
 836
 837   cad_t *c = cad_new(ctx);
 838
 839   blsurf_session_t *bls = blsurf_session_new(ctx);
 840
 841   FacesWithSizeMap.Clear();
 842   FaceId2SizeMap.clear();
 843   EdgesWithSizeMap.Clear();
 844   EdgeId2SizeMap.clear();
 845   VertecesWithSizeMap.Clear();
 846   VertexId2SizeMap.clear();
 847
 848   MESSAGE("BEGIN SetParameters");
 849   SetParameters(_hypothesis, bls);
 850   MESSAGE("END SetParameters");
 851
 852   TopTools_IndexedMapOfShape fmap;
 853   TopTools_IndexedMapOfShape emap;
 854   TopTools_IndexedMapOfShape pmap;
 855   vector<Handle(Geom2d_Curve)> curves;
 856   vector<Handle(Geom_Surface)> surfaces;
 857
 858   fmap.Clear();
 859   FaceId2PythonSmp.clear();
 860   emap.Clear();
 861   EdgeId2PythonSmp.clear();
 862   pmap.Clear();
 863   VertexId2PythonSmp.clear();
 864   surfaces.resize(0);
 865   curves.resize(0);
 866
 867   assert(Py_IsInitialized());
 868   PyGILState_STATE gstate;
 869   gstate = PyGILState_Ensure();
 870
 871   string theSizeMapStr;
 872
 873   /****************************************************************************************
 874                                   FACES
 875   *****************************************************************************************/
 876   int iface = 0;
 877   string bad_end = "return";
 878   int faceKey = -1;
 879   for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next()) {
 880     TopoDS_Face f=TopoDS::Face(face_iter.Current());
 881
 882     if (fmap.FindIndex(f) > 0)
 883       continue;
 884
 885     fmap.Add(f);
 886     iface++;
 887     surfaces.push_back(BRep_Tool::Surface(f));
 888
 889     cad_face_t *fce = cad_face_new(c, iface, surf_fun, surfaces.back());
 890     cad_face_set_tag(fce, iface);
 891     if(f.Orientation() != TopAbs_FORWARD){
 892       cad_face_set_orientation(fce, CAD_ORIENTATION_REVERSED);
 893     } else {
 894       cad_face_set_orientation(fce, CAD_ORIENTATION_FORWARD);
 895     }
 896
 897     if (HasSizeMapOnFace){
 898       std::cout << "A size map is defined on a face" << std::endl;
 899       // Classic size map
 900       faceKey = FacesWithSizeMap.FindIndex(f);
 901
 902       if (FaceId2SizeMap.find(faceKey)!=FaceId2SizeMap.end()){
 903         theSizeMapStr = FaceId2SizeMap[faceKey];
 904         // check if function ends with "return"
 905         if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
 906           continue;
 907         // Expr To Python function, verification is performed at validation in GUI
 908         PyObject * obj = NULL;
 909         obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
 910         Py_DECREF(obj);
 911         PyObject * func = NULL;
 912         func = PyObject_GetAttrString(main_mod, "f");
 913         FaceId2PythonSmp[iface]=func;
 914         FaceId2SizeMap.erase(faceKey);
 915       }
 916
 917       // Specific size map = Attractor
 918       std::map<int,std::vector<double> >::iterator attractor_iter = FaceId2AttractorCoords.begin();
 919       int iatt=0;
 920       for (; attractor_iter != FaceId2AttractorCoords.end(); ++attractor_iter) {
 921         if (attractor_iter->first == faceKey) {
 922           MESSAGE("Face indice: " << iface);
 923           MESSAGE("Adding attractor");
 924
 925           double xyzCoords[3]  = {attractor_iter->second[2],
 926                                   attractor_iter->second[3],
 927                                   attractor_iter->second[4]};
 928
 929           MESSAGE("Check position of vertex =(" << xyzCoords[0] << "," << xyzCoords[1] << "," << xyzCoords[2] << ")");
 930           gp_Pnt P(xyzCoords[0],xyzCoords[1],xyzCoords[2]);
 931           BRepClass_FaceClassifier scl(f,P,1e-7);
 932           // scl.Perform() is bugged. The function was rewritten
 933 //          scl.Perform();
 934           BRepClass_FaceClassifierPerform(&scl,f,P,1e-7);
 935           TopAbs_State result = scl.State();
 936           MESSAGE("Position of point on face: "<<result);
 937           if ( result == TopAbs_OUT )
 938               MESSAGE("Point is out of face: node is not created");
 939           if ( result == TopAbs_UNKNOWN )
 940               MESSAGE("Point position on face is unknown: node is not created");
 941           if ( result == TopAbs_ON )
 942               MESSAGE("Point is on border of face: node is not created");
 943           if ( result == TopAbs_IN )
 944           {
 945             // Point is inside face and not on border
 946             MESSAGE("Point is in face: node is created");
 947             double uvCoords[2]   = {attractor_iter->second[0],attractor_iter->second[1]};
 948             iatt++;
 949             MESSAGE("Add cad point on (u,v)=(" << uvCoords[0] << "," << uvCoords[1] << ") with id = " << iatt);
 950             cad_point_t* point_p = cad_point_new(fce, iatt, uvCoords);
 951             cad_point_set_tag(point_p, iatt);
 952           }
 953           FaceId2AttractorCoords.erase(faceKey);
 954         }
 955       }
 956
 957       // Enforced Verteces
 958       faceKey = FacesWithEnforcedVerteces.FindIndex(f);
 959       std::map<int,std::set<std::vector<double> > >::const_iterator evmIt = FaceId2EnforcedVertexCoords.find(faceKey);
 960       if (evmIt != FaceId2EnforcedVertexCoords.end()) {
 961         std::cout << "Some enforced verteces are defined" << std::endl;
 962         int ienf = 0;
 963         std::set<std::vector<double> > evl;
 964 //         std::vector<double> ev;
 965         MESSAGE("Face indice: " << iface);
 966         MESSAGE("Adding enforced verteces");
 967         evl = evmIt->second;
 968         MESSAGE("Number of verteces to add: "<< evl.size())
 969         std::set< std::vector<double> >::const_iterator evlIt = evl.begin();
 970         for (; evlIt != evl.end(); ++evlIt) {
 971 //           ev = *evlIt;
 972 //         for (int i=0; i<evl.size() ; i++) {
 973 //           ev = evl[i];
 974
 975 //           double xyzCoords[3]  = {ev[2], ev[3], ev[4]};
 976           double xyzCoords[3]  = {evlIt->at(0), evlIt->at(3), evlIt->at(4)};
 977           MESSAGE("Check position of vertex =(" << xyzCoords[0] << "," << xyzCoords[1] << "," << xyzCoords[2] << ")");
 978           gp_Pnt P(xyzCoords[0],xyzCoords[1],xyzCoords[2]);
 979           BRepClass_FaceClassifier scl(f,P,1e-7);
 980           // scl.Perform() is bugged. The function was rewritten
 981 //          scl.Perform();
 982           BRepClass_FaceClassifierPerform(&scl,f,P,1e-7);
 983           TopAbs_State result = scl.State();
 984           MESSAGE("Position of point on face: "<<result);
 985           if ( result == TopAbs_OUT )
 986               MESSAGE("Point is out of face: node is not created");
 987           if ( result == TopAbs_UNKNOWN )
 988               MESSAGE("Point position on face is unknown: node is not created");
 989           if ( result == TopAbs_ON )
 990               MESSAGE("Point is on border of face: node is not created");
 991           if ( result == TopAbs_IN )
 992           {
 993             // Point is inside face and not on border
 994             MESSAGE("Point is in face: node is created");
 995 //             double uvCoords[2]   = {ev[0],ev[1]};
 996             double uvCoords[2]   = {evlIt->at(0),evlIt->at(1)};
 997             ienf++;
 998             MESSAGE("Add cad point on (u,v)=(" << uvCoords[0] << "," << uvCoords[1] << ") with id = " << ienf);
 999             cad_point_t* point_p = cad_point_new(fce, ienf, uvCoords);
1000             cad_point_set_tag(point_p, ienf);
1001           }
1002         }
1003         FaceId2EnforcedVertexCoords.erase(faceKey);
1004       }
1005       else
1006         std::cout << "No enforced vertex defined" << std::endl;
1007     }
1008
1009
1010     /****************************************************************************************
1011                                     EDGES
1012     *****************************************************************************************/
1013     int edgeKey = -1;
1014     for (TopExp_Explorer edge_iter(f,TopAbs_EDGE);edge_iter.More();edge_iter.Next()) {
1015       TopoDS_Edge e = TopoDS::Edge(edge_iter.Current());
1016       int ic = emap.FindIndex(e);
1017       if (ic <= 0)
1018         ic = emap.Add(e);
1019
1020       double tmin,tmax;
1021       curves.push_back(BRep_Tool::CurveOnSurface(e, f, tmin, tmax));
1022
1023       if (HasSizeMapOnEdge){
1024         edgeKey = EdgesWithSizeMap.FindIndex(e);
1025         if (EdgeId2SizeMap.find(edgeKey)!=EdgeId2SizeMap.end()) {
1026           theSizeMapStr = EdgeId2SizeMap[faceKey];
1027           if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
1028             continue;
1029           // Expr To Python function, verification is performed at validation in GUI
1030           PyObject * obj = NULL;
1031           obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
1032           Py_DECREF(obj);
1033           PyObject * func = NULL;
1034           func = PyObject_GetAttrString(main_mod, "f");
1035           EdgeId2PythonSmp[ic]=func;
1036           EdgeId2SizeMap.erase(edgeKey);
1037         }
1038       }
1039       cad_edge_t *edg = cad_edge_new(fce, ic, tmin, tmax, curv_fun, curves.back());
1040       cad_edge_set_tag(edg, ic);
1041       cad_edge_set_property(edg, EDGE_PROPERTY_SOFT_REQUIRED);
1042       if (e.Orientation() == TopAbs_INTERNAL)
1043         cad_edge_set_property(edg, EDGE_PROPERTY_INTERNAL);
1044
1045       int npts = 0;
1046       int ip1, ip2, *ip;
1047       gp_Pnt2d e0 = curves.back()->Value(tmin);
1048       gp_Pnt ee0 = surfaces.back()->Value(e0.X(), e0.Y());
1049       Standard_Real d1=0,d2=0;
1050
1051
1052       /****************************************************************************************
1053                                       VERTECES
1054       *****************************************************************************************/
1055       int vertexKey = -1;
1056       for (TopExp_Explorer ex_edge(e ,TopAbs_VERTEX); ex_edge.More(); ex_edge.Next()) {
1057         TopoDS_Vertex v = TopoDS::Vertex(ex_edge.Current());
1058         ++npts;
1059         if (npts == 1){
1060           ip = &ip1;
1061           d1 = ee0.SquareDistance(BRep_Tool::Pnt(v));
1062         } else {
1063           ip = &ip2;
1064           d2 = ee0.SquareDistance(BRep_Tool::Pnt(v));
1065         }
1066         *ip = pmap.FindIndex(v);
1067         if(*ip <= 0)
1068           *ip = pmap.Add(v);
1069
1070         vertexKey = VertecesWithSizeMap.FindIndex(v);
1071         if (HasSizeMapOnVertex){
1072           vertexKey = VertecesWithSizeMap.FindIndex(v);
1073           if (VertexId2SizeMap.find(vertexKey)!=VertexId2SizeMap.end()){
1074             theSizeMapStr = VertexId2SizeMap[faceKey];
1075             if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
1076               continue;
1077             // Expr To Python function, verification is performed at validation in GUI
1078             PyObject * obj = NULL;
1079             obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
1080             Py_DECREF(obj);
1081             PyObject * func = NULL;
1082             func = PyObject_GetAttrString(main_mod, "f");
1083             VertexId2PythonSmp[*ip]=func;
1084 //             VertexId2SizeMap.erase(vertexKey);   // do not erase if using a vector
1085           }
1086         }
1087       }
1088       if (npts != 2) {
1089         // should not happen
1090         MESSAGE("An edge does not have 2 extremities.");
1091       } else {
1092         if (d1 < d2)
1093           cad_edge_set_extremities(edg, ip1, ip2);
1094         else
1095           cad_edge_set_extremities(edg, ip2, ip1);
1096       }
1097     } // for edge
1098   } //for face
1099
1100
1101   PyGILState_Release(gstate);
1102
1103   blsurf_data_set_cad(bls, c);
1104
1105   std::cout << std::endl;
1106   std::cout << "Beginning of Surface Mesh generation" << std::endl;
1107   std::cout << std::endl;
1108
1109   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
1110 #ifndef WNT
1111   feclearexcept( FE_ALL_EXCEPT );
1112   int oldFEFlags = fedisableexcept( FE_ALL_EXCEPT );
1113 #endif
1114
1115   status_t status = STATUS_ERROR;
1116
1117   try {
1118     OCC_CATCH_SIGNALS;
1119
1120     status = blsurf_compute_mesh(bls);
1121
1122   }
1123   catch ( std::exception& exc ) {
1124     _comment += exc.what();
1125   }
1126   catch (Standard_Failure& ex) {
1127     _comment += ex.DynamicType()->Name();
1128     if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) {
1129       _comment += ": ";
1130       _comment += ex.GetMessageString();
1131     }
1132   }
1133   catch (...) {
1134     if ( _comment.empty() )
1135       _comment = "Exception in blsurf_compute_mesh()";
1136   }
1137   if ( status != STATUS_OK) {
1138     blsurf_session_delete(bls);
1139     cad_delete(c);
1140     context_delete(ctx);
1141
1142     return error(_comment);
1143   }
1144
1145   std::cout << std::endl;
1146   std::cout << "End of Surface Mesh generation" << std::endl;
1147   std::cout << std::endl;
1148
1149   mesh_t *msh;
1150   blsurf_data_get_mesh(bls, &msh);
1151   if(!msh){
1152     blsurf_session_delete(bls);
1153     cad_delete(c);
1154     context_delete(ctx);
1155
1156     return error(_comment);
1157     //return false;
1158   }
1159
1160   integer nv, ne, nt, nq, vtx[4], tag;
1161   real xyz[3];
1162
1163   mesh_get_vertex_count(msh, &nv);
1164   mesh_get_edge_count(msh, &ne);
1165   mesh_get_triangle_count(msh, &nt);
1166   mesh_get_quadrangle_count(msh, &nq);
1167
1168
1169   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
1170   SMDS_MeshNode** nodes = new SMDS_MeshNode*[nv+1];
1171   bool* tags = new bool[nv+1];
1172
1173   for(int iv=1;iv<=nv;iv++) {
1174     mesh_get_vertex_coordinates(msh, iv, xyz);
1175     mesh_get_vertex_tag(msh, iv, &tag);
1176     nodes[iv] = meshDS->AddNode(xyz[0], xyz[1], xyz[2]);
1177     // internal point are tagged to zero
1178     if(tag){
1179       meshDS->SetNodeOnVertex(nodes[iv], TopoDS::Vertex(pmap(tag)));
1180       tags[iv] = false;
1181     } else {
1182       tags[iv] = true;
1183     }
1184   }
1185
1186   for(int it=1;it<=ne;it++) {
1187     mesh_get_edge_vertices(msh, it, vtx);
1188     SMDS_MeshEdge* edg = meshDS->AddEdge(nodes[vtx[0]], nodes[vtx[1]]);
1189     mesh_get_edge_tag(msh, it, &tag);
1190
1191     if (tags[vtx[0]]) {
1192       Set_NodeOnEdge(meshDS, nodes[vtx[0]], emap(tag));
1193       tags[vtx[0]] = false;
1194     };
1195     if (tags[vtx[1]]) {
1196       Set_NodeOnEdge(meshDS, nodes[vtx[1]], emap(tag));
1197       tags[vtx[1]] = false;
1198     };
1199     meshDS->SetMeshElementOnShape(edg, TopoDS::Edge(emap(tag)));
1200
1201   }
1202
1203   for(int it=1;it<=nt;it++) {
1204     mesh_get_triangle_vertices(msh, it, vtx);
1205     SMDS_MeshFace* tri = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]]);
1206     mesh_get_triangle_tag(msh, it, &tag);
1207     meshDS->SetMeshElementOnShape(tri, TopoDS::Face(fmap(tag)));
1208     if (tags[vtx[0]]) {
1209       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
1210       tags[vtx[0]] = false;
1211     };
1212     if (tags[vtx[1]]) {
1213       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
1214       tags[vtx[1]] = false;
1215     };
1216     if (tags[vtx[2]]) {
1217       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
1218       tags[vtx[2]] = false;
1219     };
1220   }
1221
1222   for(int it=1;it<=nq;it++) {
1223     mesh_get_quadrangle_vertices(msh, it, vtx);
1224     SMDS_MeshFace* quad = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]], nodes[vtx[3]]);
1225     mesh_get_quadrangle_tag(msh, it, &tag);
1226     meshDS->SetMeshElementOnShape(quad, TopoDS::Face(fmap(tag)));
1227     if (tags[vtx[0]]) {
1228       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
1229       tags[vtx[0]] = false;
1230     };
1231     if (tags[vtx[1]]) {
1232       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
1233       tags[vtx[1]] = false;
1234     };
1235     if (tags[vtx[2]]) {
1236       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
1237       tags[vtx[2]] = false;
1238     };
1239     if (tags[vtx[3]]) {
1240       meshDS->SetNodeOnFace(nodes[vtx[3]], TopoDS::Face(fmap(tag)));
1241       tags[vtx[3]] = false;
1242     };
1243   }
1244
1245   delete nodes;
1246
1247   /* release the mesh object */
1248   blsurf_data_regain_mesh(bls, msh);
1249
1250   /* clean up everything */
1251   blsurf_session_delete(bls);
1252   cad_delete(c);
1253
1254   context_delete(ctx);
1255
1256   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
1257 #ifndef WNT
1258   if ( oldFEFlags > 0 )
1259     feenableexcept( oldFEFlags );
1260   feclearexcept( FE_ALL_EXCEPT );
1261 #endif
1262
1263   std::cout << "FacesWithSizeMap" << std::endl;
1264   FacesWithSizeMap.Statistics(std::cout);
1265   std::cout << "EdgesWithSizeMap" << std::endl;
1266   EdgesWithSizeMap.Statistics(std::cout);
1267   std::cout << "VertecesWithSizeMap" << std::endl;
1268   VertecesWithSizeMap.Statistics(std::cout);
1269   std::cout << "FacesWithEnforcedVerteces" << std::endl;
1270   FacesWithEnforcedVerteces.Statistics(std::cout);
1271
1272   return true;
1273 }
1274
1275 //=============================================================================
1276 /*!
1277  *  SetNodeOnEdge
1278  */
1279 //=============================================================================
1280
1281 void BLSURFPlugin_BLSURF::Set_NodeOnEdge(SMESHDS_Mesh* meshDS, SMDS_MeshNode* node, const TopoDS_Shape& ed) {
1282   const TopoDS_Edge edge = TopoDS::Edge(ed);
1283
1284   gp_Pnt pnt(node->X(), node->Y(), node->Z());
1285
1286   Standard_Real p0 = 0.0;
1287   Standard_Real p1 = 1.0;
1288   Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, p0, p1);
1289
1290   GeomAPI_ProjectPointOnCurve proj(pnt, curve);
1291
1292   double pa = (double)proj.Parameter(1);
1293
1294   GProp_GProps LProps;
1295   BRepGProp::LinearProperties(ed, LProps);
1296   double lg = (double)LProps.Mass();
1297
1298   meshDS->SetNodeOnEdge(node, edge, pa);
1299 }
1300
1301 //=============================================================================
1302 /*!
1303  *
1304  */
1305 //=============================================================================
1306
1307 ostream & BLSURFPlugin_BLSURF::SaveTo(ostream & save)
1308 {
1309   return save;
1310 }
1311
1312 //=============================================================================
1313 /*!
1314  *
1315  */
1316 //=============================================================================
1317
1318 istream & BLSURFPlugin_BLSURF::LoadFrom(istream & load)
1319 {
1320   return load;
1321 }
1322
1323 //=============================================================================
1324 /*!
1325  *
1326  */
1327 //=============================================================================
1328
1329 ostream & operator << (ostream & save, BLSURFPlugin_BLSURF & hyp)
1330 {
1331   return hyp.SaveTo( save );
1332 }
1333
1334 //=============================================================================
1335 /*!
1336  *
1337  */
1338 //=============================================================================
1339
1340 istream & operator >> (istream & load, BLSURFPlugin_BLSURF & hyp)
1341 {
1342   return hyp.LoadFrom( load );
1343 }
1344
1345 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data)
1346 {
1347   const Geom2d_Curve*pargeo = (const Geom2d_Curve*) user_data;
1348
1349   if (uv){
1350     gp_Pnt2d P;
1351     P=pargeo->Value(t);
1352     uv[0]=P.X(); uv[1]=P.Y();
1353   }
1354
1355   if(dt) {
1356     gp_Vec2d V1;
1357     V1=pargeo->DN(t,1);
1358     dt[0]=V1.X(); dt[1]=V1.Y();
1359   }
1360
1361   if(dtt){
1362     gp_Vec2d V2;
1363     V2=pargeo->DN(t,2);
1364     dtt[0]=V2.X(); dtt[1]=V2.Y();
1365   }
1366
1367   return 0;
1368 }
1369
1370 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
1371                   real *duu, real *duv, real *dvv, void *user_data)
1372 {
1373   const Geom_Surface* geometry = (const Geom_Surface*) user_data;
1374
1375   if(xyz){
1376    gp_Pnt P;
1377    P=geometry->Value(uv[0],uv[1]);   // S.D0(U,V,P);
1378    xyz[0]=P.X(); xyz[1]=P.Y(); xyz[2]=P.Z();
1379   }
1380
1381   if(du && dv){
1382     gp_Pnt P;
1383     gp_Vec D1U,D1V;
1384
1385     geometry->D1(uv[0],uv[1],P,D1U,D1V);
1386     du[0]=D1U.X(); du[1]=D1U.Y(); du[2]=D1U.Z();
1387     dv[0]=D1V.X(); dv[1]=D1V.Y(); dv[2]=D1V.Z();
1388   }
1389
1390   if(duu && duv && dvv){
1391     gp_Pnt P;
1392     gp_Vec D1U,D1V;
1393     gp_Vec D2U,D2V,D2UV;
1394
1395     geometry->D2(uv[0],uv[1],P,D1U,D1V,D2U,D2V,D2UV);
1396     duu[0]=D2U.X(); duu[1]=D2U.Y(); duu[2]=D2U.Z();
1397     duv[0]=D2UV.X(); duv[1]=D2UV.Y(); duv[2]=D2UV.Z();
1398     dvv[0]=D2V.X(); dvv[1]=D2V.Y(); dvv[2]=D2V.Z();
1399   }
1400
1401   return 0;
1402 }
1403
1404
1405 status_t size_on_surface(integer face_id, real *uv, real *size, void *user_data)
1406 {
1407   if (face_id == 1) {
1408     if (my_u_min > uv[0]) {
1409       my_u_min = uv[0];
1410     }
1411     if (my_v_min > uv[1]) {
1412       my_v_min = uv[1];
1413     }
1414     if (my_u_max < uv[0]) {
1415       my_u_max = uv[0];
1416     }
1417     if (my_v_max < uv[1]) {
1418       my_v_max = uv[1];
1419     }
1420   }
1421
1422   if (FaceId2PythonSmp.count(face_id) != 0){
1423     PyObject * pyresult = NULL;
1424     PyObject* new_stderr = NULL;
1425     assert(Py_IsInitialized());
1426     PyGILState_STATE gstate;
1427     gstate = PyGILState_Ensure();
1428     pyresult = PyObject_CallFunction(FaceId2PythonSmp[face_id],"(f,f)",uv[0],uv[1]);
1429     double result;
1430     if ( pyresult == NULL){
1431       fflush(stderr);
1432       string err_description="";
1433       new_stderr = newPyStdOut(err_description);
1434       PySys_SetObject("stderr", new_stderr);
1435       PyErr_Print();
1436       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1437       Py_DECREF(new_stderr);
1438       MESSAGE("Can't evaluate f(" << uv[0] << "," << uv[1] << ")" << " error is " << err_description);
1439       result = *((double*)user_data);
1440       }
1441     else {
1442       result = PyFloat_AsDouble(pyresult);
1443       Py_DECREF(pyresult);
1444     }
1445     *size = result;
1446     //MESSAGE("f(" << uv[0] << "," << uv[1] << ")" << " = " << result);
1447     PyGILState_Release(gstate);
1448   }
1449   else {
1450     *size = *((double*)user_data);
1451   }
1452   return STATUS_OK;
1453 }
1454
1455 status_t size_on_edge(integer edge_id, real t, real *size, void *user_data)
1456 {
1457   if (EdgeId2PythonSmp.count(edge_id) != 0){
1458     PyObject * pyresult = NULL;
1459     PyObject* new_stderr = NULL;
1460     assert(Py_IsInitialized());
1461     PyGILState_STATE gstate;
1462     gstate = PyGILState_Ensure();
1463     pyresult = PyObject_CallFunction(EdgeId2PythonSmp[edge_id],"(f)",t);
1464     double result;
1465     if ( pyresult == NULL){
1466       fflush(stderr);
1467       string err_description="";
1468       new_stderr = newPyStdOut(err_description);
1469       PySys_SetObject("stderr", new_stderr);
1470       PyErr_Print();
1471       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1472       Py_DECREF(new_stderr);
1473       MESSAGE("Can't evaluate f(" << t << ")" << " error is " << err_description);
1474       result = *((double*)user_data);
1475       }
1476     else {
1477       result = PyFloat_AsDouble(pyresult);
1478       Py_DECREF(pyresult);
1479     }
1480     *size = result;
1481     PyGILState_Release(gstate);
1482   }
1483   else {
1484     *size = *((double*)user_data);
1485   }
1486   return STATUS_OK;
1487 }
1488
1489 status_t size_on_vertex(integer point_id, real *size, void *user_data)
1490 {
1491   if (VertexId2PythonSmp.count(point_id) != 0){
1492     PyObject * pyresult = NULL;
1493     PyObject* new_stderr = NULL;
1494     assert(Py_IsInitialized());
1495     PyGILState_STATE gstate;
1496     gstate = PyGILState_Ensure();
1497     pyresult = PyObject_CallFunction(VertexId2PythonSmp[point_id],"");
1498     double result;
1499     if ( pyresult == NULL){
1500       fflush(stderr);
1501       string err_description="";
1502       new_stderr = newPyStdOut(err_description);
1503       PySys_SetObject("stderr", new_stderr);
1504       PyErr_Print();
1505       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1506       Py_DECREF(new_stderr);
1507       MESSAGE("Can't evaluate f()" << " error is " << err_description);
1508       result = *((double*)user_data);
1509       }
1510     else {
1511       result = PyFloat_AsDouble(pyresult);
1512       Py_DECREF(pyresult);
1513     }
1514     *size = result;
1515     PyGILState_Release(gstate);
1516   }
1517   else {
1518     *size = *((double*)user_data);
1519   }
1520  return STATUS_OK;
1521 }
1522
1523 status_t message_callback(message_t *msg, void *user_data)
1524 {
1525   integer errnumber = 0;
1526   char *desc;
1527   message_get_number(msg, &errnumber);
1528   message_get_description(msg, &desc);
1529   if ( errnumber < 0 ) {
1530     string * error = (string*)user_data;
1531 //   if ( !error->empty() )
1532 //     *error += "\n";
1533     // remove ^A from the tail
1534     int len = strlen( desc );
1535     while (len > 0 && desc[len-1] != '\n')
1536       len--;
1537     error->append( desc, len );
1538   }
1539   else {
1540       std::cout << desc << std::endl;
1541   }
1542   return STATUS_OK;
1543 }
1544
1545
1546 //=============================================================================
1547 /*!
1548  *
1549  */
1550 //=============================================================================
1551 bool BLSURFPlugin_BLSURF::Evaluate(SMESH_Mesh& aMesh,
1552                                    const TopoDS_Shape& aShape,
1553                                    MapShapeNbElems& aResMap)
1554 {
1555   int    _physicalMesh  = BLSURFPlugin_Hypothesis::GetDefaultPhysicalMesh();
1556   double _phySize       = BLSURFPlugin_Hypothesis::GetDefaultPhySize();
1557   //int    _geometricMesh = BLSURFPlugin_Hypothesis::GetDefaultGeometricMesh();
1558   //double _angleMeshS    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshS();
1559   double _angleMeshC    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshC();
1560   if(_hypothesis) {
1561     _physicalMesh  = (int) _hypothesis->GetPhysicalMesh();
1562     _phySize       = _hypothesis->GetPhySize();
1563     //_geometricMesh = (int) hyp->GetGeometricMesh();
1564     //_angleMeshS    = hyp->GetAngleMeshS();
1565     _angleMeshC    = _hypothesis->GetAngleMeshC();
1566   }
1567
1568   bool IsQuadratic = false;
1569
1570   // ----------------
1571   // evaluate 1D
1572   // ----------------
1573   TopTools_DataMapOfShapeInteger EdgesMap;
1574   double fullLen = 0.0;
1575   double fullNbSeg = 0;
1576   for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
1577     TopoDS_Edge E = TopoDS::Edge( exp.Current() );
1578     if( EdgesMap.IsBound(E) )
1579       continue;
1580     SMESH_subMesh *sm = aMesh.GetSubMesh(E);
1581     double aLen = SMESH_Algo::EdgeLength(E);
1582     fullLen += aLen;
1583     int nb1d = 0;
1584     if(_physicalMesh==1) {
1585        nb1d = (int)( aLen/_phySize + 1 );
1586     }
1587     else {
1588       // use geometry
1589       double f,l;
1590       Handle(Geom_Curve) C = BRep_Tool::Curve(E,f,l);
1591       double fullAng = 0.0;
1592       double dp = (l-f)/200;
1593       gp_Pnt P1,P2,P3;
1594       C->D0(f,P1);
1595       C->D0(f+dp,P2);
1596       gp_Vec V1(P1,P2);
1597       for(int j=2; j<=200; j++) {
1598         C->D0(f+dp*j,P3);
1599         gp_Vec V2(P2,P3);
1600         fullAng += fabs(V1.Angle(V2));
1601         V1 = V2;
1602         P2 = P3;
1603       }
1604       nb1d = (int)( fullAng/_angleMeshC + 1 );
1605     }
1606     fullNbSeg += nb1d;
1607     std::vector<int> aVec(17);
1608     for(int i=0; i<17; i++) aVec[i]=0;
1609     if( IsQuadratic > 0 ) {
1610       aVec[0] = 2*nb1d - 1;
1611       aVec[2] = nb1d;
1612     }
1613     else {
1614       aVec[0] = nb1d - 1;
1615       aVec[1] = nb1d;
1616     }
1617     aResMap.insert(std::make_pair(sm,aVec));
1618     EdgesMap.Bind(E,nb1d);
1619   }
1620   double ELen = fullLen/fullNbSeg;
1621   // ----------------
1622   // evaluate 2D
1623   // ----------------
1624   // try to evaluate as in MEFISTO
1625   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
1626     TopoDS_Face F = TopoDS::Face( exp.Current() );
1627     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
1628     GProp_GProps G;
1629     BRepGProp::SurfaceProperties(F,G);
1630     double anArea = G.Mass();
1631     int nb1d = 0;
1632     for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next()) {
1633       nb1d += EdgesMap.Find(exp1.Current());
1634     }
1635     int nbFaces = (int) ( anArea/( ELen*ELen*sqrt(3.) / 4 ) );
1636     int nbNodes = (int) ( ( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
1637     std::vector<int> aVec(17);
1638     for(int i=0; i<17; i++) aVec[i]=0;
1639     if( IsQuadratic ) {
1640       int nb1d_in = (nbFaces*3 - nb1d) / 2;
1641       aVec[0] = nbNodes + nb1d_in;
1642       aVec[4] = nbFaces;
1643     }
1644     else {
1645       aVec[0] = nbNodes;
1646       aVec[3] = nbFaces;
1647     }
1648     aResMap.insert(std::make_pair(sm,aVec));
1649   }
1650
1651   // ----------------
1652   // evaluate 3D
1653   // ----------------
1654   GProp_GProps G;
1655   BRepGProp::VolumeProperties(aShape,G);
1656   double aVolume = G.Mass();
1657   double tetrVol = 0.1179*ELen*ELen*ELen;
1658   int nbVols = (int)aVolume/tetrVol;
1659   int nb1d_in = (int) ( ( nbVols*6 - fullNbSeg ) / 6 );
1660   std::vector<int> aVec(17);
1661   for(int i=0; i<17; i++) aVec[i]=0;
1662   if( IsQuadratic ) {
1663     aVec[0] = nb1d_in/3 + 1 + nb1d_in;
1664     aVec[9] = nbVols;
1665   }
1666   else {
1667     aVec[0] = nb1d_in/3 + 1;
1668     aVec[8] = nbVols;
1669   }
1670   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
1671   aResMap.insert(std::make_pair(sm,aVec));
1672
1673   return true;
1674 }
1675
1676 //=============================================================================
1677 /*!
1678  *  Rewritting of the BRepClass_FaceClassifier::Perform function which is bugged (CAS 6.3sp6)
1679  *  Following line was added:
1680  *        myExtrem.Perform(P);
1681  */
1682 //=============================================================================
1683 void  BLSURFPlugin_BLSURF::BRepClass_FaceClassifierPerform(BRepClass_FaceClassifier* fc,
1684                     const TopoDS_Face& face,
1685                     const gp_Pnt& P,
1686                     const Standard_Real Tol)
1687 {
1688   //-- Voir BRepExtrema_ExtPF.cxx
1689   BRepAdaptor_Surface Surf(face);
1690   Standard_Real U1, U2, V1, V2;
1691   BRepTools::UVBounds(face, U1, U2, V1, V2);
1692   Extrema_ExtPS myExtrem;
1693   myExtrem.Initialize(Surf, U1, U2, V1, V2, Tol, Tol);
1694   myExtrem.Perform(P);
1695   //----------------------------------------------------------
1696   //-- On cherche le point le plus proche , PUIS
1697   //-- On le classifie.
1698   Standard_Integer nbv    = 0; // xpu
1699   Standard_Real MaxDist   =  RealLast();
1700   Standard_Integer indice = 0;
1701   if(myExtrem.IsDone()) {
1702     nbv = myExtrem.NbExt();
1703     for (Standard_Integer i = 1; i <= nbv; i++) {
1704       Standard_Real d = myExtrem.Value(i);
1705       d = Abs(d);
1706       if(d <= MaxDist) {
1707     MaxDist = d;
1708     indice = i;
1709       }
1710     }
1711   }
1712   if(indice) {
1713     gp_Pnt2d Puv;
1714     Standard_Real U1,U2;
1715     myExtrem.Point(indice).Parameter(U1, U2);
1716     Puv.SetCoord(U1, U2);
1717     fc->Perform(face, Puv, Tol);
1718   }
1719   else {
1720     fc->Perform(face, gp_Pnt2d(U1-1.0,V1 - 1.0), Tol); //-- NYI etc BUG PAS BEAU En attendant l acces a rejected
1721     //-- le resultat est TopAbs_OUT;
1722   }
1723 }