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 "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("Vertex with key " << key << " already in map");
 677               }
 678               VertexId2SizeMap[key] = smIt->second;
 679             }
 680           }
 681         }
 682         // Single face
 683         if (GeomType == TopAbs_FACE){
 684           HasSizeMapOnFace = true;
 685           if (! FacesWithSizeMap.Contains(TopoDS::Face(GeomShape))) {
 686             key = FacesWithSizeMap.Add(TopoDS::Face(GeomShape));
 687           }
 688           else {
 689             key = FacesWithSizeMap.FindIndex(TopoDS::Face(GeomShape));
 690 //             MESSAGE("Face with key " << key << " already in map");
 691           }
 692           FaceId2SizeMap[key] = smIt->second;
 693         }
 694         // Single edge
 695         if (GeomType == TopAbs_EDGE){
 696           HasSizeMapOnEdge = true;
 697           HasSizeMapOnFace = true;
 698           if (! EdgesWithSizeMap.Contains(TopoDS::Edge(GeomShape))) {
 699             key = EdgesWithSizeMap.Add(TopoDS::Edge(GeomShape));
 700           }
 701           else {
 702             key = EdgesWithSizeMap.FindIndex(TopoDS::Edge(GeomShape));
 703 //             MESSAGE("Edge with key " << key << " already in map");
 704           }
 705           EdgeId2SizeMap[key] = smIt->second;
 706         }
 707         // Single vertex
 708         if (GeomType == TopAbs_VERTEX){
 709           HasSizeMapOnVertex = true;
 710           HasSizeMapOnEdge   = true;
 711           HasSizeMapOnFace   = true;
 712           if (! VerticesWithSizeMap.Contains(TopoDS::Vertex(GeomShape))) {
 713             key = VerticesWithSizeMap.Add(TopoDS::Vertex(GeomShape));
 714           }
 715           else {
 716             key = VerticesWithSizeMap.FindIndex(TopoDS::Vertex(GeomShape));
 717 //             MESSAGE("Vertex with key " << key << " already in map");
 718           }
 719           VertexId2SizeMap[key] = smIt->second;
 720         }
 721       }
 722     }
 723
 724     //
 725     // Attractors
 726     //
 727     MESSAGE("Setting Attractors");
 728     const BLSURFPlugin_Hypothesis::TSizeMap attractors = BLSURFPlugin_Hypothesis::GetAttractorEntries(hyp);
 729     BLSURFPlugin_Hypothesis::TSizeMap::const_iterator atIt = attractors.begin();
 730     for ( ; atIt != attractors.end(); ++atIt ) {
 731       if ( !atIt->second.empty() ) {
 732         MESSAGE("blsurf_set_attractor(): " << atIt->first << " = " << atIt->second);
 733         GeomShape = entryToShape(atIt->first);
 734         GeomType  = GeomShape.ShapeType();
 735         // Group Management
 736         if (GeomType == TopAbs_COMPOUND){
 737           for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
 738             if (it.Value().ShapeType() == TopAbs_FACE){
 739               HasSizeMapOnFace = true;
 740               createAttractorOnFace(it.Value(), atIt->second);
 741             }
 742           }
 743         }
 744
 745         if (GeomType == TopAbs_FACE){
 746           HasSizeMapOnFace = true;
 747           createAttractorOnFace(GeomShape, atIt->second);
 748         }
 749 /*
 750         if (GeomType == TopAbs_EDGE){
 751           HasSizeMapOnEdge = true;
 752           HasSizeMapOnFace = true;
 753         EdgeId2SizeMap[TopoDS::Edge(GeomShape).HashCode(IntegerLast())] = atIt->second;
 754         }
 755         if (GeomType == TopAbs_VERTEX){
 756           HasSizeMapOnVertex = true;
 757           HasSizeMapOnEdge   = true;
 758           HasSizeMapOnFace   = true;
 759         VertexId2SizeMap[TopoDS::Vertex(GeomShape).HashCode(IntegerLast())] = atIt->second;
 760         }
 761 */
 762       }
 763     }
 764
 765
 766     //
 767     // Enforced Vertices
 768     //
 769     MESSAGE("Setting Enforced Vertices");
 770     const BLSURFPlugin_Hypothesis::TEnforcedVertexMap enforcedVertexMap = BLSURFPlugin_Hypothesis::GetAllEnforcedVertices(hyp);
 771     BLSURFPlugin_Hypothesis::TEnforcedVertexMap::const_iterator enfIt = enforcedVertexMap.begin();
 772     for ( ; enfIt != enforcedVertexMap.end(); ++enfIt ) {
 773       if ( !enfIt->second.empty() ) {
 774         GeomShape = entryToShape(enfIt->first);
 775         GeomType  = GeomShape.ShapeType();
 776         // Group Management
 777         if (GeomType == TopAbs_COMPOUND){
 778           for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
 779             if (it.Value().ShapeType() == TopAbs_FACE){
 780               HasSizeMapOnFace = true;
 781               createEnforcedVertexOnFace(it.Value(), enfIt->second);
 782             }
 783           }
 784         }
 785
 786         if (GeomType == TopAbs_FACE){
 787           HasSizeMapOnFace = true;
 788           createEnforcedVertexOnFace(GeomShape, enfIt->second);
 789         }
 790       }
 791     }
 792
 793 //    if (HasSizeMapOnFace){
 794     // In all size map cases (hphy_flag = 2), at least map on face must be defined
 795     MESSAGE("Setting Size Map on FACES ");
 796     blsurf_data_set_sizemap_iso_cad_face(bls, size_on_surface, &_smp_phy_size);
 797 //    }
 798
 799     if (HasSizeMapOnEdge){
 800       MESSAGE("Setting Size Map on EDGES ");
 801       blsurf_data_set_sizemap_iso_cad_edge(bls, size_on_edge, &_smp_phy_size);
 802     }
 803     if (HasSizeMapOnVertex){
 804       MESSAGE("Setting Size Map on VERTICES ");
 805       blsurf_data_set_sizemap_iso_cad_point(bls, size_on_vertex, &_smp_phy_size);
 806     }
 807   }
 808   blsurf_set_param(bls, "hphydef",           to_string(_phySize).c_str());
 809   blsurf_set_param(bls, "hgeo_flag",         to_string(_geometricMesh).c_str());
 810   blsurf_set_param(bls, "relax_size",        _decimesh ? "0": to_string(_geometricMesh).c_str());
 811   blsurf_set_param(bls, "angle_meshs",       to_string(_angleMeshS).c_str());
 812   blsurf_set_param(bls, "angle_meshc",       to_string(_angleMeshC).c_str());
 813   blsurf_set_param(bls, "gradation",         to_string(_gradation).c_str());
 814   blsurf_set_param(bls, "patch_independent", _decimesh ? "1" : "0");
 815   blsurf_set_param(bls, "element",           _quadAllowed ? "q1.0" : "p1");
 816   blsurf_set_param(bls, "verb",              to_string(_verb).c_str());
 817 }
 818
 819 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data);
 820 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
 821                   real *duu, real *duv, real *dvv, void *user_data);
 822 status_t message_callback(message_t *msg, void *user_data);
 823
 824 //=============================================================================
 825 /*!
 826  *
 827  */
 828 //=============================================================================
 829
 830 bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) {
 831
 832   MESSAGE("BLSURFPlugin_BLSURF::Compute");
 833
 834   if (aShape.ShapeType() == TopAbs_COMPOUND) {
 835     MESSAGE("  the shape is a COMPOUND");
 836   }
 837   else {
 838     MESSAGE("  the shape is UNKNOWN");
 839   };
 840
 841   context_t *ctx =  context_new();
 842   context_set_message_callback(ctx, message_callback, &_comment);
 843
 844   cad_t *c = cad_new(ctx);
 845
 846   blsurf_session_t *bls = blsurf_session_new(ctx);
 847
 848   FacesWithSizeMap.Clear();
 849   FaceId2SizeMap.clear();
 850   EdgesWithSizeMap.Clear();
 851   EdgeId2SizeMap.clear();
 852   VerticesWithSizeMap.Clear();
 853   VertexId2SizeMap.clear();
 854
 855   MESSAGE("BEGIN SetParameters");
 856   SetParameters(_hypothesis, bls);
 857   MESSAGE("END SetParameters");
 858
 859   TopTools_IndexedMapOfShape fmap;
 860   TopTools_IndexedMapOfShape emap;
 861   TopTools_IndexedMapOfShape pmap;
 862   vector<Handle(Geom2d_Curve)> curves;
 863   vector<Handle(Geom_Surface)> surfaces;
 864
 865   fmap.Clear();
 866   FaceId2PythonSmp.clear();
 867   emap.Clear();
 868   EdgeId2PythonSmp.clear();
 869   pmap.Clear();
 870   VertexId2PythonSmp.clear();
 871   surfaces.resize(0);
 872   curves.resize(0);
 873
 874   assert(Py_IsInitialized());
 875   PyGILState_STATE gstate;
 876   gstate = PyGILState_Ensure();
 877
 878   string theSizeMapStr;
 879
 880   /****************************************************************************************
 881                                   FACES
 882   *****************************************************************************************/
 883   int iface = 0;
 884   string bad_end = "return";
 885   int faceKey = -1;
 886   for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next()) {
 887     TopoDS_Face f=TopoDS::Face(face_iter.Current());
 888
 889     if (fmap.FindIndex(f) > 0)
 890       continue;
 891
 892     fmap.Add(f);
 893     iface++;
 894     surfaces.push_back(BRep_Tool::Surface(f));
 895
 896     cad_face_t *fce = cad_face_new(c, iface, surf_fun, surfaces.back());
 897     cad_face_set_tag(fce, iface);
 898     if(f.Orientation() != TopAbs_FORWARD){
 899       cad_face_set_orientation(fce, CAD_ORIENTATION_REVERSED);
 900     } else {
 901       cad_face_set_orientation(fce, CAD_ORIENTATION_FORWARD);
 902     }
 903
 904     if (HasSizeMapOnFace){
 905       std::cout << "A size map is defined on a face" << std::endl;
 906       // Classic size map
 907       faceKey = FacesWithSizeMap.FindIndex(f);
 908
 909       if (FaceId2SizeMap.find(faceKey)!=FaceId2SizeMap.end()){
 910         theSizeMapStr = FaceId2SizeMap[faceKey];
 911         // check if function ends with "return"
 912         if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
 913           continue;
 914         // Expr To Python function, verification is performed at validation in GUI
 915         PyObject * obj = NULL;
 916         obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
 917         Py_DECREF(obj);
 918         PyObject * func = NULL;
 919         func = PyObject_GetAttrString(main_mod, "f");
 920         FaceId2PythonSmp[iface]=func;
 921         FaceId2SizeMap.erase(faceKey);
 922       }
 923
 924       // Specific size map = Attractor
 925       std::map<int,std::vector<double> >::iterator attractor_iter = FaceId2AttractorCoords.begin();
 926       int iatt=0;
 927       for (; attractor_iter != FaceId2AttractorCoords.end(); ++attractor_iter) {
 928         if (attractor_iter->first == faceKey) {
 929           MESSAGE("Face indice: " << iface);
 930           MESSAGE("Adding attractor");
 931
 932           double xyzCoords[3]  = {attractor_iter->second[2],
 933                                   attractor_iter->second[3],
 934                                   attractor_iter->second[4]};
 935
 936           MESSAGE("Check position of vertex =(" << xyzCoords[0] << "," << xyzCoords[1] << "," << xyzCoords[2] << ")");
 937           gp_Pnt P(xyzCoords[0],xyzCoords[1],xyzCoords[2]);
 938           BRepClass_FaceClassifier scl(f,P,1e-7);
 939           // scl.Perform() is bugged. The function was rewritten
 940 //          scl.Perform();
 941           BRepClass_FaceClassifierPerform(&scl,f,P,1e-7);
 942           TopAbs_State result = scl.State();
 943           MESSAGE("Position of point on face: "<<result);
 944           if ( result == TopAbs_OUT )
 945               MESSAGE("Point is out of face: node is not created");
 946           if ( result == TopAbs_UNKNOWN )
 947               MESSAGE("Point position on face is unknown: node is not created");
 948           if ( result == TopAbs_ON )
 949               MESSAGE("Point is on border of face: node is not created");
 950           if ( result == TopAbs_IN )
 951           {
 952             // Point is inside face and not on border
 953             MESSAGE("Point is in face: node is created");
 954             double uvCoords[2]   = {attractor_iter->second[0],attractor_iter->second[1]};
 955             iatt++;
 956             MESSAGE("Add cad point on (u,v)=(" << uvCoords[0] << "," << uvCoords[1] << ") with id = " << iatt);
 957             cad_point_t* point_p = cad_point_new(fce, iatt, uvCoords);
 958             cad_point_set_tag(point_p, iatt);
 959           }
 960           FaceId2AttractorCoords.erase(faceKey);
 961         }
 962       }
 963
 964       // Enforced Vertices
 965       faceKey = FacesWithEnforcedVertices.FindIndex(f);
 966       std::map<int,std::set<std::vector<double> > >::const_iterator evmIt = FaceId2EnforcedVertexCoords.find(faceKey);
 967       if (evmIt != FaceId2EnforcedVertexCoords.end()) {
 968         std::cout << "Some enforced vertices are defined" << std::endl;
 969         int ienf = 0;
 970         std::set<std::vector<double> > evl;
 971 //         std::vector<double> ev;
 972         MESSAGE("Face indice: " << iface);
 973         MESSAGE("Adding enforced vertices");
 974         evl = evmIt->second;
 975         MESSAGE("Number of vertices to add: "<< evl.size())
 976         std::set< std::vector<double> >::const_iterator evlIt = evl.begin();
 977         for (; evlIt != evl.end(); ++evlIt) {
 978 //           ev = *evlIt;
 979 //         for (int i=0; i<evl.size() ; i++) {
 980 //           ev = evl[i];
 981
 982 //           double xyzCoords[3]  = {ev[2], ev[3], ev[4]};
 983           double xyzCoords[3]  = {evlIt->at(0), evlIt->at(3), evlIt->at(4)};
 984           MESSAGE("Check position of vertex =(" << xyzCoords[0] << "," << xyzCoords[1] << "," << xyzCoords[2] << ")");
 985           gp_Pnt P(xyzCoords[0],xyzCoords[1],xyzCoords[2]);
 986           BRepClass_FaceClassifier scl(f,P,1e-7);
 987           // scl.Perform() is bugged. The function was rewritten
 988 //          scl.Perform();
 989           BRepClass_FaceClassifierPerform(&scl,f,P,1e-7);
 990           TopAbs_State result = scl.State();
 991           MESSAGE("Position of point on face: "<<result);
 992           if ( result == TopAbs_OUT )
 993               MESSAGE("Point is out of face: node is not created");
 994           if ( result == TopAbs_UNKNOWN )
 995               MESSAGE("Point position on face is unknown: node is not created");
 996           if ( result == TopAbs_ON )
 997               MESSAGE("Point is on border of face: node is not created");
 998           if ( result == TopAbs_IN )
 999           {
1000             // Point is inside face and not on border
1001             MESSAGE("Point is in face: node is created");
1002 //             double uvCoords[2]   = {ev[0],ev[1]};
1003             double uvCoords[2]   = {evlIt->at(0),evlIt->at(1)};
1004             ienf++;
1005             MESSAGE("Add cad point on (u,v)=(" << uvCoords[0] << "," << uvCoords[1] << ") with id = " << ienf);
1006             cad_point_t* point_p = cad_point_new(fce, ienf, uvCoords);
1007             cad_point_set_tag(point_p, ienf);
1008           }
1009         }
1010         FaceId2EnforcedVertexCoords.erase(faceKey);
1011       }
1012       else
1013         std::cout << "No enforced vertex defined" << std::endl;
1014     }
1015
1016
1017     /****************************************************************************************
1018                                     EDGES
1019     *****************************************************************************************/
1020     int edgeKey = -1;
1021     for (TopExp_Explorer edge_iter(f,TopAbs_EDGE);edge_iter.More();edge_iter.Next()) {
1022       TopoDS_Edge e = TopoDS::Edge(edge_iter.Current());
1023       int ic = emap.FindIndex(e);
1024       if (ic <= 0)
1025         ic = emap.Add(e);
1026
1027       double tmin,tmax;
1028       curves.push_back(BRep_Tool::CurveOnSurface(e, f, tmin, tmax));
1029
1030       if (HasSizeMapOnEdge){
1031         edgeKey = EdgesWithSizeMap.FindIndex(e);
1032         if (EdgeId2SizeMap.find(edgeKey)!=EdgeId2SizeMap.end()) {
1033           theSizeMapStr = EdgeId2SizeMap[faceKey];
1034           if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
1035             continue;
1036           // Expr To Python function, verification is performed at validation in GUI
1037           PyObject * obj = NULL;
1038           obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
1039           Py_DECREF(obj);
1040           PyObject * func = NULL;
1041           func = PyObject_GetAttrString(main_mod, "f");
1042           EdgeId2PythonSmp[ic]=func;
1043           EdgeId2SizeMap.erase(edgeKey);
1044         }
1045       }
1046       cad_edge_t *edg = cad_edge_new(fce, ic, tmin, tmax, curv_fun, curves.back());
1047       cad_edge_set_tag(edg, ic);
1048       cad_edge_set_property(edg, EDGE_PROPERTY_SOFT_REQUIRED);
1049       if (e.Orientation() == TopAbs_INTERNAL)
1050         cad_edge_set_property(edg, EDGE_PROPERTY_INTERNAL);
1051
1052       int npts = 0;
1053       int ip1, ip2, *ip;
1054       gp_Pnt2d e0 = curves.back()->Value(tmin);
1055       gp_Pnt ee0 = surfaces.back()->Value(e0.X(), e0.Y());
1056       Standard_Real d1=0,d2=0;
1057
1058
1059       /****************************************************************************************
1060                                       VERTICES
1061       *****************************************************************************************/
1062       int vertexKey = -1;
1063       for (TopExp_Explorer ex_edge(e ,TopAbs_VERTEX); ex_edge.More(); ex_edge.Next()) {
1064         TopoDS_Vertex v = TopoDS::Vertex(ex_edge.Current());
1065         ++npts;
1066         if (npts == 1){
1067           ip = &ip1;
1068           d1 = ee0.SquareDistance(BRep_Tool::Pnt(v));
1069         } else {
1070           ip = &ip2;
1071           d2 = ee0.SquareDistance(BRep_Tool::Pnt(v));
1072         }
1073         *ip = pmap.FindIndex(v);
1074         if(*ip <= 0)
1075           *ip = pmap.Add(v);
1076
1077         vertexKey = VerticesWithSizeMap.FindIndex(v);
1078         if (HasSizeMapOnVertex){
1079           vertexKey = VerticesWithSizeMap.FindIndex(v);
1080           if (VertexId2SizeMap.find(vertexKey)!=VertexId2SizeMap.end()){
1081             theSizeMapStr = VertexId2SizeMap[faceKey];
1082             if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
1083               continue;
1084             // Expr To Python function, verification is performed at validation in GUI
1085             PyObject * obj = NULL;
1086             obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
1087             Py_DECREF(obj);
1088             PyObject * func = NULL;
1089             func = PyObject_GetAttrString(main_mod, "f");
1090             VertexId2PythonSmp[*ip]=func;
1091 //             VertexId2SizeMap.erase(vertexKey);   // do not erase if using a vector
1092           }
1093         }
1094       }
1095       if (npts != 2) {
1096         // should not happen
1097         MESSAGE("An edge does not have 2 extremities.");
1098       } else {
1099         if (d1 < d2)
1100           cad_edge_set_extremities(edg, ip1, ip2);
1101         else
1102           cad_edge_set_extremities(edg, ip2, ip1);
1103       }
1104     } // for edge
1105   } //for face
1106
1107
1108   PyGILState_Release(gstate);
1109
1110   blsurf_data_set_cad(bls, c);
1111
1112   std::cout << std::endl;
1113   std::cout << "Beginning of Surface Mesh generation" << std::endl;
1114   std::cout << std::endl;
1115
1116   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
1117 #ifndef WNT
1118   feclearexcept( FE_ALL_EXCEPT );
1119   int oldFEFlags = fedisableexcept( FE_ALL_EXCEPT );
1120 #endif
1121
1122   status_t status = STATUS_ERROR;
1123
1124   try {
1125     OCC_CATCH_SIGNALS;
1126
1127     status = blsurf_compute_mesh(bls);
1128
1129   }
1130   catch ( std::exception& exc ) {
1131     _comment += exc.what();
1132   }
1133   catch (Standard_Failure& ex) {
1134     _comment += ex.DynamicType()->Name();
1135     if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) {
1136       _comment += ": ";
1137       _comment += ex.GetMessageString();
1138     }
1139   }
1140   catch (...) {
1141     if ( _comment.empty() )
1142       _comment = "Exception in blsurf_compute_mesh()";
1143   }
1144   if ( status != STATUS_OK) {
1145     blsurf_session_delete(bls);
1146     cad_delete(c);
1147     context_delete(ctx);
1148
1149     return error(_comment);
1150   }
1151
1152   std::cout << std::endl;
1153   std::cout << "End of Surface Mesh generation" << std::endl;
1154   std::cout << std::endl;
1155
1156   mesh_t *msh;
1157   blsurf_data_get_mesh(bls, &msh);
1158   if(!msh){
1159     blsurf_session_delete(bls);
1160     cad_delete(c);
1161     context_delete(ctx);
1162
1163     return error(_comment);
1164     //return false;
1165   }
1166
1167   integer nv, ne, nt, nq, vtx[4], tag;
1168   real xyz[3];
1169
1170   mesh_get_vertex_count(msh, &nv);
1171   mesh_get_edge_count(msh, &ne);
1172   mesh_get_triangle_count(msh, &nt);
1173   mesh_get_quadrangle_count(msh, &nq);
1174
1175
1176   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
1177   SMDS_MeshNode** nodes = new SMDS_MeshNode*[nv+1];
1178   bool* tags = new bool[nv+1];
1179
1180   for(int iv=1;iv<=nv;iv++) {
1181     mesh_get_vertex_coordinates(msh, iv, xyz);
1182     mesh_get_vertex_tag(msh, iv, &tag);
1183     nodes[iv] = meshDS->AddNode(xyz[0], xyz[1], xyz[2]);
1184     // internal point are tagged to zero
1185     if(tag){
1186       meshDS->SetNodeOnVertex(nodes[iv], TopoDS::Vertex(pmap(tag)));
1187       tags[iv] = false;
1188     } else {
1189       tags[iv] = true;
1190     }
1191   }
1192
1193   for(int it=1;it<=ne;it++) {
1194     mesh_get_edge_vertices(msh, it, vtx);
1195     SMDS_MeshEdge* edg = meshDS->AddEdge(nodes[vtx[0]], nodes[vtx[1]]);
1196     mesh_get_edge_tag(msh, it, &tag);
1197
1198     if (tags[vtx[0]]) {
1199       Set_NodeOnEdge(meshDS, nodes[vtx[0]], emap(tag));
1200       tags[vtx[0]] = false;
1201     };
1202     if (tags[vtx[1]]) {
1203       Set_NodeOnEdge(meshDS, nodes[vtx[1]], emap(tag));
1204       tags[vtx[1]] = false;
1205     };
1206     meshDS->SetMeshElementOnShape(edg, TopoDS::Edge(emap(tag)));
1207
1208   }
1209
1210   for(int it=1;it<=nt;it++) {
1211     mesh_get_triangle_vertices(msh, it, vtx);
1212     SMDS_MeshFace* tri = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]]);
1213     mesh_get_triangle_tag(msh, it, &tag);
1214     meshDS->SetMeshElementOnShape(tri, TopoDS::Face(fmap(tag)));
1215     if (tags[vtx[0]]) {
1216       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
1217       tags[vtx[0]] = false;
1218     };
1219     if (tags[vtx[1]]) {
1220       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
1221       tags[vtx[1]] = false;
1222     };
1223     if (tags[vtx[2]]) {
1224       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
1225       tags[vtx[2]] = false;
1226     };
1227   }
1228
1229   for(int it=1;it<=nq;it++) {
1230     mesh_get_quadrangle_vertices(msh, it, vtx);
1231     SMDS_MeshFace* quad = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]], nodes[vtx[3]]);
1232     mesh_get_quadrangle_tag(msh, it, &tag);
1233     meshDS->SetMeshElementOnShape(quad, TopoDS::Face(fmap(tag)));
1234     if (tags[vtx[0]]) {
1235       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
1236       tags[vtx[0]] = false;
1237     };
1238     if (tags[vtx[1]]) {
1239       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
1240       tags[vtx[1]] = false;
1241     };
1242     if (tags[vtx[2]]) {
1243       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
1244       tags[vtx[2]] = false;
1245     };
1246     if (tags[vtx[3]]) {
1247       meshDS->SetNodeOnFace(nodes[vtx[3]], TopoDS::Face(fmap(tag)));
1248       tags[vtx[3]] = false;
1249     };
1250   }
1251
1252   delete nodes;
1253
1254   /* release the mesh object */
1255   blsurf_data_regain_mesh(bls, msh);
1256
1257   /* clean up everything */
1258   blsurf_session_delete(bls);
1259   cad_delete(c);
1260
1261   context_delete(ctx);
1262
1263   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
1264 #ifndef WNT
1265   if ( oldFEFlags > 0 )
1266     feenableexcept( oldFEFlags );
1267   feclearexcept( FE_ALL_EXCEPT );
1268 #endif
1269
1270   std::cout << "FacesWithSizeMap" << std::endl;
1271   FacesWithSizeMap.Statistics(std::cout);
1272   std::cout << "EdgesWithSizeMap" << std::endl;
1273   EdgesWithSizeMap.Statistics(std::cout);
1274   std::cout << "VerticesWithSizeMap" << std::endl;
1275   VerticesWithSizeMap.Statistics(std::cout);
1276   std::cout << "FacesWithEnforcedVertices" << std::endl;
1277   FacesWithEnforcedVertices.Statistics(std::cout);
1278
1279   return true;
1280 }
1281
1282 //=============================================================================
1283 /*!
1284  *  SetNodeOnEdge
1285  */
1286 //=============================================================================
1287
1288 void BLSURFPlugin_BLSURF::Set_NodeOnEdge(SMESHDS_Mesh* meshDS, SMDS_MeshNode* node, const TopoDS_Shape& ed) {
1289   const TopoDS_Edge edge = TopoDS::Edge(ed);
1290
1291   gp_Pnt pnt(node->X(), node->Y(), node->Z());
1292
1293   Standard_Real p0 = 0.0;
1294   Standard_Real p1 = 1.0;
1295   Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, p0, p1);
1296
1297   GeomAPI_ProjectPointOnCurve proj(pnt, curve);
1298
1299   double pa = (double)proj.Parameter(1);
1300
1301   GProp_GProps LProps;
1302   BRepGProp::LinearProperties(ed, LProps);
1303   double lg = (double)LProps.Mass();
1304
1305   meshDS->SetNodeOnEdge(node, edge, pa);
1306 }
1307
1308 //=============================================================================
1309 /*!
1310  *
1311  */
1312 //=============================================================================
1313
1314 ostream & BLSURFPlugin_BLSURF::SaveTo(ostream & save)
1315 {
1316   return save;
1317 }
1318
1319 //=============================================================================
1320 /*!
1321  *
1322  */
1323 //=============================================================================
1324
1325 istream & BLSURFPlugin_BLSURF::LoadFrom(istream & load)
1326 {
1327   return load;
1328 }
1329
1330 //=============================================================================
1331 /*!
1332  *
1333  */
1334 //=============================================================================
1335
1336 ostream & operator << (ostream & save, BLSURFPlugin_BLSURF & hyp)
1337 {
1338   return hyp.SaveTo( save );
1339 }
1340
1341 //=============================================================================
1342 /*!
1343  *
1344  */
1345 //=============================================================================
1346
1347 istream & operator >> (istream & load, BLSURFPlugin_BLSURF & hyp)
1348 {
1349   return hyp.LoadFrom( load );
1350 }
1351
1352 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data)
1353 {
1354   const Geom2d_Curve*pargeo = (const Geom2d_Curve*) user_data;
1355
1356   if (uv){
1357     gp_Pnt2d P;
1358     P=pargeo->Value(t);
1359     uv[0]=P.X(); uv[1]=P.Y();
1360   }
1361
1362   if(dt) {
1363     gp_Vec2d V1;
1364     V1=pargeo->DN(t,1);
1365     dt[0]=V1.X(); dt[1]=V1.Y();
1366   }
1367
1368   if(dtt){
1369     gp_Vec2d V2;
1370     V2=pargeo->DN(t,2);
1371     dtt[0]=V2.X(); dtt[1]=V2.Y();
1372   }
1373
1374   return 0;
1375 }
1376
1377 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
1378                   real *duu, real *duv, real *dvv, void *user_data)
1379 {
1380   const Geom_Surface* geometry = (const Geom_Surface*) user_data;
1381
1382   if(xyz){
1383    gp_Pnt P;
1384    P=geometry->Value(uv[0],uv[1]);   // S.D0(U,V,P);
1385    xyz[0]=P.X(); xyz[1]=P.Y(); xyz[2]=P.Z();
1386   }
1387
1388   if(du && dv){
1389     gp_Pnt P;
1390     gp_Vec D1U,D1V;
1391
1392     geometry->D1(uv[0],uv[1],P,D1U,D1V);
1393     du[0]=D1U.X(); du[1]=D1U.Y(); du[2]=D1U.Z();
1394     dv[0]=D1V.X(); dv[1]=D1V.Y(); dv[2]=D1V.Z();
1395   }
1396
1397   if(duu && duv && dvv){
1398     gp_Pnt P;
1399     gp_Vec D1U,D1V;
1400     gp_Vec D2U,D2V,D2UV;
1401
1402     geometry->D2(uv[0],uv[1],P,D1U,D1V,D2U,D2V,D2UV);
1403     duu[0]=D2U.X(); duu[1]=D2U.Y(); duu[2]=D2U.Z();
1404     duv[0]=D2UV.X(); duv[1]=D2UV.Y(); duv[2]=D2UV.Z();
1405     dvv[0]=D2V.X(); dvv[1]=D2V.Y(); dvv[2]=D2V.Z();
1406   }
1407
1408   return 0;
1409 }
1410
1411
1412 status_t size_on_surface(integer face_id, real *uv, real *size, void *user_data)
1413 {
1414   if (face_id == 1) {
1415     if (my_u_min > uv[0]) {
1416       my_u_min = uv[0];
1417     }
1418     if (my_v_min > uv[1]) {
1419       my_v_min = uv[1];
1420     }
1421     if (my_u_max < uv[0]) {
1422       my_u_max = uv[0];
1423     }
1424     if (my_v_max < uv[1]) {
1425       my_v_max = uv[1];
1426     }
1427   }
1428
1429   if (FaceId2PythonSmp.count(face_id) != 0){
1430     PyObject * pyresult = NULL;
1431     PyObject* new_stderr = NULL;
1432     assert(Py_IsInitialized());
1433     PyGILState_STATE gstate;
1434     gstate = PyGILState_Ensure();
1435     pyresult = PyObject_CallFunction(FaceId2PythonSmp[face_id],"(f,f)",uv[0],uv[1]);
1436     double result;
1437     if ( pyresult == NULL){
1438       fflush(stderr);
1439       string err_description="";
1440       new_stderr = newPyStdOut(err_description);
1441       PySys_SetObject("stderr", new_stderr);
1442       PyErr_Print();
1443       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1444       Py_DECREF(new_stderr);
1445       MESSAGE("Can't evaluate f(" << uv[0] << "," << uv[1] << ")" << " error is " << err_description);
1446       result = *((double*)user_data);
1447       }
1448     else {
1449       result = PyFloat_AsDouble(pyresult);
1450       Py_DECREF(pyresult);
1451     }
1452     *size = result;
1453     //MESSAGE("f(" << uv[0] << "," << uv[1] << ")" << " = " << result);
1454     PyGILState_Release(gstate);
1455   }
1456   else {
1457     *size = *((double*)user_data);
1458   }
1459   return STATUS_OK;
1460 }
1461
1462 status_t size_on_edge(integer edge_id, real t, real *size, void *user_data)
1463 {
1464   if (EdgeId2PythonSmp.count(edge_id) != 0){
1465     PyObject * pyresult = NULL;
1466     PyObject* new_stderr = NULL;
1467     assert(Py_IsInitialized());
1468     PyGILState_STATE gstate;
1469     gstate = PyGILState_Ensure();
1470     pyresult = PyObject_CallFunction(EdgeId2PythonSmp[edge_id],"(f)",t);
1471     double result;
1472     if ( pyresult == NULL){
1473       fflush(stderr);
1474       string err_description="";
1475       new_stderr = newPyStdOut(err_description);
1476       PySys_SetObject("stderr", new_stderr);
1477       PyErr_Print();
1478       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1479       Py_DECREF(new_stderr);
1480       MESSAGE("Can't evaluate f(" << t << ")" << " error is " << err_description);
1481       result = *((double*)user_data);
1482       }
1483     else {
1484       result = PyFloat_AsDouble(pyresult);
1485       Py_DECREF(pyresult);
1486     }
1487     *size = result;
1488     PyGILState_Release(gstate);
1489   }
1490   else {
1491     *size = *((double*)user_data);
1492   }
1493   return STATUS_OK;
1494 }
1495
1496 status_t size_on_vertex(integer point_id, real *size, void *user_data)
1497 {
1498   if (VertexId2PythonSmp.count(point_id) != 0){
1499     PyObject * pyresult = NULL;
1500     PyObject* new_stderr = NULL;
1501     assert(Py_IsInitialized());
1502     PyGILState_STATE gstate;
1503     gstate = PyGILState_Ensure();
1504     pyresult = PyObject_CallFunction(VertexId2PythonSmp[point_id],"");
1505     double result;
1506     if ( pyresult == NULL){
1507       fflush(stderr);
1508       string err_description="";
1509       new_stderr = newPyStdOut(err_description);
1510       PySys_SetObject("stderr", new_stderr);
1511       PyErr_Print();
1512       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1513       Py_DECREF(new_stderr);
1514       MESSAGE("Can't evaluate f()" << " error is " << err_description);
1515       result = *((double*)user_data);
1516       }
1517     else {
1518       result = PyFloat_AsDouble(pyresult);
1519       Py_DECREF(pyresult);
1520     }
1521     *size = result;
1522     PyGILState_Release(gstate);
1523   }
1524   else {
1525     *size = *((double*)user_data);
1526   }
1527  return STATUS_OK;
1528 }
1529
1530 status_t message_callback(message_t *msg, void *user_data)
1531 {
1532   integer errnumber = 0;
1533   char *desc;
1534   message_get_number(msg, &errnumber);
1535   message_get_description(msg, &desc);
1536   if ( errnumber < 0 ) {
1537     string * error = (string*)user_data;
1538 //   if ( !error->empty() )
1539 //     *error += "\n";
1540     // remove ^A from the tail
1541     int len = strlen( desc );
1542     while (len > 0 && desc[len-1] != '\n')
1543       len--;
1544     error->append( desc, len );
1545   }
1546   else {
1547       std::cout << desc << std::endl;
1548   }
1549   return STATUS_OK;
1550 }
1551
1552
1553 //=============================================================================
1554 /*!
1555  *
1556  */
1557 //=============================================================================
1558 bool BLSURFPlugin_BLSURF::Evaluate(SMESH_Mesh& aMesh,
1559                                    const TopoDS_Shape& aShape,
1560                                    MapShapeNbElems& aResMap)
1561 {
1562   int    _physicalMesh  = BLSURFPlugin_Hypothesis::GetDefaultPhysicalMesh();
1563   double _phySize       = BLSURFPlugin_Hypothesis::GetDefaultPhySize();
1564   //int    _geometricMesh = BLSURFPlugin_Hypothesis::GetDefaultGeometricMesh();
1565   //double _angleMeshS    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshS();
1566   double _angleMeshC    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshC();
1567   if(_hypothesis) {
1568     _physicalMesh  = (int) _hypothesis->GetPhysicalMesh();
1569     _phySize       = _hypothesis->GetPhySize();
1570     //_geometricMesh = (int) hyp->GetGeometricMesh();
1571     //_angleMeshS    = hyp->GetAngleMeshS();
1572     _angleMeshC    = _hypothesis->GetAngleMeshC();
1573   }
1574
1575   bool IsQuadratic = false;
1576
1577   // ----------------
1578   // evaluate 1D
1579   // ----------------
1580   TopTools_DataMapOfShapeInteger EdgesMap;
1581   double fullLen = 0.0;
1582   double fullNbSeg = 0;
1583   for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
1584     TopoDS_Edge E = TopoDS::Edge( exp.Current() );
1585     if( EdgesMap.IsBound(E) )
1586       continue;
1587     SMESH_subMesh *sm = aMesh.GetSubMesh(E);
1588     double aLen = SMESH_Algo::EdgeLength(E);
1589     fullLen += aLen;
1590     int nb1d = 0;
1591     if(_physicalMesh==1) {
1592        nb1d = (int)( aLen/_phySize + 1 );
1593     }
1594     else {
1595       // use geometry
1596       double f,l;
1597       Handle(Geom_Curve) C = BRep_Tool::Curve(E,f,l);
1598       double fullAng = 0.0;
1599       double dp = (l-f)/200;
1600       gp_Pnt P1,P2,P3;
1601       C->D0(f,P1);
1602       C->D0(f+dp,P2);
1603       gp_Vec V1(P1,P2);
1604       for(int j=2; j<=200; j++) {
1605         C->D0(f+dp*j,P3);
1606         gp_Vec V2(P2,P3);
1607         fullAng += fabs(V1.Angle(V2));
1608         V1 = V2;
1609         P2 = P3;
1610       }
1611       nb1d = (int)( fullAng/_angleMeshC + 1 );
1612     }
1613     fullNbSeg += nb1d;
1614     std::vector<int> aVec(SMDSEntity_Last);
1615     for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
1616     if( IsQuadratic > 0 ) {
1617       aVec[SMDSEntity_Node] = 2*nb1d - 1;
1618       aVec[SMDSEntity_Quad_Edge] = nb1d;
1619     }
1620     else {
1621       aVec[SMDSEntity_Node] = nb1d - 1;
1622       aVec[SMDSEntity_Edge] = nb1d;
1623     }
1624     aResMap.insert(std::make_pair(sm,aVec));
1625     EdgesMap.Bind(E,nb1d);
1626   }
1627   double ELen = fullLen/fullNbSeg;
1628   // ----------------
1629   // evaluate 2D
1630   // ----------------
1631   // try to evaluate as in MEFISTO
1632   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
1633     TopoDS_Face F = TopoDS::Face( exp.Current() );
1634     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
1635     GProp_GProps G;
1636     BRepGProp::SurfaceProperties(F,G);
1637     double anArea = G.Mass();
1638     int nb1d = 0;
1639     for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next()) {
1640       nb1d += EdgesMap.Find(exp1.Current());
1641     }
1642     int nbFaces = (int) ( anArea/( ELen*ELen*sqrt(3.) / 4 ) );
1643     int nbNodes = (int) ( ( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
1644     std::vector<int> aVec(SMDSEntity_Last);
1645     for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
1646     if( IsQuadratic ) {
1647       int nb1d_in = (nbFaces*3 - nb1d) / 2;
1648       aVec[SMDSEntity_Node] = nbNodes + nb1d_in;
1649       aVec[SMDSEntity_Quad_Triangle] = nbFaces;
1650     }
1651     else {
1652       aVec[SMDSEntity_Node] = nbNodes;
1653       aVec[SMDSEntity_Triangle] = nbFaces;
1654     }
1655     aResMap.insert(std::make_pair(sm,aVec));
1656   }
1657
1658   // ----------------
1659   // evaluate 3D
1660   // ----------------
1661   GProp_GProps G;
1662   BRepGProp::VolumeProperties(aShape,G);
1663   double aVolume = G.Mass();
1664   double tetrVol = 0.1179*ELen*ELen*ELen;
1665   int nbVols = (int)aVolume/tetrVol;
1666   int nb1d_in = (int) ( ( nbVols*6 - fullNbSeg ) / 6 );
1667   std::vector<int> aVec(SMDSEntity_Last);
1668   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
1669   if( IsQuadratic ) {
1670     aVec[SMDSEntity_Node] = nb1d_in/3 + 1 + nb1d_in;
1671     aVec[SMDSEntity_Quad_Tetra] = nbVols;
1672   }
1673   else {
1674     aVec[SMDSEntity_Node] = nb1d_in/3 + 1;
1675     aVec[SMDSEntity_Tetra] = nbVols;
1676   }
1677   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
1678   aResMap.insert(std::make_pair(sm,aVec));
1679
1680   return true;
1681 }
1682
1683 //=============================================================================
1684 /*!
1685  *  Rewritting of the BRepClass_FaceClassifier::Perform function which is bugged (CAS 6.3sp6)
1686  *  Following line was added:
1687  *        myExtrem.Perform(P);
1688  */
1689 //=============================================================================
1690 void  BLSURFPlugin_BLSURF::BRepClass_FaceClassifierPerform(BRepClass_FaceClassifier* fc,
1691                     const TopoDS_Face& face,
1692                     const gp_Pnt& P,
1693                     const Standard_Real Tol)
1694 {
1695   //-- Voir BRepExtrema_ExtPF.cxx
1696   BRepAdaptor_Surface Surf(face);
1697   Standard_Real U1, U2, V1, V2;
1698   BRepTools::UVBounds(face, U1, U2, V1, V2);
1699   Extrema_ExtPS myExtrem;
1700   myExtrem.Initialize(Surf, U1, U2, V1, V2, Tol, Tol);
1701   myExtrem.Perform(P);
1702   //----------------------------------------------------------
1703   //-- On cherche le point le plus proche , PUIS
1704   //-- On le classifie.
1705   Standard_Integer nbv    = 0; // xpu
1706   Standard_Real MaxDist   =  RealLast();
1707   Standard_Integer indice = 0;
1708   if(myExtrem.IsDone()) {
1709     nbv = myExtrem.NbExt();
1710     for (Standard_Integer i = 1; i <= nbv; i++) {
1711       Standard_Real d = myExtrem.Value(i);
1712       d = Abs(d);
1713       if(d <= MaxDist) {
1714     MaxDist = d;
1715     indice = i;
1716       }
1717     }
1718   }
1719   if(indice) {
1720     gp_Pnt2d Puv;
1721     Standard_Real U1,U2;
1722     myExtrem.Point(indice).Parameter(U1, U2);
1723     Puv.SetCoord(U1, U2);
1724     fc->Perform(face, Puv, Tol);
1725   }
1726   else {
1727     fc->Perform(face, gp_Pnt2d(U1-1.0,V1 - 1.0), Tol); //-- NYI etc BUG PAS BEAU En attendant l acces a rejected
1728     //-- le resultat est TopAbs_OUT;
1729   }
1730 }
1731