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