Salome HOME
Merge from V6_main_20120808 08Aug12
[plugins/ghs3dprlplugin.git] / src / GHS3DPRLPlugin / GHS3DPRLPlugin_GHS3DPRL.cxx
index 0278eacc40d2762c10828897594a3ea4a0e7de33..cb93731912bed4d270d9656eeb64d97f15efede5 100755 (executable)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2008 OPEN CASCADE, CEA/DEN, EDF R&D
+// Copyright (C) 2007-2012  CEA/DEN, EDF R&D
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 //
 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
+
 // ---
-//
 // File   : GHS3DPRLPlugin_GHS3DPRL.cxx
 // Author : Christian VAN WAMBEKE (CEA) (from Hexotic plugin Lioka RAZAFINDRAZAKA)
-//
 // ---
-
+//
 #include "GHS3DPRLPlugin_GHS3DPRL.hxx"
 #include "GHS3DPRLPlugin_Hypothesis.hxx"
 
 #include <SMESH_ControlsDef.hxx>
 
 #include <list>
+#include <Standard_ProgramError.hxx>
 #include <TCollection_AsciiString.hxx>
+#include <TopTools_MapOfShape.hxx>
+#include <TopoDS.hxx>
+#include <BRepGProp.hxx>
+#include <GProp_GProps.hxx>
 
 /* 
 extern "C" {
@@ -62,6 +66,16 @@ extern "C" {
 }
 using namespace med_2_2;*/
 
+static void removeFile( const TCollection_AsciiString& fileName )
+{
+  try {
+    OSD_File( fileName ).Remove();
+  }
+  catch ( Standard_ProgramError ) {
+    MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied");
+  }
+}
+
 //=============================================================================
 GHS3DPRLPlugin_GHS3DPRL::GHS3DPRLPlugin_GHS3DPRL(int hypId, int studyId, SMESH_Gen* gen)
   : SMESH_3D_Algo(hypId, studyId, gen)
@@ -126,11 +140,11 @@ static bool writeGHS3DPRLFiles (const TCollection_AsciiString &  GHS3DPRL_In,
    int ifam=0;
    TCollection_AsciiString namefile(GHS3DPRL_In);
    namefile+=".points";
-   OSD_File(namefile).Remove();
+   removeFile(namefile);
    ofstream theFile;
    theFile.open(namefile.ToCString(),ios::out);
 #ifdef WIN32
-   Ok=theFile->is_open();
+   Ok = theFile.is_open();
 #else
    Ok=theFile.rdbuf()->is_open();
 #endif
@@ -169,10 +183,10 @@ static bool writeGHS3DPRLFiles (const TCollection_AsciiString &  GHS3DPRL_In,
    theFile.close();
 
    namefile=GHS3DPRL_In+".faces";
-   OSD_File(namefile).Remove();
+   removeFile(namefile);
    theFile.open(namefile.ToCString(),ios::out);
 #ifdef WIN32
-   Ok=theFile->is_open();
+   Ok=theFile.is_open();
 #else
    Ok=theFile.rdbuf()->is_open();
 #endif
@@ -281,10 +295,10 @@ static bool readResult(FILE *                           theFile,
   int shapeID = theMesh->ShapeToIndex( theShape );
 
   int line = 1, EndOfFile = 0, nbElem = 0, nField = 10, nbRef = 0, aGHS3DPRLNodeID = 0;
-  char * theField;
+  const char * theField;
 
-  char * tabField [nField];
-  int    tabRef [nField];
+  vector<const char*> tabField = vector<const char*>(nField);
+  vector<int> tabRef = vector<int>(nField);
 
   tabField[0] = "MeshVersionFormatted";    tabRef[0] = 0;
   tabField[1] = "Dimension";               tabRef[1] = 0;
@@ -331,7 +345,7 @@ static bool readResult(FILE *                           theFile,
 
                       if ( strcmp(theField, "Vertices") == 0 ) {
                           int aGHS3DPRLID;
-                          double coord[nbRef];
+                          vector<double> coord = vector<double>(nbRef);
                           SMDS_MeshNode * aGHS3DPRLNode;
 
                           for ( int iElem = 0; iElem < nbElem; iElem++ ) {
@@ -345,7 +359,7 @@ static bool readResult(FILE *                           theFile,
                           }
                       }
                       else {
-                          const SMDS_MeshNode * node[nbRef];
+                          vector<const SMDS_MeshNode*> node = vector<const SMDS_MeshNode*>(nbRef);
                           SMDS_MeshElement* aGHS3DPRLElement;
                           map <int,const SMDS_MeshNode*>::iterator itOnGHS3DPRLNode;
 
@@ -506,7 +520,7 @@ bool GHS3DPRLPlugin_GHS3DPRL::Compute(SMESH_Mesh& theMesh,
       cout<<"GHS3DPRL command : "<<run_GHS3DPRL.ToCString()<<endl;
       cout<<"Write files .faces .point ...";
       GHS3DPRL_Out = path + casenamemed;
-      OSD_File( GHS3DPRL_Outxml ).Remove(); //only the master xml file
+      removeFile( GHS3DPRL_Outxml ); //only the master xml file
       Ok=writeGHS3DPRLFiles(GHS3DPRL_In, meshDS, aSmdsToGHS3DPRLIdMap, aGHS3DPRLIdToNodeMap);
       if (Ok) {cout<<" ...done\n";}
       else {
@@ -624,3 +638,84 @@ istream & operator >> (istream & load, GHS3DPRLPlugin_GHS3DPRL & hyp)
 {
   return hyp.LoadFrom( load );
 }
+
+
+//=============================================================================
+/*!
+ *  
+ */
+//=============================================================================
+bool GHS3DPRLPlugin_GHS3DPRL::Evaluate(SMESH_Mesh& aMesh,
+                                      const TopoDS_Shape& aShape,
+                                      MapShapeNbElems& aResMap)
+{
+  int nbtri = 0, nbqua = 0;
+  double fullArea = 0.0;
+  for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
+    TopoDS_Face F = TopoDS::Face( exp.Current() );
+    SMESH_subMesh *sm = aMesh.GetSubMesh(F);
+    MapShapeNbElemsItr anIt = aResMap.find(sm);
+    if( anIt==aResMap.end() ) {
+      SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+      smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
+                                           "Submesh can not be evaluated",this));
+      return false;
+    }
+    std::vector<int> aVec = (*anIt).second;
+    nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
+    nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
+    GProp_GProps G;
+    BRepGProp::SurfaceProperties(F,G);
+    double anArea = G.Mass();
+    fullArea += anArea;
+  }
+
+  // collect info from edges
+  int nb0d_e = 0, nb1d_e = 0;
+  bool IsQuadratic = false;
+  bool IsFirst = true;
+  TopTools_MapOfShape tmpMap;
+  for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
+    TopoDS_Edge E = TopoDS::Edge(exp.Current());
+    if( tmpMap.Contains(E) )
+      continue;
+    tmpMap.Add(E);
+    SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
+    MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
+    std::vector<int> aVec = (*anIt).second;
+    nb0d_e += aVec[SMDSEntity_Node];
+    nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
+    if(IsFirst) {
+      IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
+      IsFirst = false;
+    }
+  }
+  tmpMap.Clear();
+
+  double ELen = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
+
+  GProp_GProps G;
+  BRepGProp::VolumeProperties(aShape,G);
+  double aVolume = G.Mass();
+  double tetrVol = 0.1179*ELen*ELen*ELen;
+  double CoeffQuality = 0.9;
+  int nbVols = (int)aVolume/tetrVol/CoeffQuality;
+  int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
+  int nb1d_in = (int) ( nbVols*6 - nb1d_e - nb1d_f ) / 5;
+  std::vector<int> aVec(SMDSEntity_Last);
+  for(int i=0; i<SMDSEntity_Last; i++) aVec[i]=0;
+  if( IsQuadratic ) {
+    aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
+    aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
+    aVec[SMDSEntity_Quad_Pyramid] = nbqua;
+  }
+  else {
+    aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
+    aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
+    aVec[SMDSEntity_Pyramid] = nbqua;
+  }
+  SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
+  aResMap.insert(std::make_pair(sm,aVec));
+
+  return true;
+}