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 // ---
  24 //
  25 #include "BLSURFPlugin_BLSURF.hxx"
  26 #include "BLSURFPlugin_Hypothesis.hxx"
  27
  28 #include <structmember.h>
  29
  30
  31 #include <SMESH_Gen.hxx>
  32 #include <SMESH_Mesh.hxx>
  33 #include <SMESH_ControlsDef.hxx>
  34 #include <SMESHGUI_Utils.h>
  35
  36 #include <SMESHDS_Mesh.hxx>
  37 #include <SMDS_MeshElement.hxx>
  38 #include <SMDS_MeshNode.hxx>
  39
  40 #include <utilities.h>
  41
  42 #include <limits>
  43 #include <list>
  44 #include <vector>
  45 #include <tr1/unordered_map>
  46 #include <cstdlib>
  47
  48 #include <BRep_Tool.hxx>
  49 #include <TopExp.hxx>
  50 #include <TopExp_Explorer.hxx>
  51 #include <TopoDS.hxx>
  52 #include <NCollection_Map.hxx>
  53 #include <Standard_ErrorHandler.hxx>
  54
  55 extern "C"{
  56 #include <distene/api.h>
  57 }
  58
  59 #include <Geom_Surface.hxx>
  60 #include <Handle_Geom_Surface.hxx>
  61 #include <Geom2d_Curve.hxx>
  62 #include <Handle_Geom2d_Curve.hxx>
  63 #include <Geom_Curve.hxx>
  64 #include <Handle_Geom_Curve.hxx>
  65 #include <TopoDS_Vertex.hxx>
  66 #include <TopoDS_Edge.hxx>
  67 #include <TopoDS_Wire.hxx>
  68 #include <TopoDS_Face.hxx>
  69 #include <TopoDS_Shape.hxx>
  70 #include <gp_Pnt2d.hxx>
  71 #include <TopTools_IndexedMapOfShape.hxx>
  72 #include <BRepTools.hxx>
  73
  74 #ifndef WNT
  75 #include <fenv.h>
  76 #endif
  77
  78 #include  <GeomSelectionTools.h>
  79
  80
  81 /* ===========  PYTHON ==============
  82  *
  83  *
  84  *  On est trop des brutes !
  85  *
  86  *
  87  *
  88  * ==================================*/
  89
  90 typedef struct {
  91   PyObject_HEAD
  92   int softspace;
  93   std::string *out;
  94   } PyStdOut;
  95
  96 static void
  97 PyStdOut_dealloc(PyStdOut *self)
  98 {
  99   PyObject_Del(self);
 100 }
 101
 102 static PyObject *
 103 PyStdOut_write(PyStdOut *self, PyObject *args)
 104 {
 105   char *c;
 106   int l;
 107   if (!PyArg_ParseTuple(args, "t#:write",&c, &l))
 108     return NULL;
 109
 110   //std::cerr << c ;
 111   *(self->out)=*(self->out)+c;
 112
 113   Py_INCREF(Py_None);
 114   return Py_None;
 115 }
 116
 117 static PyMethodDef PyStdOut_methods[] = {
 118   {"write",  (PyCFunction)PyStdOut_write,  METH_VARARGS,
 119     PyDoc_STR("write(string) -> None")},
 120   {NULL,    NULL}   /* sentinel */
 121 };
 122
 123 static PyMemberDef PyStdOut_memberlist[] = {
 124   {"softspace", T_INT,  offsetof(PyStdOut, softspace), 0,
 125    "flag indicating that a space needs to be printed; used by print"},
 126   {NULL} /* Sentinel */
 127 };
 128
 129 static PyTypeObject PyStdOut_Type = {
 130   /* The ob_type field must be initialized in the module init function
 131    * to be portable to Windows without using C++. */
 132   PyObject_HEAD_INIT(NULL)
 133   0,      /*ob_size*/
 134   "PyOut",   /*tp_name*/
 135   sizeof(PyStdOut),  /*tp_basicsize*/
 136   0,      /*tp_itemsize*/
 137   /* methods */
 138   (destructor)PyStdOut_dealloc, /*tp_dealloc*/
 139   0,      /*tp_print*/
 140   0, /*tp_getattr*/
 141   0, /*tp_setattr*/
 142   0,      /*tp_compare*/
 143   0,      /*tp_repr*/
 144   0,      /*tp_as_number*/
 145   0,      /*tp_as_sequence*/
 146   0,      /*tp_as_mapping*/
 147   0,      /*tp_hash*/
 148         0,                      /*tp_call*/
 149         0,                      /*tp_str*/
 150         PyObject_GenericGetAttr,                      /*tp_getattro*/
 151         /* softspace is writable:  we must supply tp_setattro */
 152         PyObject_GenericSetAttr,    /* tp_setattro */
 153         0,                      /*tp_as_buffer*/
 154         Py_TPFLAGS_DEFAULT,     /*tp_flags*/
 155         0,                      /*tp_doc*/
 156         0,                      /*tp_traverse*/
 157         0,                      /*tp_clear*/
 158         0,                      /*tp_richcompare*/
 159         0,                      /*tp_weaklistoffset*/
 160         0,                      /*tp_iter*/
 161         0,                      /*tp_iternext*/
 162         PyStdOut_methods,                      /*tp_methods*/
 163         PyStdOut_memberlist,                      /*tp_members*/
 164         0,                      /*tp_getset*/
 165         0,                      /*tp_base*/
 166         0,                      /*tp_dict*/
 167         0,                      /*tp_descr_get*/
 168         0,                      /*tp_descr_set*/
 169         0,                      /*tp_dictoffset*/
 170         0,                      /*tp_init*/
 171         0,                      /*tp_alloc*/
 172         0,                      /*tp_new*/
 173         0,                      /*tp_free*/
 174         0,                      /*tp_is_gc*/
 175 };
 176
 177 PyObject * newPyStdOut( std::string& out )
 178 {
 179   PyStdOut *self;
 180   self = PyObject_New(PyStdOut, &PyStdOut_Type);
 181   if (self == NULL)
 182     return NULL;
 183   self->softspace = 0;
 184   self->out=&out;
 185   return (PyObject*)self;
 186 }
 187
 188
 189 ////////////////////////END PYTHON///////////////////////////
 190
 191
 192 //////////////////TO WELL DEFINE UNORDERED TOPODS_MAP////////
 193
 194 namespace std {
 195 namespace tr1 {
 196 template<>
 197   struct hash<TopoDS_Face>
 198   : public std::unary_function<TopoDS_Face,std::size_t>
 199   {
 200     std::size_t operator() (const TopoDS_Face & Face) const
 201     { return Face.HashCode(471662);}
 202   //  { return Face.HashCode(std::numeric_limits<int>::max());}
 203 };
 204
 205 template<>
 206   struct hash<TopoDS_Edge>
 207   : public std::unary_function<TopoDS_Edge,std::size_t>
 208   {
 209     std::size_t operator() (const TopoDS_Edge & Edge) const
 210     { return Edge.HashCode(471662);}
 211   //  { return Face.HashCode(std::numeric_limits<int>::max());}
 212 };
 213
 214 template<>
 215   struct hash<TopoDS_Vertex>
 216   : public std::unary_function<TopoDS_Vertex,std::size_t>
 217   {
 218     std::size_t operator() (const TopoDS_Vertex & Vertex) const
 219     { return Vertex.HashCode(471662);}
 220   //  { return Face.HashCode(std::numeric_limits<int>::max());}
 221 };
 222 }
 223 }
 224
 225 //////////////////MY MAPS////////////////////////////////////////
 226 std::tr1::unordered_map<TopoDS_Face,string> Face2SizeMap;
 227 std::tr1::unordered_map<TopoDS_Edge,string> Edge2SizeMap;
 228 std::tr1::unordered_map<TopoDS_Vertex,string> Vertex2SizeMap;
 229 std::map<int,PyObject*> FaceId2PythonSmp;
 230 std::map<int,PyObject*> EdgeId2PythonSmp;
 231 std::map<int,PyObject*> VertexId2PythonSmp;
 232
 233
 234 bool HasSizeMapOnFace=false;
 235 bool HasSizeMapOnEdge=false;
 236 bool HasSizeMapOnVertex=false;
 237
 238 //=============================================================================
 239 /*!
 240  *
 241  */
 242 //=============================================================================
 243
 244 BLSURFPlugin_BLSURF::BLSURFPlugin_BLSURF(int hypId, int studyId,
 245                                                SMESH_Gen* gen)
 246   : SMESH_2D_Algo(hypId, studyId, gen)
 247 {
 248   MESSAGE("BLSURFPlugin_BLSURF::BLSURFPlugin_BLSURF");
 249
 250   _name = "BLSURF";
 251   _shapeType = (1 << TopAbs_FACE); // 1 bit /shape type
 252   _compatibleHypothesis.push_back("BLSURF_Parameters");
 253   _requireDescretBoundary = false;
 254   _onlyUnaryInput = false;
 255   _hypothesis = NULL;
 256
 257
 258   /* Initialize the Python interpreter */
 259   assert(Py_IsInitialized());
 260   PyGILState_STATE gstate;
 261   gstate = PyGILState_Ensure();
 262
 263   main_mod = NULL;
 264   main_mod = PyImport_AddModule("__main__");
 265
 266   main_dict = NULL;
 267   main_dict = PyModule_GetDict(main_mod);
 268
 269   PyRun_SimpleString("from math import *");
 270   PyGILState_Release(gstate);
 271
 272 }
 273
 274 //=============================================================================
 275 /*!
 276  *
 277  */
 278 //=============================================================================
 279
 280 BLSURFPlugin_BLSURF::~BLSURFPlugin_BLSURF()
 281 {
 282   MESSAGE("BLSURFPlugin_BLSURF::~BLSURFPlugin_BLSURF");
 283 }
 284
 285
 286 //=============================================================================
 287 /*!
 288  *
 289  */
 290 //=============================================================================
 291
 292 bool BLSURFPlugin_BLSURF::CheckHypothesis
 293                          (SMESH_Mesh&                          aMesh,
 294                           const TopoDS_Shape&                  aShape,
 295                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 296 {
 297   _hypothesis = NULL;
 298
 299   list<const SMESHDS_Hypothesis*>::const_iterator itl;
 300   const SMESHDS_Hypothesis* theHyp;
 301
 302   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
 303   int nbHyp = hyps.size();
 304   if (!nbHyp)
 305   {
 306     aStatus = SMESH_Hypothesis::HYP_OK;
 307     return true;  // can work with no hypothesis
 308   }
 309
 310   itl = hyps.begin();
 311   theHyp = (*itl); // use only the first hypothesis
 312
 313   string hypName = theHyp->GetName();
 314
 315   if (hypName == "BLSURF_Parameters")
 316   {
 317     _hypothesis = static_cast<const BLSURFPlugin_Hypothesis*> (theHyp);
 318     ASSERT(_hypothesis);
 319     if ( _hypothesis->GetPhysicalMesh() == BLSURFPlugin_Hypothesis::DefaultSize &&
 320          _hypothesis->GetGeometricMesh() == BLSURFPlugin_Hypothesis::DefaultGeom )
 321       //  hphy_flag = 0 and hgeo_flag = 0 is not allowed (spec)
 322       aStatus = SMESH_Hypothesis::HYP_BAD_PARAMETER;
 323     else
 324       aStatus = SMESH_Hypothesis::HYP_OK;
 325   }
 326   else
 327     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 328
 329   return aStatus == SMESH_Hypothesis::HYP_OK;
 330 }
 331
 332 //=============================================================================
 333 /*!
 334  * Pass parameters to BLSURF
 335  */
 336 //=============================================================================
 337
 338 inline std::string to_string(double d)
 339 {
 340    std::ostringstream o;
 341    o << d;
 342    return o.str();
 343 }
 344
 345 inline std::string to_string(int i)
 346 {
 347    std::ostringstream o;
 348    o << i;
 349    return o.str();
 350 }
 351
 352 double _smp_phy_size;
 353 status_t size_on_surface(integer face_id, real *uv, real *size, void *user_data);
 354 status_t size_on_edge(integer edge_id, real t, real *size, void *user_data);
 355 status_t size_on_vertex(integer vertex_id, real *size, void *user_data);
 356
 357 double my_u_min=1e6,my_v_min=1e6,my_u_max=-1e6,my_v_max=-1e6;
 358
 359 void BLSURFPlugin_BLSURF::SetParameters(const BLSURFPlugin_Hypothesis* hyp, blsurf_session_t *bls)
 360 {
 361   int    _topology      = BLSURFPlugin_Hypothesis::GetDefaultTopology();
 362   int    _physicalMesh  = BLSURFPlugin_Hypothesis::GetDefaultPhysicalMesh();
 363   double _phySize       = BLSURFPlugin_Hypothesis::GetDefaultPhySize();
 364   int    _geometricMesh = BLSURFPlugin_Hypothesis::GetDefaultGeometricMesh();
 365   double _angleMeshS    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshS();
 366   double _angleMeshC    = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshC();
 367   double _gradation     = BLSURFPlugin_Hypothesis::GetDefaultGradation();
 368   bool   _quadAllowed   = BLSURFPlugin_Hypothesis::GetDefaultQuadAllowed();
 369   bool   _decimesh      = BLSURFPlugin_Hypothesis::GetDefaultDecimesh();
 370   int    _verb          = BLSURFPlugin_Hypothesis::GetDefaultVerbosity();
 371
 372   if (hyp) {
 373     MESSAGE("BLSURFPlugin_BLSURF::SetParameters");
 374     _topology      = (int) hyp->GetTopology();
 375     _physicalMesh  = (int) hyp->GetPhysicalMesh();
 376     _phySize       = hyp->GetPhySize();
 377     _geometricMesh = (int) hyp->GetGeometricMesh();
 378     _angleMeshS    = hyp->GetAngleMeshS();
 379     _angleMeshC    = hyp->GetAngleMeshC();
 380     _gradation     = hyp->GetGradation();
 381     _quadAllowed   = hyp->GetQuadAllowed();
 382     _decimesh      = hyp->GetDecimesh();
 383     _verb          = hyp->GetVerbosity();
 384
 385     if ( hyp->GetPhyMin() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 386       blsurf_set_param(bls, "hphymin", to_string(hyp->GetPhyMin()).c_str());
 387     if ( hyp->GetPhyMax() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 388       blsurf_set_param(bls, "hphymax", to_string(hyp->GetPhyMax()).c_str());
 389     if ( hyp->GetGeoMin() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 390       blsurf_set_param(bls, "hgeomin", to_string(hyp->GetGeoMin()).c_str());
 391     if ( hyp->GetGeoMax() != ::BLSURFPlugin_Hypothesis::undefinedDouble() )
 392       blsurf_set_param(bls, "hgeomax", to_string(hyp->GetGeoMax()).c_str());
 393
 394     const BLSURFPlugin_Hypothesis::TOptionValues & opts = hyp->GetOptionValues();
 395     BLSURFPlugin_Hypothesis::TOptionValues::const_iterator opIt;
 396     for ( opIt = opts.begin(); opIt != opts.end(); ++opIt )
 397       if ( !opIt->second.empty() ) {
 398 #ifdef _DEBUG_
 399         cout << "blsurf_set_param(): " << opIt->first << " = " << opIt->second << endl;
 400 #endif
 401         blsurf_set_param(bls, opIt->first.c_str(), opIt->second.c_str());
 402       }
 403
 404   } else {
 405     MESSAGE("BLSURFPlugin_BLSURF::SetParameters using defaults");
 406   }
 407   _smp_phy_size = _phySize;
 408   blsurf_set_param(bls, "topo_points",       _topology > 0 ? "1" : "0");
 409   blsurf_set_param(bls, "topo_curves",       _topology > 0 ? "1" : "0");
 410   blsurf_set_param(bls, "topo_project",      _topology > 0 ? "1" : "0");
 411   blsurf_set_param(bls, "clean_boundary",    _topology > 1 ? "1" : "0");
 412   blsurf_set_param(bls, "close_boundary",    _topology > 1 ? "1" : "0");
 413   blsurf_set_param(bls, "hphy_flag",         to_string(_physicalMesh).c_str());
 414 //  blsurf_set_param(bls, "hphy_flag",         "2");
 415   if ((to_string(_physicalMesh))=="2"){
 416     MESSAGE("Setting a Size Map");
 417     const BLSURFPlugin_Hypothesis::TSizeMap & sizeMaps = hyp->GetSizeMapEntries();
 418     BLSURFPlugin_Hypothesis::TSizeMap::const_iterator smIt;
 419     GeomSelectionTools* GeomST;
 420     GeomST=new GeomSelectionTools::GeomSelectionTools( SMESH::GetActiveStudyDocument());
 421
 422     Face2SizeMap.clear();
 423     Edge2SizeMap.clear();
 424     Vertex2SizeMap.clear();
 425
 426     for ( smIt = sizeMaps.begin(); smIt != sizeMaps.end(); ++smIt ) {
 427       if ( !smIt->second.empty() ) {
 428 #ifdef _DEBUG_
 429         cout << "blsurf_set_sizeMap(): " << smIt->first << " = " << smIt->second << endl;
 430 #endif
 431         TopoDS_Shape GeomShape=GeomST->entryToShape(smIt->first);
 432         TopAbs_ShapeEnum GeomType= GeomShape.ShapeType();
 433         if (GeomType==TopAbs_FACE){
 434           HasSizeMapOnFace=true;
 435           Face2SizeMap[TopoDS::Face(GeomShape)]= smIt->second;
 436         }
 437         if (GeomType==TopAbs_EDGE){
 438           HasSizeMapOnEdge=true;
 439           HasSizeMapOnFace=true;
 440           Edge2SizeMap[TopoDS::Edge(GeomShape)]= smIt->second;
 441         }
 442         if (GeomType==TopAbs_VERTEX){
 443           HasSizeMapOnVertex=true;
 444           HasSizeMapOnEdge=true;
 445           HasSizeMapOnFace=true;
 446           Vertex2SizeMap[TopoDS::Vertex(GeomShape)]= smIt->second;
 447         }
 448       }
 449     }
 450 //    if (HasSizeMapOnFace){
 451     // In all size map cases (hphy_flag = 2), at least map on face must be defined
 452     std::cout << "Setting Size Map on FACES " << std::endl;
 453     blsurf_data_set_sizemap_iso_cad_face(bls, size_on_surface, &_smp_phy_size);
 454 //    }
 455
 456     if (HasSizeMapOnEdge){
 457       std::cout << "Setting Size Map on EDGES " << std::endl;
 458       blsurf_data_set_sizemap_iso_cad_edge(bls, size_on_edge, &_smp_phy_size);
 459     }
 460     if (HasSizeMapOnVertex){
 461       std::cout << "Setting Size Map on VERTICES " << std::endl;
 462       blsurf_data_set_sizemap_iso_cad_point(bls, size_on_vertex, &_smp_phy_size);
 463     }
 464   }
 465   blsurf_set_param(bls, "hphydef",           to_string(_phySize).c_str());
 466   blsurf_set_param(bls, "hgeo_flag",         to_string(_geometricMesh).c_str());
 467   blsurf_set_param(bls, "angle_meshs",       to_string(_angleMeshS).c_str());
 468   blsurf_set_param(bls, "angle_meshc",       to_string(_angleMeshC).c_str());
 469   blsurf_set_param(bls, "gradation",         to_string(_gradation).c_str());
 470   blsurf_set_param(bls, "patch_independent", _decimesh ? "1" : "0");
 471   blsurf_set_param(bls, "element",           _quadAllowed ? "q1.0" : "p1");
 472   blsurf_set_param(bls, "verb",              to_string(_verb).c_str());
 473 }
 474
 475 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data);
 476 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
 477                   real *duu, real *duv, real *dvv, void *user_data);
 478 status_t message_callback(message_t *msg, void *user_data);
 479
 480 //=============================================================================
 481 /*!
 482  *
 483  */
 484 //=============================================================================
 485
 486 bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) {
 487
 488   MESSAGE("BLSURFPlugin_BLSURF::Compute");
 489
 490   if (aShape.ShapeType() == TopAbs_COMPOUND) {
 491     cout << "  the shape is a COMPOUND" << endl;
 492   }
 493   else {
 494     cout << "  the shape is UNKNOWN" << endl;
 495   };
 496
 497   context_t *ctx =  context_new();
 498   context_set_message_callback(ctx, message_callback, &_comment);
 499
 500   cad_t *c = cad_new(ctx);
 501
 502   blsurf_session_t *bls = blsurf_session_new(ctx);
 503
 504
 505   SetParameters(_hypothesis, bls);
 506
 507   TopTools_IndexedMapOfShape fmap;
 508   TopTools_IndexedMapOfShape emap;
 509   TopTools_IndexedMapOfShape pmap;
 510   vector<Handle(Geom2d_Curve)> curves;
 511   vector<Handle(Geom_Surface)> surfaces;
 512
 513
 514
 515   fmap.Clear();
 516   FaceId2PythonSmp.clear();
 517   emap.Clear();
 518   EdgeId2PythonSmp.clear();
 519   pmap.Clear();
 520   VertexId2PythonSmp.clear();
 521   surfaces.resize(0);
 522   curves.resize(0);
 523
 524   assert(Py_IsInitialized());
 525   PyGILState_STATE gstate;
 526   gstate = PyGILState_Ensure();
 527 /*
 528   Standard_Real u_min;
 529   Standard_Real v_min;
 530   Standard_Real u_max;
 531   Standard_Real v_max;
 532 */
 533   int iface = 0;
 534   for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next()) {
 535     TopoDS_Face f=TopoDS::Face(face_iter.Current());
 536     if (fmap.FindIndex(f) > 0)
 537       continue;
 538
 539     fmap.Add(f);
 540     iface++;
 541     surfaces.push_back(BRep_Tool::Surface(f));
 542     // Get bound values of uv surface
 543     //BRep_Tool::Surface(f)->Bounds(u_min,u_max,v_min,v_max);
 544     //std::cout << "BRep_Tool::Surface(f)->Bounds(u_min,u_max,v_min,v_max): " << u_min << ", " << u_max << ", " << v_min << ", " << v_max << std::endl;
 545     if ((HasSizeMapOnFace) && Face2SizeMap.find(f)!=Face2SizeMap.end()){
 546       // Expr To Python function, verification is performed at validation in GUI
 547       PyObject * obj = NULL;
 548       obj= PyRun_String(Face2SizeMap[f].c_str(), Py_file_input, main_dict, NULL);
 549       Py_DECREF(obj);
 550       PyObject * func = NULL;
 551       func = PyObject_GetAttrString(main_mod, "f");
 552       FaceId2PythonSmp[iface]=func;
 553       Face2SizeMap.erase(f);
 554     }
 555     cad_face_t *fce = cad_face_new(c, iface, surf_fun, surfaces.back());
 556     cad_face_set_tag(fce, iface);
 557     if(f.Orientation() != TopAbs_FORWARD){
 558       cad_face_set_orientation(fce, CAD_ORIENTATION_REVERSED);
 559     } else {
 560       cad_face_set_orientation(fce, CAD_ORIENTATION_FORWARD);
 561     }
 562
 563     for (TopExp_Explorer edge_iter(f,TopAbs_EDGE);edge_iter.More();edge_iter.Next()) {
 564       TopoDS_Edge e = TopoDS::Edge(edge_iter.Current());
 565       int ic = emap.FindIndex(e);
 566       if (ic <= 0)
 567         ic = emap.Add(e);
 568
 569       double tmin,tmax;
 570       curves.push_back(BRep_Tool::CurveOnSurface(e, f, tmin, tmax));
 571       if ((HasSizeMapOnEdge) && Edge2SizeMap.find(e)!=Edge2SizeMap.end()){
 572         // Expr To Python function, verification is performed at validation in GUI
 573         PyObject * obj = NULL;
 574         obj= PyRun_String(Edge2SizeMap[e].c_str(), Py_file_input, main_dict, NULL);
 575         Py_DECREF(obj);
 576         PyObject * func = NULL;
 577         func = PyObject_GetAttrString(main_mod, "f");
 578         EdgeId2PythonSmp[ic]=func;
 579         Edge2SizeMap.erase(e);
 580       }
 581       cad_edge_t *edg = cad_edge_new(fce, ic, tmin, tmax, curv_fun, curves.back());
 582       cad_edge_set_tag(edg, ic);
 583       cad_edge_set_property(edg, EDGE_PROPERTY_SOFT_REQUIRED);
 584       if (e.Orientation() == TopAbs_INTERNAL)
 585         cad_edge_set_property(edg, EDGE_PROPERTY_INTERNAL);
 586
 587       int npts = 0;
 588       int ip1, ip2, *ip;
 589       gp_Pnt2d e0 = curves.back()->Value(tmin);
 590       gp_Pnt ee0 = surfaces.back()->Value(e0.X(), e0.Y());
 591       Standard_Real d1=0,d2=0;
 592       for (TopExp_Explorer ex_edge(e ,TopAbs_VERTEX); ex_edge.More(); ex_edge.Next()) {
 593         TopoDS_Vertex v = TopoDS::Vertex(ex_edge.Current());
 594         ++npts;
 595         if (npts == 1){
 596           ip = &ip1;
 597           d1 = ee0.SquareDistance(BRep_Tool::Pnt(v));
 598         } else {
 599           ip = &ip2;
 600           d2 = ee0.SquareDistance(BRep_Tool::Pnt(v));
 601         }
 602         *ip = pmap.FindIndex(v);
 603         if(*ip <= 0)
 604           *ip = pmap.Add(v);
 605         if ((HasSizeMapOnVertex) && Vertex2SizeMap.find(v)!=Vertex2SizeMap.end()){
 606           // Expr To Python function, verification is performed at validation in GUI
 607           PyObject * obj = NULL;
 608           obj= PyRun_String(Vertex2SizeMap[v].c_str(), Py_file_input, main_dict, NULL);
 609           Py_DECREF(obj);
 610           PyObject * func = NULL;
 611           func = PyObject_GetAttrString(main_mod, "f");
 612           VertexId2PythonSmp[*ip]=func;
 613           Vertex2SizeMap.erase(v);
 614         }
 615       }
 616       if (npts != 2) {
 617         // should not happen
 618         cout << "An edge does not have 2 extremities." << endl;
 619       } else {
 620         if (d1 < d2)
 621           cad_edge_set_extremities(edg, ip1, ip2);
 622         else
 623           cad_edge_set_extremities(edg, ip2, ip1);
 624       }
 625     } // for edge
 626   } //for face
 627
 628
 629   PyGILState_Release(gstate);
 630
 631   blsurf_data_set_cad(bls, c);
 632
 633   cout << endl;
 634   cout << "Beginning of Surface Mesh generation" << endl;
 635   cout << endl;
 636
 637   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
 638 #ifndef WNT
 639   feclearexcept( FE_ALL_EXCEPT );
 640   int oldFEFlags = fedisableexcept( FE_ALL_EXCEPT );
 641 #endif
 642
 643     status_t status = STATUS_ERROR;
 644
 645   try {
 646     OCC_CATCH_SIGNALS;
 647
 648     status = blsurf_compute_mesh(bls);
 649
 650   }
 651   catch ( std::exception& exc ) {
 652     _comment += exc.what();
 653   }
 654   catch (Standard_Failure& ex) {
 655     _comment += ex.DynamicType()->Name();
 656     if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) {
 657       _comment += ": ";
 658       _comment += ex.GetMessageString();
 659     }
 660   }
 661   catch (...) {
 662     if ( _comment.empty() )
 663       _comment = "Exception in blsurf_compute_mesh()";
 664   }
 665   if ( status != STATUS_OK) {
 666     blsurf_session_delete(bls);
 667     cad_delete(c);
 668     context_delete(ctx);
 669
 670     return error(_comment);
 671   }
 672
 673   cout << endl;
 674   cout << "End of Surface Mesh generation" << endl;
 675   cout << endl;
 676
 677     cout << "****************** U MIN: " << my_u_min << endl;
 678     cout << "****************** V MIN: " << my_v_min << endl;
 679     cout << "****************** U MAX: " << my_u_max << endl;
 680     cout << "****************** V MAX: " << my_v_max << endl;
 681
 682   mesh_t *msh;
 683   blsurf_data_get_mesh(bls, &msh);
 684   if(!msh){
 685     blsurf_session_delete(bls);
 686     cad_delete(c);
 687     context_delete(ctx);
 688
 689     return error(_comment);
 690     //return false;
 691   }
 692
 693   integer nv, ne, nt, nq, vtx[4], tag;
 694   real xyz[3];
 695
 696   mesh_get_vertex_count(msh, &nv);
 697   mesh_get_edge_count(msh, &ne);
 698   mesh_get_triangle_count(msh, &nt);
 699   mesh_get_quadrangle_count(msh, &nq);
 700
 701
 702   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 703   SMDS_MeshNode** nodes = new SMDS_MeshNode*[nv+1];
 704   bool* tags = new bool[nv+1];
 705
 706   for(int iv=1;iv<=nv;iv++) {
 707     mesh_get_vertex_coordinates(msh, iv, xyz);
 708     mesh_get_vertex_tag(msh, iv, &tag);
 709     nodes[iv] = meshDS->AddNode(xyz[0], xyz[1], xyz[2]);
 710     // internal point are tagged to zero
 711     if(tag){
 712       meshDS->SetNodeOnVertex(nodes[iv], TopoDS::Vertex(pmap(tag)));
 713       tags[iv] = false;
 714     } else {
 715       tags[iv] = true;
 716     }
 717   }
 718
 719   for(int it=1;it<=ne;it++) {
 720     mesh_get_edge_vertices(msh, it, vtx);
 721     SMDS_MeshEdge* edg = meshDS->AddEdge(nodes[vtx[0]], nodes[vtx[1]]);
 722     mesh_get_edge_tag(msh, it, &tag);
 723
 724     if (tags[vtx[0]]) {
 725       meshDS->SetNodeOnEdge(nodes[vtx[0]], TopoDS::Edge(emap(tag)));
 726       tags[vtx[0]] = false;
 727     };
 728     if (tags[vtx[1]]) {
 729       meshDS->SetNodeOnEdge(nodes[vtx[1]], TopoDS::Edge(emap(tag)));
 730       tags[vtx[1]] = false;
 731     };
 732     meshDS->SetMeshElementOnShape(edg, TopoDS::Edge(emap(tag)));
 733
 734   }
 735
 736   for(int it=1;it<=nt;it++) {
 737     mesh_get_triangle_vertices(msh, it, vtx);
 738     SMDS_MeshFace* tri = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]]);
 739     mesh_get_triangle_tag(msh, it, &tag);
 740     meshDS->SetMeshElementOnShape(tri, TopoDS::Face(fmap(tag)));
 741     if (tags[vtx[0]]) {
 742       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
 743       tags[vtx[0]] = false;
 744     };
 745     if (tags[vtx[1]]) {
 746       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
 747       tags[vtx[1]] = false;
 748     };
 749     if (tags[vtx[2]]) {
 750       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
 751       tags[vtx[2]] = false;
 752     };
 753   }
 754
 755   for(int it=1;it<=nq;it++) {
 756     mesh_get_quadrangle_vertices(msh, it, vtx);
 757     SMDS_MeshFace* quad = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]], nodes[vtx[3]]);
 758     mesh_get_quadrangle_tag(msh, it, &tag);
 759     meshDS->SetMeshElementOnShape(quad, TopoDS::Face(fmap(tag)));
 760     if (tags[vtx[0]]) {
 761       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
 762       tags[vtx[0]] = false;
 763     };
 764     if (tags[vtx[1]]) {
 765       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
 766       tags[vtx[1]] = false;
 767     };
 768     if (tags[vtx[2]]) {
 769       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
 770       tags[vtx[2]] = false;
 771     };
 772     if (tags[vtx[3]]) {
 773       meshDS->SetNodeOnFace(nodes[vtx[3]], TopoDS::Face(fmap(tag)));
 774       tags[vtx[3]] = false;
 775     };
 776   }
 777
 778   delete nodes;
 779
 780   /* release the mesh object */
 781   blsurf_data_regain_mesh(bls, msh);
 782
 783   /* clean up everything */
 784   blsurf_session_delete(bls);
 785   cad_delete(c);
 786
 787   context_delete(ctx);
 788
 789   // Issue 0019864. On DebianSarge, FE signals do not obey to OSD::SetSignal(false)
 790 #ifndef WNT
 791   if ( oldFEFlags > 0 )
 792     feenableexcept( oldFEFlags );
 793   feclearexcept( FE_ALL_EXCEPT );
 794 #endif
 795
 796   return true;
 797 }
 798
 799 //=============================================================================
 800 /*!
 801  *
 802  */
 803 //=============================================================================
 804
 805 ostream & BLSURFPlugin_BLSURF::SaveTo(ostream & save)
 806 {
 807   return save;
 808 }
 809
 810 //=============================================================================
 811 /*!
 812  *
 813  */
 814 //=============================================================================
 815
 816 istream & BLSURFPlugin_BLSURF::LoadFrom(istream & load)
 817 {
 818   return load;
 819 }
 820
 821 //=============================================================================
 822 /*!
 823  *
 824  */
 825 //=============================================================================
 826
 827 ostream & operator << (ostream & save, BLSURFPlugin_BLSURF & hyp)
 828 {
 829   return hyp.SaveTo( save );
 830 }
 831
 832 //=============================================================================
 833 /*!
 834  *
 835  */
 836 //=============================================================================
 837
 838 istream & operator >> (istream & load, BLSURFPlugin_BLSURF & hyp)
 839 {
 840   return hyp.LoadFrom( load );
 841 }
 842
 843 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data)
 844 {
 845   const Geom2d_Curve*pargeo = (const Geom2d_Curve*) user_data;
 846
 847   if (uv){
 848     gp_Pnt2d P;
 849     P=pargeo->Value(t);
 850     uv[0]=P.X(); uv[1]=P.Y();
 851   }
 852
 853   if(dt) {
 854     gp_Vec2d V1;
 855     V1=pargeo->DN(t,1);
 856     dt[0]=V1.X(); dt[1]=V1.Y();
 857   }
 858
 859   if(dtt){
 860     gp_Vec2d V2;
 861     V2=pargeo->DN(t,2);
 862     dtt[0]=V2.X(); dtt[1]=V2.Y();
 863   }
 864
 865   return 0;
 866 }
 867
 868 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
 869                   real *duu, real *duv, real *dvv, void *user_data)
 870 {
 871   const Geom_Surface* geometry = (const Geom_Surface*) user_data;
 872
 873   if(xyz){
 874    gp_Pnt P;
 875    P=geometry->Value(uv[0],uv[1]);   // S.D0(U,V,P);
 876    xyz[0]=P.X(); xyz[1]=P.Y(); xyz[2]=P.Z();
 877   }
 878
 879   if(du && dv){
 880     gp_Pnt P;
 881     gp_Vec D1U,D1V;
 882
 883     geometry->D1(uv[0],uv[1],P,D1U,D1V);
 884     du[0]=D1U.X(); du[1]=D1U.Y(); du[2]=D1U.Z();
 885     dv[0]=D1V.X(); dv[1]=D1V.Y(); dv[2]=D1V.Z();
 886   }
 887
 888   if(duu && duv && dvv){
 889     gp_Pnt P;
 890     gp_Vec D1U,D1V;
 891     gp_Vec D2U,D2V,D2UV;
 892
 893     geometry->D2(uv[0],uv[1],P,D1U,D1V,D2U,D2V,D2UV);
 894     duu[0]=D2U.X(); duu[1]=D2U.Y(); duu[2]=D2U.Z();
 895     duv[0]=D2UV.X(); duv[1]=D2UV.Y(); duv[2]=D2UV.Z();
 896     dvv[0]=D2V.X(); dvv[1]=D2V.Y(); dvv[2]=D2V.Z();
 897   }
 898
 899   return 0;
 900 }
 901
 902
 903 status_t size_on_surface(integer face_id, real *uv, real *size, void *user_data)
 904 {
 905   if (face_id == 1) {
 906     if (my_u_min > uv[0]) {
 907       my_u_min = uv[0];
 908 //      cout << "****************** FACE ID / U MIN: " << face_id << " / " << my_u_min << endl;
 909     }
 910     if (my_v_min > uv[1]) {
 911       my_v_min = uv[1];
 912 //      cout << "****************** FACE ID / V MIN: " << face_id << " / " << my_v_min << endl;
 913     }
 914     if (my_u_max < uv[0]) {
 915       my_u_max = uv[0];
 916 //      cout << "****************** FACE ID / U MAX: " << face_id << " / " << my_u_max << endl;
 917     }
 918     if (my_v_max < uv[1]) {
 919       my_v_max = uv[1];
 920 //      cout << "****************** FACE ID / V MAX: " << face_id << " / " << my_v_max << endl;
 921     }
 922   }
 923
 924   if (FaceId2PythonSmp.count(face_id) != 0){
 925     PyObject * pyresult = NULL;
 926     PyObject* new_stderr = NULL;
 927     assert(Py_IsInitialized());
 928     PyGILState_STATE gstate;
 929     gstate = PyGILState_Ensure();
 930     pyresult = PyObject_CallFunction(FaceId2PythonSmp[face_id],"(f,f)",uv[0],uv[1]);
 931     double result;
 932     if ( pyresult == NULL){
 933       fflush(stderr);
 934       string err_description="";
 935       new_stderr = newPyStdOut(err_description);
 936       PySys_SetObject("stderr", new_stderr);
 937       PyErr_Print();
 938       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
 939       Py_DECREF(new_stderr);
 940       std::cout << "Can't evaluate f(" << uv[0] << "," << uv[1] << ")" << " error is " << err_description << std::endl;
 941       result = *((double*)user_data);
 942       }
 943     else {
 944       result = PyFloat_AsDouble(pyresult);
 945       Py_DECREF(pyresult);
 946     }
 947     *size = result;
 948     //std::cout << "f(" << uv[0] << "," << uv[1] << ")" << " = " << result << std::endl;
 949     PyGILState_Release(gstate);
 950   }
 951   else {
 952     *size = *((double*)user_data);
 953   }
 954   return STATUS_OK;
 955 }
 956
 957 status_t size_on_edge(integer edge_id, real t, real *size, void *user_data)
 958 {
 959   if (EdgeId2PythonSmp.count(edge_id) != 0){
 960     PyObject * pyresult = NULL;
 961     PyObject* new_stderr = NULL;
 962     assert(Py_IsInitialized());
 963     PyGILState_STATE gstate;
 964     gstate = PyGILState_Ensure();
 965     pyresult = PyObject_CallFunction(EdgeId2PythonSmp[edge_id],"(f)",t);
 966     double result;
 967     if ( pyresult == NULL){
 968       fflush(stderr);
 969       string err_description="";
 970       new_stderr = newPyStdOut(err_description);
 971       PySys_SetObject("stderr", new_stderr);
 972       PyErr_Print();
 973       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
 974       Py_DECREF(new_stderr);
 975       std::cout << "Can't evaluate f(" << t << ")" << " error is " << err_description << std::endl;
 976       result = *((double*)user_data);
 977       }
 978     else {
 979       result = PyFloat_AsDouble(pyresult);
 980       Py_DECREF(pyresult);
 981     }
 982     *size = result;
 983     PyGILState_Release(gstate);
 984   }
 985   else {
 986     *size = *((double*)user_data);
 987   }
 988   return STATUS_OK;
 989 }
 990
 991 status_t size_on_vertex(integer point_id, real *size, void *user_data)
 992 {
 993   if (VertexId2PythonSmp.count(point_id) != 0){
 994     PyObject * pyresult = NULL;
 995     PyObject* new_stderr = NULL;
 996     assert(Py_IsInitialized());
 997     PyGILState_STATE gstate;
 998     gstate = PyGILState_Ensure();
 999     pyresult = PyObject_CallFunction(VertexId2PythonSmp[point_id],"");
1000     double result;
1001     if ( pyresult == NULL){
1002       fflush(stderr);
1003       string err_description="";
1004       new_stderr = newPyStdOut(err_description);
1005       PySys_SetObject("stderr", new_stderr);
1006       PyErr_Print();
1007       PySys_SetObject("stderr", PySys_GetObject("__stderr__"));
1008       Py_DECREF(new_stderr);
1009       std::cout << "Can't evaluate f()" << " error is " << err_description << std::endl;
1010       result = *((double*)user_data);
1011       }
1012     else {
1013       result = PyFloat_AsDouble(pyresult);
1014       Py_DECREF(pyresult);
1015     }
1016     *size = result;
1017     PyGILState_Release(gstate);
1018   }
1019   else {
1020     *size = *((double*)user_data);
1021   }
1022  return STATUS_OK;
1023 }
1024
1025 status_t message_callback(message_t *msg, void *user_data)
1026 {
1027   integer errnumber = 0;
1028   char *desc;
1029   message_get_number(msg, &errnumber);
1030   message_get_description(msg, &desc);
1031   if ( errnumber < 0 ) {
1032     string * error = (string*)user_data;
1033 //   if ( !error->empty() )
1034 //     *error += "\n";
1035     // remove ^A from the tail
1036     int len = strlen( desc );
1037     while (len > 0 && desc[len-1] != '\n')
1038       len--;
1039     error->append( desc, len );
1040   }
1041   else {
1042     cout << desc;
1043   }
1044   return STATUS_OK;
1045 }