src/BLSURFPlugin/BLSURFPlugin_BLSURF.cxx

   1 //  BLSURFPlugin : C++ implementation
   2 //
   3 //  Copyright (C) 2006  OPEN CASCADE, CEA/DEN, EDF R&D
   4 //
   5 //  This library is free software; you can redistribute it and/or
   6 //  modify it under the terms of the GNU Lesser General Public
   7 //  License as published by the Free Software Foundation; either
   8 //  version 2.1 of the License.
   9 //
  10 //  This library is distributed in the hope that it will be useful,
  11 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 //  Lesser General Public License for more details.
  14 //
  15 //  You should have received a copy of the GNU Lesser General Public
  16 //  License along with this library; if not, write to the Free Software
  17 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 //
  19 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  20 //
  21 //
  22 // File    : BLSURFPlugin_BLSURF.cxx
  23 // Authors : Francis KLOSS (OCC) & Patrick LAUG (INRIA) & Lioka RAZAFINDRAZAKA (CEA)
  24 //           & Aurelien ALLEAUME (DISTENE)
  25 // Date    : 20/03/2006
  26 // Project : SALOME
  27 //=============================================================================
  28 using namespace std;
  29
  30 #include "BLSURFPlugin_BLSURF.hxx"
  31 #include "BLSURFPlugin_Hypothesis.hxx"
  32
  33 #include <SMESH_Gen.hxx>
  34 #include <SMESH_Mesh.hxx>
  35 #include <SMESH_ControlsDef.hxx>
  36
  37 #include <SMESHDS_Mesh.hxx>
  38 #include <SMDS_MeshElement.hxx>
  39 #include <SMDS_MeshNode.hxx>
  40
  41 #include <utilities.h>
  42
  43 #include <list>
  44 #include <vector>
  45
  46 #include <BRep_Tool.hxx>
  47 #include <TopExp.hxx>
  48 #include <TopExp_Explorer.hxx>
  49 #include <TopoDS.hxx>
  50 #include <NCollection_Map.hxx>
  51 #include <Standard_ErrorHandler.hxx>
  52
  53 extern "C"{
  54 #include <distene/api.h>
  55 }
  56
  57 #include <Geom_Surface.hxx>
  58 #include <Handle_Geom_Surface.hxx>
  59 #include <Geom2d_Curve.hxx>
  60 #include <Handle_Geom2d_Curve.hxx>
  61 #include <Geom_Curve.hxx>
  62 #include <Handle_Geom_Curve.hxx>
  63 #include <TopoDS_Vertex.hxx>
  64 #include <TopoDS_Edge.hxx>
  65 #include <TopoDS_Wire.hxx>
  66 #include <TopoDS_Face.hxx>
  67 #include <TopoDS_Shape.hxx>
  68 #include <gp_Pnt2d.hxx>
  69 #include <TopTools_IndexedMapOfShape.hxx>
  70 #include <BRepTools.hxx>
  71
  72 //=============================================================================
  73 /*!
  74  *
  75  */
  76 //=============================================================================
  77
  78 BLSURFPlugin_BLSURF::BLSURFPlugin_BLSURF(int hypId, int studyId,
  79                                                SMESH_Gen* gen)
  80   : SMESH_2D_Algo(hypId, studyId, gen)
  81 {
  82   MESSAGE("BLSURFPlugin_BLSURF::BLSURFPlugin_BLSURF");
  83
  84   _name = "BLSURF";
  85   _shapeType = (1 << TopAbs_FACE); // 1 bit /shape type
  86   _compatibleHypothesis.push_back("BLSURF_Parameters");
  87   _requireDescretBoundary = false;
  88   _onlyUnaryInput = false;
  89   _hypothesis = NULL;
  90 }
  91
  92 //=============================================================================
  93 /*!
  94  *
  95  */
  96 //=============================================================================
  97
  98 BLSURFPlugin_BLSURF::~BLSURFPlugin_BLSURF()
  99 {
 100   MESSAGE("BLSURFPlugin_BLSURF::~BLSURFPlugin_BLSURF");
 101 }
 102
 103 //=============================================================================
 104 /*!
 105  *
 106  */
 107 //=============================================================================
 108
 109 bool BLSURFPlugin_BLSURF::CheckHypothesis
 110                          (SMESH_Mesh&                          aMesh,
 111                           const TopoDS_Shape&                  aShape,
 112                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 113 {
 114   _hypothesis = NULL;
 115
 116   list<const SMESHDS_Hypothesis*>::const_iterator itl;
 117   const SMESHDS_Hypothesis* theHyp;
 118
 119   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
 120   int nbHyp = hyps.size();
 121   if (!nbHyp)
 122   {
 123     aStatus = SMESH_Hypothesis::HYP_OK;
 124     return true;  // can work with no hypothesis
 125   }
 126
 127   itl = hyps.begin();
 128   theHyp = (*itl); // use only the first hypothesis
 129
 130   string hypName = theHyp->GetName();
 131
 132   if (hypName == "BLSURF_Parameters")
 133   {
 134     _hypothesis = static_cast<const BLSURFPlugin_Hypothesis*> (theHyp);
 135     ASSERT(_hypothesis);
 136     aStatus = SMESH_Hypothesis::HYP_OK;
 137   }
 138   else
 139     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 140
 141   return aStatus == SMESH_Hypothesis::HYP_OK;
 142 }
 143
 144 //=============================================================================
 145 /*!
 146  * Pass parameters to BLSURF
 147  */
 148 //=============================================================================
 149
 150 inline std::string to_string(double d)
 151 {
 152    std::ostringstream o;
 153    o << d;
 154    return o.str();
 155 }
 156
 157 inline std::string to_string(int i)
 158 {
 159    std::ostringstream o;
 160    o << i;
 161    return o.str();
 162 }
 163
 164 void BLSURFPlugin_BLSURF::SetParameters(const BLSURFPlugin_Hypothesis* hyp, blsurf_session_t *bls) {
 165   if (hyp) {
 166     MESSAGE("BLSURFPlugin_BLSURF::SetParameters");
 167     _topology = (int) hyp->GetTopology();
 168     _physicalMesh = (int) hyp->GetPhysicalMesh();
 169     _phySize = hyp->GetPhySize();
 170     _geometricMesh = (int) hyp->GetGeometricMesh();
 171     _angleMeshS = hyp->GetAngleMeshS();
 172     _gradation = hyp->GetGradation();
 173     _quadAllowed = hyp->GetQuadAllowed();
 174     _decimesh = hyp->GetDecimesh();
 175   } else {
 176     MESSAGE("BLSURFPlugin_BLSURF::SetParameters using defaults");
 177     _topology = BLSURFPlugin_Hypothesis::GetDefaultTopology();
 178     _physicalMesh = BLSURFPlugin_Hypothesis::GetDefaultPhysicalMesh();
 179     _phySize = BLSURFPlugin_Hypothesis::GetDefaultPhySize();
 180     _geometricMesh = BLSURFPlugin_Hypothesis::GetDefaultGeometricMesh();
 181     _angleMeshS = BLSURFPlugin_Hypothesis::GetDefaultAngleMeshS();
 182     _gradation = BLSURFPlugin_Hypothesis::GetDefaultGradation();
 183     _quadAllowed = BLSURFPlugin_Hypothesis::GetDefaultQuadAllowed();
 184     _decimesh = BLSURFPlugin_Hypothesis::GetDefaultDecimesh();
 185
 186   }
 187
 188   blsurf_set_param(bls, "topo_points", _topology > 0 ? "1" : "0");
 189   blsurf_set_param(bls, "topo_curves", _topology > 0 ? "1" : "0");
 190   blsurf_set_param(bls, "topo_project", _topology > 0 ? "1" : "0");
 191   blsurf_set_param(bls, "clean_boundary", _topology > 1 ? "1" : "0");
 192   blsurf_set_param(bls, "close_boundary", _topology > 1 ? "1" : "0");
 193   blsurf_set_param(bls, "hphy_flag", to_string(_physicalMesh).c_str());
 194   blsurf_set_param(bls, "hphydef", to_string(_phySize).c_str());
 195   blsurf_set_param(bls, "hgeo_flag", to_string(_geometricMesh).c_str());
 196   blsurf_set_param(bls, "angle_meshs", to_string(_angleMeshS).c_str());
 197   blsurf_set_param(bls, "angle_meshc", to_string(_angleMeshS).c_str());
 198   blsurf_set_param(bls, "gradation", to_string(_gradation).c_str());
 199   //  blsurf_set_param(bls, "patch_independent", to_string(_decimesh).c_str());
 200   blsurf_set_param(bls, "patch_independent", _decimesh ? "1" : "0");
 201   blsurf_set_param(bls, "element",  _quadAllowed ? "q1.0" : "p1");
 202   blsurf_set_param(bls, "verb", "10");
 203 }
 204
 205 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data);
 206 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
 207                   real *duu, real *duv, real *dvv, void *user_data);
 208 status_t message_callback(message_t *msg, void *user_data);
 209
 210 //=============================================================================
 211 /*!
 212  *
 213  */
 214 //=============================================================================
 215
 216 bool BLSURFPlugin_BLSURF::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) {
 217
 218   MESSAGE("BLSURFPlugin_BLSURF::Compute");
 219
 220   if (aShape.ShapeType() == TopAbs_COMPOUND) {
 221     cout << "  the shape is a COMPOUND" << endl;
 222   }
 223   else {
 224     cout << "  the shape is UNKNOWN" << endl;
 225   };
 226
 227   context_t *ctx =  context_new();
 228   context_set_message_callback(ctx, message_callback, &_comment);
 229
 230   cad_t *c = cad_new(ctx);
 231
 232   TopTools_IndexedMapOfShape fmap;
 233   TopTools_IndexedMapOfShape emap;
 234   TopTools_IndexedMapOfShape pmap;
 235   vector<Handle(Geom2d_Curve)> curves;
 236   vector<Handle(Geom_Surface)> surfaces;
 237
 238   fmap.Clear();
 239   emap.Clear();
 240   pmap.Clear();
 241   surfaces.resize(0);
 242   curves.resize(0);
 243
 244   int iface = 0;
 245   for (TopExp_Explorer face_iter(aShape,TopAbs_FACE);face_iter.More();face_iter.Next()) {
 246     TopoDS_Face f=TopoDS::Face(face_iter.Current());
 247     if (fmap.FindIndex(f) > 0)
 248       continue;
 249
 250     fmap.Add(f);
 251     iface++;
 252     surfaces.push_back(BRep_Tool::Surface(f));
 253     cad_face_t *fce = cad_face_new(c, iface, surf_fun, surfaces.back());
 254     cad_face_set_tag(fce, iface);
 255     if(f.Orientation() != TopAbs_FORWARD){
 256       cad_face_set_orientation(fce, CAD_ORIENTATION_REVERSED);
 257     } else {
 258       cad_face_set_orientation(fce, CAD_ORIENTATION_FORWARD);
 259     }
 260
 261     for (TopExp_Explorer edge_iter(f,TopAbs_EDGE);edge_iter.More();edge_iter.Next()) {
 262       TopoDS_Edge e = TopoDS::Edge(edge_iter.Current());
 263       int ic = emap.FindIndex(e);
 264       if (ic <= 0)
 265         ic = emap.Add(e);
 266
 267       double tmin,tmax;
 268       curves.push_back(BRep_Tool::CurveOnSurface(e, f, tmin, tmax));
 269       cad_edge_t *edg = cad_edge_new(fce, ic, tmin, tmax, curv_fun, curves.back());
 270       cad_edge_set_tag(edg, ic);
 271       cad_edge_set_property(edg, EDGE_PROPERTY_SOFT_REQUIRED);
 272
 273       int npts = 0;
 274       int ip1, ip2, *ip;
 275       gp_Pnt2d e0 = curves.back()->Value(tmin);
 276       gp_Pnt ee0 = surfaces.back()->Value(e0.X(), e0.Y());
 277       Standard_Real d1=0,d2=0;
 278       for (TopExp_Explorer ex_edge(e ,TopAbs_VERTEX); ex_edge.More(); ex_edge.Next()) {
 279         TopoDS_Vertex v = TopoDS::Vertex(ex_edge.Current());
 280
 281         ++npts;
 282         if (npts == 1){
 283           ip = &ip1;
 284           d1 = ee0.SquareDistance(BRep_Tool::Pnt(v));
 285         } else {
 286           ip = &ip2;
 287           d2 = ee0.SquareDistance(BRep_Tool::Pnt(v));
 288         }
 289         *ip = pmap.FindIndex(v);
 290         if(*ip <= 0)
 291           *ip = pmap.Add(v);
 292       }
 293       if (npts != 2) {
 294         // should not happen
 295         cout << "An edge does not have 2 extremities." << endl;
 296       } else {
 297         if (d1 < d2)
 298           cad_edge_set_extremities(edg, ip1, ip2);
 299         else
 300           cad_edge_set_extremities(edg, ip2, ip1);
 301       }
 302     } // for edge
 303   } //for face
 304
 305
 306
 307
 308   blsurf_session_t *bls = blsurf_session_new(ctx);
 309   blsurf_data_set_cad(bls, c);
 310
 311   SetParameters(_hypothesis, bls);
 312
 313   cout << endl;
 314   cout << "Beginning of Surface Mesh generation" << endl;
 315   cout << endl;
 316
 317   status_t status = STATUS_ERROR;
 318
 319   try {
 320     OCC_CATCH_SIGNALS;
 321     status = blsurf_compute_mesh(bls);
 322   }
 323   catch ( std::exception& exc ) {
 324 //     if ( !_comment.empty() )
 325 //       _comment += "\n";
 326     _comment += exc.what();
 327   }
 328   catch (Standard_Failure& ex) {
 329 //     if ( !_comment.empty() )
 330 //       _comment += "\n";
 331     _comment += ex.DynamicType()->Name();
 332     if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) {
 333       _comment += ": ";
 334       _comment += ex.GetMessageString();
 335     }
 336   }
 337   catch (...) {
 338     if ( _comment.empty() )
 339       _comment = "Exception in blsurf_compute_mesh()";
 340   }
 341   if ( status != STATUS_OK) {
 342     blsurf_session_delete(bls);
 343     cad_delete(c);
 344     context_delete(ctx);
 345
 346     return error(_comment);
 347     //return false;
 348   }
 349
 350   cout << endl;
 351   cout << "End of Surface Mesh generation" << endl;
 352   cout << endl;
 353
 354   mesh_t *msh;
 355   blsurf_data_get_mesh(bls, &msh);
 356   if(!msh){
 357     blsurf_session_delete(bls);
 358     cad_delete(c);
 359     context_delete(ctx);
 360
 361     return error(_comment);
 362     //return false;
 363   }
 364
 365   integer nv, ne, nt, nq, vtx[4], tag;
 366   real xyz[3];
 367
 368   mesh_get_vertex_count(msh, &nv);
 369   mesh_get_edge_count(msh, &ne);
 370   mesh_get_triangle_count(msh, &nt);
 371   mesh_get_quadrangle_count(msh, &nq);
 372
 373
 374   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 375   SMDS_MeshNode** nodes = new SMDS_MeshNode*[nv+1];
 376   bool* tags = new bool[nv+1];
 377
 378   for(int iv=1;iv<=nv;iv++) {
 379     mesh_get_vertex_coordinates(msh, iv, xyz);
 380     mesh_get_vertex_tag(msh, iv, &tag);
 381     nodes[iv] = meshDS->AddNode(xyz[0], xyz[1], xyz[2]);
 382     // internal point are tagged to zero
 383     if(tag){
 384       meshDS->SetNodeOnVertex(nodes[iv], TopoDS::Vertex(pmap(tag)));
 385       tags[iv] = false;
 386     } else {
 387       tags[iv] = true;
 388     }
 389   }
 390
 391   for(int it=1;it<=ne;it++) {
 392     mesh_get_edge_vertices(msh, it, vtx);
 393     SMDS_MeshEdge* edg = meshDS->AddEdge(nodes[vtx[0]], nodes[vtx[1]]);
 394     mesh_get_edge_tag(msh, it, &tag);
 395
 396     if (tags[vtx[0]]) {
 397       meshDS->SetNodeOnEdge(nodes[vtx[0]], TopoDS::Edge(emap(tag)));
 398       tags[vtx[0]] = false;
 399     };
 400     if (tags[vtx[1]]) {
 401       meshDS->SetNodeOnEdge(nodes[vtx[1]], TopoDS::Edge(emap(tag)));
 402       tags[vtx[1]] = false;
 403     };
 404     meshDS->SetMeshElementOnShape(edg, TopoDS::Edge(emap(tag)));
 405
 406   }
 407
 408   for(int it=1;it<=nt;it++) {
 409     mesh_get_triangle_vertices(msh, it, vtx);
 410     SMDS_MeshFace* tri = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]]);
 411     mesh_get_triangle_tag(msh, it, &tag);
 412     meshDS->SetMeshElementOnShape(tri, TopoDS::Face(fmap(tag)));
 413     if (tags[vtx[0]]) {
 414       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
 415       tags[vtx[0]] = false;
 416     };
 417     if (tags[vtx[1]]) {
 418       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
 419       tags[vtx[1]] = false;
 420     };
 421     if (tags[vtx[2]]) {
 422       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
 423       tags[vtx[2]] = false;
 424     };
 425   }
 426
 427   for(int it=1;it<=nq;it++) {
 428     mesh_get_quadrangle_vertices(msh, it, vtx);
 429     SMDS_MeshFace* quad = meshDS->AddFace(nodes[vtx[0]], nodes[vtx[1]], nodes[vtx[2]], nodes[vtx[3]]);
 430     mesh_get_quadrangle_tag(msh, it, &tag);
 431     meshDS->SetMeshElementOnShape(quad, TopoDS::Face(fmap(tag)));
 432     if (tags[vtx[0]]) {
 433       meshDS->SetNodeOnFace(nodes[vtx[0]], TopoDS::Face(fmap(tag)));
 434       tags[vtx[0]] = false;
 435     };
 436     if (tags[vtx[1]]) {
 437       meshDS->SetNodeOnFace(nodes[vtx[1]], TopoDS::Face(fmap(tag)));
 438       tags[vtx[1]] = false;
 439     };
 440     if (tags[vtx[2]]) {
 441       meshDS->SetNodeOnFace(nodes[vtx[2]], TopoDS::Face(fmap(tag)));
 442       tags[vtx[2]] = false;
 443     };
 444     if (tags[vtx[3]]) {
 445       meshDS->SetNodeOnFace(nodes[vtx[3]], TopoDS::Face(fmap(tag)));
 446       tags[vtx[3]] = false;
 447     };
 448   }
 449
 450   delete nodes;
 451
 452   /* release the mesh object */
 453   blsurf_data_regain_mesh(bls, msh);
 454
 455   /* clean up everything */
 456   blsurf_session_delete(bls);
 457   cad_delete(c);
 458
 459   context_delete(ctx);
 460
 461   return true;
 462 }
 463
 464 //=============================================================================
 465 /*!
 466  *
 467  */
 468 //=============================================================================
 469
 470 ostream & BLSURFPlugin_BLSURF::SaveTo(ostream & save)
 471 {
 472   return save;
 473 }
 474
 475 //=============================================================================
 476 /*!
 477  *
 478  */
 479 //=============================================================================
 480
 481 istream & BLSURFPlugin_BLSURF::LoadFrom(istream & load)
 482 {
 483   return load;
 484 }
 485
 486 //=============================================================================
 487 /*!
 488  *
 489  */
 490 //=============================================================================
 491
 492 ostream & operator << (ostream & save, BLSURFPlugin_BLSURF & hyp)
 493 {
 494   return hyp.SaveTo( save );
 495 }
 496
 497 //=============================================================================
 498 /*!
 499  *
 500  */
 501 //=============================================================================
 502
 503 istream & operator >> (istream & load, BLSURFPlugin_BLSURF & hyp)
 504 {
 505   return hyp.LoadFrom( load );
 506 }
 507
 508 status_t curv_fun(real t, real *uv, real *dt, real *dtt, void *user_data)
 509 {
 510   const Geom2d_Curve*pargeo = (const Geom2d_Curve*) user_data;
 511
 512   if (uv){
 513     gp_Pnt2d P;
 514     P=pargeo->Value(t);
 515     uv[0]=P.X(); uv[1]=P.Y();
 516   }
 517
 518   if(dt) {
 519     gp_Vec2d V1;
 520     V1=pargeo->DN(t,1);
 521     dt[0]=V1.X(); dt[1]=V1.Y();
 522   }
 523
 524   if(dtt){
 525     gp_Vec2d V2;
 526     V2=pargeo->DN(t,2);
 527     dtt[0]=V2.X(); dtt[1]=V2.Y();
 528   }
 529
 530   return 0;
 531 }
 532
 533 status_t surf_fun(real *uv, real *xyz, real*du, real *dv,
 534                   real *duu, real *duv, real *dvv, void *user_data)
 535 {
 536   const Geom_Surface* geometry = (const Geom_Surface*) user_data;
 537
 538   if(xyz){
 539    gp_Pnt P;
 540    P=geometry->Value(uv[0],uv[1]);   // S.D0(U,V,P);
 541    xyz[0]=P.X(); xyz[1]=P.Y(); xyz[2]=P.Z();
 542   }
 543
 544   if(du && dv){
 545     gp_Pnt P;
 546     gp_Vec D1U,D1V;
 547
 548     geometry->D1(uv[0],uv[1],P,D1U,D1V);
 549     du[0]=D1U.X(); du[1]=D1U.Y(); du[2]=D1U.Z();
 550     dv[0]=D1V.X(); dv[1]=D1V.Y(); dv[2]=D1V.Z();
 551   }
 552
 553   if(duu && duv && dvv){
 554     gp_Pnt P;
 555     gp_Vec D1U,D1V;
 556     gp_Vec D2U,D2V,D2UV;
 557
 558     geometry->D2(uv[0],uv[1],P,D1U,D1V,D2U,D2V,D2UV);
 559     duu[0]=D2U.X(); duu[1]=D2U.Y(); duu[2]=D2U.Z();
 560     duv[0]=D2UV.X(); duv[1]=D2UV.Y(); duv[2]=D2UV.Z();
 561     dvv[0]=D2V.X(); dvv[1]=D2V.Y(); dvv[2]=D2V.Z();
 562   }
 563
 564   return 0;
 565 }
 566
 567 status_t message_callback(message_t *msg, void *user_data)
 568 {
 569   integer errnumber = 0;
 570   message_get_number(msg, &errnumber);
 571   if ( errnumber < 0 ) {
 572     char *desc;
 573     string * error = (string*)user_data;
 574     message_get_description(msg, &desc);
 575 //   if ( !error->empty() )
 576 //     *error += "\n";
 577     // remove ^A from the tail
 578     int len = strlen( desc );
 579     while (len > 0 && desc[len-1] != '\n')
 580       len--;
 581     error->append( desc, len );
 582   }
 583   return STATUS_OK;
 584 }