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