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     TopoDS_Shape S = TopoDS_Shape();
 392     SALOMEDS::SObject_var aSO = myStudy->FindObjectID(entry.c_str());
 393     SALOMEDS::GenericAttribute_var anAttr;
 394     if (!aSO->_is_nil()){
 395       SALOMEDS::SObject_var aRefSObj;
 396       GEOM::GEOM_Object_var aShape;
 397       SALOMEDS::AttributeIOR_var myAttribute;
 398       CORBA::String_var myAttrValue;
 399       CORBA::Object_var myCorbaObj;
 400       // If selected object is a reference
 401       if ( aSO->ReferencedObject( aRefSObj ))
 402         aSO = aRefSObj;
 403       SALOMEDS::SComponent_var myFatherCpnt = aSO->GetFatherComponent();
 404       CORBA::String_var myFatherCpntDataType = myFatherCpnt->ComponentDataType();
 405       if (  strcmp(myFatherCpntDataType,"GEOM")==0) {
 406         MESSAGE("aSO father component is GEOM");
 407         if (!aSO->FindAttribute(anAttr, "AttributeIOR")) return S;
 408         myAttribute=SALOMEDS::AttributeIOR::_narrow(anAttr);
 409         myAttrValue=myAttribute->Value();
 410         MESSAGE("aSO IOR: "<< myAttrValue);
 411         myCorbaObj=smeshGen_i->GetORB()->string_to_object(myAttrValue);
 412         aShape = GEOM::GEOM_Object::_narrow(myCorbaObj);
 413       }
 414       if ( !aShape->_is_nil() )
 415         S=smeshGen_i->GeomObjectToShape( aShape.in() );
 416     }
 417     return S;
 418 }
 419
 420 /////////////////////////////////////////////////////////
 421 TopoDS_Shape BLSURFPlugin_BLSURF::entryToShape(std::string entry)
 422 {
 423   MESSAGE("BLSURFPlugin_BLSURF::entryToShape"<<entry );
 424   GEOM::GEOM_Object_var aGeomObj;
 425   TopoDS_Shape S = TopoDS_Shape();
 426   SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() );
 427   SALOMEDS::GenericAttribute_var anAttr;
 428
 429   if (!aSObj->_is_nil() && aSObj->FindAttribute(anAttr, "AttributeIOR")) {
 430     SALOMEDS::AttributeIOR_var anIOR = SALOMEDS::AttributeIOR::_narrow(anAttr);
 431     CORBA::String_var aVal = anIOR->Value();
 432     CORBA::Object_var obj = myStudy->ConvertIORToObject(aVal);
 433     aGeomObj = GEOM::GEOM_Object::_narrow(obj);
 434   }
 435   if ( !aGeomObj->_is_nil() )
 436     S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
 437   return S;
 438 }
 439
 440 /////////////////////////////////////////////////////////
 441 void createEnforcedVertexOnFace(TopoDS_Shape GeomShape, BLSURFPlugin_Hypothesis::TEnforcedVertexList enforcedVertexList)
 442 {
 443   double xe, ye, ze;
 444   std::vector<double> coords;
 445   BLSURFPlugin_Hypothesis::TEnforcedVertex enforcedVertex;
 446   BLSURFPlugin_Hypothesis::TEnforcedVertexList::const_iterator evlIt = enforcedVertexList.begin();
 447
 448   for( ; evlIt != enforcedVertexList.end() ; ++evlIt ) {
 449     coords.clear();
 450     enforcedVertex = *evlIt;
 451     xe = enforcedVertex[0];
 452     ye = enforcedVertex[1];
 453     ze = enforcedVertex[2];
 454     MESSAGE("Enforced Vertex: " << xe << ", " << ye << ", " << ze);
 455     // Get the (u,v) values of the enforced vertex on the face
 456     projectionPoint myPoint = getProjectionPoint(TopoDS::Face(GeomShape),gp_XYZ(xe,ye,ze));
 457     gp_XY uvPoint = myPoint.uv;
 458     gp_XYZ xyzPoint = myPoint.xyz;
 459     Standard_Real u0 = uvPoint.X();
 460     Standard_Real v0 = uvPoint.Y();
 461     Standard_Real x0 = xyzPoint.X();
 462     Standard_Real y0 = xyzPoint.Y();
 463     Standard_Real z0 = xyzPoint.Z();
 464     MESSAGE("Projected Vertex: " << x0 << ", " << y0 << ", " << z0);
 465     coords.push_back(u0);
 466     coords.push_back(v0);
 467     coords.push_back(x0);
 468     coords.push_back(y0);
 469     coords.push_back(z0);
 470
 471     int key = 0;
 472     if (! FacesWithEnforcedVertices.Contains(TopoDS::Face(GeomShape))) {
 473       key = FacesWithEnforcedVertices.Add(TopoDS::Face(GeomShape));
 474     }
 475     else {
 476       key = FacesWithEnforcedVertices.FindIndex(TopoDS::Face(GeomShape));
 477     }
 478
 479     // If a node is already created by an attractor, do not create enforced vertex
 480     int attractorKey = FacesWithSizeMap.FindIndex(TopoDS::Face(GeomShape));
 481     bool sameAttractor = false;
 482     if (attractorKey >= 0)
 483       if (FaceId2AttractorCoords.count(attractorKey) > 0)
 484         if (FaceId2AttractorCoords[attractorKey] == coords)
 485           sameAttractor = true;
 486
 487     if (FaceId2EnforcedVertexCoords.find(key) != FaceId2EnforcedVertexCoords.end()) {
 488       MESSAGE("Map of enf. vertex has key " << key)
 489       MESSAGE("Enf. vertex list size is: " << FaceId2EnforcedVertexCoords[key].size())
 490       if (not sameAttractor)
 491         FaceId2EnforcedVertexCoords[key].insert(coords); // there should be no redondant coords here (see std::set management)
 492       else
 493         MESSAGE("An attractor node is already defined: I don't add the enforced vertex");
 494       MESSAGE("New Enf. vertex list size is: " << FaceId2EnforcedVertexCoords[key].size())
 495     }
 496     else {
 497       MESSAGE("Map of enf. vertex has not key " << key << ": creating it")
 498       if (not sameAttractor) {
 499         std::set< std::vector<double> > ens;
 500         ens.insert(coords);
 501         FaceId2EnforcedVertexCoords[key] = ens;
 502       }
 503       else
 504         MESSAGE("An attractor node is already defined: I don't add the enforced vertex");
 505     }
 506   }
 507 }
 508
 509 /////////////////////////////////////////////////////////
 510 void createAttractorOnFace(TopoDS_Shape GeomShape, std::string AttractorFunction)
 511 {
 512   MESSAGE("Attractor function: "<< AttractorFunction);
 513   double xa, ya, za; // Coordinates of attractor point
 514   double a, b;       // Attractor parameter
 515   bool createNode=false; // To create a node on attractor projection
 516   int pos1, pos2;
 517   const char *sep = ";";
 518   // atIt->second has the following pattern:
 519   // ATTRACTOR(xa;ya;za;a;b)
 520   // where:
 521   // xa;ya;za : coordinates of  attractor
 522   // a        : desired size on attractor
 523   // b        : distance of influence of attractor
 524   //
 525   // We search the parameters in the string
 526   // xa
 527   pos1 = AttractorFunction.find(sep);
 528   if (pos1!=string::npos)
 529   xa = atof(AttractorFunction.substr(10, pos1-10).c_str());
 530   // ya
 531   pos2 = AttractorFunction.find(sep, pos1+1);
 532   if (pos2!=string::npos) {
 533   ya = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
 534   pos1 = pos2;
 535     }
 536   // za
 537   pos2 = AttractorFunction.find(sep, pos1+1);
 538   if (pos2!=string::npos) {
 539   za = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
 540   pos1 = pos2;
 541   }
 542   // a
 543   pos2 = AttractorFunction.find(sep, pos1+1);
 544   if (pos2!=string::npos) {
 545   a = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
 546   pos1 = pos2;
 547   }
 548   // b
 549   pos2 = AttractorFunction.find(sep, pos1+1);
 550   if (pos2!=string::npos) {
 551   b = atof(AttractorFunction.substr(pos1+1, pos2-pos1-1).c_str());
 552     pos1 = pos2;
 553   }
 554   // createNode
 555   pos2 = AttractorFunction.find(")");
 556   if (pos2!=string::npos) {
 557     string createNodeStr = AttractorFunction.substr(pos1+1, pos2-pos1-1);
 558     MESSAGE("createNode: " << createNodeStr);
 559     createNode = (AttractorFunction.substr(pos1+1, pos2-pos1-1) == "True");
 560   }
 561
 562   // Get the (u,v) values of the attractor on the face
 563   projectionPoint myPoint = getProjectionPoint(TopoDS::Face(GeomShape),gp_XYZ(xa,ya,za));
 564   gp_XY uvPoint = myPoint.uv;
 565   gp_XYZ xyzPoint = myPoint.xyz;
 566   Standard_Real u0 = uvPoint.X();
 567   Standard_Real v0 = uvPoint.Y();
 568   Standard_Real x0 = xyzPoint.X();
 569   Standard_Real y0 = xyzPoint.Y();
 570   Standard_Real z0 = xyzPoint.Z();
 571   std::vector<double> coords;
 572   coords.push_back(u0);
 573   coords.push_back(v0);
 574   coords.push_back(x0);
 575   coords.push_back(y0);
 576   coords.push_back(z0);
 577   // We construct the python function
 578   ostringstream attractorFunctionStream;
 579   attractorFunctionStream << "def f(u,v): return ";
 580   attractorFunctionStream << _smp_phy_size << "-(" << _smp_phy_size <<"-" << a << ")";
 581   attractorFunctionStream << "*exp(-((u-("<<u0<<"))*(u-("<<u0<<"))+(v-("<<v0<<"))*(v-("<<v0<<")))/(" << b << "*" << b <<"))";
 582
 583   MESSAGE("Python function for attractor:" << std::endl << attractorFunctionStream.str());
 584
 585   int key;
 586   if (! FacesWithSizeMap.Contains(TopoDS::Face(GeomShape))) {
 587     key = FacesWithSizeMap.Add(TopoDS::Face(GeomShape));
 588   }
 589   else {
 590     key = FacesWithSizeMap.FindIndex(TopoDS::Face(GeomShape));
 591   }
 592   FaceId2SizeMap[key] =attractorFunctionStream.str();
 593   if (createNode) {
 594     MESSAGE("Creating node on ("<<x0<<","<<y0<<","<<z0<<")");
 595     FaceId2AttractorCoords[key] = coords;
 596   }
 597 }
 598
 599 /////////////////////////////////////////////////////////
 600
 601 void BLSURFPlugin_BLSURF::SetParameters(const BLSURFPlugin_Hypothesis* hyp, blsurf_session_t *bls)
 602 {
 603   int    _topology      = BLSURFPlugin_Hypothesis::GetDefaultTopology();
 604   int    _physicalMesh  = BLSURFPlugin_Hypothesis::GetDefaultPhysicalMesh();
 605   double _phySize       = BLSURFPlugin_Hypothesis::GetDefaultPhySize();
 606   int    _geometricMesh = BLSURFPlugin_Hypothesis::GetDefaultGeometricMesh();
 607   double _angleMeshS    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshS();
 608   double _angleMeshC    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshC();
 609   double _gradation     = BLSURFPlugin_Hypothesis::GetDefaultGradation();
 610   bool   _quadAllowed   = BLSURFPlugin_Hypothesis::GetDefaultQuadAllowed();
 611   bool   _decimesh      = BLSURFPlugin_Hypothesis::GetDefaultDecimesh();
 612   int    _verb          = BLSURFPlugin_Hypothesis::GetDefaultVerbosity();
 613
 614   if (hyp) {
 615     MESSAGE("BLSURFPlugin_BLSURF::SetParameters");
 616     _topology      = (int) hyp->GetTopology();
 617     _physicalMesh  = (int) hyp->GetPhysicalMesh();
 618     _phySize       = hyp->GetPhySize();
 619     _geometricMesh = (int) hyp->GetGeometricMesh();
 620     _angleMeshS    = hyp->GetAngleMeshS();
 621     _angleMeshC    = hyp->GetAngleMeshC();
 622     _gradation     = hyp->GetGradation();
 623     _quadAllowed   = hyp->GetQuadAllowed();
 624     _decimesh      = hyp->GetDecimesh();
 625     _verb          = hyp->GetVerbosity();
 626
 627     if ( hyp->GetPhyMin() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 628       blsurf_set_param(bls, "hphymin", to_string(hyp->GetPhyMin()).c_str());
 629     if ( hyp->GetPhyMax() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 630       blsurf_set_param(bls, "hphymax", to_string(hyp->GetPhyMax()).c_str());
 631     if ( hyp->GetGeoMin() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 632       blsurf_set_param(bls, "hgeomin", to_string(hyp->GetGeoMin()).c_str());
 633     if ( hyp->GetGeoMax() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 634       blsurf_set_param(bls, "hgeomax", to_string(hyp->GetGeoMax()).c_str());
 635
 636     const BLSURFPlugin_Hypothesis::TOptionValues & opts = hyp->GetOptionValues();
 637     BLSURFPlugin_Hypothesis::TOptionValues::const_iterator opIt;
 638     for ( opIt = opts.begin(); opIt != opts.end(); ++opIt )
 639       if ( !opIt->second.empty() ) {
 640         MESSAGE("blsurf_set_param(): " << opIt->first << " = " << opIt->second);
 641         blsurf_set_param(bls, opIt->first.c_str(), opIt->second.c_str());
 642       }
 643
 644   } else {
 645     MESSAGE("BLSURFPlugin_BLSURF::SetParameters using defaults");
 646   }
 647   _smp_phy_size = _phySize;
 648   blsurf_set_param(bls, "topo_points",       _topology > 0 ? "1" : "0");
 649   blsurf_set_param(bls, "topo_curves",       _topology > 0 ? "1" : "0");
 650   blsurf_set_param(bls, "topo_project",      _topology > 0 ? "1" : "0");
 651   blsurf_set_param(bls, "clean_boundary",    _topology > 1 ? "1" : "0");
 652   blsurf_set_param(bls, "close_boundary",    _topology > 1 ? "1" : "0");
 653   blsurf_set_param(bls, "hphy_flag",         to_string(_physicalMesh).c_str());
 654 //  blsurf_set_param(bls, "hphy_flag",         "2");
 655   if ((to_string(_physicalMesh))=="2"){
 656     TopoDS_Shape GeomShape;
 657     TopAbs_ShapeEnum GeomType;
 658     //
 659     // Standard Size Maps
 660     //
 661     MESSAGE("Setting a Size Map");
 662     const BLSURFPlugin_Hypothesis::TSizeMap & sizeMaps = hyp->GetSizeMapEntries();
 663     BLSURFPlugin_Hypothesis::TSizeMap::const_iterator smIt;
 664     for ( smIt = sizeMaps.begin(); smIt != sizeMaps.end(); ++smIt ) {
 665       if ( !smIt->second.empty() ) {
 666         MESSAGE("blsurf_set_sizeMap(): " << smIt->first << " = " << smIt->second);
 667         GeomShape = entryToShape(smIt->first);
 668         GeomType  = GeomShape.ShapeType();
 669         MESSAGE("Geomtype is " << GeomType);
 670         int key = -1;
 671         // Group Management
 672         if (GeomType == TopAbs_COMPOUND){
 673           for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
 674             // Group of faces
 675             if (it.Value().ShapeType() == TopAbs_FACE){
 676               HasSizeMapOnFace = true;
 677               if (! FacesWithSizeMap.Contains(TopoDS::Face(it.Value()))) {
 678                 key = FacesWithSizeMap.Add(TopoDS::Face(it.Value()));
 679               }
 680               else {
 681                 key = FacesWithSizeMap.FindIndex(TopoDS::Face(it.Value()));
 682 //                 MESSAGE("Face with key " << key << " already in map");
 683               }
 684               FaceId2SizeMap[key] = smIt->second;
 685             }
 686             // Group of edges
 687             if (it.Value().ShapeType() == TopAbs_EDGE){
 688               HasSizeMapOnEdge = true;
 689               HasSizeMapOnFace = true;
 690               if (! EdgesWithSizeMap.Contains(TopoDS::Edge(it.Value()))) {
 691                 key = EdgesWithSizeMap.Add(TopoDS::Edge(it.Value()));
 692               }
 693               else {
 694                 key = EdgesWithSizeMap.FindIndex(TopoDS::Edge(it.Value()));
 695 //                 MESSAGE("Edge with key " << key << " already in map");
 696               }
 697               EdgeId2SizeMap[key] = smIt->second;
 698             }
 699             // Group of vertices
 700             if (it.Value().ShapeType() == TopAbs_VERTEX){
 701               HasSizeMapOnVertex = true;
 702               HasSizeMapOnEdge = true;
 703               HasSizeMapOnFace = true;
 704               if (! VerticesWithSizeMap.Contains(TopoDS::Vertex(it.Value()))) {
 705                 key = VerticesWithSizeMap.Add(TopoDS::Vertex(it.Value()));
 706               }
 707               else {
 708                 key = VerticesWithSizeMap.FindIndex(TopoDS::Vertex(it.Value()));
 709 //                 MESSAGE("Vertex with key " << key << " already in map");
 710               }
 711               VertexId2SizeMap[key] = smIt->second;
 712             }
 713           }
 714         }
 715         // Single face
 716         if (GeomType == TopAbs_FACE){
 717           HasSizeMapOnFace = true;
 718           if (! FacesWithSizeMap.Contains(TopoDS::Face(GeomShape))) {
 719             key = FacesWithSizeMap.Add(TopoDS::Face(GeomShape));
 720           }
 721           else {
 722             key = FacesWithSizeMap.FindIndex(TopoDS::Face(GeomShape));
 723 //             MESSAGE("Face with key " << key << " already in map");
 724           }
 725           FaceId2SizeMap[key] = smIt->second;
 726         }
 727         // Single edge
 728         if (GeomType == TopAbs_EDGE){
 729           HasSizeMapOnEdge = true;
 730           HasSizeMapOnFace = true;
 731           if (! EdgesWithSizeMap.Contains(TopoDS::Edge(GeomShape))) {
 732             key = EdgesWithSizeMap.Add(TopoDS::Edge(GeomShape));
 733           }
 734           else {
 735             key = EdgesWithSizeMap.FindIndex(TopoDS::Edge(GeomShape));
 736 //             MESSAGE("Edge with key " << key << " already in map");
 737           }
 738           EdgeId2SizeMap[key] = smIt->second;
 739         }
 740         // Single vertex
 741         if (GeomType == TopAbs_VERTEX){
 742           HasSizeMapOnVertex = true;
 743           HasSizeMapOnEdge   = true;
 744           HasSizeMapOnFace   = true;
 745           if (! VerticesWithSizeMap.Contains(TopoDS::Vertex(GeomShape))) {
 746             key = VerticesWithSizeMap.Add(TopoDS::Vertex(GeomShape));
 747           }
 748           else {
 749             key = VerticesWithSizeMap.FindIndex(TopoDS::Vertex(GeomShape));
 750 //             MESSAGE("Vertex with key " << key << " already in map");
 751           }
 752           VertexId2SizeMap[key] = smIt->second;
 753         }
 754       }
 755     }
 756
 757     //
 758     // Attractors
 759     //
 760     MESSAGE("Setting Attractors");
 761     const BLSURFPlugin_Hypothesis::TSizeMap attractors = BLSURFPlugin_Hypothesis::GetAttractorEntries(hyp);
 762     BLSURFPlugin_Hypothesis::TSizeMap::const_iterator atIt = attractors.begin();
 763     for ( ; atIt != attractors.end(); ++atIt ) {
 764       if ( !atIt->second.empty() ) {
 765         MESSAGE("blsurf_set_attractor(): " << atIt->first << " = " << atIt->second);
 766         GeomShape = entryToShape(atIt->first);
 767         GeomType  = GeomShape.ShapeType();
 768         // Group Management
 769         if (GeomType == TopAbs_COMPOUND){
 770           for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
 771             if (it.Value().ShapeType() == TopAbs_FACE){
 772               HasSizeMapOnFace = true;
 773               createAttractorOnFace(it.Value(), atIt->second);
 774             }
 775           }
 776         }
 777
 778         if (GeomType == TopAbs_FACE){
 779           HasSizeMapOnFace = true;
 780           createAttractorOnFace(GeomShape, atIt->second);
 781         }
 782 /*
 783         if (GeomType == TopAbs_EDGE){
 784           HasSizeMapOnEdge = true;
 785           HasSizeMapOnFace = true;
 786         EdgeId2SizeMap[TopoDS::Edge(GeomShape).HashCode(IntegerLast())] = atIt->second;
 787         }
 788         if (GeomType == TopAbs_VERTEX){
 789           HasSizeMapOnVertex = true;
 790           HasSizeMapOnEdge   = true;
 791           HasSizeMapOnFace   = true;
 792         VertexId2SizeMap[TopoDS::Vertex(GeomShape).HashCode(IntegerLast())] = atIt->second;
 793         }
 794 */
 795       }
 796     }
 797
 798
 799     //
 800     // Enforced Vertices
 801     //
 802     MESSAGE("Setting Enforced Vertices");
 803     const BLSURFPlugin_Hypothesis::TEnforcedVertexMap enforcedVertexMap = BLSURFPlugin_Hypothesis::GetAllEnforcedVertices(hyp);
 804     BLSURFPlugin_Hypothesis::TEnforcedVertexMap::const_iterator enfIt = enforcedVertexMap.begin();
 805     for ( ; enfIt != enforcedVertexMap.end(); ++enfIt ) {
 806       if ( !enfIt->second.empty() ) {
 807         GeomShape = entryToShape(enfIt->first);
 808         GeomType  = GeomShape.ShapeType();
 809         // Group Management
 810         if (GeomType == TopAbs_COMPOUND){
 811           for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
 812             if (it.Value().ShapeType() == TopAbs_FACE){
 813               HasSizeMapOnFace = true;
 814               createEnforcedVertexOnFace(it.Value(), enfIt->second);
 815             }
 816           }
 817         }
 818
 819         if (GeomType == TopAbs_FACE){
 820           HasSizeMapOnFace = true;
 821           createEnforcedVertexOnFace(GeomShape, enfIt->second);
 822         }
 823       }
 824     }
 825
 826 //    if (HasSizeMapOnFace){
 827     // In all size map cases (hphy_flag = 2), at least map on face must be defined
 828     MESSAGE("Setting Size Map on FACES ");
 829     blsurf_data_set_sizemap_iso_cad_face(bls, size_on_surface, &_smp_phy_size);
 830 //    }
 831
 832     if (HasSizeMapOnEdge){
 833       MESSAGE("Setting Size Map on EDGES ");
 834       blsurf_data_set_sizemap_iso_cad_edge(bls, size_on_edge, &_smp_phy_size);
 835     }
 836     if (HasSizeMapOnVertex){
 837       MESSAGE("Setting Size Map on VERTICES ");
 838       blsurf_data_set_sizemap_iso_cad_point(bls, size_on_vertex, &_smp_phy_size);
 839     }
 840   }
 841   blsurf_set_param(bls, "hphydef",           to_string(_phySize).c_str());
 842   blsurf_set_param(bls, "hgeo_flag",         to_string(_geometricMesh).c_str());
 843   blsurf_set_param(bls, "relax_size",        _decimesh ? "0": to_string(_geometricMesh).c_str());
 844   blsurf_set_param(bls, "angle_meshs",       to_string(_angleMeshS).c_str());
 845   blsurf_set_param(bls, "angle_meshc",       to_string(_angleMeshC).c_str());
 846   blsurf_set_param(bls, "gradation",         to_string(_gradation).c_str());
 847   blsurf_set_param(bls, "patch_independent", _decimesh ? "1" : "0");
 848   blsurf_set_param(bls, "element",           _quadAllowed ? "q1.0" : "p1");
 849   blsurf_set_param(bls, "verb",              to_string(_verb).c_str());
 850 }
 851
 852 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data);
 853 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
 854                   real *duu, real *duv, real *dvv, void *user_data);
 855 status_t message_callback(message_t *msg, void *user_data);
 856
 857 //=============================================================================
 858 /*!
 859  *
 860  */
 861 //=============================================================================
 862
 863 bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) {
 864
 865   MESSAGE("BLSURFPlugin_BLSURF::Compute");
 866
 867   if (aShape.ShapeType() == TopAbs_COMPOUND) {
 868     MESSAGE("  the shape is a COMPOUND");
 869   }
 870   else {
 871     MESSAGE("  the shape is UNKNOWN");
 872   };
 873
 874   context_t *ctx =  context_new();
 875   context_set_message_callback(ctx, message_callback, &_comment);
 876
 877   cad_t *c = cad_new(ctx);
 878
 879   blsurf_session_t *bls = blsurf_session_new(ctx);
 880
 881   FacesWithSizeMap.Clear();
 882   FaceId2SizeMap.clear();
 883   EdgesWithSizeMap.Clear();
 884   EdgeId2SizeMap.clear();
 885   VerticesWithSizeMap.Clear();
 886   VertexId2SizeMap.clear();
 887
 888   MESSAGE("BEGIN SetParameters");
 889   SetParameters(_hypothesis, bls);
 890   MESSAGE("END SetParameters");
 891
 892   TopTools_IndexedMapOfShape fmap;
 893   TopTools_IndexedMapOfShape emap;
 894   TopTools_IndexedMapOfShape pmap;
 895   vector<Handle(Geom2d_Curve)> curves;
 896   vector<Handle(Geom_Surface)> surfaces;
 897
 898   fmap.Clear();
 899   FaceId2PythonSmp.clear();
 900   emap.Clear();
 901   EdgeId2PythonSmp.clear();
 902   pmap.Clear();
 903   VertexId2PythonSmp.clear();
 904   surfaces.resize(0);
 905   curves.resize(0);
 906
 907   assert(Py_IsInitialized());
 908   PyGILState_STATE gstate;
 909   gstate = PyGILState_Ensure();
 910
 911   string theSizeMapStr;
 912
 913   /****************************************************************************************
 914                                   FACES
 915   *****************************************************************************************/
 916   int iface = 0;
 917   string bad_end = "return";
 918   int faceKey = -1;
 919   for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next()) {
 920     TopoDS_Face f=TopoDS::Face(face_iter.Current());
 921
 922     if (fmap.FindIndex(f) > 0)
 923       continue;
 924
 925     fmap.Add(f);
 926     iface++;
 927     surfaces.push_back(BRep_Tool::Surface(f));
 928
 929     cad_face_t *fce = cad_face_new(c, iface, surf_fun, surfaces.back());
 930     cad_face_set_tag(fce, iface);
 931     if(f.Orientation() != TopAbs_FORWARD){
 932       cad_face_set_orientation(fce, CAD_ORIENTATION_REVERSED);
 933     } else {
 934       cad_face_set_orientation(fce, CAD_ORIENTATION_FORWARD);
 935     }
 936
 937     if (HasSizeMapOnFace){
 938       std::cout << "A size map is defined on a face" << std::endl;
 939       // Classic size map
 940       faceKey = FacesWithSizeMap.FindIndex(f);
 941
 942       if (FaceId2SizeMap.find(faceKey)!=FaceId2SizeMap.end()){
 943         theSizeMapStr = FaceId2SizeMap[faceKey];
 944         // check if function ends with "return"
 945         if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
 946           continue;
 947         // Expr To Python function, verification is performed at validation in GUI
 948         PyObject * obj = NULL;
 949         obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
 950         Py_DECREF(obj);
 951         PyObject * func = NULL;
 952         func = PyObject_GetAttrString(main_mod, "f");
 953         FaceId2PythonSmp[iface]=func;
 954         FaceId2SizeMap.erase(faceKey);
 955       }
 956
 957       // Specific size map = Attractor
 958       std::map<int,std::vector<double> >::iterator attractor_iter = FaceId2AttractorCoords.begin();
 959       int iatt=0;
 960       for (; attractor_iter != FaceId2AttractorCoords.end(); ++attractor_iter) {
 961         if (attractor_iter->first == faceKey) {
 962           MESSAGE("Face indice: " << iface);
 963           MESSAGE("Adding attractor");
 964
 965           double xyzCoords[3]  = {attractor_iter->second[2],
 966                                   attractor_iter->second[3],
 967                                   attractor_iter->second[4]};
 968
 969           MESSAGE("Check position of vertex =(" << xyzCoords[0] << "," << xyzCoords[1] << "," << xyzCoords[2] << ")");
 970           gp_Pnt P(xyzCoords[0],xyzCoords[1],xyzCoords[2]);
 971           BRepClass_FaceClassifier scl(f,P,1e-7);
 972           // scl.Perform() is bugged. The function was rewritten
 973 //          scl.Perform();
 974           BRepClass_FaceClassifierPerform(&scl,f,P,1e-7);
 975           TopAbs_State result = scl.State();
 976           MESSAGE("Position of point on face: "<<result);
 977           if ( result == TopAbs_OUT )
 978               MESSAGE("Point is out of face: node is not created");
 979           if ( result == TopAbs_UNKNOWN )
 980               MESSAGE("Point position on face is unknown: node is not created");
 981           if ( result == TopAbs_ON )
 982               MESSAGE("Point is on border of face: node is not created");
 983           if ( result == TopAbs_IN )
 984           {
 985             // Point is inside face and not on border
 986             MESSAGE("Point is in face: node is created");
 987             double uvCoords[2]   = {attractor_iter->second[0],attractor_iter->second[1]};
 988             iatt++;
 989             MESSAGE("Add cad point on (u,v)=(" << uvCoords[0] << "," << uvCoords[1] << ") with id = " << iatt);
 990             cad_point_t* point_p = cad_point_new(fce, iatt, uvCoords);
 991             cad_point_set_tag(point_p, iatt);
 992           }
 993           FaceId2AttractorCoords.erase(faceKey);
 994         }
 995       }
 996
 997       // Enforced Vertices
 998       faceKey = FacesWithEnforcedVertices.FindIndex(f);
 999       std::map<int,std::set<std::vector<double> > >::const_iterator evmIt = FaceId2EnforcedVertexCoords.find(faceKey);
1000       if (evmIt != FaceId2EnforcedVertexCoords.end()) {
1001         std::cout << "Some enforced vertices are defined" << std::endl;
1002         int ienf = 0;
1003         std::set<std::vector<double> > evl;
1004 //         std::vector<double> ev;
1005         MESSAGE("Face indice: " << iface);
1006         MESSAGE("Adding enforced vertices");
1007         evl = evmIt->second;
1008         MESSAGE("Number of vertices to add: "<< evl.size())
1009         std::set< std::vector<double> >::const_iterator evlIt = evl.begin();
1010         for (; evlIt != evl.end(); ++evlIt) {
1011 //           ev = *evlIt;
1012 //         for (int i=0; i<evl.size() ; i++) {
1013 //           ev = evl[i];
1014
1015 //           double xyzCoords[3]  = {ev[2], ev[3], ev[4]};
1016           double xyzCoords[3]  = {evlIt->at(0), evlIt->at(3), evlIt->at(4)};
1017           MESSAGE("Check position of vertex =(" << xyzCoords[0] << "," << xyzCoords[1] << "," << xyzCoords[2] << ")");
1018           gp_Pnt P(xyzCoords[0],xyzCoords[1],xyzCoords[2]);
1019           BRepClass_FaceClassifier scl(f,P,1e-7);
1020           // scl.Perform() is bugged. The function was rewritten
1021 //          scl.Perform();
1022           BRepClass_FaceClassifierPerform(&scl,f,P,1e-7);
1023           TopAbs_State result = scl.State();
1024           MESSAGE("Position of point on face: "<<result);
1025           if ( result == TopAbs_OUT )
1026               MESSAGE("Point is out of face: node is not created");
1027           if ( result == TopAbs_UNKNOWN )
1028               MESSAGE("Point position on face is unknown: node is not created");
1029           if ( result == TopAbs_ON )
1030               MESSAGE("Point is on border of face: node is not created");
1031           if ( result == TopAbs_IN )
1032           {
1033             // Point is inside face and not on border
1034             MESSAGE("Point is in face: node is created");
1035 //             double uvCoords[2]   = {ev[0],ev[1]};
1036             double uvCoords[2]   = {evlIt->at(0),evlIt->at(1)};
1037             ienf++;
1038             MESSAGE("Add cad point on (u,v)=(" << uvCoords[0] << "," << uvCoords[1] << ") with id = " << ienf);
1039             cad_point_t* point_p = cad_point_new(fce, ienf, uvCoords);
1040             cad_point_set_tag(point_p, ienf);
1041           }
1042         }
1043         FaceId2EnforcedVertexCoords.erase(faceKey);
1044       }
1045       else
1046         std::cout << "No enforced vertex defined" << std::endl;
1047     }
1048
1049
1050     /****************************************************************************************
1051                                     EDGES
1052     *****************************************************************************************/
1053     int edgeKey = -1;
1054     for (TopExp_Explorer edge_iter(f,TopAbs_EDGE);edge_iter.More();edge_iter.Next()) {
1055       TopoDS_Edge e = TopoDS::Edge(edge_iter.Current());
1056       int ic = emap.FindIndex(e);
1057       if (ic <= 0)
1058         ic = emap.Add(e);
1059
1060       double tmin,tmax;
1061       curves.push_back(BRep_Tool::CurveOnSurface(e, f, tmin, tmax));
1062
1063       if (HasSizeMapOnEdge){
1064         edgeKey = EdgesWithSizeMap.FindIndex(e);
1065         if (EdgeId2SizeMap.find(edgeKey)!=EdgeId2SizeMap.end()) {
1066           theSizeMapStr = EdgeId2SizeMap[faceKey];
1067           if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
1068             continue;
1069           // Expr To Python function, verification is performed at validation in GUI
1070           PyObject * obj = NULL;
1071           obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
1072           Py_DECREF(obj);
1073           PyObject * func = NULL;
1074           func = PyObject_GetAttrString(main_mod, "f");
1075           EdgeId2PythonSmp[ic]=func;
1076           EdgeId2SizeMap.erase(edgeKey);
1077         }
1078       }
1079       cad_edge_t *edg = cad_edge_new(fce, ic, tmin, tmax, curv_fun, curves.back());
1080       cad_edge_set_tag(edg, ic);
1081       cad_edge_set_property(edg, EDGE_PROPERTY_SOFT_REQUIRED);
1082       if (e.Orientation() == TopAbs_INTERNAL)
1083         cad_edge_set_property(edg, EDGE_PROPERTY_INTERNAL);
1084
1085       int npts = 0;
1086       int ip1, ip2, *ip;
1087       gp_Pnt2d e0 = curves.back()->Value(tmin);
1088       gp_Pnt ee0 = surfaces.back()->Value(e0.X(), e0.Y());
1089       Standard_Real d1=0,d2=0;
1090
1091
1092       /****************************************************************************************
1093                                       VERTICES
1094       *****************************************************************************************/
1095       int vertexKey = -1;
1096       for (TopExp_Explorer ex_edge(e ,TopAbs_VERTEX); ex_edge.More(); ex_edge.Next()) {
1097         TopoDS_Vertex v = TopoDS::Vertex(ex_edge.Current());
1098         ++npts;
1099         if (npts == 1){
1100           ip = &ip1;
1101           d1 = ee0.SquareDistance(BRep_Tool::Pnt(v));
1102         } else {
1103           ip = &ip2;
1104           d2 = ee0.SquareDistance(BRep_Tool::Pnt(v));
1105         }
1106         *ip = pmap.FindIndex(v);
1107         if(*ip <= 0)
1108           *ip = pmap.Add(v);
1109
1110         vertexKey = VerticesWithSizeMap.FindIndex(v);
1111         if (HasSizeMapOnVertex){
1112           vertexKey = VerticesWithSizeMap.FindIndex(v);
1113           if (VertexId2SizeMap.find(vertexKey)!=VertexId2SizeMap.end()){
1114             theSizeMapStr = VertexId2SizeMap[faceKey];
1115             if (theSizeMapStr.find(bad_end) == (theSizeMapStr.size()-bad_end.size()-1))
1116               continue;
1117             // Expr To Python function, verification is performed at validation in GUI
1118             PyObject * obj = NULL;
1119             obj= PyRun_String(theSizeMapStr.c_str(), Py_file_input, main_dict, NULL);
1120             Py_DECREF(obj);
1121             PyObject * func = NULL;
1122             func = PyObject_GetAttrString(main_mod, "f");
1123             VertexId2PythonSmp[*ip]=func;
1124 //             VertexId2SizeMap.erase(vertexKey);   // do not erase if using a vector
1125           }
1126         }
1127       }
1128       if (npts != 2) {
1129         // should not happen
1130         MESSAGE("An edge does not have 2 extremities.");
1131       } else {
1132         if (d1 < d2)
1133           cad_edge_set_extremities(edg, ip1, ip2);
1134         else
1135           cad_edge_set_extremities(edg, ip2, ip1);
1136       }
1137     } // for edge
1138   } //for face
1139
1140
1141   PyGILState_Release(gstate);
1142
1143   blsurf_data_set_cad(bls, c);
1144
1145   std::cout << std::endl;
1146   std::cout << "Beginning of Surface Mesh generation" << std::endl;
1147   std::cout << std::endl;
1148
1149   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
1150 #ifndef WNT
1151   feclearexcept( FE_ALL_EXCEPT );
1152   int oldFEFlags = fedisableexcept( FE_ALL_EXCEPT );
1153 #endif
1154
1155   status_t status = STATUS_ERROR;
1156
1157   try {
1158     OCC_CATCH_SIGNALS;
1159
1160     status = blsurf_compute_mesh(bls);
1161
1162   }
1163   catch ( std::exception& exc ) {
1164     _comment += exc.what();
1165   }
1166   catch (Standard_Failure& ex) {
1167     _comment += ex.DynamicType()->Name();
1168     if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) {
1169       _comment += ": ";
1170       _comment += ex.GetMessageString();
1171     }
1172   }
1173   catch (...) {
1174     if ( _comment.empty() )
1175       _comment = "Exception in blsurf_compute_mesh()";
1176   }
1177   if ( status != STATUS_OK) {
1178     blsurf_session_delete(bls);
1179     cad_delete(c);
1180     context_delete(ctx);
1181
1182     return error(_comment);
1183   }
1184
1185   std::cout << std::endl;
1186   std::cout << "End of Surface Mesh generation" << std::endl;
1187   std::cout << std::endl;
1188
1189   mesh_t *msh;
1190   blsurf_data_get_mesh(bls, &msh);
1191   if(!msh){
1192     blsurf_session_delete(bls);
1193     cad_delete(c);
1194     context_delete(ctx);
1195
1196     return error(_comment);
1197     //return false;
1198   }
1199
1200   integer nv, ne, nt, nq, vtx[4], tag;
1201   real xyz[3];
1202
1203   mesh_get_vertex_count(msh, &nv);
1204   mesh_get_edge_count(msh, &ne);
1205   mesh_get_triangle_count(msh, &nt);
1206   mesh_get_quadrangle_count(msh, &nq);
1207
1208
1209   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
1210   SMDS_MeshNode** nodes = new SMDS_MeshNode*[nv+1];
1211   bool* tags = new bool[nv+1];
1212
1213   for(int iv=1;iv<=nv;iv++) {
1214     mesh_get_vertex_coordinates(msh, iv, xyz);
1215     mesh_get_vertex_tag(msh, iv, &tag);
1216     nodes[iv] = meshDS->AddNode(xyz[0], xyz[1], xyz[2]);
1217     // internal point are tagged to zero
1218     if(tag){
1219       meshDS->SetNodeOnVertex(nodes[iv], TopoDS::Vertex(pmap(tag)));
1220       tags[iv] = false;
1221     } else {
1222       tags[iv] = true;
1223     }
1224   }
1225
1226   for(int it=1;it<=ne;it++) {
1227     mesh_get_edge_vertices(msh, it, vtx);
1228     SMDS_MeshEdge* edg = meshDS->AddEdge(nodes[vtx[0]], nodes[vtx[1]]);
1229     mesh_get_edge_tag(msh, it, &tag);
1230
1231     if (tags[vtx[0]]) {
1232       Set_NodeOnEdge(meshDS, nodes[vtx[0]], emap(tag));
1233       tags[vtx[0]] = false;
1234     };
1235     if (tags[vtx[1]]) {
1236       Set_NodeOnEdge(meshDS, nodes[vtx[1]], emap(tag));
1237       tags[vtx[1]] = false;
1238     };
1239     meshDS->SetMeshElementOnShape(edg, TopoDS::Edge(emap(tag)));
1240
1241   }
1242
1243   for(int it=1;it<=nt;it++) {
1244     mesh_get_triangle_vertices(msh, it, vtx);
1245     SMDS_MeshFace* tri = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]]);
1246     mesh_get_triangle_tag(msh, it, &tag);
1247     meshDS->SetMeshElementOnShape(tri, TopoDS::Face(fmap(tag)));
1248     if (tags[vtx[0]]) {
1249       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
1250       tags[vtx[0]] = false;
1251     };
1252     if (tags[vtx[1]]) {
1253       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
1254       tags[vtx[1]] = false;
1255     };
1256     if (tags[vtx[2]]) {
1257       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
1258       tags[vtx[2]] = false;
1259     };
1260   }
1261
1262   for(int it=1;it<=nq;it++) {
1263     mesh_get_quadrangle_vertices(msh, it, vtx);
1264     SMDS_MeshFace* quad = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]], nodes[vtx[3]]);
1265     mesh_get_quadrangle_tag(msh, it, &tag);
1266     meshDS->SetMeshElementOnShape(quad, TopoDS::Face(fmap(tag)));
1267     if (tags[vtx[0]]) {
1268       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
1269       tags[vtx[0]] = false;
1270     };
1271     if (tags[vtx[1]]) {
1272       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
1273       tags[vtx[1]] = false;
1274     };
1275     if (tags[vtx[2]]) {
1276       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
1277       tags[vtx[2]] = false;
1278     };
1279     if (tags[vtx[3]]) {
1280       meshDS->SetNodeOnFace(nodes[vtx[3]], TopoDS::Face(fmap(tag)));
1281       tags[vtx[3]] = false;
1282     };
1283   }
1284
1285   delete nodes;
1286
1287   /* release the mesh object */
1288   blsurf_data_regain_mesh(bls, msh);
1289
1290   /* clean up everything */
1291   blsurf_session_delete(bls);
1292   cad_delete(c);
1293
1294   context_delete(ctx);
1295
1296   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
1297 #ifndef WNT
1298   if ( oldFEFlags > 0 )
1299     feenableexcept( oldFEFlags );
1300   feclearexcept( FE_ALL_EXCEPT );
1301 #endif
1302
1303   std::cout << "FacesWithSizeMap" << std::endl;
1304   FacesWithSizeMap.Statistics(std::cout);
1305   std::cout << "EdgesWithSizeMap" << std::endl;
1306   EdgesWithSizeMap.Statistics(std::cout);
1307   std::cout << "VerticesWithSizeMap" << std::endl;
1308   VerticesWithSizeMap.Statistics(std::cout);
1309   std::cout << "FacesWithEnforcedVertices" << std::endl;
1310   FacesWithEnforcedVertices.Statistics(std::cout);
1311
1312   return true;
1313 }
1314
1315 //=============================================================================
1316 /*!
1317  *  SetNodeOnEdge
1318  */
1319 //=============================================================================
1320
1321 void BLSURFPlugin_BLSURF::Set_NodeOnEdge(SMESHDS_Mesh* meshDS, SMDS_MeshNode* node, const TopoDS_Shape& ed) {
1322   const TopoDS_Edge edge = TopoDS::Edge(ed);
1323
1324   gp_Pnt pnt(node->X(), node->Y(), node->Z());
1325
1326   Standard_Real p0 = 0.0;
1327   Standard_Real p1 = 1.0;
1328   Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, p0, p1);
1329
1330   GeomAPI_ProjectPointOnCurve proj(pnt, curve);
1331
1332   double pa = (double)proj.Parameter(1);
1333
1334   GProp_GProps LProps;
1335   BRepGProp::LinearProperties(ed, LProps);
1336   double lg = (double)LProps.Mass();
1337
1338   meshDS->SetNodeOnEdge(node, edge, pa);
1339 }
1340
1341 //=============================================================================
1342 /*!
1343  *
1344  */
1345 //=============================================================================
1346
1347 ostream & BLSURFPlugin_BLSURF::SaveTo(ostream & save)
1348 {
1349   return save;
1350 }
1351
1352 //=============================================================================
1353 /*!
1354  *
1355  */
1356 //=============================================================================
1357
1358 istream & BLSURFPlugin_BLSURF::LoadFrom(istream & load)
1359 {
1360   return load;
1361 }
1362
1363 //=============================================================================
1364 /*!
1365  *
1366  */
1367 //=============================================================================
1368
1369 ostream & operator << (ostream & save, BLSURFPlugin_BLSURF & hyp)
1370 {
1371   return hyp.SaveTo( save );
1372 }
1373
1374 //=============================================================================
1375 /*!
1376  *
1377  */
1378 //=============================================================================
1379
1380 istream & operator >> (istream & load, BLSURFPlugin_BLSURF & hyp)
1381 {
1382   return hyp.LoadFrom( load );
1383 }
1384
1385 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data)
1386 {
1387   const Geom2d_Curve*pargeo = (const Geom2d_Curve*) user_data;
1388
1389   if (uv){
1390     gp_Pnt2d P;
1391     P=pargeo->Value(t);
1392     uv[0]=P.X(); uv[1]=P.Y();
1393   }
1394
1395   if(dt) {
1396     gp_Vec2d V1;
1397     V1=pargeo->DN(t,1);
1398     dt[0]=V1.X(); dt[1]=V1.Y();
1399   }
1400
1401   if(dtt){
1402     gp_Vec2d V2;
1403     V2=pargeo->DN(t,2);
1404     dtt[0]=V2.X(); dtt[1]=V2.Y();
1405   }
1406
1407   return 0;
1408 }
1409
1410 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
1411                   real *duu, real *duv, real *dvv, void *user_data)
1412 {
1413   const Geom_Surface* geometry = (const Geom_Surface*) user_data;
1414
1415   if(xyz){
1416    gp_Pnt P;
1417    P=geometry->Value(uv[0],uv[1]);   // S.D0(U,V,P);
1418    xyz[0]=P.X(); xyz[1]=P.Y(); xyz[2]=P.Z();
1419   }
1420
1421   if(du && dv){
1422     gp_Pnt P;
1423     gp_Vec D1U,D1V;
1424
1425     geometry->D1(uv[0],uv[1],P,D1U,D1V);
1426     du[0]=D1U.X(); du[1]=D1U.Y(); du[2]=D1U.Z();
1427     dv[0]=D1V.X(); dv[1]=D1V.Y(); dv[2]=D1V.Z();
1428   }
1429
1430   if(duu && duv && dvv){
1431     gp_Pnt P;
1432     gp_Vec D1U,D1V;
1433     gp_Vec D2U,D2V,D2UV;
1434
1435     geometry->D2(uv[0],uv[1],P,D1U,D1V,D2U,D2V,D2UV);
1436     duu[0]=D2U.X(); duu[1]=D2U.Y(); duu[2]=D2U.Z();
1437     duv[0]=D2UV.X(); duv[1]=D2UV.Y(); duv[2]=D2UV.Z();
1438     dvv[0]=D2V.X(); dvv[1]=D2V.Y(); dvv[2]=D2V.Z();
1439   }
1440
1441   return 0;
1442 }
1443
1444
1445 status_t size_on_surface(integer face_id, real *uv, real *size, void *user_data)
1446 {
1447   if (face_id == 1) {
1448     if (my_u_min > uv[0]) {
1449       my_u_min = uv[0];
1450     }
1451     if (my_v_min > uv[1]) {
1452       my_v_min = uv[1];
1453     }
1454     if (my_u_max < uv[0]) {
1455       my_u_max = uv[0];
1456     }
1457     if (my_v_max < uv[1]) {
1458       my_v_max = uv[1];
1459     }
1460   }
1461
1462   if (FaceId2PythonSmp.count(face_id) != 0){
1463     PyObject * pyresult = NULL;
1464     PyObject* new_stderr = NULL;
1465     assert(Py_IsInitialized());
1466     PyGILState_STATE gstate;
1467     gstate = PyGILState_Ensure();
1468     pyresult = PyObject_CallFunction(FaceId2PythonSmp[face_id],"(f,f)",uv[0],uv[1]);
1469     double result;
1470     if ( pyresult == NULL){
1471       fflush(stderr);
1472       string err_description="";
1473       new_stderr = newPyStdOut(err_description);
1474       PySys_SetObject("stderr", new_stderr);
1475       PyErr_Print();
1476       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1477       Py_DECREF(new_stderr);
1478       MESSAGE("Can't evaluate f(" << uv[0] << "," << uv[1] << ")" << " error is " << err_description);
1479       result = *((double*)user_data);
1480       }
1481     else {
1482       result = PyFloat_AsDouble(pyresult);
1483       Py_DECREF(pyresult);
1484     }
1485     *size = result;
1486     //MESSAGE("f(" << uv[0] << "," << uv[1] << ")" << " = " << result);
1487     PyGILState_Release(gstate);
1488   }
1489   else {
1490     *size = *((double*)user_data);
1491   }
1492   return STATUS_OK;
1493 }
1494
1495 status_t size_on_edge(integer edge_id, real t, real *size, void *user_data)
1496 {
1497   if (EdgeId2PythonSmp.count(edge_id) != 0){
1498     PyObject * pyresult = NULL;
1499     PyObject* new_stderr = NULL;
1500     assert(Py_IsInitialized());
1501     PyGILState_STATE gstate;
1502     gstate = PyGILState_Ensure();
1503     pyresult = PyObject_CallFunction(EdgeId2PythonSmp[edge_id],"(f)",t);
1504     double result;
1505     if ( pyresult == NULL){
1506       fflush(stderr);
1507       string err_description="";
1508       new_stderr = newPyStdOut(err_description);
1509       PySys_SetObject("stderr", new_stderr);
1510       PyErr_Print();
1511       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1512       Py_DECREF(new_stderr);
1513       MESSAGE("Can't evaluate f(" << t << ")" << " error is " << err_description);
1514       result = *((double*)user_data);
1515       }
1516     else {
1517       result = PyFloat_AsDouble(pyresult);
1518       Py_DECREF(pyresult);
1519     }
1520     *size = result;
1521     PyGILState_Release(gstate);
1522   }
1523   else {
1524     *size = *((double*)user_data);
1525   }
1526   return STATUS_OK;
1527 }
1528
1529 status_t size_on_vertex(integer point_id, real *size, void *user_data)
1530 {
1531   if (VertexId2PythonSmp.count(point_id) != 0){
1532     PyObject * pyresult = NULL;
1533     PyObject* new_stderr = NULL;
1534     assert(Py_IsInitialized());
1535     PyGILState_STATE gstate;
1536     gstate = PyGILState_Ensure();
1537     pyresult = PyObject_CallFunction(VertexId2PythonSmp[point_id],"");
1538     double result;
1539     if ( pyresult == NULL){
1540       fflush(stderr);
1541       string err_description="";
1542       new_stderr = newPyStdOut(err_description);
1543       PySys_SetObject("stderr", new_stderr);
1544       PyErr_Print();
1545       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1546       Py_DECREF(new_stderr);
1547       MESSAGE("Can't evaluate f()" << " error is " << err_description);
1548       result = *((double*)user_data);
1549       }
1550     else {
1551       result = PyFloat_AsDouble(pyresult);
1552       Py_DECREF(pyresult);
1553     }
1554     *size = result;
1555     PyGILState_Release(gstate);
1556   }
1557   else {
1558     *size = *((double*)user_data);
1559   }
1560  return STATUS_OK;
1561 }
1562
1563 status_t message_callback(message_t *msg, void *user_data)
1564 {
1565   integer errnumber = 0;
1566   char *desc;
1567   message_get_number(msg, &errnumber);
1568   message_get_description(msg, &desc);
1569   if ( errnumber < 0 ) {
1570     string * error = (string*)user_data;
1571 //   if ( !error->empty() )
1572 //     *error += "\n";
1573     // remove ^A from the tail
1574     int len = strlen( desc );
1575     while (len > 0 && desc[len-1] != '\n')
1576       len--;
1577     error->append( desc, len );
1578   }
1579   else {
1580       std::cout << desc << std::endl;
1581   }
1582   return STATUS_OK;
1583 }
1584
1585
1586 //=============================================================================
1587 /*!
1588  *
1589  */
1590 //=============================================================================
1591 bool BLSURFPlugin_BLSURF::Evaluate(SMESH_Mesh& aMesh,
1592                                    const TopoDS_Shape& aShape,
1593                                    MapShapeNbElems& aResMap)
1594 {
1595   int    _physicalMesh  = BLSURFPlugin_Hypothesis::GetDefaultPhysicalMesh();
1596   double _phySize       = BLSURFPlugin_Hypothesis::GetDefaultPhySize();
1597   //int    _geometricMesh = BLSURFPlugin_Hypothesis::GetDefaultGeometricMesh();
1598   //double _angleMeshS    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshS();
1599   double _angleMeshC    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshC();
1600   if(_hypothesis) {
1601     _physicalMesh  = (int) _hypothesis->GetPhysicalMesh();
1602     _phySize       = _hypothesis->GetPhySize();
1603     //_geometricMesh = (int) hyp->GetGeometricMesh();
1604     //_angleMeshS    = hyp->GetAngleMeshS();
1605     _angleMeshC    = _hypothesis->GetAngleMeshC();
1606   }
1607
1608   bool IsQuadratic = false;
1609
1610   // ----------------
1611   // evaluate 1D
1612   // ----------------
1613   TopTools_DataMapOfShapeInteger EdgesMap;
1614   double fullLen = 0.0;
1615   double fullNbSeg = 0;
1616   for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
1617     TopoDS_Edge E = TopoDS::Edge( exp.Current() );
1618     if( EdgesMap.IsBound(E) )
1619       continue;
1620     SMESH_subMesh *sm = aMesh.GetSubMesh(E);
1621     double aLen = SMESH_Algo::EdgeLength(E);
1622     fullLen += aLen;
1623     int nb1d = 0;
1624     if(_physicalMesh==1) {
1625        nb1d = (int)( aLen/_phySize + 1 );
1626     }
1627     else {
1628       // use geometry
1629       double f,l;
1630       Handle(Geom_Curve) C = BRep_Tool::Curve(E,f,l);
1631       double fullAng = 0.0;
1632       double dp = (l-f)/200;
1633       gp_Pnt P1,P2,P3;
1634       C->D0(f,P1);
1635       C->D0(f+dp,P2);
1636       gp_Vec V1(P1,P2);
1637       for(int j=2; j<=200; j++) {
1638         C->D0(f+dp*j,P3);
1639         gp_Vec V2(P2,P3);
1640         fullAng += fabs(V1.Angle(V2));
1641         V1 = V2;
1642         P2 = P3;
1643       }
1644       nb1d = (int)( fullAng/_angleMeshC + 1 );
1645     }
1646     fullNbSeg += nb1d;
1647     std::vector<int> aVec(SMDSEntity_Last);
1648     for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
1649     if( IsQuadratic > 0 ) {
1650       aVec[SMDSEntity_Node] = 2*nb1d - 1;
1651       aVec[SMDSEntity_Quad_Edge] = nb1d;
1652     }
1653     else {
1654       aVec[SMDSEntity_Node] = nb1d - 1;
1655       aVec[SMDSEntity_Edge] = nb1d;
1656     }
1657     aResMap.insert(std::make_pair(sm,aVec));
1658     EdgesMap.Bind(E,nb1d);
1659   }
1660   double ELen = fullLen/fullNbSeg;
1661   // ----------------
1662   // evaluate 2D
1663   // ----------------
1664   // try to evaluate as in MEFISTO
1665   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
1666     TopoDS_Face F = TopoDS::Face( exp.Current() );
1667     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
1668     GProp_GProps G;
1669     BRepGProp::SurfaceProperties(F,G);
1670     double anArea = G.Mass();
1671     int nb1d = 0;
1672     for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next()) {
1673       nb1d += EdgesMap.Find(exp1.Current());
1674     }
1675     int nbFaces = (int) ( anArea/( ELen*ELen*sqrt(3.) / 4 ) );
1676     int nbNodes = (int) ( ( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
1677     std::vector<int> aVec(SMDSEntity_Last);
1678     for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
1679     if( IsQuadratic ) {
1680       int nb1d_in = (nbFaces*3 - nb1d) / 2;
1681       aVec[SMDSEntity_Node] = nbNodes + nb1d_in;
1682       aVec[SMDSEntity_Quad_Triangle] = nbFaces;
1683     }
1684     else {
1685       aVec[SMDSEntity_Node] = nbNodes;
1686       aVec[SMDSEntity_Triangle] = nbFaces;
1687     }
1688     aResMap.insert(std::make_pair(sm,aVec));
1689   }
1690
1691   // ----------------
1692   // evaluate 3D
1693   // ----------------
1694   GProp_GProps G;
1695   BRepGProp::VolumeProperties(aShape,G);
1696   double aVolume = G.Mass();
1697   double tetrVol = 0.1179*ELen*ELen*ELen;
1698   int nbVols = (int)aVolume/tetrVol;
1699   int nb1d_in = (int) ( ( nbVols*6 - fullNbSeg ) / 6 );
1700   std::vector<int> aVec(SMDSEntity_Last);
1701   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
1702   if( IsQuadratic ) {
1703     aVec[SMDSEntity_Node] = nb1d_in/3 + 1 + nb1d_in;
1704     aVec[SMDSEntity_Quad_Tetra] = nbVols;
1705   }
1706   else {
1707     aVec[SMDSEntity_Node] = nb1d_in/3 + 1;
1708     aVec[SMDSEntity_Tetra] = nbVols;
1709   }
1710   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
1711   aResMap.insert(std::make_pair(sm,aVec));
1712
1713   return true;
1714 }
1715
1716 //=============================================================================
1717 /*!
1718  *  Rewritting of the BRepClass_FaceClassifier::Perform function which is bugged (CAS 6.3sp6)
1719  *  Following line was added:
1720  *        myExtrem.Perform(P);
1721  */
1722 //=============================================================================
1723 void  BLSURFPlugin_BLSURF::BRepClass_FaceClassifierPerform(BRepClass_FaceClassifier* fc,
1724                     const TopoDS_Face& face,
1725                     const gp_Pnt& P,
1726                     const Standard_Real Tol)
1727 {
1728   //-- Voir BRepExtrema_ExtPF.cxx
1729   BRepAdaptor_Surface Surf(face);
1730   Standard_Real U1, U2, V1, V2;
1731   BRepTools::UVBounds(face, U1, U2, V1, V2);
1732   Extrema_ExtPS myExtrem;
1733   myExtrem.Initialize(Surf, U1, U2, V1, V2, Tol, Tol);
1734   myExtrem.Perform(P);
1735   //----------------------------------------------------------
1736   //-- On cherche le point le plus proche , PUIS
1737   //-- On le classifie.
1738   Standard_Integer nbv    = 0; // xpu
1739   Standard_Real MaxDist   =  RealLast();
1740   Standard_Integer indice = 0;
1741   if(myExtrem.IsDone()) {
1742     nbv = myExtrem.NbExt();
1743     for (Standard_Integer i = 1; i <= nbv; i++) {
1744       Standard_Real d = myExtrem.Value(i);
1745       d = Abs(d);
1746       if(d <= MaxDist) {
1747     MaxDist = d;
1748     indice = i;
1749       }
1750     }
1751   }
1752   if(indice) {
1753     gp_Pnt2d Puv;
1754     Standard_Real U1,U2;
1755     myExtrem.Point(indice).Parameter(U1, U2);
1756     Puv.SetCoord(U1, U2);
1757     fc->Perform(face, Puv, Tol);
1758   }
1759   else {
1760     fc->Perform(face, gp_Pnt2d(U1-1.0,V1 - 1.0), Tol); //-- NYI etc BUG PAS BEAU En attendant l acces a rejected
1761     //-- le resultat est TopAbs_OUT;
1762   }
1763 }
1764