src/SMESH/SMESH_Homard.cxx

   1 // SMESH HOMARD : implementation of SMESHHOMARD idl descriptions
   2 //
   3 // Copyright (C) 2011-2021  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, or (at your option) any later version.
   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.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 //
  21
  22 #include "SMESH_Homard.hxx"
  23
  24 #include <Utils_SALOME_Exception.hxx>
  25 #include <utilities.h>
  26
  27 #include <iostream>
  28 #include <sstream>
  29 #include <cstdlib>
  30 #include <sys/stat.h>
  31
  32 #ifndef WIN32
  33 #include <unistd.h>
  34 #else
  35 #include <direct.h>
  36 #endif
  37
  38 ////
  39
  40 #include <MCAuto.hxx>
  41 #include <MEDCouplingMemArray.hxx>
  42 #include <MEDFileMesh.hxx>
  43
  44 #include <XAO_Xao.hxx>
  45 #include <XAO_BrepGeometry.hxx>
  46 #include <XAO_Group.hxx>
  47
  48 #include <stdexcept>
  49 #include <cstdio>
  50 #include <cstdlib>
  51 #include <list>
  52 #include <limits>
  53
  54 #include <fcntl.h>
  55 #include <boost/filesystem.hpp>
  56
  57 #include <OSD_Parallel.hxx>
  58 #include <BRepAdaptor_Curve.hxx>
  59 #include <BRepBndLib.hxx>
  60 #include <BRepTopAdaptor_FClass2d.hxx>
  61 #include <BRep_Tool.hxx>
  62 #include <Bnd_Box.hxx>
  63 #include <TopExp.hxx>
  64 #include <TopExp_Explorer.hxx>
  65 #include <TopTools_IndexedMapOfShape.hxx>
  66 #include <TopoDS.hxx>
  67 #include <TopoDS_Edge.hxx>
  68 #include <TopoDS_Face.hxx>
  69 #include <TopoDS_Vertex.hxx>
  70 #include <TopoDS_Shape.hxx>
  71 #include <ElCLib.hxx>
  72 #include <ElSLib.hxx>
  73 #include <GCPnts_UniformDeflection.hxx>
  74 #include <GeomAdaptor_Curve.hxx>
  75 #include <GeomLib_IsPlanarSurface.hxx>
  76 #include <ShapeAnalysis_Curve.hxx>
  77 #include <ShapeAnalysis_Surface.hxx>
  78 #include <gp_Circ.hxx>
  79 #include <gp_Cylinder.hxx>
  80 #include <gp_Dir.hxx>
  81 #include <gp_Pln.hxx>
  82 #include <gp_Pnt.hxx>
  83 #include <gp_Sphere.hxx>
  84 #include <gp_Vec.hxx>
  85
  86 namespace boofs = boost::filesystem;
  87
  88 ////
  89
  90 namespace SMESHHOMARDImpl
  91 {
  92
  93   std::string SEPARATOR = "|" ;
  94
  95   /*!
  96     \brief Read next chunk of data from the string
  97     \internal
  98
  99     The function tries to read next chunk of the data from the input string \a str.
 100     The parameter \a start specifies the start position of next chunk. If the operation
 101     read the chunk successfully, after its completion this parameter will refer to the
 102     start position of the next chunk. The function returns resulting chunk as a string.
 103     The status of the operation is returned via \a ok parameter.
 104
 105     \param str source data stream string
 106     \param start start position to get next chunk
 107     \param ok in this variable the status of the chunk reading operation is returned
 108     \return next chunk read from the string
 109   */
 110   static std::string getNextChunk( const std::string& str, std::string::size_type& start, bool& ok )
 111   {
 112     std::string chunk = "";
 113     ok = false;
 114     if ( start <= str.size() ) {
 115       std::string::size_type end = str.find( separator(), start );
 116       chunk = str.substr( start, end == std::string::npos ? std::string::npos : end-start );
 117       start = end == std::string::npos ? str.size()+1 : end + separator().size();
 118       ok = true;
 119     }
 120     return chunk;
 121   }
 122
 123   /*!
 124     \brief Get persistence signature
 125     \param type persistence entity type
 126     \return persistence signature
 127   */
 128   std::string GetSignature( SignatureType type )
 129   {
 130     std::string signature = "";
 131     switch ( type ) {
 132     case Case:       signature = "CASE"; break;
 133     case Zone:       signature = "ZONE"; break;
 134     case Hypothesis: signature = "HYPO"; break;
 135     case Iteration:  signature = "ITER"; break;
 136     case Boundary:   signature = "BOUNDARY"; break;
 137     default: break;
 138     }
 139     signature += separator();
 140     return signature;
 141   }
 142
 143   /*!
 144     \brief Get data separator
 145     \return string that is used to separate data entities in the stream
 146   */
 147   std::string separator()
 148   {
 149     return SEPARATOR ;
 150   }
 151
 152 // =======================
 153 // 1.1. Case
 154 // =======================
 155   /*!
 156     \brief Dump case to the string
 157     \param cas case being dumped
 158     \return string representation of the case
 159   */
 160   std::string Dump( const HOMARD_Cas& cas )
 161   {
 162     std::stringstream os;
 163     std::string saux ;
 164     // ...
 165     MESSAGE( ". Sauvegarde du cas "<<cas.GetName());
 166     os << cas.GetName();
 167     os << separator() << cas.GetDirName();
 168     os << separator() << cas.GetConfType();
 169
 170     std::vector<double> coor = cas.GetBoundingBox();
 171     os << separator() << coor.size();
 172     for ( unsigned int i = 0; i < coor.size(); i++ )
 173           os << separator() << coor[i];
 174
 175     std::list<std::string> ListString = cas.GetIterations();
 176     os << separator() << ListString.size();
 177     std::list<std::string>::const_iterator it;
 178     for ( it = ListString.begin(); it != ListString.end(); ++it )
 179           os << separator() << *it;
 180
 181     ListString = cas.GetGroups();
 182     os << separator() << ListString.size();
 183     for ( it = ListString.begin(); it != ListString.end(); ++it )
 184          os << separator() << *it;
 185     ListString = cas.GetBoundaryGroup();
 186     os << separator() << ListString.size();
 187     for ( it = ListString.begin(); it != ListString.end(); ++it )
 188          os << separator() << *it;
 189
 190     os << separator() << 0; //cas.GetPyram()
 191
 192     saux = os.str();
 193 //     MESSAGE( ". Fin avec "<<saux);
 194     return saux ;
 195   }
 196 //
 197 // ==============
 198 // 1.2. Iteration
 199 // ==============
 200 //
 201   /*!
 202     \brief Dump iteration to the string
 203     \param iteration iteration being dumped
 204     \return string representation of the iteration
 205   */
 206   std::string Dump( const HOMARD_Iteration& iteration )
 207   {
 208     std::stringstream os;
 209     std::string saux ;
 210     // ...
 211     MESSAGE( ". Sauvegarde de l'iteration "<<iteration.GetName());
 212     os << iteration.GetName();
 213     os << separator() << iteration.GetState();
 214     os << separator() << iteration.GetNumber();
 215     os << separator() << iteration.GetMeshFile();
 216     os << separator() << iteration.GetLogFile();
 217     os << separator() << iteration.GetMeshName();
 218     os << separator() << iteration.GetFieldFile();
 219     os << separator() << iteration.GetTimeStep();
 220     os << separator() << iteration.GetRank();
 221     os << separator() << iteration.GetIterParentName();
 222     //
 223     std::list<std::string> ListString = iteration.GetIterations();
 224     os << separator() << ListString.size();
 225     std::list<std::string>::const_iterator it;
 226     for ( it = ListString.begin(); it != ListString.end(); ++it )
 227       os << separator() << *it;
 228
 229     os << separator() << iteration.GetHypoName();
 230     os << separator() << iteration.GetCaseName();
 231     os << separator() << iteration.GetDirNameLoc();
 232
 233     saux = os.str();
 234 //     MESSAGE( ". Fin avec "<<saux);
 235     return saux ;
 236   }
 237 //
 238 // ==============
 239 // 1.3. hypothese
 240 // ==============
 241   /*!
 242     \brief Dump hypothesis to the string
 243     \param hypothesis hypothesis being dumped
 244     \return string representation of the hypothesis
 245   */
 246   std::string Dump( const HOMARD_Hypothesis& hypothesis )
 247   {
 248     std::stringstream os;
 249     std::string saux ;
 250     // ...
 251     MESSAGE( ". Sauvegarde de l'hypothese "<<hypothesis.GetName());
 252     os << hypothesis.GetName();
 253     os << separator() << hypothesis.GetCaseCreation();
 254     os << separator() << hypothesis.GetAdapType();
 255     os << separator() << hypothesis.GetRefinType();
 256     os << separator() << hypothesis.GetUnRefType();
 257     os << separator() << hypothesis.GetFieldName();
 258     os << separator() << hypothesis.GetRefinThrType();
 259     os << separator() << hypothesis.GetThreshR();
 260     os << separator() << hypothesis.GetUnRefThrType();
 261     os << separator() << hypothesis.GetThreshC();
 262     os << separator() << hypothesis.GetUseField();
 263     os << separator() << hypothesis.GetUseComp();
 264     os << separator() << hypothesis.GetTypeFieldInterp();
 265
 266     std::list<std::string> ListString = hypothesis.GetIterations();
 267     std::list<std::string>::const_iterator it;
 268     os << separator() << ListString.size();
 269     for ( it = ListString.begin(); it != ListString.end(); ++it )
 270          os << separator() << *it;
 271
 272     ListString = hypothesis.GetZones();
 273     os << separator() << ListString.size();
 274     for ( it = ListString.begin(); it != ListString.end(); ++it )
 275           os << separator() << *it;
 276
 277     ListString = hypothesis.GetComps();
 278     os << separator() << ListString.size();
 279     for ( it = ListString.begin(); it != ListString.end(); ++it )
 280          os << separator() << *it;
 281
 282     ListString = hypothesis.GetGroups();
 283     os << separator() << ListString.size();
 284     for ( it = ListString.begin(); it != ListString.end(); ++it )
 285           os << separator() << *it;
 286
 287     ListString = hypothesis.GetFieldInterps();
 288     os << separator() << ListString.size();
 289     for ( it = ListString.begin(); it != ListString.end(); ++it )
 290           os << separator() << *it;
 291
 292     os << separator() << hypothesis.GetNivMax();
 293     os << separator() << hypothesis.GetDiamMin();
 294     os << separator() << hypothesis.GetAdapInit();
 295     os << separator() << hypothesis.GetExtraOutput();
 296
 297     saux = os.str();
 298 //     MESSAGE( ". Fin avec "<<saux);
 299     return saux ;
 300   }
 301 //
 302 // ==============================
 303 // 1.5. Archivage d'une frontiere
 304 // ==============================
 305
 306   /*!
 307     \brief Dump boundary to the string
 308     \param boundary boundary being dumped
 309     \return string representation of the boundary
 310   */
 311   std::string Dump( const HOMARD_Boundary& boundary )
 312   {
 313     std::stringstream os;
 314     std::string saux ;
 315     MESSAGE( ". Sauvegarde de la frontiere "<<boundary.GetName());
 316
 317     int BoundaryType = boundary.GetType() ;
 318
 319     os << boundary.GetName() ;
 320     os << separator() << BoundaryType ;
 321     os << separator() << boundary.GetCaseCreation() ;
 322
 323     if ( BoundaryType == -1 )
 324     {
 325       os << separator() << boundary.GetDataFile();
 326     }
 327     else if ( BoundaryType == 0 )
 328     {
 329       os << separator() << boundary.GetMeshName();
 330       os << separator() << boundary.GetDataFile();
 331     }
 332     else {
 333       std::vector<double> coor = boundary.GetCoords() ;
 334       for ( unsigned int i = 0; i < coor.size(); i++ )
 335             os << separator() << coor[i];
 336       std::vector<double> limit = boundary.GetLimit();
 337       for ( unsigned int i = 0; i < limit.size(); i++ )
 338             os << separator() << limit[i];
 339     }
 340
 341     std::list<std::string> ListString = boundary.GetGroups();
 342     std::list<std::string>::const_iterator it;
 343     os << separator() << ListString.size();
 344     for ( it = ListString.begin(); it != ListString.end(); ++it )
 345           os << separator() << *it;
 346
 347     saux = os.str();
 348 //     MESSAGE( ". Fin avec "<<saux);
 349     return saux ;
 350   }
 351
 352 //
 353 // 2. Restauration des objets
 354 // ==========================
 355 // 2.1. Case
 356 // ==========================
 357 //
 358   /*!
 359     \brief Restore case from the string
 360     \param cas case being restored
 361     \param stream string representation of the case
 362     \return \c true if case is correctly restored or \c false otherwise
 363   */
 364   bool Restore( HOMARD_Cas& cas, const std::string& stream )
 365   {
 366     MESSAGE( ". Restoration du cas ");
 367     std::string::size_type start = 0;
 368     std::string chunk, chunkNext;
 369     bool ok;
 370     // ...
 371     chunk = getNextChunk( stream, start, ok );
 372     if ( !ok ) return false;
 373     cas.SetName( chunk.c_str() );
 374
 375     chunk = getNextChunk( stream, start, ok );
 376     if ( !ok ) return false;
 377     cas.SetDirName( chunk.c_str() );
 378
 379     chunk = getNextChunk( stream, start, ok );
 380     if ( !ok ) return false;
 381     cas.SetConfType( atoi( chunk.c_str() ) );
 382
 383     chunk = getNextChunk( stream, start, ok );
 384     if ( !ok ) return false;
 385
 386     int size = atoi( chunk.c_str() );
 387     std::vector<double> boite;
 388     boite.resize( size );
 389     for ( int i = 0; i < size; i++ ) {
 390       chunk = getNextChunk( stream, start, ok );
 391       if ( !ok ) return false;
 392       boite[i] = strtod( chunk.c_str(), 0 );
 393     }
 394     cas.SetBoundingBox( boite );
 395
 396     chunk = getNextChunk( stream, start, ok );
 397     if ( !ok ) return false;
 398
 399     size = atoi( chunk.c_str() );
 400     for ( int i = 0; i < size; i++ ) {
 401       chunk = getNextChunk( stream, start, ok );
 402       if ( !ok ) return false;
 403       cas.AddIteration( chunk.c_str() );
 404     }
 405
 406     chunk = getNextChunk( stream, start, ok );
 407     if ( !ok ) return false;
 408     size = atoi( chunk.c_str() );
 409     for ( int i = 0; i < size; i++ )
 410     {
 411       chunk = getNextChunk( stream, start, ok );
 412       if ( !ok ) return false;
 413       cas.AddGroup( chunk.c_str() );
 414     }
 415
 416     chunk = getNextChunk( stream, start, ok );
 417     if ( !ok ) return false;
 418     size = atoi( chunk.c_str() );
 419     for ( int i = 0; i < size; i++ ) {
 420       chunk = getNextChunk( stream, start, ok );
 421       if ( !ok ) return false;
 422       i++;
 423       chunkNext = getNextChunk( stream, start, ok );
 424       if ( !ok ) return false;
 425       cas.AddBoundaryGroup( chunk.c_str(), chunkNext.c_str() );
 426     }
 427
 428     chunk = getNextChunk( stream, start, ok );
 429     if ( !ok ) return false;
 430     //cas.SetPyram( atoi( chunk.c_str() ) );
 431
 432     return true;
 433   }
 434 //
 435 // ==============
 436 // 2.2. Iteration
 437 // ==============
 438   /*!
 439     \brief Restore iteration from the string
 440     \param iteration iteration being restored
 441     \param stream string representation of the iteration
 442     \return \c true if iteration is correctly restored or \c false otherwise
 443   */
 444   bool Restore( HOMARD_Iteration& iteration, const std::string& stream )
 445   {
 446     std::string::size_type start = 0;
 447     std::string chunk;
 448     bool ok;
 449     chunk = getNextChunk( stream, start, ok );
 450     if ( !ok ) return false;
 451
 452     iteration.SetName( chunk.c_str() );
 453     chunk = getNextChunk( stream, start, ok );
 454     if ( !ok ) return false;
 455     iteration.SetState( atoi( chunk.c_str() ) );
 456     chunk = getNextChunk( stream, start, ok );
 457     if ( !ok ) return false;
 458     iteration.SetNumber( atoi( chunk.c_str() ) );
 459     chunk = getNextChunk( stream, start, ok );
 460     if ( !ok ) return false;
 461     iteration.SetMeshFile( chunk.c_str() );
 462     chunk = getNextChunk( stream, start, ok );
 463     if ( !ok ) return false;
 464     iteration.SetLogFile( chunk.c_str() );
 465     chunk = getNextChunk( stream, start, ok );
 466     if ( !ok ) return false;
 467     iteration.SetMeshName( chunk.c_str() );
 468     chunk = getNextChunk( stream, start, ok );
 469     if ( !ok ) return false;
 470     iteration.SetFieldFile( chunk.c_str() );
 471     // .
 472     int timestep, rank;
 473     chunk = getNextChunk( stream, start, ok );
 474     if ( !ok ) return false;
 475     timestep = atoi( chunk.c_str() );
 476     chunk = getNextChunk( stream, start, ok );
 477     if ( !ok ) return false;
 478     rank = atoi( chunk.c_str() );
 479     iteration.SetTimeStepRank( timestep, rank );
 480     chunk = getNextChunk( stream, start, ok );
 481     if ( !ok ) return false;
 482     iteration.SetIterParentName( chunk.c_str() );
 483     //
 484     chunk = getNextChunk( stream, start, ok );
 485     if ( !ok ) return false;
 486     int size = atoi( chunk.c_str() );
 487     for ( int i = 0; i < size; i++ ) {
 488       chunk = getNextChunk( stream, start, ok );
 489       if ( !ok ) return false;
 490       iteration.LinkNextIteration( chunk.c_str() );
 491     }
 492     //
 493     chunk = getNextChunk( stream, start, ok );
 494     if ( !ok ) return false;
 495     iteration.SetHypoName( chunk.c_str() );
 496     chunk = getNextChunk( stream, start, ok );
 497     if ( !ok ) return false;
 498     iteration.SetCaseName( chunk.c_str() );
 499     chunk = getNextChunk( stream, start, ok );
 500     if ( !ok ) return false;
 501     iteration.SetDirNameLoc( chunk.c_str() );
 502     return true;
 503   }
 504
 505 //
 506 // ==============
 507 // 2.3. hypothese
 508 // ==============
 509   /*!
 510     \brief Restore hypothesis from the string
 511     \param hypothesis hypothesis being restored
 512     \param stream string representation of the hypothesis
 513     \return \c true if hypothesis is correctly restored or \c false otherwise
 514   */
 515   bool Restore( HOMARD_Hypothesis& hypothesis, const std::string& stream )
 516   {
 517     std::string::size_type start = 0;
 518     std::string chunk, chunkNext;
 519     bool ok;
 520
 521     chunk = getNextChunk( stream, start, ok );
 522     if ( !ok ) return false;
 523     hypothesis.SetName( chunk.c_str() );
 524
 525     chunk = getNextChunk( stream, start, ok );
 526     if ( !ok ) return false;
 527     hypothesis.SetCaseCreation( chunk.c_str() );
 528
 529     chunk = getNextChunk( stream, start, ok );
 530     if ( !ok ) return false;
 531     hypothesis.SetAdapType( atoi( chunk.c_str() ) );
 532
 533     chunk = getNextChunk( stream, start, ok );
 534     if ( !ok ) return false;
 535     int typeraff = atoi( chunk.c_str() );
 536     chunk = getNextChunk( stream, start, ok );
 537     if ( !ok ) return false;
 538     int typedera = atoi( chunk.c_str() );
 539     hypothesis.SetRefinTypeDera( typeraff, typedera );
 540
 541     chunk = getNextChunk( stream, start, ok );
 542     if ( !ok ) return false;
 543     hypothesis.SetField( chunk.c_str() );
 544
 545     chunk = getNextChunk( stream, start, ok );
 546     if ( !ok ) return false;
 547     int typethr = atoi( chunk.c_str() );
 548     chunk = getNextChunk( stream, start, ok );
 549     if ( !ok ) return false;
 550     double threshr = strtod( chunk.c_str(), 0 );
 551     hypothesis.SetRefinThr( typethr, threshr );
 552
 553     chunk = getNextChunk( stream, start, ok );
 554     if ( !ok ) return false;
 555     int typethc = atoi( chunk.c_str() );
 556     chunk = getNextChunk( stream, start, ok );
 557     if ( !ok ) return false;
 558     double threshc = strtod( chunk.c_str(), 0 );
 559     hypothesis.SetUnRefThr( typethc, threshc );
 560
 561     chunk = getNextChunk( stream, start, ok );
 562     if ( !ok ) return false;
 563     hypothesis.SetUseField(atoi(chunk.c_str()));
 564
 565     chunk = getNextChunk( stream, start, ok );
 566     if ( !ok ) return false;
 567     hypothesis.SetUseComp(atoi(chunk.c_str()));
 568
 569     chunk = getNextChunk( stream, start, ok );
 570     if ( !ok ) return false;
 571     hypothesis.SetTypeFieldInterp(atoi(chunk.c_str()));
 572
 573     chunk = getNextChunk( stream, start, ok );
 574     if ( !ok ) return false;
 575     int size = atoi( chunk.c_str() );
 576     for ( int i = 0; i < size; i++ ) {
 577       chunk = getNextChunk( stream, start, ok );
 578       if ( !ok ) return false;
 579       hypothesis.LinkIteration( chunk.c_str() );
 580     }
 581
 582     chunk = getNextChunk( stream, start, ok );
 583     if ( !ok ) return false;
 584     size = atoi( chunk.c_str() );
 585     for ( int i = 0; i < size; i++ ) {
 586       chunk = getNextChunk( stream, start, ok );
 587       if ( !ok ) return false;
 588       i++;
 589       chunkNext = getNextChunk( stream, start, ok );
 590       int typeuse = atoi( chunkNext.c_str() );
 591       if ( !ok ) return false;
 592       hypothesis.AddZone( chunk.c_str(), typeuse );
 593     }
 594
 595     chunk = getNextChunk( stream, start, ok );
 596     if ( !ok ) return false;
 597     size = atoi( chunk.c_str() );
 598     for ( int i = 0; i < size; i++ ) {
 599       chunk = getNextChunk( stream, start, ok );
 600       if ( !ok ) return false;
 601       hypothesis.AddComp( chunk.c_str() );
 602     }
 603
 604     chunk = getNextChunk( stream, start, ok );
 605     if ( !ok ) return false;
 606     size = atoi( chunk.c_str() );
 607     for ( int i = 0; i < size; i++ ) {
 608       chunk = getNextChunk( stream, start, ok );
 609       if ( !ok ) return false;
 610       hypothesis.AddGroup( chunk.c_str() );
 611     }
 612
 613     chunk = getNextChunk( stream, start, ok );
 614     if ( !ok ) return false;
 615     size = atoi( chunk.c_str() );
 616     for ( int i = 0; i < size; i++ ) {
 617       chunk = getNextChunk( stream, start, ok );
 618       if ( !ok ) return false;
 619       i++;
 620       chunkNext = getNextChunk( stream, start, ok );
 621       int TypeInterp = atoi( chunkNext.c_str() );
 622       if ( !ok ) return false;
 623       hypothesis.AddFieldInterpType( chunk.c_str(), TypeInterp );
 624     }
 625
 626     chunk = getNextChunk( stream, start, ok );
 627     if ( !ok ) return false;
 628     hypothesis.SetNivMax( atoi( chunk.c_str() ) );
 629
 630     chunk = getNextChunk( stream, start, ok );
 631     if ( !ok ) return false;
 632     hypothesis.SetDiamMin( strtod( chunk.c_str(), 0 ) );
 633
 634     chunk = getNextChunk( stream, start, ok );
 635     if ( !ok ) return false;
 636     hypothesis.SetAdapInit( strtod( chunk.c_str(), 0 ) );
 637
 638     chunk = getNextChunk( stream, start, ok );
 639     if ( !ok ) return false;
 640     hypothesis.SetExtraOutput( strtod( chunk.c_str(), 0 ) );
 641
 642     return true;
 643   }
 644
 645 //
 646 // =================================
 647 // 2.5. Restauration d'une frontiere
 648 // =================================
 649
 650   /*!
 651     \brief Restore boundary from the string
 652     \param boundary boundary being restored
 653     \param stream string representation of the boundary
 654     \return \c true if the boundary is correctly restored or \c false otherwise
 655   */
 656   bool Restore( HOMARD_Boundary& boundary, const std::string& stream )
 657   {
 658     std::string::size_type start = 0;
 659     std::string chunk;
 660     bool ok;
 661
 662     chunk = getNextChunk( stream, start, ok );
 663     if ( !ok ) return false;
 664     boundary.SetName( chunk.c_str() );
 665
 666     chunk = getNextChunk( stream, start, ok );
 667     if ( !ok ) return false;
 668     int BoundaryType = atoi( chunk.c_str() ) ;
 669     boundary.SetType( BoundaryType );
 670
 671     chunk = getNextChunk( stream, start, ok );
 672     if ( !ok ) return false;
 673     boundary.SetCaseCreation( chunk.c_str() );
 674
 675     // Si analytique, les coordonnees des frontieres : le nombre depend du type
 676     // Si discret, le maillage
 677     // Si CAO, la géométrie
 678     int lgcoords ;
 679     if ( BoundaryType == -1 ) { lgcoords = -1 ; }
 680     else if ( BoundaryType == 1 ) { lgcoords = 7 ; }
 681     else if ( BoundaryType == 2 ) { lgcoords = 4 ; }
 682     else { lgcoords = 0 ; }
 683 //
 684     if ( lgcoords == -1 )
 685     {
 686       chunk = getNextChunk( stream, start, ok );
 687       if ( !ok ) return false;
 688       boundary.SetDataFile( chunk.c_str() );
 689     }
 690     else if ( lgcoords == 0 )
 691     {
 692       chunk = getNextChunk( stream, start, ok );
 693       if ( !ok ) return false;
 694       boundary.SetMeshName( chunk.c_str() );
 695
 696       chunk = getNextChunk( stream, start, ok );
 697       if ( !ok ) return false;
 698       boundary.SetDataFile( chunk.c_str() );
 699     }
 700     else
 701     { std::vector<double> coords;
 702       coords.resize( lgcoords );
 703       for ( int i = 0; i < lgcoords; i++ ) {
 704         chunk = getNextChunk( stream, start, ok );
 705         if ( !ok ) return false;
 706         coords[i] = strtod( chunk.c_str(), 0 );
 707       }
 708       if ( BoundaryType == 1 )
 709       { boundary.SetCylinder(coords[0],coords[1],coords[2],coords[3],coords[4],coords[5],coords[6]); }
 710       else if ( BoundaryType == 2 )
 711       { boundary.SetSphere( coords[0], coords[1], coords[2], coords[3]); }
 712       else if ( BoundaryType == 3 )
 713       { boundary.SetConeA( coords[0], coords[1], coords[2], coords[3], coords[4], coords[5], coords[6]); }
 714       else if ( BoundaryType == 4 )
 715       { boundary.SetConeR( coords[0], coords[1], coords[2], coords[3], coords[4], coords[5], coords[6], coords[7]); }
 716       // Remarque : la taille de coords est suffisante pour les limites
 717       for ( int i = 0; i < 3; i++ ) {
 718         chunk = getNextChunk( stream, start, ok );
 719         if ( !ok ) return false;
 720         coords[i] = strtod( chunk.c_str(), 0 );
 721       }
 722       boundary.SetLimit( coords[0], coords[1], coords[2]);
 723     }
 724     // Les groupes
 725     chunk = getNextChunk( stream, start, ok );
 726     if ( !ok ) return false;
 727     int size = atoi( chunk.c_str() );
 728     for ( int i = 0; i < size; i++ ) {
 729       chunk = getNextChunk( stream, start, ok );
 730       if ( !ok ) return false;
 731       boundary.AddGroup( chunk.c_str() );
 732     }
 733
 734     return true;
 735   }
 736
 737 //=============================================================================
 738 /*!
 739  *  default constructor:
 740  */
 741 //=============================================================================
 742 HOMARD_Boundary::HOMARD_Boundary():
 743   _Name( "" ),_Type( 1 ),
 744   _Xmin( 0 ), _Xmax( 0 ), _Ymin( 0 ), _Ymax( 0 ), _Zmin( 0 ), _Zmax( 0 ),
 745   _Xaxe( 0 ), _Yaxe( 0 ), _Zaxe( 0 ),
 746   _Xcentre( 0 ), _Ycentre( 0 ), _Zcentre( 0 ), _rayon( 0 ),
 747   _Xincr( 0 ), _Yincr( 0 ), _Zincr( 0 )
 748 {
 749   MESSAGE("HOMARD_Boundary");
 750 }
 751
 752 //=============================================================================
 753 HOMARD_Boundary::~HOMARD_Boundary()
 754 {
 755   MESSAGE("~HOMARD_Boundary");
 756 }
 757 //=============================================================================
 758 //=============================================================================
 759 // Generalites
 760 //=============================================================================
 761 //=============================================================================
 762 void HOMARD_Boundary::SetName( const char* Name )
 763 {
 764   _Name = std::string( Name );
 765 }
 766 //=============================================================================
 767 std::string HOMARD_Boundary::GetName() const
 768 {
 769   return _Name;
 770 }
 771 //=============================================================================
 772 std::string HOMARD_Boundary::GetDumpPython() const
 773 {
 774   std::ostringstream aScript;
 775   switch (_Type) {
 776     case -1:
 777     {
 778       aScript << _Name << " = smeshhomard.CreateBoundaryCAO(\"" << _Name << "\", ";
 779       aScript << "\"" << _DataFile << "\")\n";
 780       break ;
 781     }
 782     case 0:
 783     {
 784       aScript << _Name << " = smeshhomard.CreateBoundaryDi(\"" << _Name << "\", ";
 785       aScript << "\"" << _MeshName << "\", ";
 786       aScript << "\"" << _DataFile << "\")\n";
 787       break ;
 788     }
 789     case 1:
 790     {
 791       aScript << _Name << " = smeshhomard.CreateBoundaryCylinder(\"" << _Name << "\", ";
 792       aScript << _Xcentre << ", " << _Ycentre << ", " << _Zcentre << ", " << _Xaxe << ", " << _Yaxe << ", " << _Zaxe << ", " << _rayon << ")\n";
 793       break ;
 794     }
 795     case 2:
 796     {
 797       aScript << _Name << " = smeshhomard.CreateBoundarySphere(\"" << _Name << "\", ";
 798       aScript << _Xcentre << ", " << _Ycentre << ", " << _Zcentre << ", " << _rayon << ")\n";
 799       break ;
 800     }
 801     case 3:
 802     {
 803       aScript << _Name << " = smeshhomard.CreateBoundaryConeA(\"" << _Name << "\", ";
 804       aScript << _Xaxe << ", " << _Yaxe << ", " << _Zaxe << ", " << _Angle << ", " << _Xcentre << ", " << _Ycentre << ", " << _Zcentre << ")\n";
 805       break ;
 806     }
 807     case 4:
 808     {
 809       aScript << _Name << " = smeshhomard.CreateBoundaryConeR(\"" << _Name << "\", ";
 810       aScript << _Xcentre1 << ", " << _Ycentre1 << ", " << _Zcentre1 << ", " << _Rayon1 << ", " << _Xcentre2 << ", " << _Ycentre2 << ", " << _Zcentre2 << ", " << _Rayon2 << ")\n";
 811       break ;
 812     }
 813     case 5:
 814     {
 815       aScript << _Name << " = smeshhomard.CreateBoundaryTorus(\"" << _Name << "\", ";
 816       aScript << _Xcentre << ", " << _Ycentre << ", " << _Zcentre << ", " << _Xaxe << ", " << _Yaxe << ", " << _Zaxe << ", " << _Rayon1 << ", " << _Rayon2 << ")\n";
 817       break ;
 818     }
 819   }
 820
 821   return aScript.str();
 822 }
 823 //=============================================================================
 824 //=============================================================================
 825 // Caracteristiques
 826 //=============================================================================
 827 //=============================================================================
 828 void HOMARD_Boundary::SetType( int Type )
 829 {
 830   _Type = Type;
 831 }
 832 //=============================================================================
 833 int HOMARD_Boundary::GetType() const
 834 {
 835   return _Type;
 836 }
 837 //=============================================================================
 838 void HOMARD_Boundary::SetMeshName( const char* MeshName )
 839 {
 840   _MeshName = std::string( MeshName );
 841 }
 842 //=============================================================================
 843 std::string HOMARD_Boundary::GetMeshName() const
 844 {
 845   return _MeshName;
 846 }
 847 //=============================================================================
 848 void HOMARD_Boundary::SetDataFile( const char* DataFile )
 849 {
 850   _DataFile = std::string( DataFile );
 851 }
 852 //=============================================================================
 853 std::string HOMARD_Boundary::GetDataFile() const
 854 {
 855   return _DataFile;
 856 }
 857 //=======================================================================================
 858 void HOMARD_Boundary::SetCylinder( double X0, double X1, double X2,
 859                                    double X3, double X4, double X5, double X6 )
 860 {
 861   _Xcentre = X0; _Ycentre = X1; _Zcentre = X2;
 862   _Xaxe = X3; _Yaxe = X4; _Zaxe = X5;
 863   _rayon = X6;
 864 }
 865 //======================================================================
 866 void HOMARD_Boundary::SetSphere( double X0, double X1, double X2, double X3 )
 867 {
 868   _Xcentre = X0; _Ycentre = X1; _Zcentre = X2;
 869   _rayon = X3;
 870 }
 871 //======================================================================
 872 void HOMARD_Boundary::SetConeR( double Xcentre1, double Ycentre1, double Zcentre1, double Rayon1,
 873                                 double Xcentre2, double Ycentre2, double Zcentre2, double Rayon2)
 874 {
 875   _Xcentre1 = Xcentre1; _Ycentre1 = Ycentre1; _Zcentre1 = Zcentre1;
 876   _Rayon1 = Rayon1;
 877   _Xcentre2 = Xcentre2; _Ycentre2 = Ycentre2; _Zcentre2 = Zcentre2;
 878   _Rayon2 = Rayon2;
 879 }
 880 //======================================================================
 881 void HOMARD_Boundary::SetConeA( double Xaxe, double Yaxe, double Zaxe, double Angle,
 882                                 double Xcentre, double Ycentre, double Zcentre)
 883 {
 884   _Xaxe = Xaxe; _Yaxe = Yaxe; _Zaxe = Zaxe;
 885   _Angle = Angle;
 886   _Xcentre = Xcentre; _Ycentre = Ycentre; _Zcentre = Zcentre;
 887 }
 888 //=======================================================================================
 889 void HOMARD_Boundary::SetTorus( double X0, double X1, double X2,
 890                                 double X3, double X4, double X5, double X6, double X7 )
 891 {
 892   _Xcentre = X0; _Ycentre = X1; _Zcentre = X2;
 893   _Xaxe = X3; _Yaxe = X4; _Zaxe = X5;
 894   _Rayon1 = X6;
 895   _Rayon2 = X7;
 896 }
 897 //=======================================================================================
 898 std::vector<double> HOMARD_Boundary::GetCoords() const
 899 {
 900   std::vector<double> mesCoor;
 901 //
 902   switch (_Type)
 903   {
 904 //  Cylindre
 905     case 1:
 906     {
 907       mesCoor.push_back( _Xcentre );
 908       mesCoor.push_back( _Ycentre );
 909       mesCoor.push_back( _Zcentre );
 910       mesCoor.push_back( _Xaxe );
 911       mesCoor.push_back( _Yaxe );
 912       mesCoor.push_back( _Zaxe );
 913       mesCoor.push_back( _rayon );
 914       break ;
 915     }
 916 //  Sphere
 917     case 2:
 918     {
 919       mesCoor.push_back( _Xcentre );
 920       mesCoor.push_back( _Ycentre );
 921       mesCoor.push_back( _Zcentre );
 922       mesCoor.push_back( _rayon );
 923       break ;
 924     }
 925 //  Cone defini par un axe et un angle
 926     case 3:
 927     {
 928       mesCoor.push_back( _Xaxe );
 929       mesCoor.push_back( _Yaxe );
 930       mesCoor.push_back( _Zaxe );
 931       mesCoor.push_back( _Angle );
 932       mesCoor.push_back( _Xcentre );
 933       mesCoor.push_back( _Ycentre );
 934       mesCoor.push_back( _Zcentre );
 935       break ;
 936     }
 937 //  Cone defini par les 2 rayons
 938     case 4:
 939     {
 940       mesCoor.push_back( _Xcentre1 );
 941       mesCoor.push_back( _Ycentre1 );
 942       mesCoor.push_back( _Zcentre1 );
 943       mesCoor.push_back( _Rayon1 );
 944       mesCoor.push_back( _Xcentre2 );
 945       mesCoor.push_back( _Ycentre2 );
 946       mesCoor.push_back( _Zcentre2 );
 947       mesCoor.push_back( _Rayon2 );
 948       break ;
 949     }
 950 //  Tore
 951     case 5:
 952     {
 953       mesCoor.push_back( _Xcentre );
 954       mesCoor.push_back( _Ycentre );
 955       mesCoor.push_back( _Zcentre );
 956       mesCoor.push_back( _Xaxe );
 957       mesCoor.push_back( _Yaxe );
 958       mesCoor.push_back( _Zaxe );
 959       mesCoor.push_back( _Rayon1 );
 960       mesCoor.push_back( _Rayon2 );
 961       break ;
 962     }
 963     VERIFICATION( (_Type>=1) && (_Type<=5) ) ;
 964   }
 965   return mesCoor;
 966 }
 967 //======================================================================
 968 void HOMARD_Boundary::SetLimit( double X0, double X1, double X2 )
 969 {
 970   _Xincr = X0; _Yincr = X1; _Zincr = X2;
 971 }
 972 //=======================================================================================
 973 std::vector<double> HOMARD_Boundary::GetLimit() const
 974 {
 975   std::vector<double> mesLimit;
 976   mesLimit.push_back( _Xincr );
 977   mesLimit.push_back( _Yincr );
 978   mesLimit.push_back( _Zincr );
 979   return mesLimit;
 980 }
 981 //=============================================================================
 982 void HOMARD_Boundary::AddGroup( const char* Group)
 983 {
 984   _ListGroupSelected.push_back(Group);
 985 }
 986 //=============================================================================
 987 void HOMARD_Boundary::SetGroups( const std::list<std::string>& ListGroup )
 988 {
 989   _ListGroupSelected.clear();
 990   std::list<std::string>::const_iterator it = ListGroup.begin();
 991   while(it != ListGroup.end())
 992     _ListGroupSelected.push_back((*it++));
 993 }
 994 //=============================================================================
 995 const std::list<std::string>& HOMARD_Boundary::GetGroups() const
 996 {
 997   return _ListGroupSelected;
 998 }
 999 //=============================================================================
1000 //=============================================================================
1001 // Liens avec les autres structures
1002 //=============================================================================
1003 //=============================================================================
1004 void HOMARD_Boundary::SetCaseCreation( const char* NomCasCreation )
1005 {
1006   _NomCasCreation = std::string( NomCasCreation );
1007 }
1008 //=============================================================================
1009 std::string HOMARD_Boundary::GetCaseCreation() const
1010 {
1011   return _NomCasCreation;
1012 }
1013 //=============================================================================
1014
1015 //=============================================================================
1016 /*!
1017  *  default constructor:
1018  *  Par defaut, l'adaptation est conforme, sans suivi de frontiere
1019  */
1020 //=============================================================================
1021 HOMARD_Cas::HOMARD_Cas():
1022   _Name(""), _NomDir("/tmp"), _ConfType(0)
1023 {
1024   MESSAGE("HOMARD_Cas");
1025 }
1026 //=============================================================================
1027 HOMARD_Cas::~HOMARD_Cas()
1028 //=============================================================================
1029 {
1030   MESSAGE("~HOMARD_Cas");
1031 }
1032 //=============================================================================
1033 //=============================================================================
1034 // Generalites
1035 //=============================================================================
1036 //=============================================================================
1037 void HOMARD_Cas::SetName( const char* Name )
1038 {
1039   _Name = std::string( Name );
1040 }
1041 //=============================================================================
1042 std::string HOMARD_Cas::GetName() const
1043 {
1044   return _Name;
1045 }
1046 //=============================================================================
1047 std::string HOMARD_Cas::GetDumpPython() const
1048 {
1049   std::ostringstream aScript;
1050   //aScript << _Name << ".SetDirName(\"" << _NomDir << "\")\n";
1051   aScript << _Name << ".SetConfType(" << _ConfType << ")\n";
1052   // Suivi de frontieres
1053   std::list<std::string>::const_iterator it = _ListBoundaryGroup.begin();
1054   while (it != _ListBoundaryGroup.end()) {
1055     aScript << _Name << ".AddBoundaryGroup(\"" << *it << "\", \"";
1056     it++;
1057     aScript << *it << "\")\n";
1058     it++;
1059   }
1060
1061   return aScript.str();
1062 }
1063 //=============================================================================
1064 //=============================================================================
1065 // Caracteristiques
1066 //=============================================================================
1067 //=============================================================================
1068 int HOMARD_Cas::SetDirName( const char* NomDir )
1069 {
1070 //   MESSAGE("SetDirName,  NomDir : "<<NomDir);
1071 //   MESSAGE("SetDirName, _NomDir : "<<_NomDir);
1072   int erreur = 0 ;
1073   // On vérifie qu'aucun calcul n'a eu lieu pour ce cas
1074 //   MESSAGE("SetDirName, _ListIter.size() : "<<_ListIter.size());
1075   if ( _ListIter.size() > 1 ) { erreur = 1 ; }
1076   // Creation
1077   if ( CHDIR(NomDir) == 0 )
1078   { _NomDir = std::string( NomDir ); }
1079   else
1080   {
1081
1082 #ifndef WIN32
1083     if ( mkdir(NomDir, S_IRWXU|S_IRGRP|S_IXGRP) == 0 )
1084 #else
1085     if ( _mkdir(NomDir) == 0 )
1086 #endif
1087     {
1088       if ( CHDIR(NomDir) == 0 ) { _NomDir = std::string( NomDir ); }
1089       else                      { erreur = 2 ; }
1090     }
1091     else { erreur = 2 ; }
1092   };
1093   return erreur ;
1094 }
1095 //=============================================================================
1096 std::string HOMARD_Cas::GetDirName() const
1097 {
1098   return _NomDir;
1099 }
1100 //=============================================================================
1101 int HOMARD_Cas::GetNumberofIter()
1102 {
1103   return _ListIter.size();
1104 }
1105 //
1106 // Le type de conformite ou non conformite
1107 //
1108 //=============================================================================
1109 void HOMARD_Cas::SetConfType( int Conftype )
1110 {
1111 //   VERIFICATION( (Conftype>=-2) && (Conftype<=3) );
1112   _ConfType = Conftype;
1113 }
1114 //=============================================================================
1115 const int HOMARD_Cas::GetConfType() const
1116 {
1117   return _ConfType;
1118 }
1119 //
1120 // La boite englobante
1121 //
1122 //=============================================================================
1123 void HOMARD_Cas::SetBoundingBox( const std::vector<double>& extremas )
1124 {
1125   _Boite.clear();
1126   _Boite.resize( extremas.size() );
1127   for ( unsigned int i = 0; i < extremas.size(); i++ )
1128     _Boite[i] = extremas[i];
1129 }
1130 //=============================================================================
1131 const std::vector<double>& HOMARD_Cas::GetBoundingBox() const
1132 {
1133   return _Boite;
1134 }
1135 //
1136 // Les groupes
1137 //
1138 //=============================================================================
1139 void HOMARD_Cas::AddGroup( const char* Group )
1140 {
1141   _ListGroup.push_back(Group);
1142 }
1143 //=============================================================================
1144 void HOMARD_Cas::SetGroups( const std::list<std::string>& ListGroup )
1145 {
1146   _ListGroup.clear();
1147   std::list<std::string>::const_iterator it = ListGroup.begin();
1148   while(it != ListGroup.end())
1149   {
1150     _ListGroup.push_back((*it++));
1151   }
1152 }
1153 //=============================================================================
1154 const std::list<std::string>& HOMARD_Cas::GetGroups() const
1155 {
1156   return _ListGroup;
1157 }
1158 //=============================================================================
1159 void HOMARD_Cas::SupprGroups()
1160 {
1161   _ListGroup.clear();
1162 }
1163 //
1164 // Les frontieres
1165 //
1166 //=============================================================================
1167 void HOMARD_Cas::AddBoundary( const char* Boundary )
1168 {
1169 //   MESSAGE ( ". HOMARD_Cas::AddBoundary : Boundary = " << Boundary );
1170   const char* Group = "";
1171   AddBoundaryGroup( Boundary, Group );
1172 }
1173 //=============================================================================
1174 void HOMARD_Cas::AddBoundaryGroup( const char* Boundary, const char* Group )
1175 {
1176 //   MESSAGE ( ". HOMARD_Cas::AddBoundaryGroup : Boundary = " << Boundary );
1177 //   MESSAGE ( ". HOMARD_Cas::AddBoundaryGroup : Group = " << Group );
1178   _ListBoundaryGroup.push_back( Boundary );
1179   _ListBoundaryGroup.push_back( Group    );
1180 }
1181 //=============================================================================
1182 const std::list<std::string>& HOMARD_Cas::GetBoundaryGroup() const
1183 {
1184   return _ListBoundaryGroup;
1185 }
1186 //=============================================================================
1187 void HOMARD_Cas::SupprBoundaryGroup()
1188 {
1189   _ListBoundaryGroup.clear();
1190 }
1191 //=============================================================================
1192 //=============================================================================
1193 // Liens avec les autres structures
1194 //=============================================================================
1195 //=============================================================================
1196 std::string HOMARD_Cas::GetIter0Name() const
1197 {
1198 // Par construction de la liste, l'iteration a ete mise en tete.
1199   return (*(_ListIter.begin()));
1200 }
1201 //=============================================================================
1202 void HOMARD_Cas::AddIteration( const char* NomIteration )
1203 {
1204   _ListIter.push_back( std::string( NomIteration ) );
1205 }
1206 //=============================================================================
1207 const std::list<std::string>& HOMARD_Cas::GetIterations() const
1208 {
1209   return _ListIter;
1210 }
1211 //=============================================================================
1212 void HOMARD_Cas::SupprIterations()
1213 {
1214   _ListIter.clear();
1215 }
1216
1217 //=============================================================================
1218 //=============================================================================
1219 HomardDriver::HomardDriver(const std::string siter, const std::string siterp1):
1220   _HOMARD_Exec( "" ), _NomDir( "" ), _NomFichierConfBase( "HOMARD.Configuration" ),
1221   _NomFichierConf( "" ), _NomFichierDonn( "" ), _siter( "" ), _siterp1( "" ),
1222   _Texte( "" ), _bLu( false )
1223 {
1224   MESSAGE("siter = "<<siter<<", siterp1 = "<<siterp1);
1225   // Le repertoire ou se trouve l'executable HOMARD
1226   std::string dir ;
1227   // TODO?
1228   if ( getenv("HOMARD_ROOT_DIR") != NULL ) { dir = getenv("HOMARD_ROOT_DIR") ; }
1229   dir += "/bin/salome";
1230   MESSAGE("dir ="<<dir);
1231   // L'executable HOMARD
1232   std::string executable = "homard";
1233   MESSAGE("executable ="<<executable);
1234   // Memorisation du nom complet de l'executable HOMARD
1235   _HOMARD_Exec = dir + "/" + executable ;
1236   MESSAGE("==> _HOMARD_Exec ="<<_HOMARD_Exec) ;
1237   //
1238   _siter = siter ;
1239   _siterp1 = siterp1 ;
1240 }
1241 //=============================================================================
1242 //=============================================================================
1243 HomardDriver::~HomardDriver()
1244 {
1245 }
1246 //===============================================================================
1247 // A. Generalites
1248 //===============================================================================
1249 void HomardDriver::TexteInit( const std::string DirCompute, const std::string LogFile, const std::string Langue )
1250 {
1251   MESSAGE("TexteInit, DirCompute ="<<DirCompute<<", LogFile ="<<LogFile);
1252 //
1253   _Texte  = "ListeStd \"" + LogFile + "\"\n" ;
1254   _Texte += "RepeTrav \"" + DirCompute + "\"\n" ;
1255   _Texte += "RepeInfo \"" + DirCompute + "\"\n" ;
1256   _Texte += "Langue \"" + Langue + "\"\n" ;
1257 //
1258 }
1259 //===============================================================================
1260 void HomardDriver::TexteAdap()
1261 {
1262   MESSAGE("TexteAdap");
1263
1264   _Texte += "Action   homa\n";
1265   _Texte += "CCAssoci med\n";
1266   _Texte += "ModeHOMA 1\n";
1267   _Texte += "NumeIter " + _siter + "\n";
1268   _modeHOMARD = 1;
1269 }
1270 //===============================================================================
1271 void HomardDriver::TexteInfo( int TypeBila, int NumeIter )
1272 {
1273   MESSAGE("TexteInfo: TypeBila ="<<TypeBila<<", NumeIter ="<<NumeIter);
1274 //
1275   _Texte += "ModeHOMA 2\n" ;
1276   std::stringstream saux1 ;
1277   saux1 << TypeBila ;
1278   std::string saux2 = saux1.str() ;
1279   _Texte += "TypeBila " + saux2 + "\n" ;
1280   if ( NumeIter ==  0 )
1281   {
1282     _Texte += "NumeIter 0\n" ;
1283     _Texte += "Action   info_av\n" ;
1284     _Texte += "CCAssoci med\n" ;
1285   }
1286   else
1287   {
1288     _Texte += "NumeIter " + _siter + "\n" ;
1289     _Texte += "Action   info_ap\n" ;
1290     _Texte += "CCAssoci homard\n" ;
1291   }
1292   _modeHOMARD = 2 ;
1293 //
1294 }
1295 //===============================================================================
1296 void HomardDriver::TexteMajCoords( int NumeIter )
1297 {
1298   MESSAGE("TexteMajCoords: NumeIter ="<<NumeIter);
1299 //
1300   _Texte += "ModeHOMA 5\n" ;
1301   _Texte += "NumeIter " + _siterp1 + "\n" ;
1302   _Texte += "Action   homa\n" ;
1303   _Texte += "CCAssoci med\n" ;
1304   _Texte += "EcriFiHO N_SANS_FRONTIERE\n" ;
1305   _modeHOMARD = 5 ;
1306 //
1307 }
1308 //===============================================================================
1309 // B. Les maillages en entree et en sortie
1310 //===============================================================================
1311 void HomardDriver::TexteMaillage( const std::string NomMesh, const std::string MeshFile, int apres )
1312 {
1313   MESSAGE("TexteMaillage, NomMesh  = "<<NomMesh);
1314   MESSAGE("TexteMaillage, MeshFile = "<<MeshFile);
1315   MESSAGE("TexteMaillage, apres = "<<apres);
1316   std::string saux ;
1317   saux = "P1" ;
1318   if ( apres < 1 ) { saux = "__" ; }
1319
1320   _Texte += "# Maillages Med " + saux + "\n" ;
1321   _Texte += "CCNoMN" + saux + " \"" + NomMesh  + "\"\n" ;
1322   _Texte += "CCMaiN" + saux + " \"" + MeshFile + "\"\n" ;
1323 }
1324
1325 //===============================================================================
1326 void HomardDriver::TexteMaillageHOMARD( const std::string Dir, const std::string liter, int apres )
1327 {
1328   MESSAGE("TexteMaillageHOMARD, Dir ="<<Dir<<", liter ="<<liter<<", apres ="<<apres);
1329   std::string saux ;
1330   if ( apres < 1 ) { saux = "__" ; }
1331   else             { saux = "P1" ; }
1332
1333   _Texte += "# Maillage HOMARD " + liter + "\n" ;
1334   _Texte += "HOMaiN" + saux + " Mai" + liter   + " \"" + Dir + "/maill." + liter   + ".hom.med\"\n" ;
1335 }
1336
1337 //===============================================================================
1338 // C. Le pilotage de l'adaptation
1339 //===============================================================================
1340 void HomardDriver::TexteConfRaffDera( int ConfType, int TypeAdap, int TypeRaff, int TypeDera )
1341 {
1342   MESSAGE("TexteConfRaffDera, ConfType ="<<ConfType);
1343   MESSAGE("TexteConfRaffDera, TypeAdap ="<<TypeAdap<<", TypeRaff ="<<TypeRaff<<", TypeDera ="<<TypeDera);
1344 //
1345 // Type de conformite
1346 //
1347   std::string saux ;
1348   switch (ConfType)
1349   {
1350     case -2: //
1351     {
1352       saux = "NON_CONFORME_1_ARETE" ;
1353       break;
1354     }
1355     case -1: //
1356     {
1357       saux = "CONFORME_BOITES" ;
1358       break;
1359     }
1360     case 0: //
1361     {
1362       saux = "CONFORME" ;
1363       break;
1364     }
1365     case 1: //
1366     {
1367       saux = "NON_CONFORME" ;
1368       break;
1369     }
1370     case 2: //
1371     {
1372       saux = "NON_CONFORME_1_NOEUD" ;
1373       break;
1374     }
1375     case 3: //
1376     {
1377       saux = "NON_CONFORME_INDICATEUR" ;
1378       break;
1379     }
1380   }
1381   _Texte += "# Type de conformite\nTypeConf " + saux + "\n" ;
1382 //
1383 // Type de raffinement/deraffinement
1384 //
1385   if ( TypeAdap == -1 )
1386   {
1387     if ( TypeRaff == 1 )
1388     {
1389       saux = "TypeRaff uniforme\n" ;
1390     }
1391     else
1392     {
1393       saux = "TypeRaff non\n" ;
1394     }
1395     if ( TypeDera == 1 )
1396     {
1397       saux += "TypeDera uniforme" ;
1398     }
1399     else
1400     {
1401       saux += "TypeDera non" ;
1402     }
1403   }
1404   else
1405   {
1406     if ( TypeRaff == 1 )
1407     {
1408       saux = "TypeRaff libre\n" ;
1409     }
1410     else
1411     {
1412       saux = "TypeRaff non\n" ;
1413     }
1414     if ( TypeDera == 1 )
1415     {
1416       saux += "TypeDera libre" ;
1417     }
1418     else
1419     {
1420       saux += "TypeDera non" ;
1421     }
1422   }
1423   _Texte += "# Type de raffinement/deraffinement\n" + saux + "\n" ;
1424 //
1425 //   MESSAGE("A la fin de HomardDriver::TexteConfRaffDera, _Texte ="<<_Texte);
1426 }
1427 //===============================================================================
1428 void HomardDriver::TexteCompo( int NumeComp, const std::string NomCompo)
1429 {
1430   MESSAGE("TexteCompo, NumeComp = "<<NumeComp<<", NomCompo = "<<NomCompo);
1431   _Texte +="CCCoChaI \"" + NomCompo + "\"\n" ;
1432 }
1433 //===============================================================================
1434 void HomardDriver::TexteZone( int NumeZone, int ZoneType, int TypeUse, double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8 )
1435 {
1436   MESSAGE("TexteZone, NumeZone = "<<NumeZone<<", ZoneType = "<<ZoneType<<", TypeUse = "<<TypeUse);
1437   MESSAGE("TexteZone, coor = "<< x0<<","<<x1<< ","<< x2<< ","<< x3<<","<<x4<<","<<x5<<","<<x6<<","<<x7<<","<<x8);
1438 //
1439   std::string saux, saux2 ;
1440 //
1441 // Type de zones
1442 // On convertit le type de zone au sens du module HOMARD dans Salome, ZoneType, dans le
1443 // type au sens de l'executable HOMARD, ZoneTypeHOMARD
1444 // Attention a mettre le bon signe a ZoneTypeHOMARD :
1445 //    >0 signifie que l'on raffinera les mailles contenues dans la zone,
1446 //    <0 signifie que l'on deraffinera
1447 //
1448   int ZoneTypeHOMARD ;
1449   if ( ZoneType >= 11 && ZoneType <= 13 ) { ZoneTypeHOMARD = 1 ; }
1450   else if ( ZoneType >= 31 && ZoneType <= 33 ) { ZoneTypeHOMARD = 3 ; }
1451   else if ( ZoneType >= 61 && ZoneType <= 63 ) { ZoneTypeHOMARD = 6 ; }
1452   else { ZoneTypeHOMARD = ZoneType ; }
1453 //
1454   if ( TypeUse < 0 ) { ZoneTypeHOMARD = -ZoneTypeHOMARD ; }
1455 //
1456   std::stringstream saux1 ;
1457   saux1 << NumeZone ;
1458   saux = "#\n# Zone numero " + saux1.str() + "\n" ;
1459 //
1460   { std::stringstream saux1 ;
1461     saux1 << NumeZone << " " << ZoneTypeHOMARD ;
1462     saux += "ZoRaType " + saux1.str() + "\n" ;
1463   }
1464 //
1465 // Cas du rectangle
1466 //
1467   if ( ZoneType == 11 ) // Z est constant X Homard <=> X Salome
1468 //                                        Y Homard <=> Y Salome
1469   {
1470     saux += "#Rectangle\n" ;
1471     { std::stringstream saux1 ;
1472       saux1 << NumeZone << " " << x0 ;
1473       saux += "ZoRaXmin " + saux1.str() + "\n" ;
1474     }
1475     { std::stringstream saux1 ;
1476       saux1 << NumeZone << " " << x1 ;
1477       saux += "ZoRaXmax " + saux1.str() + "\n" ;
1478     }
1479     { std::stringstream saux1 ;
1480       saux1 << NumeZone << " " << x2 ;
1481       saux += "ZoRaYmin " + saux1.str() + "\n" ;
1482     }
1483     { std::stringstream saux1 ;
1484       saux1 << NumeZone << " " << x3 ;
1485       saux += "ZoRaYmax " + saux1.str() + "\n" ;
1486     }
1487   }
1488 //
1489   else if ( ZoneType == 12 ) // X est constant X Homard <=> Y Salome
1490 //                                             Y Homard <=> Z Salome
1491   {
1492     saux += "#Rectangle\n" ;
1493     { std::stringstream saux1 ;
1494       saux1 << NumeZone << " " << x2 ;
1495       saux += "ZoRaXmin " + saux1.str() + "\n" ;
1496     }
1497     { std::stringstream saux1 ;
1498       saux1 << NumeZone << " " << x3 ;
1499       saux += "ZoRaXmax " + saux1.str() + "\n" ;
1500     }
1501     { std::stringstream saux1 ;
1502       saux1 << NumeZone << " " << x4 ;
1503       saux += "ZoRaYmin " + saux1.str() + "\n" ;
1504     }
1505     { std::stringstream saux1 ;
1506       saux1 << NumeZone << " " << x5 ;
1507       saux += "ZoRaYmax " + saux1.str() + "\n" ;
1508     }
1509   }
1510 //
1511   else if ( ZoneType == 13 ) // Y est constant X Homard <=> X Salome
1512 //                                             Y Homard <=> Z Salome
1513   {
1514     saux += "#Rectangle\n" ;
1515     { std::stringstream saux1 ;
1516       saux1 << NumeZone << " " << x0 ;
1517       saux += "ZoRaXmin " + saux1.str() + "\n" ;
1518     }
1519     { std::stringstream saux1 ;
1520       saux1 << NumeZone << " " << x1 ;
1521       saux += "ZoRaXmax " + saux1.str() + "\n" ;
1522     }
1523     { std::stringstream saux1 ;
1524       saux1 << NumeZone << " " << x4 ;
1525       saux += "ZoRaYmin " + saux1.str() + "\n" ;
1526     }
1527     { std::stringstream saux1 ;
1528       saux1 << NumeZone << " " << x5 ;
1529       saux += "ZoRaYmax " + saux1.str() + "\n" ;
1530     }
1531   }
1532 //
1533 // Cas du parallelepipede
1534 //
1535   else if ( ZoneType == 2 )
1536   {
1537     saux += "# Boite\n" ;
1538     { std::stringstream saux1 ;
1539       saux1 << NumeZone << " " << x0 ;
1540       saux += "ZoRaXmin " + saux1.str() + "\n" ;
1541     }
1542     { std::stringstream saux1 ;
1543       saux1 << NumeZone << " " << x1 ;
1544       saux += "ZoRaXmax " + saux1.str() + "\n" ;
1545     }
1546     { std::stringstream saux1 ;
1547       saux1 << NumeZone << " " << x2 ;
1548       saux += "ZoRaYmin " + saux1.str() + "\n" ;
1549     }
1550     { std::stringstream saux1 ;
1551       saux1 << NumeZone << " " << x3 ;
1552       saux += "ZoRaYmax " + saux1.str() + "\n" ;
1553     }
1554     { std::stringstream saux1 ;
1555       saux1 << NumeZone << " " << x4 ;
1556       saux += "ZoRaZmin " + saux1.str() + "\n" ;
1557     }
1558     { std::stringstream saux1 ;
1559       saux1 << NumeZone << " " << x5 ;
1560       saux += "ZoRaZmax " + saux1.str() + "\n" ;
1561     }
1562   }
1563 //
1564 // Cas du disque
1565 //
1566   else if ( ZoneType == 31 || ZoneType == 61 )
1567   {
1568     saux += "# Sphere\n" ;
1569     { std::stringstream saux1 ;
1570       saux1 << NumeZone << " " << x0 ;
1571       saux += "ZoRaXCen " + saux1.str() + "\n" ;
1572     }
1573     { std::stringstream saux1 ;
1574       saux1 << NumeZone << " " << x1 ;
1575       saux += "ZoRaYCen " + saux1.str() + "\n" ;
1576     }
1577     { std::stringstream saux1 ;
1578       saux1 << NumeZone << " " << x6 ;
1579       saux2 = saux1.str() ;
1580       if ( ZoneType == 61 ) { saux += "ZoRaRayE " + saux2 + "\n" ; }
1581       else                  { saux += "ZoRaRayo " + saux2 + "\n" ; }
1582     }
1583     if ( ZoneType == 61 )
1584     { std::stringstream saux1 ;
1585       saux1 << NumeZone << " " << x8 ;
1586       saux += "ZoRaRayI " + saux1.str() + "\n" ;
1587     }
1588   }
1589   else if ( ZoneType == 32 || ZoneType == 62 )
1590   {
1591     saux += "# Sphere\n" ;
1592     { std::stringstream saux1 ;
1593       saux1 << NumeZone << " " << x1 ;
1594       saux += "ZoRaXCen " + saux1.str() + "\n" ;
1595     }
1596     { std::stringstream saux1 ;
1597       saux1 << NumeZone << " " << x2 ;
1598       saux += "ZoRaYCen " + saux1.str() + "\n" ;
1599     }
1600     { std::stringstream saux1 ;
1601       saux1 << NumeZone << " " << x6 ;
1602       saux2 = saux1.str() ;
1603       if ( ZoneType == 62 ) { saux += "ZoRaRayE " + saux2 + "\n" ; }
1604       else                  { saux += "ZoRaRayo " + saux2 + "\n" ; }
1605     }
1606     if ( ZoneType == 62 )
1607     { std::stringstream saux1 ;
1608       saux1 << NumeZone << " " << x8 ;
1609       saux += "ZoRaRayI " + saux1.str() + "\n" ;
1610     }
1611   }
1612   else if ( ZoneType == 33 || ZoneType == 63 )
1613   {
1614     saux += "# Sphere\n" ;
1615     { std::stringstream saux1 ;
1616       saux1 << NumeZone << " " << x0 ;
1617       saux += "ZoRaXCen " + saux1.str() + "\n" ;
1618     }
1619     { std::stringstream saux1 ;
1620       saux1 << NumeZone << " " << x2 ;
1621       saux += "ZoRaYCen " + saux1.str() + "\n" ;
1622     }
1623     { std::stringstream saux1 ;
1624       saux1 << NumeZone << " " << x6 ;
1625       saux2 = saux1.str() ;
1626       if ( ZoneType == 63 ) { saux += "ZoRaRayE " + saux2 + "\n" ; }
1627       else                  { saux += "ZoRaRayo " + saux2 + "\n" ; }
1628     }
1629     if ( ZoneType == 63 )
1630     { std::stringstream saux1 ;
1631       saux1 << NumeZone << " " << x8 ;
1632       saux += "ZoRaRayI " + saux1.str() + "\n" ;
1633     }
1634   }
1635 //
1636 // Cas de la sphere
1637 //
1638   else if ( ZoneType == 4 )
1639   {
1640     saux += "# Sphere\n" ;
1641     { std::stringstream saux1 ;
1642       saux1 << NumeZone << " " << x0 ;
1643       saux += "ZoRaXCen " + saux1.str() + "\n" ;
1644     }
1645     { std::stringstream saux1 ;
1646       saux1 << NumeZone << " " << x1 ;
1647       saux += "ZoRaYCen " + saux1.str() + "\n" ;
1648     }
1649     { std::stringstream saux1 ;
1650       saux1 << NumeZone << " " << x2 ;
1651       saux += "ZoRaZCen " + saux1.str() + "\n" ;
1652     }
1653     { std::stringstream saux1 ;
1654       saux1 << NumeZone << " " << x3 ;
1655       saux += "ZoRaRayo " + saux1.str() + "\n" ;
1656     }
1657   }
1658 //
1659 // Cas du cylindre ou du tuyau
1660 //
1661   else if ( ZoneType == 5 || ZoneType == 7 )
1662   {
1663     if ( ZoneType == 5 ) { saux += "# Cylindre\n" ; }
1664     else                 { saux += "# Tuyau\n" ; }
1665     { std::stringstream saux1 ;
1666       saux1 << NumeZone << " " << x0 ;
1667       saux += "ZoRaXBas " + saux1.str() + "\n" ;
1668     }
1669     { std::stringstream saux1 ;
1670       saux1 << NumeZone << " " << x1 ;
1671       saux += "ZoRaYBas " + saux1.str() + "\n" ;
1672     }
1673     { std::stringstream saux1 ;
1674       saux1 << NumeZone << " " << x2 ;
1675       saux += "ZoRaZBas " + saux1.str() + "\n" ;
1676     }
1677     { std::stringstream saux1 ;
1678       saux1 << NumeZone << " " << x3 ;
1679       saux += "ZoRaXAxe " + saux1.str() + "\n" ;
1680     }
1681     { std::stringstream saux1 ;
1682       saux1 << NumeZone << " " << x4 ;
1683       saux += "ZoRaYAxe " + saux1.str() + "\n" ;
1684     }
1685     { std::stringstream saux1 ;
1686       saux1 << NumeZone << " " << x5 ;
1687       saux += "ZoRaZAxe " + saux1.str() + "\n" ;
1688     }
1689     { std::stringstream saux1 ;
1690       saux1 << NumeZone << " " << x6 ;
1691       saux2 = saux1.str() ;
1692      if ( ZoneType == 5 ) { saux += "ZoRaRayo " + saux2 + "\n" ; }
1693      else                 { saux += "ZoRaRayE " + saux2 + "\n" ; }
1694     }
1695     { std::stringstream saux1 ;
1696       saux1 << NumeZone << " " << x7 ;
1697       saux += "ZoRaHaut " + saux1.str() + "\n" ;
1698     }
1699     if ( ZoneType == 7 )
1700     { std::stringstream saux1 ;
1701       saux1 << NumeZone << " " << x8 ;
1702       saux += "ZoRaRayI " + saux1.str() + "\n" ;
1703     }
1704   }
1705 //
1706   _Texte += saux + "#\n" ;
1707 //
1708 //   MESSAGE("A la fin de HomardDriver::TexteZone, _Texte ="<<_Texte);
1709 }
1710 //===============================================================================
1711 void HomardDriver::TexteField( const std::string FieldName, const std::string FieldFile, int TimeStep, int Rank,
1712                int TypeThR, double ThreshR, int TypeThC, double ThreshC,
1713                int UsField, int UsCmpI )
1714 {
1715   MESSAGE("TexteField, FieldName = "<<FieldName<<", FieldFile = "<<FieldFile);
1716   MESSAGE("TexteField, TimeStep = "<<TimeStep<<", Rank = "<<Rank);
1717
1718   std::string saux, saux2 ;
1719 //
1720 //
1721   _Texte += "# Champ d'indicateurs\n" ;
1722   _Texte += "CCIndica \"" + FieldFile  + "\"\n" ;
1723   _Texte += "CCNoChaI \"" + FieldName  + "\"\n" ;
1724
1725 // Cas ou on prend le dernier pas de temps
1726   if ( TimeStep == -2 )
1727   { _Texte += "CCNumPTI Last\n" ; }
1728 // Cas avec pas de temps
1729   else if ( TimeStep >= 0 )
1730   {
1731     {
1732       std::stringstream saux1 ;
1733       saux1 << TimeStep ;
1734       saux2 = saux1.str() ;
1735       _Texte += "CCNumPTI " + saux2  + "\n" ;
1736     }
1737     if ( Rank >= 0 )
1738     {
1739       std::stringstream saux1 ;
1740       saux1 << Rank ;
1741       saux2 = saux1.str() ;
1742       _Texte += "CCNumOrI " + saux2  + "\n" ;
1743     }
1744   }
1745 //
1746   saux = " " ;
1747   if ( TypeThR == 1 )
1748   { saux = "Hau" ; }
1749   if ( TypeThR == 2 )
1750   { saux = "HRe" ; }
1751   if ( TypeThR == 3 )
1752   { saux = "HPE" ; }
1753   if ( TypeThR == 4 )
1754   { saux = "HMS" ; }
1755   if ( saux != " " )
1756   {
1757     std::stringstream saux1 ;
1758     saux1 << ThreshR ;
1759     _Texte += "Seuil" + saux + " " + saux1.str()  + "\n" ;
1760   }
1761 //
1762   saux = " " ;
1763   if ( TypeThC == 1 )
1764   { saux = "Bas" ; }
1765   if ( TypeThC == 2 )
1766   { saux = "BRe" ; }
1767   if ( TypeThC == 3 )
1768   { saux = "BPE" ; }
1769   if ( TypeThC == 4 )
1770   { saux = "BMS" ; }
1771   if ( saux != " " )
1772   {
1773     std::stringstream saux1 ;
1774     saux1 << ThreshC ;
1775     _Texte += "Seuil" + saux + " " + saux1.str()  + "\n" ;
1776   }
1777 //
1778   saux = " " ;
1779   if ( UsField == 0 )
1780   { saux = "MAILLE" ; }
1781   if ( UsField == 1 )
1782   { saux = "SAUT" ; }
1783   if ( saux != " " )
1784   {
1785     _Texte += "CCModeFI " + saux  + "\n" ;
1786   }
1787 //
1788   saux = " " ;
1789   if ( UsCmpI == 0 )
1790   { saux = "L2" ; }
1791   if ( UsCmpI == 1 )
1792   { saux = "INFINI" ; }
1793   if ( UsCmpI == 2 )
1794   { saux = "RELATIF" ; }
1795   if ( saux != " " )
1796   {
1797     _Texte += "CCUsCmpI " + saux  + "\n" ;
1798   }
1799 }
1800 //===============================================================================
1801 void HomardDriver::TexteGroup( const std::string GroupName )
1802 {
1803   MESSAGE("TexteGroup, GroupName = "<<GroupName);
1804 //
1805   _Texte += "CCGroAda \"" + GroupName  + "\"\n" ;
1806 //
1807 }
1808 //===============================================================================
1809 // D. Les frontieres
1810 //===============================================================================
1811 void HomardDriver::TexteBoundaryOption( int BoundaryOption )
1812 {
1813   MESSAGE("TexteBoundaryOption, BoundaryOption = "<<BoundaryOption);
1814 //
1815 // Type de suivi de frontiere
1816 //
1817   std::stringstream saux1 ;
1818   saux1 << BoundaryOption ;
1819   std::string saux = saux1.str() ;
1820   _Texte += "SuivFron " + saux + "\n" ;
1821 //
1822 }//===============================================================================
1823 void HomardDriver::TexteBoundaryCAOGr(  const std::string GroupName )
1824 {
1825   MESSAGE("TexteBoundaryCAOGr, GroupName  = "<<GroupName);
1826 //
1827   _Texte += "GrFroCAO \"" + GroupName + "\"\n" ;
1828 //
1829 }
1830
1831 //===============================================================================
1832 void HomardDriver::TexteBoundaryDi(  const std::string MeshName, const std::string MeshFile )
1833 {
1834   MESSAGE("TexteBoundaryDi, MeshName  = "<<MeshName);
1835   MESSAGE("TexteBoundaryDi, MeshFile  = "<<MeshFile);
1836 //
1837   _Texte += "#\n# Frontiere discrete\n" ;
1838   _Texte += "CCNoMFro \"" + MeshName + "\"\n" ;
1839   _Texte += "CCFronti \"" + MeshFile + "\"\n" ;
1840 //
1841 }
1842 //===============================================================================
1843 void HomardDriver::TexteBoundaryDiGr(  const std::string GroupName )
1844 {
1845   MESSAGE("TexteBoundaryDiGr, GroupName  = "<<GroupName);
1846 //
1847   _Texte += "CCGroFro \"" + GroupName + "\"\n" ;
1848 //
1849 }
1850 //===============================================================================
1851 void HomardDriver::TexteBoundaryAn( const std::string NameBoundary, int NumeBoundary, int BoundaryType, double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7 )
1852 {
1853   MESSAGE("TexteBoundaryAn, NameBoundary = "<<NameBoundary);
1854 //   MESSAGE("TexteBoundaryAn, NumeBoundary = "<<NumeBoundary);
1855   MESSAGE("TexteBoundaryAn, BoundaryType = "<<BoundaryType);
1856 //   MESSAGE("TexteBoundaryAn, coor         = "<< x0<<","<<x1<< ","<< x2<< ","<< x3<<","<<x4<<","<<x5<<","<<x6","<<x7);
1857 //
1858   std::string saux, saux2 ;
1859 //
1860 // Commentaires
1861 //
1862   std::stringstream saux1 ;
1863   saux1 << NumeBoundary ;
1864   saux2 = saux1.str() ;
1865   saux = "#\n# Frontiere numero " + saux2 + "\n" ;
1866   if ( BoundaryType == 1 )
1867   { saux += "# Cylindre\n" ; }
1868   if ( BoundaryType == 2 )
1869   { saux += "# Sphere\n" ; }
1870   if ( BoundaryType == 3 || BoundaryType == 4 )
1871   { saux += "# Cone\n" ; }
1872   if ( BoundaryType == 5 )
1873   { saux += "# Tore\n" ; }
1874 //
1875 // Le nom de la frontiere
1876 //
1877   { std::stringstream saux1 ;
1878     saux1 << NumeBoundary ;
1879     saux += "FANom " + saux1.str() + " \"" + NameBoundary + "\"\n" ;
1880   }
1881 //
1882 // Type de frontiere
1883 //
1884   { std::stringstream saux1 ;
1885     saux1 << NumeBoundary << " " << BoundaryType ;
1886     saux += "FAType " + saux1.str() + "\n" ;
1887   }
1888 //
1889 // Cas du cylindre
1890 //
1891   if ( BoundaryType == 1 )
1892   {
1893     { std::stringstream saux1 ;
1894       saux1 << NumeBoundary << " " << x0 ;
1895       saux2 = saux1.str() ;
1896       saux += "FAXCen " + saux1.str() + "\n" ;
1897     }
1898     { std::stringstream saux1 ;
1899       saux1 << NumeBoundary << " " << x1 ;
1900       saux += "FAYCen " + saux1.str() + "\n" ;
1901     }
1902     { std::stringstream saux1 ;
1903       saux1 << NumeBoundary << " " << x2 ;
1904       saux += "FAZCen " + saux1.str() + "\n" ;
1905     }
1906     { std::stringstream saux1 ;
1907       saux1 << NumeBoundary << " " << x3 ;
1908       saux += "FAXAxe " + saux1.str() + "\n" ;
1909     }
1910     { std::stringstream saux1 ;
1911       saux1 << NumeBoundary << " " << x4 ;
1912       saux += "FAYAxe " + saux1.str() + "\n" ;
1913     }
1914     { std::stringstream saux1 ;
1915       saux1 << NumeBoundary << " " << x5 ;
1916       saux += "FAZAxe " + saux1.str() + "\n" ;
1917     }
1918     { std::stringstream saux1 ;
1919       saux1 << NumeBoundary << " " << x6 ;
1920       saux += "FARayon " + saux1.str() + "\n" ;
1921     }
1922  }
1923 //
1924 // Cas de la sphere
1925 //
1926   else if ( BoundaryType == 2 )
1927   {
1928     { std::stringstream saux1 ;
1929       saux1 << NumeBoundary << " " << x0 ;
1930       saux += "FAXCen " + saux1.str() + "\n" ;
1931     }
1932     { std::stringstream saux1 ;
1933       saux1 << NumeBoundary << " " << x1 ;
1934       saux += "FAYCen " + saux1.str() + "\n" ;
1935     }
1936     { std::stringstream saux1 ;
1937       saux1 << NumeBoundary << " " << x2 ;
1938       saux += "FAZCen " + saux1.str() + "\n" ;
1939     }
1940     { std::stringstream saux1 ;
1941       saux1 << NumeBoundary << " " << x3 ;
1942       saux += "FARayon " + saux1.str() + "\n" ;
1943     }
1944   }
1945 //
1946 // Cas du cone defini par un axe et un angle
1947 //
1948   if ( BoundaryType == 3 )
1949   {
1950     { std::stringstream saux1 ;
1951       saux1 << NumeBoundary << " " << x0 ;
1952       saux += "FAXAxe " + saux1.str() + "\n" ;
1953     }
1954     { std::stringstream saux1 ;
1955       saux1 << NumeBoundary << " " << x1 ;
1956       saux += "FAYAxe " + saux1.str() + "\n" ;
1957     }
1958     { std::stringstream saux1 ;
1959       saux1 << NumeBoundary << " " << x2 ;
1960       saux += "FAZAxe " + saux1.str() + "\n" ;
1961     }
1962     { std::stringstream saux1 ;
1963       saux1 << NumeBoundary << " " << x3 ;
1964       saux += "FAAngle " + saux1.str() + "\n" ;
1965     }
1966     { std::stringstream saux1 ;
1967       saux1 << NumeBoundary << " " << x4 ;
1968       saux += "FAXCen " + saux1.str() + "\n" ;
1969     }
1970     { std::stringstream saux1 ;
1971       saux1 << NumeBoundary << " " << x5 ;
1972       saux += "FAYCen " + saux1.str() + "\n" ;
1973     }
1974     { std::stringstream saux1 ;
1975       saux1 << NumeBoundary << " " << x6 ;
1976       saux += "FAZCen " + saux1.str() + "\n" ;
1977     }
1978  }
1979 //
1980 // Cas du cone defini par les 2 rayons
1981 //
1982   if ( BoundaryType == 4 )
1983   {
1984     { std::stringstream saux1 ;
1985       saux1 << NumeBoundary << " " << x0 ;
1986       saux += "FAXCen " + saux1.str() + "\n" ;
1987     }
1988     { std::stringstream saux1 ;
1989       saux1 << NumeBoundary << " " << x1 ;
1990       saux += "FAYCen " + saux1.str() + "\n" ;
1991     }
1992     { std::stringstream saux1 ;
1993       saux1 << NumeBoundary << " " << x2 ;
1994       saux += "FAZCen " + saux1.str() + "\n" ;
1995     }
1996     { std::stringstream saux1 ;
1997       saux1 << NumeBoundary << " " << x3 ;
1998       saux += "FARayon " + saux1.str() + "\n" ;
1999     }
2000     { std::stringstream saux1 ;
2001       saux1 << NumeBoundary << " " << x4 ;
2002       saux += "FAXCen2 " + saux1.str() + "\n" ;
2003     }
2004     { std::stringstream saux1 ;
2005       saux1 << NumeBoundary << " " << x5 ;
2006       saux += "FAYCen2 " + saux1.str() + "\n" ;
2007     }
2008     { std::stringstream saux1 ;
2009       saux1 << NumeBoundary << " " << x6 ;
2010       saux += "FAZCen2 " + saux1.str() + "\n" ;
2011     }
2012     { std::stringstream saux1 ;
2013       saux1 << NumeBoundary << " " << x7 ;
2014       saux += "FARayon2 " + saux1.str() + "\n" ;
2015     }
2016  }
2017 //
2018 // Cas du tore
2019 //
2020   if ( BoundaryType == 5 )
2021   {
2022     { std::stringstream saux1 ;
2023       saux1 << NumeBoundary << " " << x0 ;
2024       saux2 = saux1.str() ;
2025       saux += "FAXCen " + saux1.str() + "\n" ;
2026     }
2027     { std::stringstream saux1 ;
2028       saux1 << NumeBoundary << " " << x1 ;
2029       saux += "FAYCen " + saux1.str() + "\n" ;
2030     }
2031     { std::stringstream saux1 ;
2032       saux1 << NumeBoundary << " " << x2 ;
2033       saux += "FAZCen " + saux1.str() + "\n" ;
2034     }
2035     { std::stringstream saux1 ;
2036       saux1 << NumeBoundary << " " << x3 ;
2037       saux += "FAXAxe " + saux1.str() + "\n" ;
2038     }
2039     { std::stringstream saux1 ;
2040       saux1 << NumeBoundary << " " << x4 ;
2041       saux += "FAYAxe " + saux1.str() + "\n" ;
2042     }
2043     { std::stringstream saux1 ;
2044       saux1 << NumeBoundary << " " << x5 ;
2045       saux += "FAZAxe " + saux1.str() + "\n" ;
2046     }
2047     { std::stringstream saux1 ;
2048       saux1 << NumeBoundary << " " << x6 ;
2049       saux += "FARayon  " + saux1.str() + "\n" ;
2050     }
2051     { std::stringstream saux1 ;
2052       saux1 << NumeBoundary << " " << x7 ;
2053       saux += "FARayon2 " + saux1.str() + "\n" ;
2054     }
2055  }
2056 //
2057   _Texte += saux + "#\n" ;
2058 //
2059 }
2060 //===============================================================================
2061 void HomardDriver::TexteBoundaryAnGr( const std::string NameBoundary, int NumeBoundary, const std::string GroupName )
2062 {
2063   MESSAGE("TexteBoundaryAnGr, NameBoundary  = "<<NameBoundary);
2064 //   MESSAGE("TexteBoundaryAnGr, NumeBoundary  = "<<NumeBoundary);
2065 //   MESSAGE("TexteBoundaryAnGr, GroupName  = "<<GroupName);
2066 //
2067 // Commentaires
2068 //
2069   std::string saux, saux2 ;
2070   std::stringstream saux1 ;
2071   saux1 << NumeBoundary ;
2072   saux2 = saux1.str() ;
2073   saux = "#\n# Lien Frontiere/Groupe numero " + saux2 + "\n" ;
2074 //
2075   saux += "FGNomFro " + saux2 + " \"" + NameBoundary + "\"\n" ;
2076   saux += "FGNomGro " + saux2 + " \"" + GroupName + "\"\n" ;
2077 //
2078   _Texte += saux + "#\n" ;
2079 //
2080 }
2081 //===============================================================================
2082 // E. Les interpolations
2083 //===============================================================================
2084 // Les fichiers d'entree et de sortie des champs a interpoler
2085 void HomardDriver::TexteFieldInterp( const std::string FieldFile, const std::string MeshFile )
2086 {
2087   MESSAGE("TexteFieldInterp, FieldFile = "<<FieldFile<<", MeshFile = "<<MeshFile);
2088 //
2089   _Texte += "#\n# Interpolations des champs\n" ;
2090 //
2091 // Fichier en entree
2092   _Texte += "CCSolN__ \"" + FieldFile + "\"\n" ;
2093 // Fichier en sortie
2094   _Texte += "CCSolNP1 \"" + MeshFile  + "\"\n" ;
2095 //
2096 //  std::cerr << "A la fin de TexteFieldInterp _Texte ="<<_Texte << std::endl;
2097 }
2098 //===============================================================================
2099 // Tous les champs sont a interpoler
2100 void HomardDriver::TexteFieldInterpAll( )
2101 {
2102   MESSAGE("TexteFieldInterpAll");
2103 //
2104   _Texte += "CCChaTou oui\n" ;
2105 }
2106 //===============================================================================
2107 // Ecrit les caracteristiques de chaque interpolation sous la forme :
2108 //   CCChaNom 1 "DEPL"     ! Nom du 1er champ a interpoler
2109 //   CCChaTIn 1 0          ! Mode d'interpolation : automatique
2110 //   CCChaNom 2 "VOLUME"   ! Nom du 2nd champ a interpoler
2111 //   CCChaTIn 2 1          ! Mode d'interpolation : une variable extensive
2112 //   CCChaPdT 2 14         ! Pas de temps 14
2113 //   CCChaNuO 2 14         ! Numero d'ordre 14
2114 //   etc.
2115 //
2116 // NumeChamp : numero d'ordre du champ a interpoler
2117 // FieldName : nom du champ
2118 // TypeInterp : type d'interpolation
2119 // TimeStep : pas de temps retenu (>0 si pas de precision)
2120 // Rank : numero d'ordre retenu
2121 //
2122 void HomardDriver::TexteFieldInterpNameType( int NumeChamp, const std::string FieldName, const std::string TypeInterp, int TimeStep, int Rank)
2123 {
2124   MESSAGE("TexteFieldInterpNameType, NumeChamp = "<<NumeChamp<<", FieldName = "<<FieldName<<", TypeInterp = "<<TypeInterp);
2125   MESSAGE("TexteFieldInterpNameType, TimeStep = "<<TimeStep<<", Rank = "<<Rank);
2126 // Numero d'ordre du champ a interpoler
2127   std::stringstream saux1 ;
2128   saux1 << NumeChamp ;
2129   std::string saux = saux1.str() ;
2130 // Nom du champ
2131   _Texte +="CCChaNom " + saux + " \"" + FieldName + "\"\n" ;
2132 // Type d'interpolation pour le champ
2133   _Texte +="CCChaTIn " + saux + " " + TypeInterp + "\n" ;
2134 //
2135   if ( TimeStep >= 0 )
2136   {
2137     {
2138       std::stringstream saux1 ;
2139       saux1 << TimeStep ;
2140       _Texte += "CCChaPdT " + saux + " " + saux1.str()  + "\n" ;
2141     }
2142     {
2143       std::stringstream saux1 ;
2144       saux1 << Rank ;
2145       _Texte += "CCChaNuO " + saux + " " + saux1.str()  + "\n" ;
2146     }
2147   }
2148 }
2149 //===============================================================================
2150 // F. Les options avancees
2151 //===============================================================================
2152 void HomardDriver::TexteAdvanced( int NivMax, double DiamMin, int AdapInit, int ExtraOutput )
2153 {
2154   MESSAGE("TexteAdvanced, NivMax ="<<NivMax<<", DiamMin ="<<DiamMin<<", AdapInit ="<<AdapInit<<", ExtraOutput ="<<ExtraOutput);
2155
2156   if ( NivMax > 0 )
2157   {
2158     _Texte += "# Niveaux extremes\n" ;
2159     { std::stringstream saux1 ;
2160       saux1 << NivMax ;
2161       _Texte += "NiveauMa " + saux1.str() + "\n" ;
2162     }
2163   }
2164   if ( DiamMin > 0 )
2165   {
2166     _Texte += "# Diametre minimal\n" ;
2167     { std::stringstream saux1 ;
2168       saux1 << DiamMin ;
2169       _Texte += "DiametMi " + saux1.str()  + "\n" ;
2170     }
2171   }
2172   if ( AdapInit != 0 )
2173   {
2174     if ( AdapInit > 0 )
2175     { _Texte += "# Raffinement" ; }
2176     else
2177     { _Texte += "# Deraffinement" ; }
2178     _Texte += " des regions sans indicateur\n" ;
2179     { std::stringstream saux1 ;
2180       saux1 << AdapInit ;
2181       _Texte += "AdapInit " + saux1.str() + "\n" ;
2182     }
2183   }
2184   if ( ExtraOutput % 2 == 0 )
2185   {
2186     _Texte += "# Sortie des niveaux de raffinement\n" ;
2187     _Texte += "NCNiveau NIVEAU\n" ;
2188   }
2189   if ( ExtraOutput % 3 == 0 )
2190   {
2191     _Texte += "# Sortie des qualités des mailles\n" ;
2192     _Texte += "NCQualit QUAL\n" ;
2193   }
2194   if ( ExtraOutput % 5 == 0 )
2195   {
2196     _Texte += "# Sortie des diamètres des mailles\n" ;
2197     _Texte += "NCDiamet DIAM\n" ;
2198   }
2199   if ( ExtraOutput % 7 == 0 )
2200   {
2201     _Texte += "# Sortie des parents des mailles\n" ;
2202     _Texte += "NCParent PARENT\n" ;
2203   }
2204   if ( ExtraOutput % 11 == 0 )
2205   {
2206     _Texte += "# Volumes voisins par recollement\n" ;
2207     _Texte += "NCVoisRc Voisin-Recollement\n" ;
2208   }
2209 }
2210 //===============================================================================
2211 // G. Les messages
2212 //===============================================================================
2213 void HomardDriver::TexteInfoCompute( int MessInfo )
2214 {
2215   MESSAGE("TexteAdvanced, MessInfo ="<<MessInfo);
2216
2217   if ( MessInfo != 0 )
2218   {
2219      _Texte += "# Messages d'informations\n" ;
2220     { std::stringstream saux1 ;
2221       saux1 << MessInfo ;
2222       _Texte += "MessInfo " + saux1.str()  + "\n" ;
2223     }
2224    }
2225 }
2226 //===============================================================================
2227 void HomardDriver::CreeFichier( )
2228 {
2229 //
2230   if ( _modeHOMARD == 1 )
2231   { _NomFichierConf = _NomFichierConfBase + "." + _siter + ".vers." + _siterp1 ; }
2232   else if ( _modeHOMARD == 2 )
2233   { _NomFichierConf = _NomFichierConfBase + "." + _siter + ".info" ; }
2234   else if ( _modeHOMARD == 5 )
2235   { _NomFichierConf = _NomFichierConfBase + ".majc" ; }
2236 //
2237   std::ofstream Fic(_NomFichierConf.c_str(), std::ios::out ) ;
2238   if (Fic.is_open() == true) { Fic << _Texte << std::endl ; }
2239   Fic.close() ;
2240 //
2241 }
2242 //===============================================================================
2243 // Creation du fichier de donnees pour l'information
2244 //===============================================================================
2245 void HomardDriver::CreeFichierDonn( )
2246 {
2247 //
2248   MESSAGE("CreeFichierDonn");
2249   _NomFichierDonn = "info.donn" ;
2250 //
2251   std::string data ;
2252   data  = "0\n" ;
2253   data += "0\n" ;
2254   data += "q\n" ;
2255   std::ofstream Fic(_NomFichierDonn.c_str(), std::ios::out ) ;
2256   if (Fic.is_open() == true) { Fic << data << std::endl ; }
2257   Fic.close() ;
2258 //
2259 }
2260 //===============================================================================
2261 int HomardDriver::ExecuteHomard()
2262 {
2263   MESSAGE("ExecuteHomard");
2264   std::string commande ;
2265   int codret ;
2266   // Copie des Fichiers HOMARD
2267   commande = "cp " + _NomFichierConf + " " + _NomFichierConfBase ;
2268   codret = system(commande.c_str()) ;
2269
2270 // Execution de HOMARD
2271   if ( codret == 0)
2272   {
2273     commande = _HOMARD_Exec.c_str() ;
2274     if ( _NomFichierDonn != "" ) { commande += " < " + _NomFichierDonn ; }
2275     codret = system(commande.c_str());
2276     if ( codret != 0) { MESSAGE ( "Erreur en executant HOMARD : " << codret ); };
2277     _NomFichierDonn = "" ;
2278   };
2279   return codret ;
2280 }
2281
2282 //=============================================================================
2283 //=============================================================================
2284 HOMARD_Gen::HOMARD_Gen()
2285 {
2286   MESSAGE("HOMARD_Gen");
2287 }
2288
2289 //=============================================================================
2290 //=============================================================================
2291 HOMARD_Gen::~HOMARD_Gen()
2292 {
2293   MESSAGE("~HOMARD_Gen");
2294 }
2295 //=============================================================================
2296
2297 //=============================================================================
2298 /*!
2299  *  default constructor:
2300  */
2301 //=============================================================================
2302 HOMARD_Hypothesis::HOMARD_Hypothesis():
2303   _Name(""), _NomCasCreation(""),
2304   _TypeAdap(-1), _TypeRaff(0), _TypeDera(0),
2305   _Field(""),
2306   _TypeThR(0), _TypeThC(0),
2307   _ThreshR(0), _ThreshC(0),
2308   _UsField(0), _UsCmpI(0), _TypeFieldInterp(0),
2309
2310   _NivMax(-1), _DiamMin(-1.0), _AdapInit(0), _ExtraOutput(1)
2311 {
2312   MESSAGE("HOMARD_Hypothesis");
2313 }
2314
2315 //=============================================================================
2316 /*!
2317  */
2318 //=============================================================================
2319 HOMARD_Hypothesis::~HOMARD_Hypothesis()
2320 {
2321   MESSAGE("~HOMARD_Hypothesis");
2322 }
2323 //=============================================================================
2324 //=============================================================================
2325 // Generalites
2326 //=============================================================================
2327 //=============================================================================
2328 void HOMARD_Hypothesis::SetName( const char* Name )
2329 {
2330   _Name = std::string( Name );
2331 }
2332 //=============================================================================
2333 std::string HOMARD_Hypothesis::GetName() const
2334 {
2335   return _Name;
2336 }
2337 //=============================================================================
2338 //=============================================================================
2339 // Caracteristiques
2340 //=============================================================================
2341 //=============================================================================
2342 void HOMARD_Hypothesis::SetAdapType( int TypeAdap )
2343 {
2344   VERIFICATION( (TypeAdap>=-1) && (TypeAdap<=1) );
2345   _TypeAdap = TypeAdap;
2346 }
2347 //=============================================================================
2348 int HOMARD_Hypothesis::GetAdapType() const
2349 {
2350   return _TypeAdap;
2351 }
2352 //=============================================================================
2353 void HOMARD_Hypothesis::SetRefinTypeDera( int TypeRaff, int TypeDera )
2354 {
2355   VERIFICATION( (TypeRaff>=-1) && (TypeRaff<=1) );
2356   _TypeRaff = TypeRaff;
2357   VERIFICATION( (TypeDera>=-1) && (TypeDera<=1) );
2358   _TypeDera = TypeDera;
2359 }
2360 //=============================================================================
2361 int HOMARD_Hypothesis::GetRefinType() const
2362 {
2363   return _TypeRaff;
2364 }
2365 //=============================================================================
2366 int HOMARD_Hypothesis::GetUnRefType() const
2367 {
2368   return _TypeDera;
2369 }
2370 //=============================================================================
2371 void HOMARD_Hypothesis::SetField( const char* FieldName )
2372 {
2373   _Field = std::string( FieldName );
2374   MESSAGE( "SetField : FieldName = " << FieldName );
2375 }
2376 //=============================================================================
2377 std::string HOMARD_Hypothesis::GetFieldName() const
2378 {
2379   return _Field;
2380 }
2381 //=============================================================================
2382 void HOMARD_Hypothesis::SetUseField( int UsField )
2383 {
2384   VERIFICATION( (UsField>=0) && (UsField<=1) );
2385   _UsField = UsField;
2386 }
2387 //=============================================================================
2388 int HOMARD_Hypothesis::GetUseField() const
2389 {
2390   return _UsField;
2391 }
2392 //=============================================================================
2393 void HOMARD_Hypothesis::SetUseComp( int UsCmpI )
2394 {
2395   MESSAGE ("SetUseComp pour UsCmpI = "<<UsCmpI) ;
2396   VERIFICATION( (UsCmpI>=0) && (UsCmpI<=2) );
2397   _UsCmpI = UsCmpI;
2398 }
2399 //=============================================================================
2400 int HOMARD_Hypothesis::GetUseComp() const
2401 {
2402   return _UsCmpI;
2403 }
2404 //=============================================================================
2405 void HOMARD_Hypothesis::AddComp( const char* NomComp )
2406 {
2407 // On commence par supprimer la composante au cas ou elle aurait deja ete inseree
2408 // Cela peut se produire dans un schema YACS quand on repasse plusieurs fois par la
2409 // definition de l'hypothese
2410   SupprComp( NomComp ) ;
2411 // Insertion veritable
2412   _ListComp.push_back( std::string( NomComp ) );
2413 }
2414 //=============================================================================
2415 void HOMARD_Hypothesis::SupprComp( const char* NomComp )
2416 {
2417   MESSAGE ("SupprComp pour "<<NomComp) ;
2418   std::list<std::string>::iterator it = find( _ListComp.begin(), _ListComp.end(), NomComp );
2419   if ( it != _ListComp.end() ) { it = _ListComp.erase( it ); }
2420 }
2421 //=============================================================================
2422 void HOMARD_Hypothesis::SupprComps()
2423 {
2424   _ListComp.clear();
2425 }
2426 //=============================================================================
2427 const std::list<std::string>& HOMARD_Hypothesis::GetComps() const
2428 {
2429   return _ListComp;
2430 }
2431 //=============================================================================
2432 void HOMARD_Hypothesis::SetRefinThr( int TypeThR, double ThreshR )
2433 {
2434   MESSAGE( "SetRefinThr : TypeThR = " << TypeThR << ", ThreshR = " << ThreshR );
2435   VERIFICATION( (TypeThR>=0) && (TypeThR<=4) );
2436   _TypeThR = TypeThR;
2437   _ThreshR = ThreshR;
2438 }
2439 //=============================================================================
2440 int HOMARD_Hypothesis::GetRefinThrType() const
2441 {
2442   return _TypeThR;
2443 }
2444 //=============================================================================
2445 double HOMARD_Hypothesis::GetThreshR() const
2446 {
2447   return _ThreshR;
2448 }
2449 //=============================================================================
2450 void HOMARD_Hypothesis::SetUnRefThr( int TypeThC, double ThreshC )
2451 {
2452   VERIFICATION( (TypeThC>=0) && (TypeThC<=4) );
2453   _TypeThC = TypeThC;
2454   _ThreshC = ThreshC;
2455 }
2456 //=============================================================================
2457 int HOMARD_Hypothesis::GetUnRefThrType() const
2458 {
2459   return _TypeThC;
2460 }
2461 //=============================================================================
2462 double HOMARD_Hypothesis::GetThreshC() const
2463 {
2464   return _ThreshC;
2465 }
2466 //=============================================================================
2467 void HOMARD_Hypothesis::SetNivMax( int NivMax )
2468 //=============================================================================
2469 {
2470   _NivMax = NivMax;
2471 }
2472 //=============================================================================
2473 const int HOMARD_Hypothesis::GetNivMax() const
2474 //=============================================================================
2475 {
2476   return _NivMax;
2477 }
2478 //=============================================================================
2479 void HOMARD_Hypothesis::SetDiamMin( double DiamMin )
2480 //=============================================================================
2481 {
2482   _DiamMin = DiamMin;
2483 }
2484 //=============================================================================
2485 const double HOMARD_Hypothesis::GetDiamMin() const
2486 //=============================================================================
2487 {
2488   return _DiamMin;
2489 }
2490 //=============================================================================
2491 void HOMARD_Hypothesis::SetAdapInit( int AdapInit )
2492 //=============================================================================
2493 {
2494   _AdapInit = AdapInit;
2495 }
2496 //=============================================================================
2497 const int HOMARD_Hypothesis::GetAdapInit() const
2498 //=============================================================================
2499 {
2500   return _AdapInit;
2501 }
2502 //=============================================================================
2503 void HOMARD_Hypothesis::SetExtraOutput( int ExtraOutput )
2504 //=============================================================================
2505 {
2506   _ExtraOutput = ExtraOutput;
2507 }
2508 //=============================================================================
2509 const int HOMARD_Hypothesis::GetExtraOutput() const
2510 //=============================================================================
2511 {
2512   return _ExtraOutput;
2513 }
2514 //=============================================================================
2515 void HOMARD_Hypothesis::AddGroup( const char* Group)
2516 {
2517 // On commence par supprimer le groupe au cas ou il aurait deja ete insere
2518 // Cela peut se produire dans un schema YACS quand on repasse plusieurs fois par la
2519 // definition de l'hypothese
2520   SupprGroup( Group ) ;
2521 // Insertion veritable
2522   _ListGroupSelected.push_back(Group);
2523 }
2524 //=============================================================================
2525 void HOMARD_Hypothesis::SupprGroup( const char* Group )
2526 {
2527   MESSAGE ("SupprGroup pour "<<Group) ;
2528   std::list<std::string>::iterator it = find( _ListGroupSelected.begin(), _ListGroupSelected.end(), Group );
2529   if ( it != _ListGroupSelected.end() ) { it = _ListGroupSelected.erase( it ); }
2530 }
2531 //=============================================================================
2532 void HOMARD_Hypothesis::SupprGroups()
2533 {
2534   _ListGroupSelected.clear();
2535 }
2536 //=============================================================================
2537 void HOMARD_Hypothesis::SetGroups( const std::list<std::string>& ListGroup )
2538 {
2539   _ListGroupSelected.clear();
2540   std::list<std::string>::const_iterator it = ListGroup.begin();
2541   while(it != ListGroup.end())
2542     _ListGroupSelected.push_back((*it++));
2543 }
2544 //=============================================================================
2545 const std::list<std::string>& HOMARD_Hypothesis::GetGroups() const
2546 {
2547   return _ListGroupSelected;
2548 }
2549 //=============================================================================
2550 // Type d'interpolation des champs :
2551 //   0 : aucun champ n'est interpole
2552 //   1 : tous les champs sont interpoles
2553 //   2 : certains champs sont interpoles
2554 void HOMARD_Hypothesis::SetTypeFieldInterp( int TypeFieldInterp )
2555 {
2556   VERIFICATION( (TypeFieldInterp>=0) && (TypeFieldInterp<=2) );
2557   _TypeFieldInterp = TypeFieldInterp;
2558 }
2559 //=============================================================================
2560 int HOMARD_Hypothesis::GetTypeFieldInterp() const
2561 {
2562   return _TypeFieldInterp;
2563 }
2564 //=============================================================================
2565 void HOMARD_Hypothesis::AddFieldInterpType( const char* FieldInterp, int TypeInterp )
2566 {
2567   MESSAGE ("Dans AddFieldInterpType pour " << FieldInterp << " et TypeInterp = " << TypeInterp) ;
2568 // On commence par supprimer le champ au cas ou il aurait deja ete insere
2569 // Cela peut se produire dans un schema YACS quand on repasse plusieurs fois par la
2570 // definition de l'hypothese
2571   SupprFieldInterp( FieldInterp ) ;
2572 // Insertion veritable
2573 // . Nom du champ
2574   _ListFieldInterp.push_back( std::string( FieldInterp ) );
2575 // . Usage du champ
2576   std::stringstream saux1 ;
2577   saux1 << TypeInterp ;
2578   _ListFieldInterp.push_back( saux1.str() );
2579 // . Indication generale : certains champs sont a interpoler
2580   SetTypeFieldInterp ( 2 ) ;
2581 }
2582 //=============================================================================
2583 void HOMARD_Hypothesis::SupprFieldInterp( const char* FieldInterp )
2584 {
2585   MESSAGE ("Dans SupprFieldInterp pour " << FieldInterp) ;
2586   std::list<std::string>::iterator it = find( _ListFieldInterp.begin(), _ListFieldInterp.end(), FieldInterp ) ;
2587 // Attention a supprimer le nom du champ et le type d'usage
2588   if ( it != _ListFieldInterp.end() )
2589   {
2590     it = _ListFieldInterp.erase( it ) ;
2591     it = _ListFieldInterp.erase( it ) ;
2592   }
2593 // Decompte du nombre de champs restant a interpoler
2594   it = _ListFieldInterp.begin() ;
2595   int cpt = 0 ;
2596   while(it != _ListFieldInterp.end())
2597   {
2598     cpt += 1 ;
2599     (*it++);
2600   }
2601   MESSAGE("Nombre de champ restants = "<<cpt/2);
2602 // . Indication generale : aucun champ ne reste a interpoler
2603   if ( cpt == 0 )
2604   {
2605     SetTypeFieldInterp ( 0 ) ;
2606   }
2607 }
2608 //=============================================================================
2609 void HOMARD_Hypothesis::SupprFieldInterps()
2610 {
2611   MESSAGE ("SupprFieldInterps") ;
2612   _ListFieldInterp.clear();
2613 // . Indication generale : aucun champ ne reste a interpoler
2614   SetTypeFieldInterp ( 0 ) ;
2615 }
2616 //=============================================================================
2617 const std::list<std::string>& HOMARD_Hypothesis::GetFieldInterps() const
2618 {
2619   return _ListFieldInterp;
2620 }
2621 //=============================================================================
2622 //=============================================================================
2623 // Liens avec les autres structures
2624 //=============================================================================
2625 //=============================================================================
2626 void HOMARD_Hypothesis::SetCaseCreation( const char* NomCasCreation )
2627 {
2628   _NomCasCreation = std::string( NomCasCreation );
2629 }
2630 //=============================================================================
2631 std::string HOMARD_Hypothesis::GetCaseCreation() const
2632 {
2633   return _NomCasCreation;
2634 }
2635 //=============================================================================
2636 void HOMARD_Hypothesis::LinkIteration( const char* NomIteration )
2637 {
2638   _ListIter.push_back( std::string( NomIteration ) );
2639 }
2640 //=============================================================================
2641 void HOMARD_Hypothesis::UnLinkIteration( const char* NomIteration )
2642 {
2643   std::list<std::string>::iterator it = find( _ListIter.begin(), _ListIter.end(), NomIteration ) ;
2644   if ( it != _ListIter.end() )
2645   {
2646     MESSAGE ("Dans UnLinkIteration pour " << NomIteration) ;
2647     it = _ListIter.erase( it ) ;
2648   }
2649 }
2650 //=============================================================================
2651 void HOMARD_Hypothesis::UnLinkIterations()
2652 {
2653   _ListIter.clear();
2654 }
2655 //=============================================================================
2656 const std::list<std::string>& HOMARD_Hypothesis::GetIterations() const
2657 {
2658   return _ListIter;
2659 }
2660 //=============================================================================
2661 void HOMARD_Hypothesis::AddZone( const char* NomZone, int TypeUse )
2662 {
2663   MESSAGE ("Dans AddZone pour " << NomZone << " et TypeUse = " << TypeUse) ;
2664 // On commence par supprimer la zone au cas ou elle aurait deja ete inseree
2665 // Cela peut se produire dans un schema YACS quand on repasse plusieurs fois par la
2666 // definition de l'hypothese
2667   SupprZone( NomZone ) ;
2668 // Insertion veritable
2669 // . Nom de la zone
2670   _ListZone.push_back( std::string( NomZone ) );
2671 // . Usage de la zone
2672   std::stringstream saux1 ;
2673   saux1 << TypeUse ;
2674   _ListZone.push_back( saux1.str() );
2675 }
2676 //=============================================================================
2677 void HOMARD_Hypothesis::SupprZone( const char* NomZone )
2678 {
2679   MESSAGE ("Dans SupprZone pour " << NomZone) ;
2680   std::list<std::string>::iterator it = find( _ListZone.begin(), _ListZone.end(), NomZone );
2681 // Attention a supprimer le nom de zone et le type d'usage
2682   if ( it != _ListZone.end() )
2683   {
2684     it = _ListZone.erase( it );
2685     it = _ListZone.erase( it );
2686   }
2687 }
2688 //=============================================================================
2689 void HOMARD_Hypothesis::SupprZones()
2690 {
2691   _ListZone.clear();
2692 }
2693 //=============================================================================
2694 const std::list<std::string>& HOMARD_Hypothesis::GetZones() const
2695 {
2696   return _ListZone;
2697 }
2698
2699 //=============================================================================
2700 /*!
2701  *  default constructor:
2702  */
2703 //=============================================================================
2704 HOMARD_Iteration::HOMARD_Iteration():
2705   _Name( "" ), _Etat( 0 ),
2706  _NumIter( -1 ),
2707   _NomMesh( "" ), _MeshFile( "" ),
2708   _FieldFile( "" ), _TimeStep( -1 ), _Rank( -1 ),
2709   _LogFile( "" ),
2710   _IterParent( "" ),
2711   _NomHypo( "" ), _NomCas( "" ), _NomDir( "" ),
2712   _FileInfo( "" ),
2713  _MessInfo( 1 )
2714 {
2715   MESSAGE("HOMARD_Iteration");
2716 }
2717 //=============================================================================
2718 /*!
2719  *
2720  */
2721 //=============================================================================
2722 HOMARD_Iteration::~HOMARD_Iteration()
2723 {
2724   MESSAGE("~HOMARD_Iteration");
2725 }
2726 //=============================================================================
2727 //=============================================================================
2728 // Generalites
2729 //=============================================================================
2730 //=============================================================================
2731 void HOMARD_Iteration::SetName( const char* Name )
2732 {
2733   _Name = std::string( Name );
2734 }
2735 //=============================================================================
2736 std::string HOMARD_Iteration::GetName() const
2737 {
2738   return _Name;
2739 }
2740 //=============================================================================
2741 //=============================================================================
2742 // Caracteristiques
2743 //=============================================================================
2744 //=============================================================================
2745 void HOMARD_Iteration::SetDirNameLoc( const char* NomDir )
2746 {
2747   _NomDir = std::string( NomDir );
2748 }
2749 //=============================================================================
2750 std::string HOMARD_Iteration::GetDirNameLoc() const
2751 {
2752    return _NomDir;
2753 }
2754 //=============================================================================
2755 void HOMARD_Iteration::SetNumber( int NumIter )
2756 {
2757   _NumIter = NumIter;
2758 }
2759 //=============================================================================
2760 int HOMARD_Iteration::GetNumber() const
2761 {
2762   return _NumIter;
2763 }
2764 //=============================================================================
2765 void HOMARD_Iteration::SetState( int etat )
2766 {
2767   _Etat = etat;
2768 }
2769 //=============================================================================
2770 int HOMARD_Iteration::GetState() const
2771 {
2772   return _Etat;
2773 }
2774 //=============================================================================
2775 void HOMARD_Iteration::SetMeshName( const char* NomMesh )
2776 {
2777   _NomMesh = std::string( NomMesh );
2778 }
2779 //=============================================================================
2780 std::string HOMARD_Iteration::GetMeshName() const
2781 {
2782   return _NomMesh;
2783 }
2784 //=============================================================================
2785 void HOMARD_Iteration::SetMeshFile( const char* MeshFile )
2786 {
2787   _MeshFile = std::string( MeshFile );
2788 }
2789 //=============================================================================
2790 std::string HOMARD_Iteration::GetMeshFile() const
2791 {
2792   return _MeshFile;
2793 }
2794 //=============================================================================
2795 void HOMARD_Iteration::SetFieldFile( const char* FieldFile )
2796 {
2797   _FieldFile = std::string( FieldFile );
2798 }
2799 //=============================================================================
2800 std::string HOMARD_Iteration::GetFieldFile() const
2801 {
2802   return _FieldFile;
2803 }
2804 //=============================================================================
2805 // Instants pour le champ de pilotage
2806 //=============================================================================
2807 void HOMARD_Iteration::SetTimeStep( int TimeStep )
2808 {
2809   _TimeStep = TimeStep;
2810 }
2811 //=============================================================================
2812 void HOMARD_Iteration::SetTimeStepRank( int TimeStep, int Rank )
2813 {
2814   _TimeStep = TimeStep;
2815   _Rank = Rank;
2816 }
2817 //=============================================================================
2818 void HOMARD_Iteration::SetTimeStepRankLast()
2819 {
2820   _TimeStep = -2;
2821 }
2822 //=============================================================================
2823 int HOMARD_Iteration::GetTimeStep() const
2824 {
2825   return _TimeStep;
2826 }
2827 //=============================================================================
2828 int HOMARD_Iteration::GetRank() const
2829 {
2830   return _Rank;
2831 }
2832 //=============================================================================
2833 // Instants pour un champ a interpoler
2834 //=============================================================================
2835 void HOMARD_Iteration::SetFieldInterpTimeStep( const char* FieldInterp, int TimeStep )
2836 {
2837   SetFieldInterpTimeStepRank( FieldInterp, TimeStep, TimeStep ) ;
2838 }
2839 //=============================================================================
2840 void HOMARD_Iteration::SetFieldInterpTimeStepRank( const char* FieldInterp, int TimeStep, int Rank )
2841 {
2842   MESSAGE("Champ " << FieldInterp << ", hypothese " << _NomHypo )
2843 // Verification de la presence du champ dans l'hypothese
2844   std::list<std::string>::iterator it = find( _ListFieldInterp.begin(), _ListFieldInterp.end(), FieldInterp );
2845   if ( it == _ListFieldInterp.end() )
2846   {
2847     INFOS("Champ " << FieldInterp << " ; hypothese " << _NomHypo )
2848     VERIFICATION("Le champ est inconnu dans l'hypothese associee a cette iteration." == 0);
2849   }
2850
2851 // . Nom du champ
2852   _ListFieldInterpTSR.push_back( std::string( FieldInterp ) );
2853 // . Pas de temps
2854   std::stringstream saux1 ;
2855   saux1 << TimeStep ;
2856   _ListFieldInterpTSR.push_back( saux1.str() );
2857 // . Numero d'ordre
2858   std::stringstream saux2 ;
2859   saux2 << Rank ;
2860   _ListFieldInterpTSR.push_back( saux2.str() );
2861 }
2862 //=============================================================================
2863 const std::list<std::string>& HOMARD_Iteration::GetFieldInterpsTimeStepRank() const
2864 {
2865   return _ListFieldInterpTSR;
2866 }
2867 //=============================================================================
2868 void HOMARD_Iteration::SetFieldInterp( const char* FieldInterp )
2869 {
2870   _ListFieldInterp.push_back( std::string( FieldInterp ) );
2871 }
2872 //=============================================================================
2873 const std::list<std::string>& HOMARD_Iteration::GetFieldInterps() const
2874 {
2875   return _ListFieldInterp;
2876 }
2877 //=============================================================================
2878 void HOMARD_Iteration::SupprFieldInterps()
2879 {
2880   _ListFieldInterp.clear();
2881 }
2882 //=============================================================================
2883 void HOMARD_Iteration::SetLogFile( const char* LogFile )
2884 {
2885   _LogFile = std::string( LogFile );
2886 }
2887 //=============================================================================
2888 std::string HOMARD_Iteration::GetLogFile() const
2889 {
2890   return _LogFile;
2891 }
2892 //=============================================================================
2893 void HOMARD_Iteration::SetFileInfo( const char* FileInfo )
2894 {
2895   _FileInfo = std::string( FileInfo );
2896 }
2897 //=============================================================================
2898 std::string HOMARD_Iteration::GetFileInfo() const
2899 {
2900   return _FileInfo;
2901 }
2902 //=============================================================================
2903 //=============================================================================
2904 // Liens avec les autres iterations
2905 //=============================================================================
2906 //=============================================================================
2907 void HOMARD_Iteration::LinkNextIteration( const char* NomIteration )
2908 {
2909   _mesIterFilles.push_back( std::string( NomIteration ) );
2910 }
2911 //=============================================================================
2912 void HOMARD_Iteration::UnLinkNextIteration( const char* NomIteration )
2913 {
2914   std::list<std::string>::iterator it = find( _mesIterFilles.begin(), _mesIterFilles.end(), NomIteration ) ;
2915   if ( it != _mesIterFilles.end() )
2916   {
2917     MESSAGE ("Dans UnLinkNextIteration pour " << NomIteration) ;
2918     it = _mesIterFilles.erase( it ) ;
2919   }
2920 }
2921 //=============================================================================
2922 void HOMARD_Iteration::UnLinkNextIterations()
2923 {
2924   _mesIterFilles.clear();
2925 }
2926 //=============================================================================
2927 const std::list<std::string>& HOMARD_Iteration::GetIterations() const
2928 {
2929   return _mesIterFilles;
2930 }
2931 //=============================================================================
2932 void HOMARD_Iteration::SetIterParentName( const char* IterParent )
2933 {
2934   _IterParent = IterParent;
2935 }
2936 //=============================================================================
2937 std::string HOMARD_Iteration::GetIterParentName() const
2938 {
2939   return _IterParent;
2940 }
2941 //=============================================================================
2942 //=============================================================================
2943 // Liens avec les autres structures
2944 //=============================================================================
2945 //=============================================================================
2946 void HOMARD_Iteration::SetCaseName( const char* NomCas )
2947 {
2948   _NomCas = std::string( NomCas );
2949 }
2950 //=============================================================================
2951 std::string HOMARD_Iteration::GetCaseName() const
2952 {
2953   return _NomCas;
2954 }
2955 //=============================================================================
2956 void HOMARD_Iteration::SetHypoName( const char* NomHypo )
2957 {
2958   _NomHypo = std::string( NomHypo );
2959 }
2960 //=============================================================================
2961 std::string HOMARD_Iteration::GetHypoName() const
2962 {
2963   return _NomHypo;
2964 }
2965 //=============================================================================
2966 //=============================================================================
2967 // Divers
2968 //=============================================================================
2969 //=============================================================================
2970 void HOMARD_Iteration::SetInfoCompute( int MessInfo )
2971 {
2972   _MessInfo = MessInfo;
2973 }
2974 //=============================================================================
2975 int HOMARD_Iteration::GetInfoCompute() const
2976 {
2977   return _MessInfo;
2978 }
2979
2980
2981 /*!
2982  * \brief Relocate nodes to lie on geometry
2983  *  \param [in] theInputMedFile - a MED file holding a mesh including nodes that will be
2984  *         moved onto the geometry
2985  *  \param [in] theOutputMedFile - a MED file to create, that will hold a modified mesh
2986  *  \param [in] theInputNodeFiles - an array of names of files describing groups of nodes that
2987  *         will be moved onto the geometry
2988  *  \param [in] theXaoFileName - a path to a file in XAO format holding the geometry and
2989  *         the geometrical groups.
2990  *  \param [in] theIsParallel - if \c true, all processors are used to treat boundary shapes
2991  *          in parallel.
2992  */
2993 void FrontTrack::track( const std::string&                 theInputMedFile,
2994                         const std::string&                 theOutputMedFile,
2995                         const std::vector< std::string > & theInputNodeFiles,
2996                         const std::string&                 theXaoFileName,
2997                         bool                               theIsParallel )
2998 {
2999   // check arguments
3000 #ifdef _DEBUG_
3001   std::cout << "FrontTrack::track" << std::endl;
3002 #endif
3003
3004   if ( theInputNodeFiles.empty() )
3005     return;
3006
3007 #ifdef _DEBUG_
3008   std::cout << "Input MED file: " << theInputMedFile << std::endl;
3009 #endif
3010   if ( !FT_Utils::fileExists( theInputMedFile ))
3011     throw std::invalid_argument( "Input MED file does not exist: " + theInputMedFile );
3012
3013 #ifdef _DEBUG_
3014   std::cout << "Output MED file: " << theOutputMedFile << std::endl;
3015 #endif
3016   if ( !FT_Utils::canWrite( theOutputMedFile ))
3017     throw std::invalid_argument( "Can't create the output MED file: " + theOutputMedFile );
3018
3019   std::vector< std::string > theNodeFiles ;
3020   for ( size_t i = 0; i < theInputNodeFiles.size(); ++i )
3021   {
3022 #ifdef _DEBUG_
3023     std::cout << "Initial input node file #"<<i<<": " << theInputNodeFiles[i] << std::endl;
3024 #endif
3025     if ( !FT_Utils::fileExists( theInputNodeFiles[i] ))
3026       throw std::invalid_argument( "Input node file does not exist: " + theInputNodeFiles[i] );
3027     // the name of the groupe on line #1, then the numbers of nodes on line #>1
3028     // keep only files with more than 1 line:
3029     std::ifstream fichier(theInputNodeFiles[i].c_str());
3030     std::string s;
3031     unsigned int nb_lines = 0;
3032     while(std::getline(fichier,s)) ++nb_lines;
3033 //     std::cout << ". nb_lines: " << nb_lines << std::endl;
3034     if ( nb_lines >= 2 ) { theNodeFiles.push_back( theInputNodeFiles[i] ); }
3035   }
3036 #ifdef _DEBUG_
3037   for ( size_t i = 0; i < theNodeFiles.size(); ++i )
3038   { std::cout << "Valid input node file #"<<i<<": " << theNodeFiles[i] << std::endl; }
3039 #endif
3040
3041 #ifdef _DEBUG_
3042   std::cout << "XAO file: " << theXaoFileName << std::endl;
3043 #endif
3044   if ( !FT_Utils::fileExists( theXaoFileName ))
3045     throw std::invalid_argument( "Input XAO file does not exist: " + theXaoFileName );
3046
3047   // read a mesh
3048
3049 #ifdef _DEBUG_
3050   std::cout << "Lecture du maillage" << std::endl;
3051 #endif
3052   MEDCoupling::MCAuto< MEDCoupling::MEDFileUMesh >
3053     mfMesh( MEDCoupling::MEDFileUMesh::New( theInputMedFile ));
3054   if ( mfMesh.isNull() )
3055     throw std::invalid_argument( "Failed to read the input MED file: " + theInputMedFile );
3056
3057   MEDCoupling::DataArrayDouble * nodeCoords = mfMesh->getCoords();
3058   if ( !nodeCoords || nodeCoords->empty() )
3059     throw std::invalid_argument( "No nodes in the input mesh" );
3060
3061
3062   // read a geometry
3063
3064 #ifdef _DEBUG_
3065   std::cout << "Lecture de la geometrie" << std::endl;
3066 #endif
3067   XAO::Xao xao;
3068   if ( !xao.importXAO( theXaoFileName ) || !xao.getGeometry() )
3069     throw std::invalid_argument( "Failed to read the XAO input file: " + theXaoFileName );
3070
3071 #ifdef _DEBUG_
3072   std::cout << "Conversion en BREP" << std::endl;
3073 #endif
3074   XAO::BrepGeometry* xaoGeom = dynamic_cast<XAO::BrepGeometry*>( xao.getGeometry() );
3075   if ( !xaoGeom || xaoGeom->getTopoDS_Shape().IsNull() )
3076     throw std::invalid_argument( "Failed to get a BREP shape from the XAO input file" );
3077
3078
3079   // read groups of nodes and associate them with boundary shapes using names (no projection so far)
3080
3081 #ifdef _DEBUG_
3082   std::cout << "Lecture des groupes" << std::endl;
3083 #endif
3084   FT_NodeGroups nodeGroups;
3085   nodeGroups.read( theNodeFiles, &xao, nodeCoords );
3086 #ifdef _DEBUG_
3087   std::cout << "Nombre de groupes : " << nodeGroups.nbOfGroups() << std::endl;
3088 #endif
3089
3090   // project nodes to the boundary shapes and change their coordinates
3091
3092 #ifdef _DEBUG_
3093   std::cout << "Projection des noeuds, theIsParallel=" << theIsParallel << std::endl;
3094 #endif
3095   OSD_Parallel::For( 0, nodeGroups.nbOfGroups(), nodeGroups, !theIsParallel );
3096
3097   // save the modified mesh
3098
3099 #ifdef _DEBUG_
3100   std::cout << "Ecriture du maillage" << std::endl;
3101 #endif
3102   const int erase = 2;
3103   mfMesh->write( theOutputMedFile, /*mode=*/erase );
3104
3105   if ( !nodeGroups.isOK() )
3106     throw std::runtime_error("Unable to project some nodes");
3107 }
3108
3109   //================================================================================
3110   /*!
3111    * \brief Initialize FT_Projector's with all sub-shapes of given type
3112    *  \param [in] theMainShape - the shape to explore
3113    *  \param [in] theSubType - the type of sub-shapes
3114    *  \param [out] theProjectors - the projectors
3115    */
3116   //================================================================================
3117
3118   void getProjectors( const TopoDS_Shape&           theMainShape,
3119                       const TopAbs_ShapeEnum        theSubType,
3120                       std::vector< FT_Projector > & theProjectors )
3121   {
3122     TopTools_IndexedMapOfShape subShapes;
3123     TopExp::MapShapes( theMainShape, theSubType, subShapes );
3124 #ifdef _DEBUG_
3125     std::cout << ". Nombre de subShapes : " << subShapes.Size() << std::endl;
3126 #endif
3127
3128     theProjectors.resize( subShapes.Size() );
3129     for ( int i = 1; i <= subShapes.Size(); ++i )
3130       theProjectors[ i-1 ].setBoundaryShape( subShapes( i ));
3131   }
3132
3133 //================================================================================
3134 /*!
3135  * \brief Load node groups from files
3136  *  \param [in] theNodeFiles - an array of names of files describing groups of nodes that
3137  *         will be moved onto geometry
3138  *  \param [in] theXaoGeom - the whole geometry to project on
3139  *  \param [inout] theNodeCoords - array of node coordinates
3140  */
3141 //================================================================================
3142
3143 void FT_NodeGroups::read( const std::vector< std::string >& theNodeFiles,
3144                           const XAO::Xao*                   theXao,
3145                           MEDCoupling::DataArrayDouble*     theNodeCoords )
3146 {
3147   // get projectors for all boundary sub-shapes;
3148   // index of a projector in the vector corresponds to a XAO index of a sub-shape
3149   XAO::BrepGeometry* xaoGeom = dynamic_cast<XAO::BrepGeometry*>( theXao->getGeometry() );
3150   getProjectors( xaoGeom->getTopoDS_Shape(), TopAbs_EDGE, _projectors[0] );
3151   getProjectors( xaoGeom->getTopoDS_Shape(), TopAbs_FACE, _projectors[1] );
3152
3153   _nodesOnGeom.resize( theNodeFiles.size() );
3154
3155   // read node IDs and look for projectors to boundary sub-shapes by group name
3156   FT_Utils::XaoGroups xaoGroups( theXao );
3157   for ( size_t i = 0; i < theNodeFiles.size(); ++i )
3158   {
3159     _nodesOnGeom[i].read( theNodeFiles[i], xaoGroups, theNodeCoords, _projectors );
3160   }
3161 }
3162
3163 //================================================================================
3164 /*!
3165  * \brief Project and move nodes of a given group of nodes
3166  */
3167 //================================================================================
3168
3169 void FT_NodeGroups::projectAndMove( const int groupIndex )
3170 {
3171   _nodesOnGeom[ groupIndex ].projectAndMove();
3172 }
3173
3174 //================================================================================
3175 /*!
3176  * \brief Return true if all nodes were successfully relocated
3177  */
3178 //================================================================================
3179
3180 bool FT_NodeGroups::isOK() const
3181 {
3182   for ( size_t i = 0; i < _nodesOnGeom.size(); ++i )
3183     if ( ! _nodesOnGeom[ i ].isOK() )
3184       return false;
3185
3186   return true;
3187 }
3188
3189 //================================================================================
3190 /*!
3191  * \brief Print some statistics on node groups
3192  */
3193 //================================================================================
3194
3195 void FT_NodeGroups::dumpStat() const
3196 {
3197   for ( size_t i = 0; i < _nodesOnGeom.size(); ++i )
3198   {
3199     std::cout << _nodesOnGeom[i].getShapeDim() << "D "
3200               << _nodesOnGeom[i].nbNodes() << " nodes" << std::endl;
3201   }
3202 }
3203
3204   /*!
3205    * \brief Close a file at destruction
3206    */
3207   struct FileCloser
3208   {
3209     FILE * _file;
3210
3211     FileCloser( FILE * file ): _file( file ) {}
3212     ~FileCloser() { if ( _file ) ::fclose( _file ); }
3213   };
3214
3215 //================================================================================
3216 /*!
3217  * \brief Read node ids from a file and find shapes for projection
3218  *  \param [in] theNodeFile - a name of file holding IDs of nodes that
3219  *         will be moved onto geometry
3220  *  \param [in] theXaoGroups - a tool returning FT_Projector's by XAO group name
3221  *  \param [inout] theNodeCoords - array of node coordinates
3222  *  \param [in] theAllProjectorsByDim - all projectors of 2 dimensions, ordered so that
3223  *         a vector index corresponds to a XAO sub-shape ID
3224  */
3225 //================================================================================
3226
3227 void FT_NodesOnGeom::read( const std::string&            theNodeFile,
3228                            const FT_Utils::XaoGroups&    theXaoGroups,
3229                            MEDCoupling::DataArrayDouble* theNodeCoords,
3230                            std::vector< FT_Projector > * theAllProjectorsByDim )
3231 {
3232   _nodeCoords = theNodeCoords;
3233
3234   FILE * file = ::fopen( theNodeFile.c_str(), "r" );
3235   if ( !file )
3236     throw std::invalid_argument( "Can't open an input node file: " + theNodeFile );
3237
3238   FileCloser fileCloser( file );
3239
3240   // -------------------------------------
3241   // get shape dimension by the file name
3242   // -------------------------------------
3243
3244   // hope the file name is something like "frnD.**" with n in (1,2)
3245   int dimPos = theNodeFile.size() - 5;
3246   if ( theNodeFile[ dimPos ] == '2' )
3247     _shapeDim = 2;
3248   else if ( theNodeFile[ dimPos ] == '1' )
3249     _shapeDim = 1;
3250   else
3251     throw std::invalid_argument( "Can't define dimension by node file name " + theNodeFile );
3252 #ifdef _DEBUG_
3253   std::cout << ". Dimension of the file " << theNodeFile << ": " << _shapeDim << std::endl;
3254 #endif
3255
3256   // -------------------------------------
3257   // read geom group names; several lines
3258   // -------------------------------------
3259
3260   std::vector< std::string > geomNames;
3261
3262   const int maxLineLen = 256;
3263   char line[ maxLineLen ];
3264
3265   long int pos = ::ftell( file );
3266   while ( ::fgets( line, maxLineLen, file )) // read a line
3267   {
3268     if ( ::feof( file ))
3269     {
3270       return; // no nodes in the file
3271     }
3272
3273     // check if the line describes node ids in format 3I10 (e.g. "       120         1        43\n")
3274     size_t lineLen = strlen( line );
3275     if ( lineLen  >= 31        &&
3276          ::isdigit( line[9] )  &&
3277          line[10] == ' '       &&
3278          ::isdigit( line[19] ) &&
3279          line[20] == ' '       &&
3280          ::isdigit( line[29] ) &&
3281          ::isspace( line[30] ))
3282       break;
3283
3284     geomNames.push_back( line + 1 ); // skip the 1st white space
3285
3286     pos = ::ftell( file ); // remember the position to return if the next line holds node ids
3287   }
3288
3289   ::fseek( file, pos, SEEK_SET ); // return to the 1st line holding nodes ids
3290
3291
3292   // --------------
3293   // read node ids
3294   // --------------
3295
3296   FT_NodeToMove nodeIds;
3297   std::vector< int > ids;
3298
3299   const int nbNodes = theNodeCoords->getNumberOfTuples(); // to check validity of node IDs
3300
3301   while ( ::fgets( line, maxLineLen, file )) // read a line
3302   {
3303     // find node ids in the line
3304
3305     char *beg = line, *end = 0;
3306     long int id;
3307
3308     ids.clear();
3309     while (( id = ::strtol( beg, &end, 10 )) &&
3310            ( beg != end ))
3311     {
3312       ids.push_back( id );
3313       if ( id > nbNodes )
3314         throw std::invalid_argument( "Too large node ID: " + FT_Utils::toStr( id ));
3315       beg = end;
3316     }
3317
3318     if ( ids.size() >= 3 )
3319     {
3320       std::vector< int >::iterator i = ids.begin();
3321       nodeIds._nodeToMove = *i;
3322       nodeIds._neighborNodes.assign( ++i, ids.end() );
3323
3324       _nodes.push_back( nodeIds );
3325     }
3326
3327     if ( ::feof( file ))
3328       break;
3329   }
3330
3331   // -----------------------------------------------------------------
3332   // try to find FT_Projector's to boundary sub-shapes by group names
3333   // -----------------------------------------------------------------
3334
3335   _allProjectors = & theAllProjectorsByDim[ _shapeDim - 1 ];
3336
3337   _projectors.reserve( geomNames.size() );
3338   std::vector< const FT_Projector* >  projectors;
3339
3340   for ( size_t i = 0; i < geomNames.size(); ++i )
3341   {
3342     std::string & groupName = geomNames[i];
3343 #ifdef _DEBUG_
3344     std::cout << ". Group name: " << groupName << std::endl;
3345 #endif
3346
3347     // remove trailing white spaces
3348     for ( int iC = groupName.size() - 1; iC >= 0; --iC )
3349     {
3350       if ( ::isspace( groupName[iC] ) )
3351         groupName.resize( iC );
3352       else
3353         break;
3354     }
3355     if ( groupName.empty() )
3356       continue;
3357
3358     _groupNames.push_back( groupName ); // keep _groupNames for easier debug :)
3359
3360     // get projectors by group name
3361     theXaoGroups.getProjectors( groupName, _shapeDim,
3362                                 theAllProjectorsByDim[ _shapeDim-1 ], projectors );
3363   }
3364
3365   // ------------------------------
3366   // check the found FT_Projector's
3367   // ------------------------------
3368
3369   if ( projectors.size() == 1 )
3370   {
3371     _projectors.push_back( *projectors[ 0 ]);
3372   }
3373   else
3374   {
3375     Bnd_Box nodesBox;
3376     for ( size_t i = 0; i < _nodes.size(); ++i )
3377       nodesBox.Add( getPoint( _nodes[i]._nodeToMove ));
3378
3379     if ( projectors.size() > 1 )
3380     {
3381       // more than one boundary shape;
3382       // try to filter off unnecessary projectors using a bounding box of nodes
3383       for ( size_t i = 0; i < projectors.size(); ++i )
3384         if ( !nodesBox.IsOut( projectors[ i ]->getBoundingBox() ))
3385           _projectors.push_back( *projectors[ i ]);
3386     }
3387
3388     if ( _projectors.empty() )
3389     {
3390       // select projectors using a bounding box of nodes
3391       std::vector< FT_Projector > & allProjectors = *_allProjectors;
3392       for ( size_t i = 0; i < allProjectors.size(); ++i )
3393         if ( !nodesBox.IsOut( allProjectors[ i ].getBoundingBox() ))
3394           _projectors.push_back( allProjectors[ i ]);
3395
3396       if ( _projectors.empty() && !_nodes.empty() )
3397         throw std::runtime_error("No boundary shape found for nodes in file " + theNodeFile );
3398     }
3399   }
3400
3401   // prepare for projection - create real projectors
3402   for ( size_t i = 0; i < _projectors.size(); ++i )
3403     _projectors[ i ].prepareForProjection();
3404
3405 }
3406
3407 //================================================================================
3408 /*!
3409  * \brief Project nodes to the shapes and move them to new positions
3410  */
3411 //================================================================================
3412
3413 void FT_NodesOnGeom::projectAndMove()
3414 {
3415   _OK = true;
3416 //
3417 // 1. Préalables
3418 //
3419   // check if all the shapes are planar
3420   bool isAllPlanar = true;
3421   for ( size_t i = 0; i < _projectors.size() &&  isAllPlanar; ++i )
3422     isAllPlanar = _projectors[i].isPlanarBoundary();
3423   if ( isAllPlanar )
3424     return;
3425
3426   // set nodes in the order suitable for optimal projection
3427   putNodesInOrder();
3428
3429   // project and move nodes
3430
3431   std::vector< FT_NodeToMove* > notProjectedNodes;
3432   size_t iP, iProjector;
3433   gp_Pnt newXyz;
3434
3435 #ifdef _DEBUG_
3436     std::cout << ".. _projectors.size() = " << _projectors.size() << std::endl;
3437     std::cout << ".. _nodesOrder.size() = " << _nodesOrder.size() << std::endl;
3438 #endif
3439 //
3440 // 2. Calculs
3441 // 2.1. Avec plusieurs shapes
3442 //
3443   if ( _projectors.size() > 1 )
3444   {
3445     // the nodes are to be projected onto several boundary shapes;
3446     // in addition to the projecting, classification on a shape is necessary
3447     // in order to find out on which of the shapes a node is to be projected
3448
3449     iProjector = 0;
3450     for ( size_t i = 0; i < _nodesOrder.size(); ++i )
3451     {
3452       FT_NodeToMove& nn = _nodes[ _nodesOrder[ i ]];
3453       gp_Pnt        xyz = getPoint( nn._nodeToMove );
3454       gp_Pnt       xyz1 = getPoint( nn._neighborNodes[0] );
3455       gp_Pnt       xyz2 = getPoint( nn._neighborNodes[1] );
3456       double   maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.;
3457       if ( _projectors[ iProjector ].projectAndClassify( xyz, maxDist2, newXyz,
3458                                                          nn._params, nn._nearParams ))
3459       {
3460         moveNode( nn._nodeToMove, newXyz );
3461       }
3462       else // a node is not on iProjector-th shape, find the shape it is on
3463       {
3464         for ( iP = 1; iP < _projectors.size(); ++iP ) // check _projectors other than iProjector
3465         {
3466           iProjector = ( iProjector + 1 ) % _projectors.size();
3467           if ( _projectors[ iProjector ].projectAndClassify( xyz, maxDist2, newXyz,
3468                                                              nn._params, nn._nearParams ))
3469           {
3470             moveNode( nn._nodeToMove, newXyz );
3471             break;
3472           }
3473         }
3474         if ( iP == _projectors.size() )
3475         {
3476           notProjectedNodes.push_back( &nn );
3477
3478 #ifdef _DEBUG_
3479           std::cerr << "Warning: no shape found for node " << nn._nodeToMove << std::endl;
3480           if ( !_groupNames.empty() )
3481             std::cerr << "Warning:    group -- " << _groupNames[0] << std::endl;
3482 #endif
3483         }
3484       }
3485     }
3486   }
3487 //
3488 // 2.2. Avec une seule shape
3489 //
3490   else // one shape
3491   {
3492     for ( size_t i = 0; i < _nodesOrder.size(); ++i )
3493     {
3494       FT_NodeToMove& nn = _nodes[ _nodesOrder[ i ]];
3495       gp_Pnt        xyz = getPoint( nn._nodeToMove );
3496       gp_Pnt       xyz1 = getPoint( nn._neighborNodes[0] );
3497       gp_Pnt       xyz2 = getPoint( nn._neighborNodes[1] );
3498
3499 // maxDist2 : le quart du carré de la distance entre les deux voisins du noeud à bouger
3500       double   maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.;
3501 #ifdef _DEBUG_
3502     std::cout << "\n.. maxDist2 = " << maxDist2 << " entre " << nn._neighborNodes[0] << " et " << nn._neighborNodes[1] << " - milieu " << nn._nodeToMove << " - d/2 = " << sqrt(maxDist2) << " - d = " << sqrt(xyz1.SquareDistance( xyz2 )) << std::endl;
3503 #endif
3504       if ( _projectors[ 0 ].project( xyz, maxDist2, newXyz,
3505                                      nn._params, nn._nearParams ))
3506         moveNode( nn._nodeToMove, newXyz );
3507       else
3508         notProjectedNodes.push_back( &nn );
3509     }
3510   }
3511 //
3512 // 3. Bilan
3513 //
3514   if ( !notProjectedNodes.empty() )
3515   {
3516     // project nodes that are not projected by any of _projectors;
3517     // a proper projector is selected by evaluation of a distance between neighbor nodes
3518     // and a shape
3519
3520     std::vector< FT_Projector > & projectors = *_allProjectors;
3521
3522     iProjector = 0;
3523     for ( size_t i = 0; i < notProjectedNodes.size(); ++i )
3524     {
3525       FT_NodeToMove& nn = *notProjectedNodes[ i ];
3526       gp_Pnt        xyz = getPoint( nn._nodeToMove );
3527       gp_Pnt       xyz1 = getPoint( nn._neighborNodes[0] );
3528       gp_Pnt       xyz2 = getPoint( nn._neighborNodes[1] );
3529       double   maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.;
3530       double       tol2 = 1e-6 * maxDist2;
3531
3532       bool ok;
3533       for ( iP = 0; iP < projectors.size(); ++iP )
3534       {
3535         projectors[ iProjector ].prepareForProjection();
3536         projectors[ iProjector ].tryWithoutPrevSolution( true );
3537
3538         if (( ok = projectors[ iProjector ].isOnShape( xyz1, tol2, nn._params, nn._nearParams )) &&
3539             ( ok = projectors[ iProjector ].isOnShape( xyz2, tol2, nn._params, nn._params )))
3540         {
3541           if ( nn._neighborNodes.size() == 4 )
3542           {
3543             gp_Pnt xyz1 = getPoint( nn._neighborNodes[2] );
3544             gp_Pnt xyz2 = getPoint( nn._neighborNodes[3] );
3545             if (( ok = projectors[ iProjector ].isOnShape( xyz1, tol2, nn._params, nn._params )))
3546               ok     = projectors[ iProjector ].isOnShape( xyz2, tol2, nn._params, nn._params );
3547           }
3548         }
3549
3550         if ( ok && projectors[iProjector].project( xyz, maxDist2, newXyz, nn._params, nn._params ))
3551         {
3552           moveNode( nn._nodeToMove, newXyz );
3553           break;
3554         }
3555         iProjector = ( iProjector + 1 ) % projectors.size();
3556       }
3557       if ( iP == projectors.size() )
3558       {
3559         _OK = false;
3560
3561         std::cerr << "Error: not projected node " << nn._nodeToMove << std::endl;
3562       }
3563     }
3564   }
3565 }
3566
3567 //================================================================================
3568 /*!
3569  * \brief Put nodes in the order for optimal projection and set FT_NodeToMove::_nearParams
3570  *        to point to a FT_NodeToMove::_params of a node that will be projected earlier
3571  */
3572 //================================================================================
3573
3574 void FT_NodesOnGeom::putNodesInOrder()
3575 {
3576   if ( !_nodesOrder.empty() )
3577     return;
3578
3579   // check if any of projectors can use parameters of a previously projected node on a shape
3580   // to speed up projection
3581
3582   bool isPrevSolutionUsed = false;
3583   for ( size_t i = 0; i < _projectors.size() &&  !isPrevSolutionUsed; ++i )
3584     isPrevSolutionUsed = _projectors[i].canUsePrevSolution();
3585
3586   if ( !isPrevSolutionUsed )
3587   {
3588     _nodesOrder.resize( _nodes.size() );
3589     for ( size_t i = 0; i < _nodesOrder.size(); ++i )
3590       _nodesOrder[ i ] = i;
3591     return;
3592   }
3593
3594   // make a map to find a neighbor projected node
3595
3596   // map of { FT_NodeToMove::_neighborNodes[i] } to { FT_NodeToMove* };
3597   // here we call FT_NodeToMove a 'link' as this data links a _neighborNodes[i] node to other nodes
3598   typedef NCollection_DataMap< int, std::vector< FT_NodeToMove* > > TNodeIDToLinksMap;
3599   TNodeIDToLinksMap neigborsMap;
3600
3601   int mapSize = ( _shapeDim == 1 ) ? _nodes.size() + 1 : _nodes.size() * 3;
3602   neigborsMap.Clear();
3603   neigborsMap.ReSize( mapSize );
3604
3605   std::vector< FT_NodeToMove* > linkVec, *linkVecPtr;
3606   const int maxNbLinks = ( _shapeDim == 1 ) ? 2 : 6; // usual nb of links
3607
3608   for ( size_t i = 0; i < _nodes.size(); ++i )
3609   {
3610     FT_NodeToMove& nn = _nodes[i];
3611     for ( size_t iN = 0; iN < nn._neighborNodes.size(); ++iN )
3612     {
3613       if ( !( linkVecPtr = neigborsMap.ChangeSeek( nn._neighborNodes[ iN ] )))
3614       {
3615         linkVecPtr = neigborsMap.Bound( nn._neighborNodes[ iN ], linkVec );
3616         linkVecPtr->reserve( maxNbLinks );
3617       }
3618       linkVecPtr->push_back( & nn );
3619     }
3620   }
3621
3622   // fill in _nodesOrder
3623
3624   _nodesOrder.reserve( _nodes.size() );
3625
3626   std::list< FT_NodeToMove* > queue;
3627   queue.push_back( &_nodes[0] );
3628   _nodes[0]._nearParams = _nodes[0]._params; // to avoid re-adding to the queue
3629
3630   while ( !queue.empty() )
3631   {
3632     FT_NodeToMove* nn = queue.front();
3633     queue.pop_front();
3634
3635     _nodesOrder.push_back( nn - & _nodes[0] );
3636
3637     // add neighbors to the queue and set their _nearParams = nn->_params
3638     for ( size_t iN = 0; iN < nn->_neighborNodes.size(); ++iN )
3639     {
3640       std::vector< FT_NodeToMove* >& linkVec = neigborsMap( nn->_neighborNodes[ iN ]);
3641       for ( size_t iL = 0; iL < linkVec.size(); ++iL )
3642       {
3643         FT_NodeToMove* nnn = linkVec[ iL ];
3644         if ( nnn != nn && nnn->_nearParams == 0 )
3645         {
3646           nnn->_nearParams = nn->_params;
3647           queue.push_back( nnn );
3648         }
3649       }
3650     }
3651   }
3652   _nodes[0]._nearParams = 0; // reset
3653 }
3654
3655 //================================================================================
3656 /*!
3657  * \brief Get node coordinates. Node IDs count from a unit
3658  */
3659 //================================================================================
3660
3661 gp_Pnt FT_NodesOnGeom::getPoint( const int nodeID )
3662 {
3663   const size_t dim = _nodeCoords->getNumberOfComponents();
3664   const double * xyz = _nodeCoords->getConstPointer() + ( dim * ( nodeID - 1 ));
3665   return gp_Pnt( xyz[0], xyz[1], dim == 2 ? 0 : xyz[2] );
3666 }
3667
3668 //================================================================================
3669 /*!
3670  * \brief change node coordinates
3671  */
3672 //================================================================================
3673
3674 void FT_NodesOnGeom::moveNode( const int nodeID, const gp_Pnt& newXyz )
3675 {
3676   const size_t dim = _nodeCoords->getNumberOfComponents();
3677   double z, *xyz = _nodeCoords->getPointer() + ( dim * ( nodeID - 1 ));
3678   newXyz.Coord( xyz[0], xyz[1], dim == 2 ? z : xyz[2] );
3679 }
3680
3681 //-----------------------------------------------------------------------------
3682 /*!
3683  * \brief Root class of a projector of a point to a boundary shape
3684  */
3685 struct FT_RealProjector
3686 {
3687   virtual ~FT_RealProjector() {}
3688
3689   /*!
3690    * \brief Project a point to a boundary shape
3691    *  \param [in] point - the point to project
3692    *  \param [out] newSolution - position on the shape (U or UV) found during the projection
3693    *  \param [in] prevSolution - position already found during the projection of a neighbor point
3694    *  \return gp_Pnt - the projection point
3695    */
3696   virtual gp_Pnt project( const gp_Pnt& point,
3697                           double*       newSolution,
3698                           const double* prevSolution = 0) = 0;
3699
3700   /*!
3701    * \brief Project a point to a boundary shape and check if the projection is within
3702    *        the shape boundary
3703    *  \param [in] point - the point to project
3704    *  \param [in] maxDist2 - the maximal allowed square distance between point and projection
3705    *  \param [out] projection - the projection point
3706    *  \param [out] newSolution - position on the shape (U or UV) found during the projection
3707    *  \param [in] prevSolution - position already found during the projection of a neighbor point
3708    *  \return bool - false if the projection point lies out of the shape boundary or
3709    the distance the point and the projection is more than sqrt(maxDist2)
3710   */
3711   virtual bool projectAndClassify( const gp_Pnt& point,
3712                                    const double  maxDist2,
3713                                    gp_Pnt&       projection,
3714                                    double*       newSolution,
3715                                    const double* prevSolution = 0) = 0;
3716
3717   // return true if a previously found solution can be used to speed up the projection
3718
3719   virtual bool canUsePrevSolution() const { return false; }
3720
3721
3722   double _dist; // distance between the point being projected and its projection
3723 };
3724
3725 namespace // actual projection algorithms
3726 {
3727   const double theEPS = 1e-12;
3728
3729   //================================================================================
3730   /*!
3731    * \brief Projector to any curve
3732    */
3733   //================================================================================
3734
3735   struct CurveProjector : public FT_RealProjector
3736   {
3737     BRepAdaptor_Curve   _curve;
3738     double              _tol;
3739     ShapeAnalysis_Curve _projector;
3740     double              _uRange[2];
3741
3742     //-----------------------------------------------------------------------------
3743     CurveProjector( const TopoDS_Edge& e, const double tol ):
3744       _curve( e ), _tol( tol )
3745     {
3746       BRep_Tool::Range( e, _uRange[0], _uRange[1] );
3747     }
3748
3749     //-----------------------------------------------------------------------------
3750     // project a point to the curve
3751     virtual gp_Pnt project( const gp_Pnt& P,
3752                             double*       newSolution,
3753                             const double* prevSolution = 0)
3754     {
3755 #ifdef _DEBUG_
3756     std::cout << ".. project a point to the curve prevSolution = " << prevSolution << std::endl;
3757 #endif
3758       gp_Pnt         proj;
3759       Standard_Real param;
3760
3761       if ( prevSolution )
3762       {
3763         _dist = _projector.NextProject( prevSolution[0], _curve, P, _tol, proj, param );
3764       }
3765       else
3766       {
3767         _dist = _projector.Project( _curve, P, _tol, proj, param, false );
3768       }
3769 #ifdef _DEBUG_
3770     std::cout << "..    _dist : " << _dist << std::endl;
3771 #endif
3772       proj = _curve.Value( param );
3773
3774       newSolution[0] = param;
3775
3776       return proj;
3777     }
3778
3779     //-----------------------------------------------------------------------------
3780     // project a point to a curve and check if the projection is within the curve boundary
3781     virtual bool projectAndClassify( const gp_Pnt& point,
3782                                      const double  maxDist2,
3783                                      gp_Pnt&       projection,
3784                                      double*       newSolution,
3785                                      const double* prevSolution = 0)
3786     {
3787 #ifdef _DEBUG_
3788     std::cout << ".. project a point to a curve and check " << std::endl;
3789 #endif
3790       projection = project( point, newSolution, prevSolution );
3791       return ( _uRange[0] < newSolution[0] && newSolution[0] < _uRange[1] &&
3792                _dist * _dist < maxDist2 );
3793     }
3794
3795     //-----------------------------------------------------------------------------
3796     // return true if a previously found solution can be used to speed up the projection
3797     virtual bool canUsePrevSolution() const { return true; }
3798   };
3799
3800   //================================================================================
3801   /*!
3802    * \brief Projector to a straight curve. Don't project, classify only
3803    */
3804   //================================================================================
3805
3806   struct LineProjector : public FT_RealProjector
3807   {
3808     gp_Pnt _p0, _p1;
3809
3810     //-----------------------------------------------------------------------------
3811     LineProjector( TopoDS_Edge e )
3812     {
3813       e.Orientation( TopAbs_FORWARD );
3814       _p0 = BRep_Tool::Pnt( TopExp::FirstVertex( e ));
3815       _p1 = BRep_Tool::Pnt( TopExp::LastVertex ( e ));
3816     }
3817
3818     //-----------------------------------------------------------------------------
3819     // does nothing
3820     virtual gp_Pnt project( const gp_Pnt& P,
3821                             double*       newSolution,
3822                             const double* prevSolution = 0)
3823     {
3824       return P;
3825     }
3826     //-----------------------------------------------------------------------------
3827     // check if a point lies within the line segment
3828     virtual bool projectAndClassify( const gp_Pnt& point,
3829                                      const double  maxDist2,
3830                                      gp_Pnt&       projection,
3831                                      double*       newSolution,
3832                                      const double* prevSolution = 0)
3833     {
3834       gp_Vec edge( _p0, _p1 );
3835       gp_Vec p0p ( _p0, point  );
3836       double u = ( edge * p0p ) / edge.SquareMagnitude();  // param [0,1] on the edge
3837       projection = ( 1. - u ) * _p0.XYZ() + u * _p1.XYZ(); // projection of the point on the edge
3838       if ( u < 0 || 1 < u )
3839         return false;
3840
3841       // check distance
3842       return point.SquareDistance( projection ) < theEPS * theEPS;
3843     }
3844   };
3845
3846   //================================================================================
3847   /*!
3848    * \brief Projector to a circular edge
3849    */
3850   //================================================================================
3851
3852   struct CircleProjector : public FT_RealProjector
3853   {
3854     gp_Circ _circle;
3855     double _uRange[2];
3856
3857     //-----------------------------------------------------------------------------
3858     CircleProjector( const gp_Circ& c, const double f, const double l ):
3859       _circle( c )
3860     {
3861       _uRange[0] = f;
3862       _uRange[1] = l;
3863     }
3864
3865     //-----------------------------------------------------------------------------
3866     // project a point to the circle
3867     virtual gp_Pnt project( const gp_Pnt& P,
3868                             double*       newSolution,
3869                             const double* prevSolution = 0)
3870     {
3871       // assume that P is already on the the plane of circle, since
3872       // it is in the middle of two points lying on the circle
3873
3874       // move P to the circle
3875       const gp_Pnt& O = _circle.Location();
3876       gp_Vec radiusVec( O, P );
3877       double radius = radiusVec.Magnitude();
3878       if ( radius < std::numeric_limits<double>::min() )
3879         return P; // P in on the axe
3880
3881       gp_Pnt proj = O.Translated( radiusVec.Multiplied( _circle.Radius() / radius ));
3882
3883       _dist = _circle.Radius() - radius;
3884
3885       return proj;
3886     }
3887
3888     //-----------------------------------------------------------------------------
3889     // project and check if a projection lies within the circular edge
3890     virtual bool projectAndClassify( const gp_Pnt& point,
3891                                      const double  maxDist2,
3892                                      gp_Pnt&       projection,
3893                                      double*       newSolution,
3894                                      const double* prevSolution = 0)
3895     {
3896       _dist = -1;
3897       projection = project( point, newSolution );
3898       if ( _dist < 0 || // ?
3899            _dist * _dist > maxDist2 )
3900         return false;
3901
3902       newSolution[0] = ElCLib::Parameter( _circle, projection );
3903       return ( _uRange[0] < newSolution[0] && newSolution[0] < _uRange[1] );
3904     }
3905   };
3906
3907   //================================================================================
3908   /*!
3909    * \brief Projector to any surface
3910    */
3911   //================================================================================
3912
3913   struct SurfaceProjector : public FT_RealProjector
3914   {
3915     ShapeAnalysis_Surface    _projector;
3916     double                   _tol;
3917     BRepTopAdaptor_FClass2d* _classifier;
3918
3919     //-----------------------------------------------------------------------------
3920     SurfaceProjector( const TopoDS_Face& face, const double tol, BRepTopAdaptor_FClass2d* cls ):
3921       _projector( BRep_Tool::Surface( face )),
3922       _tol( tol ),
3923       _classifier( cls )
3924     {
3925     }
3926     //-----------------------------------------------------------------------------
3927     // delete _classifier
3928     ~SurfaceProjector()
3929     {
3930       delete _classifier;
3931     }
3932
3933     //-----------------------------------------------------------------------------
3934     // project a point to a surface
3935     virtual gp_Pnt project( const gp_Pnt& P,
3936                             double*       newSolution,
3937                             const double* prevSolution = 0)
3938     {
3939       gp_Pnt2d uv;
3940
3941       if ( prevSolution )
3942       {
3943         gp_Pnt2d prevUV( prevSolution[0], prevSolution[1] );
3944         uv = _projector.NextValueOfUV( prevUV, P, _tol );
3945       }
3946       else
3947       {
3948         uv = _projector.ValueOfUV( P, _tol );
3949       }
3950
3951       uv.Coord( newSolution[0], newSolution[1] );
3952
3953       gp_Pnt proj = _projector.Value( uv );
3954
3955       _dist = _projector.Gap();
3956
3957       return proj;
3958     }
3959
3960     //-----------------------------------------------------------------------------
3961     // project a point to a surface and check if the projection is within the surface boundary
3962     virtual bool projectAndClassify( const gp_Pnt& point,
3963                                      const double  maxDist2,
3964                                      gp_Pnt&       projection,
3965                                      double*       newSolution,
3966                                      const double* prevSolution = 0)
3967     {
3968       projection = project( point, newSolution, prevSolution );
3969       return ( _dist * _dist < maxDist2 )  &&  classify( newSolution );
3970     }
3971
3972     //-----------------------------------------------------------------------------
3973     // check if the projection is within the shape boundary
3974     bool classify( const double* newSolution )
3975     {
3976       TopAbs_State state = _classifier->Perform( gp_Pnt2d( newSolution[0], newSolution[1]) );
3977       return ( state != TopAbs_OUT );
3978     }
3979
3980     //-----------------------------------------------------------------------------
3981     // return true if a previously found solution can be used to speed up the projection
3982     virtual bool canUsePrevSolution() const { return true; }
3983   };
3984
3985   //================================================================================
3986   /*!
3987    * \brief Projector to a plane. Don't project, classify only
3988    */
3989   //================================================================================
3990
3991   struct PlaneProjector : public SurfaceProjector
3992   {
3993     gp_Pln _plane;
3994     bool   _isRealPlane; // false means that a surface is planar but parametrization is different
3995
3996     //-----------------------------------------------------------------------------
3997     PlaneProjector( const gp_Pln&            pln,
3998                     const TopoDS_Face&       face,
3999                     BRepTopAdaptor_FClass2d* cls,
4000                     bool                     isRealPlane=true):
4001       SurfaceProjector( face, 0, cls ),
4002       _plane( pln ),
4003       _isRealPlane( isRealPlane )
4004     {}
4005
4006     //-----------------------------------------------------------------------------
4007     // does nothing
4008     virtual gp_Pnt project( const gp_Pnt& P,
4009                             double*       newSolution,
4010                             const double* prevSolution = 0)
4011     {
4012       return P;
4013     }
4014     //-----------------------------------------------------------------------------
4015     // check if a point lies within the boundry of the planar face
4016     virtual bool projectAndClassify( const gp_Pnt& point,
4017                                      const double  maxDist2,
4018                                      gp_Pnt&       projection,
4019                                      double*       newSolution,
4020                                      const double* prevSolution = 0)
4021     {
4022       if ( _isRealPlane )
4023       {
4024         ElSLib::PlaneParameters( _plane.Position(), point, newSolution[0], newSolution[1]);
4025         projection = ElSLib::PlaneValue ( newSolution[0], newSolution[1], _plane.Position() );
4026         if ( projection.SquareDistance( point ) > theEPS * theEPS )
4027           return false;
4028
4029         return SurfaceProjector::classify( newSolution );
4030       }
4031       else
4032       {
4033         return SurfaceProjector::projectAndClassify( point, maxDist2, projection,
4034                                                      newSolution, prevSolution );
4035       }
4036     }
4037     //-----------------------------------------------------------------------------
4038     // return true if a previously found solution can be used to speed up the projection
4039     virtual bool canUsePrevSolution() const { return false; }
4040   };
4041
4042   //================================================================================
4043   /*!
4044    * \brief Projector to a cylinder
4045    */
4046   //================================================================================
4047
4048   struct CylinderProjector : public SurfaceProjector
4049   {
4050     gp_Cylinder _cylinder;
4051
4052     //-----------------------------------------------------------------------------
4053     CylinderProjector( const gp_Cylinder&       c,
4054                        const TopoDS_Face&       face,
4055                        BRepTopAdaptor_FClass2d* cls ):
4056       SurfaceProjector( face, 0, cls ),
4057       _cylinder( c )
4058     {}
4059
4060     //-----------------------------------------------------------------------------
4061     // project a point to the cylinder
4062     virtual gp_Pnt project( const gp_Pnt& P,
4063                             double*       newSolution,
4064                             const double* prevSolution = 0)
4065     {
4066       // project the point P to the cylinder axis -> Pp
4067       const gp_Pnt& O   = _cylinder.Position().Location();
4068       const gp_Dir& axe = _cylinder.Position().Direction();
4069       gp_Vec       trsl = gp_Vec( axe ).Multiplied( gp_Vec( O, P ).Dot( axe ));
4070       gp_Pnt       Pp   = O.Translated( trsl );
4071
4072       // move Pp to the cylinder
4073       gp_Vec radiusVec( Pp, P );
4074       double radius = radiusVec.Magnitude();
4075       if ( radius < std::numeric_limits<double>::min() )
4076         return P; // P in on the axe
4077
4078       gp_Pnt proj = Pp.Translated( radiusVec.Multiplied( _cylinder.Radius() / radius ));
4079
4080       _dist = _cylinder.Radius() - radius;
4081
4082       return proj;
4083     }
4084     //-----------------------------------------------------------------------------
4085     // project a point to the cylinder and check if the projection is within the surface boundary
4086     virtual bool projectAndClassify( const gp_Pnt& point,
4087                                      const double  maxDist2,
4088                                      gp_Pnt&       projection,
4089                                      double*       newSolution,
4090                                      const double* prevSolution = 0)
4091     {
4092       ElSLib::CylinderParameters( _cylinder.Position(), _cylinder.Radius(), point,
4093                                   newSolution[0], newSolution[1]);
4094       projection = ElSLib::CylinderValue( newSolution[0], newSolution[1],
4095                                           _cylinder.Position(), _cylinder.Radius() );
4096
4097       return ( _dist * _dist < maxDist2 )  &&  SurfaceProjector::classify( newSolution );
4098     }
4099     //-----------------------------------------------------------------------------
4100     // return true if a previously found solution can be used to speed up the projection
4101     virtual bool canUsePrevSolution() const { return false; }
4102   };
4103
4104   //================================================================================
4105   /*!
4106    * \brief Projector to a cone
4107    */
4108   //================================================================================
4109
4110   struct ConeProjector : public SurfaceProjector
4111   {
4112     gp_Cone _cone;
4113
4114     //-----------------------------------------------------------------------------
4115     ConeProjector( const gp_Cone&           c,
4116                    const TopoDS_Face&       face,
4117                    BRepTopAdaptor_FClass2d* cls ):
4118       SurfaceProjector( face, 0, cls ),
4119       _cone( c )
4120     {}
4121
4122     //-----------------------------------------------------------------------------
4123     // project a point to the cone
4124     virtual gp_Pnt project( const gp_Pnt& point,
4125                             double*       newSolution,
4126                             const double* prevSolution = 0)
4127     {
4128       ElSLib::ConeParameters( _cone.Position(), _cone.RefRadius(), _cone.SemiAngle(),
4129                               point, newSolution[0], newSolution[1]);
4130       gp_Pnt proj = ElSLib::ConeValue( newSolution[0], newSolution[1],
4131                                        _cone.Position(), _cone.RefRadius(), _cone.SemiAngle() );
4132       _dist = point.Distance( proj );
4133
4134       return proj;
4135     }
4136
4137     //-----------------------------------------------------------------------------
4138     // project a point to the cone and check if the projection is within the surface boundary
4139     virtual bool projectAndClassify( const gp_Pnt& point,
4140                                      const double  maxDist2,
4141                                      gp_Pnt&       projection,
4142                                      double*       newSolution,
4143                                      const double* prevSolution = 0)
4144     {
4145       projection = project( point, newSolution, prevSolution );
4146
4147       return ( _dist * _dist < maxDist2 )  &&  SurfaceProjector::classify( newSolution );
4148     }
4149     //-----------------------------------------------------------------------------
4150     // return true if a previously found solution can be used to speed up the projection
4151     virtual bool canUsePrevSolution() const { return false; }
4152   };
4153
4154   //================================================================================
4155   /*!
4156    * \brief Projector to a sphere
4157    */
4158   //================================================================================
4159
4160   struct SphereProjector : public SurfaceProjector
4161   {
4162     gp_Sphere _sphere;
4163
4164     //-----------------------------------------------------------------------------
4165     SphereProjector( const gp_Sphere&         s,
4166                      const TopoDS_Face&       face,
4167                      BRepTopAdaptor_FClass2d* cls ):
4168       SurfaceProjector( face, 0, cls ),
4169       _sphere( s )
4170     {}
4171
4172     //-----------------------------------------------------------------------------
4173     // project a point to the sphere
4174     virtual gp_Pnt project( const gp_Pnt& P,
4175                             double*       newSolution,
4176                             const double* prevSolution = 0)
4177     {
4178       // move Pp to the Sphere
4179       const gp_Pnt& O = _sphere.Location();
4180       gp_Vec radiusVec( O, P );
4181       double radius = radiusVec.Magnitude();
4182       if ( radius < std::numeric_limits<double>::min() )
4183         return P; // P is on O
4184
4185       gp_Pnt proj = O.Translated( radiusVec.Multiplied( _sphere.Radius() / radius ));
4186
4187       _dist = _sphere.Radius() - radius;
4188
4189       return proj;
4190     }
4191
4192     //-----------------------------------------------------------------------------
4193     // project a point to the sphere and check if the projection is within the surface boundary
4194     virtual bool projectAndClassify( const gp_Pnt& point,
4195                                      const double  maxDist2,
4196                                      gp_Pnt&       projection,
4197                                      double*       newSolution,
4198                                      const double* prevSolution = 0)
4199     {
4200       ElSLib::SphereParameters( _sphere.Position(), _sphere.Radius(), point,
4201                                   newSolution[0], newSolution[1]);
4202       projection = ElSLib::SphereValue( newSolution[0], newSolution[1],
4203                                         _sphere.Position(), _sphere.Radius() );
4204
4205       return ( _dist * _dist < maxDist2 )  &&  SurfaceProjector::classify( newSolution );
4206     }
4207     //-----------------------------------------------------------------------------
4208     // return true if a previously found solution can be used to speed up the projection
4209     virtual bool canUsePrevSolution() const { return false; }
4210   };
4211
4212   //================================================================================
4213   /*!
4214    * \brief Projector to a torus
4215    */
4216   //================================================================================
4217
4218   struct TorusProjector : public SurfaceProjector
4219   {
4220     gp_Torus _torus;
4221
4222     //-----------------------------------------------------------------------------
4223     TorusProjector( const gp_Torus&          t,
4224                     const TopoDS_Face&       face,
4225                     BRepTopAdaptor_FClass2d* cls ):
4226       SurfaceProjector( face, 0, cls ),
4227       _torus( t )
4228     {}
4229
4230     //-----------------------------------------------------------------------------
4231     // project a point to the torus
4232     virtual gp_Pnt project( const gp_Pnt& point,
4233                             double*       newSolution,
4234                             const double* prevSolution = 0)
4235     {
4236       ElSLib::TorusParameters( _torus.Position(), _torus.MajorRadius(), _torus.MinorRadius(),
4237                                point, newSolution[0], newSolution[1]);
4238       gp_Pnt proj = ElSLib::TorusValue( newSolution[0], newSolution[1],
4239                                         _torus.Position(), _torus.MajorRadius(), _torus.MinorRadius() );
4240       _dist = point.Distance( proj );
4241
4242       return proj;
4243     }
4244
4245     //-----------------------------------------------------------------------------
4246     // project a point to the torus and check if the projection is within the surface boundary
4247     virtual bool projectAndClassify( const gp_Pnt& point,
4248                                      const double  maxDist2,
4249                                      gp_Pnt&       projection,
4250                                      double*       newSolution,
4251                                      const double* prevSolution = 0)
4252     {
4253       projection = project( point, newSolution, prevSolution );
4254
4255       return ( _dist * _dist < maxDist2 )  &&  SurfaceProjector::classify( newSolution );
4256     }
4257     //-----------------------------------------------------------------------------
4258     // return true if a previously found solution can be used to speed up the projection
4259     virtual bool canUsePrevSolution() const { return false; }
4260   };
4261
4262   //================================================================================
4263   /*!
4264    * \brief Check if a curve can be considered straight
4265    */
4266   //================================================================================
4267
4268   bool isStraight( const GeomAdaptor_Curve& curve, const double tol )
4269   {
4270     // rough check: evaluate how far from a straight line connecting the curve ends
4271     // stand several internal points of the curve
4272
4273     const double  f = curve.FirstParameter();
4274     const double  l = curve.LastParameter();
4275     const gp_Pnt pf = curve.Value( f );
4276     const gp_Pnt pl = curve.Value( l );
4277     const gp_Vec lineVec( pf, pl );
4278     const double lineLen2 = lineVec.SquareMagnitude();
4279     if ( lineLen2 < std::numeric_limits< double >::min() )
4280       return false; // E seems closed
4281
4282     const double nbSamples = 7;
4283     for ( int i = 0; i < nbSamples; ++i )
4284     {
4285       const double  r = ( i + 1 ) / nbSamples;
4286       const gp_Pnt pi = curve.Value( f * r + l * ( 1 - r ));
4287       const gp_Vec vi( pf, pi );
4288       const double h2 = lineVec.Crossed( vi ).SquareMagnitude() / lineLen2;
4289       if ( h2 > tol * tol )
4290         return false;
4291     }
4292
4293     // thorough check
4294     GCPnts_UniformDeflection divider( curve, tol );
4295     return ( divider.IsDone() && divider.NbPoints() < 3 );
4296   }
4297 }
4298
4299 //================================================================================
4300 /*!
4301  * \brief Initialize with a boundary shape
4302  */
4303 //================================================================================
4304
4305 FT_Projector::FT_Projector(const TopoDS_Shape& shape)
4306 {
4307   _realProjector = 0;
4308   setBoundaryShape( shape );
4309   _tryWOPrevSolution = false;
4310 }
4311
4312 //================================================================================
4313 /*!
4314  * \brief Copy another projector
4315  */
4316 //================================================================================
4317
4318 FT_Projector::FT_Projector(const FT_Projector& other)
4319 {
4320   _realProjector = 0;
4321   _shape = other._shape;
4322   _bndBox = other._bndBox;
4323   _tryWOPrevSolution = false;
4324 }
4325
4326 //================================================================================
4327 /*!
4328  * \brief Destructor. Delete _realProjector
4329  */
4330 //================================================================================
4331
4332 FT_Projector::~FT_Projector()
4333 {
4334   delete _realProjector;
4335 }
4336
4337 //================================================================================
4338 /*!
4339  * \brief Initialize with a boundary shape. Compute the bounding box
4340  */
4341 //================================================================================
4342
4343 void FT_Projector::setBoundaryShape(const TopoDS_Shape& shape)
4344 {
4345   delete _realProjector; _realProjector = 0;
4346   _shape = shape;
4347   if ( shape.IsNull() )
4348     return;
4349
4350   BRepBndLib::Add( shape, _bndBox );
4351   _bndBox.Enlarge( 1e-5 * sqrt( _bndBox.SquareExtent() ));
4352 }
4353
4354 //================================================================================
4355 /*!
4356  * \brief Create a real projector
4357  */
4358 //================================================================================
4359
4360 void FT_Projector::prepareForProjection()
4361 {
4362   if ( _shape.IsNull() || _realProjector )
4363     return;
4364
4365   if ( _shape.ShapeType() == TopAbs_EDGE )
4366   {
4367     const TopoDS_Edge& edge = TopoDS::Edge( _shape );
4368
4369     double tol = 1e-6 * sqrt( _bndBox.SquareExtent() );
4370
4371     double f,l;
4372     Handle(Geom_Curve) curve = BRep_Tool::Curve( edge, f,l );
4373     if ( curve.IsNull() )
4374       return; // degenerated edge
4375
4376     GeomAdaptor_Curve acurve( curve, f, l );
4377     switch ( acurve.GetType() )
4378     {
4379     case GeomAbs_Line:
4380       _realProjector = new LineProjector( edge );
4381       break;
4382     case GeomAbs_Circle:
4383       _realProjector = new CircleProjector( acurve.Circle(), f, l );
4384       break;
4385     case GeomAbs_BezierCurve:
4386     case GeomAbs_BSplineCurve:
4387     case GeomAbs_OffsetCurve:
4388     case GeomAbs_OtherCurve:
4389       if ( isStraight( acurve, tol ))
4390       {
4391         _realProjector = new LineProjector( edge );
4392         break;
4393       }
4394     case GeomAbs_Ellipse:
4395     case GeomAbs_Hyperbola:
4396     case GeomAbs_Parabola:
4397       _realProjector = new CurveProjector( edge, tol );
4398     }
4399   }
4400   else if ( _shape.ShapeType() == TopAbs_FACE )
4401   {
4402     TopoDS_Face face = TopoDS::Face( _shape );
4403
4404     Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
4405     if ( surface.IsNull() )
4406       return;
4407
4408     GeomAdaptor_Surface asurface( surface );
4409     Standard_Real tol   = BRep_Tool::Tolerance( face );
4410     Standard_Real toluv = Min( asurface.UResolution( tol ), asurface.VResolution( tol ));
4411     BRepTopAdaptor_FClass2d* classifier = new BRepTopAdaptor_FClass2d( face, toluv );
4412
4413     switch ( asurface.GetType() )
4414     {
4415     case GeomAbs_Plane:
4416       _realProjector = new PlaneProjector( asurface.Plane(), face, classifier );
4417       break;
4418     case GeomAbs_Cylinder:
4419       _realProjector = new CylinderProjector( asurface.Cylinder(), face, classifier );
4420       break;
4421     case GeomAbs_Sphere:
4422       _realProjector = new SphereProjector( asurface.Sphere(), face, classifier );
4423       break;
4424     case GeomAbs_Cone:
4425       _realProjector = new ConeProjector( asurface.Cone(), face, classifier );
4426       break;
4427     case GeomAbs_Torus:
4428       _realProjector = new TorusProjector( asurface.Torus(), face, classifier );
4429       break;
4430     case GeomAbs_BezierSurface:
4431     case GeomAbs_BSplineSurface:
4432     case GeomAbs_SurfaceOfRevolution:
4433     case GeomAbs_SurfaceOfExtrusion:
4434     case GeomAbs_OffsetSurface:
4435     case GeomAbs_OtherSurface:
4436       GeomLib_IsPlanarSurface isPlaneCheck( surface, tol );
4437       if ( isPlaneCheck.IsPlanar() )
4438       {
4439         _realProjector = new PlaneProjector( isPlaneCheck.Plan(), face, classifier,
4440                                              /*isRealPlane=*/false);
4441       }
4442       else
4443       {
4444         _realProjector = new SurfaceProjector( face, tol, classifier );
4445       }
4446       break;
4447     }
4448
4449     if ( !_realProjector )
4450       delete classifier;
4451   }
4452 }
4453
4454 //================================================================================
4455 /*!
4456  * \brief Return true if projection is not needed
4457  */
4458 //================================================================================
4459
4460 bool FT_Projector::isPlanarBoundary() const
4461 {
4462   return ( dynamic_cast< LineProjector*  >( _realProjector ) ||
4463            dynamic_cast< PlaneProjector* >( _realProjector ) );
4464 }
4465
4466 //================================================================================
4467 /*!
4468  * \brief Check if a point lies on the boundary shape
4469  *  \param [in] point - the point to check
4470  *  \param [in] tol2 - a square tolerance allowing to decide whether a point is on the shape
4471  *  \param [in] newSolution - position on the shape (U or UV) of the point found
4472  *         during projecting
4473  *  \param [in] prevSolution - position on the shape (U or UV) of a neighbor point
4474  *  \return bool - \c true if the point lies on the boundary shape
4475  *
4476  * This method is used to select a shape by checking if all neighbor nodes of a node to move
4477  * lie on a shape.
4478  */
4479 //================================================================================
4480
4481 bool FT_Projector::isOnShape( const gp_Pnt& point,
4482                               const double  tol2,
4483                               double*       newSolution,
4484                               const double* prevSolution)
4485 {
4486   if ( _bndBox.IsOut( point ) || !_realProjector )
4487     return false;
4488
4489   gp_Pnt proj;
4490   if ( isPlanarBoundary() )
4491     return projectAndClassify( point, tol2, proj, newSolution, prevSolution );
4492
4493   return project( point, tol2, proj, newSolution, prevSolution );
4494 }
4495
4496 //================================================================================
4497 /*!
4498  * \brief Project a point to the boundary shape
4499  *  \param [in] point - the point to project
4500  *  \param [in] maxDist2 - the maximal square distance between the point and the projection
4501  *  \param [out] projection - the projection
4502  *  \param [out] newSolution - position on the shape (U or UV) of the point found
4503  *         during projecting
4504  *  \param [in] prevSolution - already found position on the shape (U or UV) of a neighbor point
4505  *  \return bool - false if the distance between the point and the projection
4506  *         is more than sqrt(maxDist2)
4507  *
4508  * This method is used to project a node in the case where only one shape is found by name
4509  */
4510 //================================================================================
4511
4512 bool FT_Projector::project( const gp_Pnt& point,
4513                             const double  maxDist2,
4514                             gp_Pnt&       projection,
4515                             double*       newSolution,
4516                             const double* prevSolution)
4517 {
4518   if ( !_realProjector )
4519     return false;
4520
4521   _realProjector->_dist = 1e100;
4522   projection = _realProjector->project( point, newSolution, prevSolution );
4523
4524   bool ok = ( _realProjector->_dist * _realProjector->_dist < maxDist2 );
4525   if ( !ok && _tryWOPrevSolution && prevSolution )
4526   {
4527     projection = _realProjector->project( point, newSolution );
4528     ok = ( _realProjector->_dist * _realProjector->_dist < maxDist2 );
4529   }
4530   return ok;
4531 }
4532
4533 //================================================================================
4534 /*!
4535  * \brief Project a point to the boundary shape and check if the projection lies within
4536  *        the shape boundary
4537  *  \param [in] point - the point to project
4538  *  \param [in] maxDist2 - the maximal square distance between the point and the projection
4539  *  \param [out] projection - the projection
4540  *  \param [out] newSolution - position on the shape (U or UV) of the point found
4541  *         during projecting
4542  *  \param [in] prevSolution - already found position on the shape (U or UV) of a neighbor point
4543  *  \return bool - false if the projection point lies out of the shape boundary or
4544  *          the distance between the point and the projection is more than sqrt(maxDist2)
4545  *
4546  * This method is used to project a node in the case where several shapes are selected for
4547  * projection of a node group
4548  */
4549 //================================================================================
4550
4551 bool FT_Projector::projectAndClassify( const gp_Pnt& point,
4552                                        const double  maxDist2,
4553                                        gp_Pnt&       projection,
4554                                        double*       newSolution,
4555                                        const double* prevSolution)
4556 {
4557   if ( _bndBox.IsOut( point ) || !_realProjector )
4558     return false;
4559
4560   bool ok = _realProjector->projectAndClassify( point, maxDist2, projection,
4561                                                 newSolution, prevSolution );
4562   if ( !ok && _tryWOPrevSolution && prevSolution )
4563     ok = _realProjector->projectAndClassify( point, maxDist2, projection, newSolution );
4564
4565   return ok;
4566 }
4567
4568 //================================================================================
4569 /*!
4570  * \brief Return true if a previously found solution can be used to speed up the projection
4571  */
4572 //================================================================================
4573
4574 bool FT_Projector::canUsePrevSolution() const
4575 {
4576   return ( _realProjector && _realProjector->canUsePrevSolution() );
4577 }
4578
4579 //================================================================================
4580 /*
4581  * \brief Check if a file exists
4582  */
4583 //================================================================================
4584
4585 bool FT_Utils::fileExists( const std::string& path )
4586 {
4587   if ( path.empty() )
4588     return false;
4589
4590   boost::system::error_code err;
4591   bool res = boofs::exists( path, err );
4592
4593   return err ? false : res;
4594 }
4595
4596 //================================================================================
4597 /*!
4598  * \brief Check if a file can be created/overwritten
4599  */
4600 //================================================================================
4601
4602 bool FT_Utils::canWrite( const std::string& path )
4603 {
4604   if ( path.empty() )
4605     return false;
4606
4607   bool can = false;
4608 #ifdef WIN32
4609
4610   HANDLE file = CreateFile( path.c_str(),           // name of the write
4611                             GENERIC_WRITE,          // open for writing
4612                             0,                      // do not share
4613                             NULL,                   // default security
4614                             OPEN_ALWAYS,            // CREATE NEW or OPEN EXISTING
4615                             FILE_ATTRIBUTE_NORMAL,  // normal file
4616                             NULL);                  // no attr. template
4617   can = ( file != INVALID_HANDLE_VALUE );
4618   CloseHandle( file );
4619
4620 #else
4621
4622   int file = ::open( path.c_str(),
4623                      O_WRONLY | O_CREAT,
4624                      S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH ); // rw-r--r--
4625   can = ( file >= 0 );
4626
4627 #endif
4628
4629   return can;
4630 }
4631
4632 //================================================================================
4633 /*!
4634  * \brief Make a map of XAO groups
4635  */
4636 //================================================================================
4637
4638 FT_Utils::XaoGroups::XaoGroups( const XAO::Xao* theXao )
4639 {
4640   XAO::Xao* xao = const_cast< XAO::Xao* >( theXao );
4641
4642   for ( int iG = 0; iG < theXao->countGroups(); ++iG )
4643   {
4644     XAO::Group* group = xao->getGroup( iG );
4645
4646     if ( group->getDimension() == 1 )
4647
4648       _xaoGroups[ 0 ].insert( std::make_pair( group->getName(), group ));
4649
4650     else if ( group->getDimension() == 2 )
4651
4652       _xaoGroups[ 1 ].insert( std::make_pair( group->getName(), group ));
4653   }
4654 }
4655
4656 //================================================================================
4657 /*!
4658  * \brief Return FT_Projector's by a group name
4659  *  \param [in] groupName - the group name
4660  *  \param [in] dim - the group dimension
4661  *  \param [in] allProjectors - the projector of all shapes of \a dim dimension
4662  *  \param [out] groupProjectors - projectors to shapes of the group
4663  *  \return int - number of found shapes
4664  */
4665 //================================================================================
4666
4667 int FT_Utils::XaoGroups::getProjectors( const std::string&                   groupName,
4668                                         const int                            dim,
4669                                         const std::vector< FT_Projector > &  allProjectors,
4670                                         std::vector< const FT_Projector* > & groupProjectors) const
4671 {
4672   // get namesake groups
4673
4674   const TGroupByNameMap* groupMap = 0;
4675   if ( dim == 1 )
4676     groupMap = &_xaoGroups[ 0 ];
4677   else if ( dim == 2 )
4678     groupMap = &_xaoGroups[ 1 ];
4679   else
4680     return 0;
4681
4682   TGroupByNameMap::const_iterator name2gr = groupMap->find( groupName );
4683   if ( name2gr == groupMap->end() )
4684     return 0;
4685
4686   std::vector< XAO::Group* > groups;
4687   groups.push_back( name2gr->second );
4688
4689   for ( ++name2gr; name2gr != groupMap->end(); ++name2gr )
4690   {
4691     if ( name2gr->second->getName() == groupName )
4692       groups.push_back( name2gr->second );
4693     else
4694       break;
4695   }
4696
4697   // get projectors
4698
4699   int nbFound = 0;
4700   for ( size_t i = 0; i < groups.size(); ++i )
4701   {
4702     // IDs in XAO correspond to indices of allProjectors
4703     std::set<int>::iterator id = groups[i]->begin(), end = groups[i]->end();
4704     for ( ; id != end; ++id, ++nbFound )
4705       if ( *id < (int) allProjectors.size() )
4706         groupProjectors.push_back ( & allProjectors[ *id ]);
4707   }
4708
4709   return nbFound;
4710 }
4711
4712 } // namespace SMESHHOMARDImpl /end/