This commit was generated by cvs2git to create branch 'IMPORT'.

author admin <salome-admin@opencascade.com>

Mon, 19 May 2003 14:07:01 +0000 (14:07 +0000)

committer admin <salome-admin@opencascade.com>

Mon, 19 May 2003 14:07:01 +0000 (14:07 +0000)
author admin <salome-admin@opencascade.com>
Mon, 19 May 2003 14:07:01 +0000 (14:07 +0000)
committer admin <salome-admin@opencascade.com>
Mon, 19 May 2003 14:07:01 +0000 (14:07 +0000)
diff --git a/idl/SMESH_BasicHypothesis.idl b/idl/SMESH_BasicHypothesis.idl

new file mode 100644 (file)

index 0000000..694b791
--- /dev/null
+++ b/idl/SMESH_BasicHypothesis.idl
@@ -0,0 +1,64 @@
+//=============================================================================
+// File      : SMESH_BasicHypothesis.idl
+// Created   : mer mai 15 13:37:18 CEST 2002
+// Author    : Paul RASCLE, EDF
+// Project   : SALOME
+// Copyright : EDF 2002
+// $Header$
+//=============================================================================
+
+#ifndef _SMESH_BASICHYPOTHESIS_IDL_
+#define _SMESH_BASICHYPOTHESIS_IDL_
+
+#include "SALOME_Exception.idl"
+#include "SMESH_Hypothesis.idl"
+
+
+module SMESH
+{
+  interface SMESH_LocalLength : SMESH_Hypothesis
+  {
+    void SetLength(in double length) 
+      raises (SALOME::SALOME_Exception);
+    double GetLength();
+  };
+
+  interface SMESH_NumberOfSegments : SMESH_Hypothesis
+  {
+    void SetNumberOfSegments(in long segmentsNumber) 
+      raises (SALOME::SALOME_Exception);
+    long GetNumberOfSegments();
+  };
+
+  interface SMESH_MaxElementArea : SMESH_Hypothesis
+  {
+    void SetMaxElementArea(in double area) 
+      raises (SALOME::SALOME_Exception);
+    double GetMaxElementArea();
+  };
+
+  interface SMESH_MaxElementVolume : SMESH_Hypothesis
+  {
+    void SetMaxElementVolume(in double volume) 
+      raises (SALOME::SALOME_Exception);
+    double GetMaxElementVolume();
+  };
+
+  interface SMESH_Regular_1D : SMESH_1D_Algo
+  {
+  };
+
+  interface SMESH_MEFISTO_2D : SMESH_2D_Algo
+  {
+  };
+
+  interface SMESH_Quadrangle_2D : SMESH_2D_Algo
+  {
+  };
+
+  interface SMESH_Hexa_3D : SMESH_3D_Algo
+  {
+  };
+};
+
+#endif
diff --git a/idl/SMESH_Gen.idl b/idl/SMESH_Gen.idl

new file mode 100644 (file)

index 0000000..5617d14
--- /dev/null
+++ b/idl/SMESH_Gen.idl
@@ -0,0 +1,96 @@
+//=============================================================================
+// File      : SMESH_Gen.idl
+// Created   : jeu avr 11 15:26:35 CEST 2002
+// Author    : Paul RASCLE, EDF
+// Project   : SALOME
+// Copyright : EDF 2002
+// $Header$
+//=============================================================================
+
+#ifndef _SMESH_GEN_IDL_
+#define _SMESH_GEN_IDL_
+
+#include "SALOME_Exception.idl"
+#include "SALOME_Component.idl"
+#include "SALOMEDS.idl"
+
+#include "GEOM_Gen.idl"
+#include "GEOM_Shape.idl"
+
+#include "SMESH_Mesh.idl"
+#include "SMESH_Hypothesis.idl"
+
+module SMESH
+{
+  typedef sequence<GEOM::GEOM_Shape> shape_array;
+
+  interface SMESH_Gen : Engines::Component, SALOMEDS::Driver
+  {
+    /*!
+     * Create an hypothesis that can be shared by differents parts of the mesh.
+     * An hypothesis is either:
+     * - a method used to generate or modify a part of the mesh (algorithm).
+     * - a parameter or a law used by an algorithm.
+     * Algorithms are 1D, 2D or 3D.
+     */
+    SMESH_Hypothesis CreateHypothesis( in string anHyp,
+                                      in long studyId)
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Create a Mesh object, given a geometry shape.
+     * Mesh is created empty (no points, no elements).
+     * Shape is explored via GEOM_Client to create local copies.
+     * of TopoDS_Shapes and bind CORBA references of shape & subshapes
+     * with TopoDS_Shapes
+     */
+     SMESH_Mesh Init(in GEOM::GEOM_Gen geomEngine,
+                    in long studyId,
+                    in GEOM::GEOM_Shape aShape)
+       raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Create a Mesh object, without a geometry shape reference
+     */
+//      SMESH_Mesh NewEmpty(in GEOM::GEOM_Gen geomEngine,
+//                      in long studyId)
+//        raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Mesh a subShape. 
+     * First, verify list of hypothesis associated with the subShape,
+     * return NOK if hypothesis are not sufficient
+     */
+    boolean Compute(in SMESH_Mesh aMesh, in GEOM::GEOM_Shape aSubShape)
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * 
+     */
+
+    boolean IsReadyToCompute(in SMESH_Mesh aMesh, in GEOM::GEOM_Shape aSubShape)
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * 
+     */
+    long_array GetSubShapesId(in GEOM::GEOM_Gen geomEngine,
+                            in long studyId,
+                            in GEOM::GEOM_Shape mainShape,
+                            in shape_array listOfSubShape)
+       raises (SALOME::SALOME_Exception);
+    
+    /*!
+     * 
+     */
+    //    long_array GetSubMeshesState(in GEOM::GEOM_Gen geomEngine,
+    //                          in long studyId,
+    //                          in shape_array listOfSubShape)
+    //  raises (SALOME::SALOME_Exception);
+    
+
+  };
+
+};
+
+#endif
diff --git a/idl/SMESH_Hypothesis.idl b/idl/SMESH_Hypothesis.idl

new file mode 100644 (file)

index 0000000..ec61223
--- /dev/null
+++ b/idl/SMESH_Hypothesis.idl
@@ -0,0 +1,71 @@
+//=============================================================================
+// File      : SMESH_Hypothesis.idl
+// Created   : jeu avr 11 19:26:16 CEST 2002
+// Author    : Paul RASCLE, EDF
+// Project   : SALOME
+// Copyright : EDF 2002
+// $Header$
+//=============================================================================
+
+#ifndef _SMESH_HYPOTHESIS_IDL_
+#define _SMESH_HYPOTHESIS_IDL_
+
+#include "SALOME_Exception.idl"
+
+module SMESH
+{
+  interface SMESH_Hypothesis;
+
+  typedef sequence<SMESH_Hypothesis> ListOfHypothesis;
+  typedef sequence<string> ListOfHypothesisName;
+
+  interface SMESH_Hypothesis
+  {
+    /*!
+     * Get the Hypothesis typeName 
+     */
+    string GetName();
+
+    /*!
+     * Get the internal Id 
+     */
+    long GetId();
+ };
+
+  interface SMESH_Algo : SMESH_Hypothesis
+  {
+    /*!
+     * Get list of hypothesis that can be used with this algorithm
+     */
+    ListOfHypothesisName GetCompatibleHypothesis();
+
+  };
+
+  interface SMESH_1D_Algo : SMESH_Algo
+  {
+    /*!
+     * 
+     */
+  };
+
+  interface SMESH_2D_Algo : SMESH_Algo
+  {
+    /*!
+     * 
+     */
+  };
+
+  interface SMESH_3D_Algo : SMESH_Algo
+  {
+    /*!
+     * 
+     */
+  };
+};
+
+  // -----------------------------------------------------------------
+  // Specific Algorithms in separate idl file
+  // -----------------------------------------------------------------
+
+
+#endif
diff --git a/idl/SMESH_Mesh.idl b/idl/SMESH_Mesh.idl

new file mode 100644 (file)

index 0000000..a25ba26
--- /dev/null
+++ b/idl/SMESH_Mesh.idl
@@ -0,0 +1,262 @@
+//=============================================================================
+// File      : SMESH_Mesh.idl
+// Created   : jeu avr 11 15:31:39 CEST 2002
+// Author    : Paul RASCLE, EDF
+// Project   : SALOME
+// Copyright : EDF 2002
+// $Header$
+//=============================================================================
+
+
+#ifndef _SMESH_MESH_IDL_
+#define _SMESH_MESH_IDL_
+
+#include "SALOME_Exception.idl"
+#include "SMESH_Hypothesis.idl"
+
+#include "GEOM_Shape.idl"
+#include "MED.idl"
+
+module SMESH
+{
+  typedef sequence<double> double_array ;
+  typedef sequence<long> long_array ;
+  typedef sequence<string> string_array ;
+
+  enum log_command
+    {
+      ADD_NODE,
+      ADD_EDGE,
+      ADD_TRIANGLE,
+      ADD_QUADRANGLE,
+      ADD_TETRAHEDRON,
+      ADD_PYRAMID,
+      ADD_PRISM,
+      ADD_HEXAHEDRON,
+      REMOVE_NODE,
+      REMOVE_ELEMENT
+    };
+
+  struct log_block
+  {
+    long commandType;
+    long number;
+    double_array coords;
+    long_array indexes;
+  };
+  
+  typedef sequence<log_block> log_array;
+
+  interface SMESH_subMesh;
+  interface SMESH_MeshEditor;
+  interface SMESH_Mesh
+  {
+    /*!
+     * Associate a Shape to a Mesh created with NewEmpty
+     */
+
+//     boolean SetMesh(in GEOM::GEOM_Shape aShape)
+//        raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Get the subMesh object associated to a subShape. The subMesh object
+     * gives access to nodes and elements IDs.
+     * SubMesh will be used instead of SubShape in a next idl version to 
+     * adress a specific subMesh...
+     */
+    SMESH_subMesh GetElementsOnShape(in GEOM::GEOM_Shape aSubShape)
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Create a subMesh without reference to a subShape
+     */
+//     SMESH_subMesh NewEmpty()
+//       raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Add hypothesis to the mesh, under a particular subShape
+     * (or the main shape itself)
+     * The Add method is only used to prepare the build of the mesh and store
+     * the algorithms and associated parameters.
+     * Actual job of mesh the shape is done by MESH_Gen.
+     * @params
+     * - aSubShape : subShape obtained by a shape explode in GEOM
+     *   (or main shape)
+     * - anHyp : hypothesis object 
+     * @return
+     * - OK if the hypothesis is compatible with the subShape 
+     *   (and all previous hypothesis on the subShape)
+     * - NOK if the hypothesis is not compatible with the subShape
+     *   (or one previous hypothesis on the subShape)
+     * raises exception if hypothesis has not been created
+     */
+    boolean AddHypothesis(in GEOM::GEOM_Shape aSubShape, in SMESH_Hypothesis anHyp)
+      raises (SALOME::SALOME_Exception);
+//     boolean AddHypothesis(in SMESH_subMesh aSubMesh, in SMESH_Hypothesis anHyp)
+//       raises (SALOME::SALOME_Exception);
+
+
+    /*!
+     * Remove an hypothesis previouly added with AddHypothesis.
+     */
+    boolean RemoveHypothesis(in GEOM::GEOM_Shape aSubShape,
+                            in SMESH_Hypothesis anHyp)
+      raises (SALOME::SALOME_Exception);
+//     boolean RemoveHypothesis(in SMESH_subMesh aSubMesh,
+//                          in SMESH_Hypothesis anHyp)
+//       raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Get the list of hypothesis added on a subShape
+     */
+    ListOfHypothesis GetHypothesisList(in GEOM::GEOM_Shape aSubShape)
+      raises (SALOME::SALOME_Exception);
+//     ListOfHypothesis GetHypothesisList(in SMESH_subMesh aSubMesh)
+//       raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Get the log of nodes and elements added or removed since previous
+     * clear of the log.
+     * @params
+     * - clearAfterGet : log is emptied after Get (safe if concurrents access)
+     */
+    //    string_array GetLog(in boolean clearAfterGet)
+    //      raises (SALOME::SALOME_Exception);
+    log_array GetLog(in boolean clearAfterGet)
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Clear the log of nodes and elements added or removed since previous
+     * clear. Must be used immediately after GetLog if clearAfterGet is false.
+     */
+    void ClearLog()
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Get the internal Id 
+     */
+    long GetId();     
+
+    /*!
+     * Get the study Id 
+     */
+    long GetStudyId();    
+    
+    SMESH_MeshEditor GetMeshEditor() 
+       raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Export Mesh with DAT and MED Formats
+     */
+    void ExportDAT( in string file )
+      raises (SALOME::SALOME_Exception);
+    void ExportMED( in string file )
+      raises (SALOME::SALOME_Exception);
+    void ExportUNV( in string file )
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Get MED Mesh
+     */
+    SALOME_MED::MESH GetMEDMesh()
+      raises (SALOME::SALOME_Exception);
+
+    long NbNodes()
+      raises (SALOME::SALOME_Exception);
+
+    long NbEdges()
+      raises (SALOME::SALOME_Exception);
+
+    long NbFaces()
+      raises (SALOME::SALOME_Exception);
+
+    long NbTriangles()
+      raises (SALOME::SALOME_Exception);
+
+    long NbQuadrangles()
+      raises (SALOME::SALOME_Exception);
+
+    long NbVolumes()
+      raises (SALOME::SALOME_Exception);
+
+    long NbTetras()
+      raises (SALOME::SALOME_Exception);
+
+    long NbHexas()
+      raises (SALOME::SALOME_Exception);
+
+    long NbSubMesh()
+      raises (SALOME::SALOME_Exception);
+  };
+
+  interface SMESH_subMesh
+  {
+    /*!
+     * 
+     */
+    long GetNumberOfElements()
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * 
+     */
+    long GetNumberOfNodes()
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * 
+     */
+    long_array GetElementsId()
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * 
+     */
+    long_array GetNodesId()
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Get SMESH_Mesh which stores nodes coordinates & elements definition
+     */
+    SMESH_Mesh GetFather()
+      raises (SALOME::SALOME_Exception);
+
+    /*!
+     * Get the internal Id 
+     */
+    long GetId();    
+
+    /*!
+     * Get MED subMesh
+     */
+    SALOME_MED::FAMILY GetFamily()
+      raises (SALOME::SALOME_Exception);
+  };
+  
+  /* 
+   * This interface makes modifications on the Mesh - removing elements and nodes
+   */     
+  interface SMESH_MeshEditor {
+
+
+    boolean RemoveElements(in long_array IDsOfElements)
+       raises (SALOME::SALOME_Exception);
+
+    boolean RemoveNodes(in long_array IDsOfNodes)
+      raises (SALOME::SALOME_Exception);
+
+    boolean AddNode(in double x, in double y, in double z)
+      raises (SALOME::SALOME_Exception);
+
+    boolean AddEdge(in long_array IDsOfNodes)
+      raises (SALOME::SALOME_Exception);
+
+    boolean AddFace(in long_array IDsOfNodes)
+      raises (SALOME::SALOME_Exception);
+
+    boolean AddVolume(in long_array IDsOfNodes)
+      raises (SALOME::SALOME_Exception);
+  };
+};
+
+#endif
diff --git a/src/Driver/Document_Reader.cxx b/src/Driver/Document_Reader.cxx

new file mode 100644 (file)

index 0000000..3a0ec88
--- /dev/null
+++ b/src/Driver/Document_Reader.cxx
@@ -0,0 +1,13 @@
+using namespace std;
+#include "Document_Reader.h"
+
+#include "utilities.h"
+
+void Document_Reader::SetFile(string aFile) {
+  myFile = aFile;
+}
+
+void Document_Reader::SetDocument(Handle(SMESHDS_Document)& aDoc) {
+  myDocument = aDoc;
+}
+
diff --git a/src/Driver/Document_Reader.h b/src/Driver/Document_Reader.h

new file mode 100644 (file)

index 0000000..46472e5
--- /dev/null
+++ b/src/Driver/Document_Reader.h
@@ -0,0 +1,18 @@
+#ifndef _INCLUDE_DOCUMENT_READER
+#define _INCLUDE_DOCUMENT_READER
+
+#include "SMESHDS_Document.hxx"
+#include <string>
+
+class Document_Reader {
+  
+  public :
+    virtual void Read() =0;
+    void SetFile(string);
+    void SetDocument(Handle(SMESHDS_Document)&);
+
+  protected :
+    Handle_SMESHDS_Document myDocument;
+    string myFile; 
+};
+#endif
diff --git a/src/Driver/Document_Writer.cxx b/src/Driver/Document_Writer.cxx

new file mode 100644 (file)

index 0000000..73ad7f7
--- /dev/null
+++ b/src/Driver/Document_Writer.cxx
@@ -0,0 +1,11 @@
+using namespace std;
+#include "Document_Writer.h"
+
+void Document_Writer::SetFile(string aFile) {
+  myFile = aFile;
+}
+
+void Document_Writer::SetDocument(Handle(SMESHDS_Document)& aDoc) {
+  myDocument = aDoc;
+}
+
diff --git a/src/Driver/Document_Writer.h b/src/Driver/Document_Writer.h

new file mode 100644 (file)

index 0000000..4d3a67f
--- /dev/null
+++ b/src/Driver/Document_Writer.h
@@ -0,0 +1,19 @@
+#ifndef _INCLUDE_DOCUMENT_WRITER
+#define _INCLUDE_DOCUMENT_WRITER
+
+#include "SMESHDS_Document.hxx"
+#include <string>
+
+class Document_Writer {
+
+  public :
+    virtual void Write() =0;
+    void SetFile(string);
+    void SetDocument(Handle(SMESHDS_Document)&);
+
+  protected :
+    Handle_SMESHDS_Document myDocument;
+    string myFile;
+
+};
+#endif
diff --git a/src/Driver/Driver_dl.cxx b/src/Driver/Driver_dl.cxx

new file mode 100644 (file)

index 0000000..81dd423
--- /dev/null
+++ b/src/Driver/Driver_dl.cxx
@@ -0,0 +1,183 @@
+using namespace std;
+#include "Driver.h"
+
+#include <dlfcn.h>
+#include <stdio.h>
+#include <utilities.h>
+
+#include "DriverDAT_R_SMESHDS_Document.h"
+#include "Test.h"
+
+/*!  extern "C"
+{
+  void test() {
+
+                  void *handle;
+                  double (*cosine)(double);
+                  char *error;
+
+                  handle = dlopen ("/usr/lib/libm.so", RTLD_LAZY);
+                  if (!handle) {
+                      fputs (dlerror(), stderr);
+                      exit(1);
+                  }
+
+                  cosine = dlsym(handle, "cos");
+                  if ((error = dlerror()) != NULL)  {
+                      fprintf (stderr, "%s\n", error);
+                      exit(1);
+                  }
+
+                  printf ("%f\n", (*cosine)(2.0));
+                  dlclose(handle);
+
+  char* error;
+  string Extension=string("DAT");
+  string Class=string("SMESHDS_Document");
+  string myLibrary = string("/home/barberou/barberou/SALOME_3105/build/lib/libMeshDriver")+Extension+string(".so");
+  SCRUTE(myLibrary);
+  //Document_Reader* myDriver;//a caster ???
+  DriverDAT_R_SMESHDS_Document* myDriver;
+
+  string myClass = string("Driver")+Extension+string("_R_")+Class;
+  SCRUTE(myClass);
+  
+  void* handle = dlopen (myLibrary.c_str(), RTLD_LAZY);
+  if (!handle) {
+    fputs (dlerror(), stderr);
+    exit(1);
+  }
+  MESSAGE("Open ok");
+
+  //int* res = (int*)dlsym(handle, "getOne");
+  //SCRUTE(res);
+  SCRUTE(dlsym(handle, "getOne"));
+  //int res2= (*res)();
+  myDriver = (DriverDAT_R_SMESHDS_Document*) dlsym(handle, myClass.c_str());
+  MESSAGE("Reading 1");
+  SCRUTE(myDriver);
+  if ((error = dlerror()) != NULL)  {
+    fprintf (stderr, "%s\n", error);
+    exit(1);
+  }
+  MESSAGE("Reading 2");
+    
+  dlclose(handle);
+  MESSAGE("after close");
+                 
+  }
+  }*/
+
+Document_Reader* Driver::GetDocumentReader(string Extension, string Class) {
+  test();
+  //p-e extern C ?
+  /*!
+  char* error;
+  string myLibrary = string("/home/barberou/barberou/SALOME_3105/build/lib/libMeshDriver")+Extension+string(".so");
+  SCRUTE(myLibrary);
+  //Document_Reader* myDriver;//a caster ???
+  DriverDAT_R_SMESHDS_Document* myDriver;
+
+  string myClass = string("Driver")+Extension+string("_R_")+Class;
+  SCRUTE(myClass);
+  
+  void* handle = dlopen (myLibrary.c_str(), RTLD_LAZY);
+  if (!handle) {
+    fputs (dlerror(), stderr);
+    exit(1);
+  }
+  MESSAGE("Open ok");
+  //myDriver = (Document_Reader*) dlsym(handle, myClass.c_str());
+  int* res = (int*) dlsym(handle, "getOne");
+  SCRUTE(res);
+  myDriver = (DriverDAT_R_SMESHDS_Document*) dlsym(handle, myClass.c_str());
+  MESSAGE("Reading 1");
+  SCRUTE(myDriver);
+  if ((error = dlerror()) != NULL)  {
+    fprintf (stderr, "%s\n", error);
+    exit(1);
+  }
+  MESSAGE("Reading 2");
+    
+  dlclose(handle);
+  MESSAGE("after close");
+  return (myDriver);
+  */
+
+}
+
+Document_Writer* Driver::GetDocumentWriter(string Extension, string Class) {
+
+  char* error;
+  string myLibrary = string("libMeshDriver")+Extension+string(".so");
+  Document_Writer* myDriver;//a caster ???
+  string myClass = string("Driver")+Extension+string("_W_")+Class;
+  
+  void* handle = dlopen (myLibrary.c_str(), RTLD_LAZY);
+  if (!handle) {
+    fputs (dlerror(), stderr);
+    exit(1);
+  }
+
+  myDriver = (Document_Writer*) dlsym(handle, myClass.c_str());
+  if ((error = dlerror()) != NULL)  {
+    fprintf (stderr, "%s\n", error);
+    exit(1);
+  }
+    
+  dlclose(handle);
+  
+  return (myDriver);
+
+}
+
+Mesh_Reader* Driver::GetMeshReader(string Extension, string Class) {
+
+  char* error;
+  string myLibrary = string("libMeshDriver")+Extension+string(".so");
+  Mesh_Reader* myDriver;//a caster ???
+  string myClass = string("Driver")+Extension+string("_R_")+Class;
+  
+  void* handle = dlopen (myLibrary.c_str(), RTLD_LAZY);
+  if (!handle) {
+    fputs (dlerror(), stderr);
+    exit(1);
+  }
+
+  myDriver = (Mesh_Reader*) dlsym(handle, myClass.c_str());
+  if ((error = dlerror()) != NULL)  {
+    fprintf (stderr, "%s\n", error);
+    exit(1);
+  }
+    
+  dlclose(handle);
+  
+  return (myDriver);
+
+}
+
+Mesh_Writer* Driver::GetMeshWriter(string Extension, string Class) {
+
+  char* error;
+  string myLibrary = string("libMeshDriver")+Extension+string(".so");
+  Mesh_Writer* myDriver;//a caster ???
+  string myClass = string("Driver")+Extension+string("_W_")+Class;
+  
+  void* handle = dlopen (myLibrary.c_str(), RTLD_LAZY);
+  if (!handle) {
+    fputs (dlerror(), stderr);
+    exit(1);
+  }
+
+  myDriver = (Mesh_Writer*) dlsym(handle, myClass.c_str());
+  if ((error = dlerror()) != NULL)  {
+    fprintf (stderr, "%s\n", error);
+    exit(1);
+  }
+    
+  dlclose(handle);
+  
+  return (myDriver);
+
+}
+
diff --git a/src/Driver/Makefile.in b/src/Driver/Makefile.in

new file mode 100644 (file)

index 0000000..c8bb273
--- /dev/null
+++ b/src/Driver/Makefile.in
@@ -0,0 +1,39 @@
+# -* Makefile *- 
+#
+# Author : Marc Tajchman (CEA)
+# Date : 5/07/2001
+# $Header$
+#
+
+# source path
+top_srcdir=@top_srcdir@
+top_builddir=../..
+srcdir=@srcdir@
+VPATH=.:@srcdir@
+
+
+@COMMENCE@
+
+# header files 
+EXPORT_HEADERS= Document_Reader.h Document_Writer.h Mesh_Reader.h Mesh_Writer.h
+
+# Libraries targets
+LIB = libMeshDriver.la
+LIB_SRC = Document_Reader.cxx Document_Writer.cxx Mesh_Reader.cxx Mesh_Writer.cxx
+
+LIB_CLIENT_IDL =
+
+LIB_SERVER_IDL =
+
+# additionnal information to compil and link file
+CPPFLAGS += $(OCC_INCLUDES) -I${KERNEL_ROOT_DIR}/include/salome
+CXXFLAGS += $(OCC_CXXFLAGS) $(MED2_INCLUDES) -rdynamic -ldl -I${KERNEL_ROOT_DIR}/include/salome
+LDFLAGS  += $(OCC_LIBS) $(MED2_LIBS) -lSMESHDS -lSMDS
+
+%_moc.cxx: %.h
+       $(MOC) $< -o $@
+
+@CONCLUDE@
+
+
+
diff --git a/src/Driver/Mesh_Reader.cxx b/src/Driver/Mesh_Reader.cxx

new file mode 100644 (file)

index 0000000..3001d36
--- /dev/null
+++ b/src/Driver/Mesh_Reader.cxx
@@ -0,0 +1,3 @@
+using namespace std;
+#include "Mesh_Reader.h"
+
diff --git a/src/Driver/Mesh_Reader.h b/src/Driver/Mesh_Reader.h

new file mode 100644 (file)

index 0000000..8671f5d
--- /dev/null
+++ b/src/Driver/Mesh_Reader.h
@@ -0,0 +1,17 @@
+#ifndef _INCLUDE_MESH_READER
+#define _INCLUDE_MESH_READER
+
+#include <string>
+#include "Handle_SMDS_Mesh.hxx"
+
+class Mesh_Reader {
+
+  public :
+    virtual void Add() =0;
+    virtual void Read() =0;
+    virtual void SetMesh(Handle(SMDS_Mesh)&) =0;
+    virtual void SetMeshId(int) =0;
+    virtual void SetFile(string) =0;
+
+};
+#endif
diff --git a/src/Driver/Mesh_Writer.cxx b/src/Driver/Mesh_Writer.cxx

new file mode 100644 (file)

index 0000000..034e68f
--- /dev/null
+++ b/src/Driver/Mesh_Writer.cxx
@@ -0,0 +1,2 @@
+using namespace std;
+#include "Mesh_Writer.h"
diff --git a/src/Driver/Mesh_Writer.h b/src/Driver/Mesh_Writer.h

new file mode 100644 (file)

index 0000000..e9d84df
--- /dev/null
+++ b/src/Driver/Mesh_Writer.h
@@ -0,0 +1,17 @@
+#ifndef _INCLUDE_MESH_WRITER
+#define _INCLUDE_MESH_WRITER
+
+#include <string>
+#include "Handle_SMDS_Mesh.hxx"
+
+class Mesh_Writer {
+
+  public :
+    virtual void Add() =0;
+    virtual void Write() =0;
+    virtual void SetMesh(Handle(SMDS_Mesh)&) =0;
+    virtual void SetFile(string) =0;
+    virtual void SetMeshId(int) =0;
+
+};
+#endif
diff --git a/src/Driver/SMESHDriver.cxx b/src/Driver/SMESHDriver.cxx

new file mode 100644 (file)

index 0000000..2bd13fa
--- /dev/null
+++ b/src/Driver/SMESHDriver.cxx
@@ -0,0 +1,113 @@
+using namespace std;
+#include "SMESHDriver.h"
+
+#include <dlfcn.h>
+#include <utilities.h>
+
+//A enlever
+#include "DriverMED_R_SMESHDS_Document.h"
+#include "DriverMED_R_SMESHDS_Mesh.h"
+#include "DriverMED_R_SMDS_Mesh.h"
+#include "DriverMED_W_SMESHDS_Document.h"
+#include "DriverMED_W_SMESHDS_Mesh.h"
+#include "DriverMED_W_SMDS_Mesh.h"
+
+#include "DriverDAT_R_SMESHDS_Document.h"
+#include "DriverDAT_R_SMESHDS_Mesh.h"
+#include "DriverDAT_R_SMDS_Mesh.h"
+#include "DriverDAT_W_SMESHDS_Document.h"
+#include "DriverDAT_W_SMESHDS_Mesh.h"
+#include "DriverDAT_W_SMDS_Mesh.h"
+//
+
+Document_Reader* SMESHDriver::GetDocumentReader(string Extension, string Class) {
+  if (Extension==string("MED")) {
+    DriverMED_R_SMESHDS_Document* myDriver = new DriverMED_R_SMESHDS_Document();
+    return (myDriver);
+  }
+  else if (Extension==string("DAT")) {
+    DriverDAT_R_SMESHDS_Document* myDriver = new DriverDAT_R_SMESHDS_Document();
+    return (myDriver);
+  }
+  else {
+    MESSAGE("No driver known for this extension");
+    return (Document_Reader*)NULL;
+  }
+
+
+}
+
+Document_Writer* SMESHDriver::GetDocumentWriter(string Extension, string Class) {
+  if (Extension==string("MED")) {
+    DriverMED_W_SMESHDS_Document* myDriver = new DriverMED_W_SMESHDS_Document();
+    return (myDriver);
+  }    
+  else if (Extension==string("DAT")) {
+    DriverDAT_W_SMESHDS_Document* myDriver = new DriverDAT_W_SMESHDS_Document();
+    return (myDriver);
+  }    
+  else {
+    MESSAGE("No driver known for this extension");
+    return (Document_Writer*)NULL;
+  }
+
+
+}
+
+Mesh_Reader* SMESHDriver::GetMeshReader(string Extension, string Class) {
+  if (Extension==string("MED")) {
+
+  if (strcmp(Class.c_str(),"SMESHDS_Mesh")==0) {
+    DriverMED_R_SMESHDS_Mesh* myDriver = new DriverMED_R_SMESHDS_Mesh();
+    return (myDriver);
+  }
+  else if (strcmp(Class.c_str(),"SMDS_Mesh")==0) {
+    DriverMED_R_SMDS_Mesh* myDriver = new DriverMED_R_SMDS_Mesh();
+    return (myDriver);
+  }
+
+  }
+  else if (Extension==string("DAT")) {
+
+  if (strcmp(Class.c_str(),"SMESHDS_Mesh")==0) {
+    DriverDAT_R_SMESHDS_Mesh* myDriver = new DriverDAT_R_SMESHDS_Mesh();
+    return (myDriver);
+  }
+  else if (strcmp(Class.c_str(),"SMDS_Mesh")==0) {
+    DriverDAT_R_SMDS_Mesh* myDriver = new DriverDAT_R_SMDS_Mesh();
+    return (myDriver);
+  }
+
+  }
+
+
+}
+
+Mesh_Writer* SMESHDriver::GetMeshWriter(string Extension, string Class) {
+  if (Extension==string("MED")) {
+
+  if (strcmp(Class.c_str(),"SMESHDS_Mesh")==0) {
+    DriverMED_W_SMESHDS_Mesh* myDriver = new DriverMED_W_SMESHDS_Mesh();
+    return (myDriver);
+  }
+  else if (strcmp(Class.c_str(),"SMDS_Mesh")==0) {
+    DriverMED_W_SMDS_Mesh* myDriver = new DriverMED_W_SMDS_Mesh();
+    return (myDriver);
+  }
+
+  }
+  else if (Extension==string("DAT")) {
+
+  if (strcmp(Class.c_str(),"SMESHDS_Mesh")==0) {
+    DriverDAT_W_SMESHDS_Mesh* myDriver = new DriverDAT_W_SMESHDS_Mesh();
+    return (myDriver);
+  }
+  else if (strcmp(Class.c_str(),"SMDS_Mesh")==0) {
+    DriverDAT_W_SMDS_Mesh* myDriver = new DriverDAT_W_SMDS_Mesh();
+    return (myDriver);
+  }
+
+  }
+
+}
+
diff --git a/src/Driver/SMESHDriver.h b/src/Driver/SMESHDriver.h

new file mode 100644 (file)

index 0000000..ba0e7ad
--- /dev/null
+++ b/src/Driver/SMESHDriver.h
@@ -0,0 +1,19 @@
+#ifndef _INCLUDE_SMESHDRIVER
+#define _INCLUDE_SMESHDRIVER
+
+#include "Document_Reader.h"
+#include "Document_Writer.h"
+#include "Mesh_Reader.h"
+#include "Mesh_Writer.h"
+
+class SMESHDriver {
+
+  public :
+    static Document_Reader* GetDocumentReader(string Extension, string Class);
+    static Document_Writer* GetDocumentWriter(string Extension, string Class);
+
+    static Mesh_Reader* GetMeshReader(string Extension, string Class);
+    static Mesh_Writer* GetMeshWriter(string Extension, string Class);
+
+};
+#endif
diff --git a/src/DriverDAT/DriverDAT_R_SMDS_Mesh.cxx b/src/DriverDAT/DriverDAT_R_SMDS_Mesh.cxx

new file mode 100644 (file)

index 0000000..72e0d30
--- /dev/null
+++ b/src/DriverDAT/DriverDAT_R_SMDS_Mesh.cxx
@@ -0,0 +1,175 @@
+using namespace std;
+#include "DriverDAT_R_SMDS_Mesh.h"
+
+#include "utilities.h"
+
+DriverDAT_R_SMDS_Mesh::DriverDAT_R_SMDS_Mesh() {
+;
+}
+
+DriverDAT_R_SMDS_Mesh::~DriverDAT_R_SMDS_Mesh() {
+;
+}
+
+void DriverDAT_R_SMDS_Mesh::SetMesh(Handle(SMDS_Mesh)& aMesh) {
+  myMesh = aMesh;
+}
+
+void DriverDAT_R_SMDS_Mesh::SetFile(string aFile) {
+  myFile = aFile;
+}
+
+void DriverDAT_R_SMDS_Mesh::SetFileId(FILE* aFileId) {
+  myFileId = aFileId;
+}
+
+void DriverDAT_R_SMDS_Mesh::SetMeshId(int aMeshId) {
+  myMeshId = aMeshId;
+}
+
+void DriverDAT_R_SMDS_Mesh::Add() {
+  ;
+}
+
+void DriverDAT_R_SMDS_Mesh::Read() {
+
+  int i,j;
+  int nbNodes,nbCells;
+  int intNumPoint;
+  float coordX, coordY, coordZ;
+  int nbNoeuds;
+
+  int intNumMaille,Degre;
+  int ValElement;
+  int ValDegre;
+  int NoeudsMaille[20];
+  int NoeudMaille;
+
+  bool ok;
+
+  MESSAGE("in DriverDAT_R_SMDS_Mesh::Read()");
+  /****************************************************************************
+  *                      OUVERTURE DU FICHIER EN LECTURE                      *
+  ****************************************************************************/
+  char* file2Read = (char*)myFile.c_str();
+  myFileId = fopen(file2Read,"r");
+  if (myFileId < 0)
+    {
+      fprintf(stderr,">> ERREUR : ouverture du fichier %s \n",file2Read);
+      exit(EXIT_FAILURE);
+    }
+
+  fscanf(myFileId,"%d %d\n",&nbNodes,&nbCells);
+
+  /****************************************************************************
+  *                       LECTURE DES NOEUDS                                  *
+  ****************************************************************************/
+  fprintf(stdout,"\n(************************)\n");
+  fprintf(stdout,"(* NOEUDS DU MAILLAGE : *)\n");
+  fprintf(stdout,"(************************)\n");
+
+  for (i=0;i<nbNodes;i++) {
+    fscanf(myFileId,"%d %e %e %e\n",&intNumPoint,&coordX,&coordY,&coordZ);
+    ok = myMesh->AddNodeWithID(coordX,coordY,coordZ,intNumPoint);
+  }
+
+  fprintf(stdout,"%d noeuds\n",myMesh->NbNodes());
+  /****************************************************************************
+  *                       LECTURE DES ELEMENTS                                *
+  ****************************************************************************/
+  fprintf(stdout,"\n(**************************)\n");
+  fprintf(stdout,"(* ELEMENTS DU MAILLAGE : *)\n");
+  fprintf(stdout,"(**************************)");
+
+  fprintf(stdout,"%d elements\n",nbCells);
+
+  for (i=0; i<nbCells; i++) {
+    fscanf(myFileId,"%d %d",&intNumMaille,&ValElement);
+    Degre=abs(ValElement/100);
+    nbNoeuds=ValElement-(Degre*100);
+    
+    // Recuperation des noeuds de la maille
+    for (j=0; j<nbNoeuds; j++) {
+      fscanf(myFileId,"%d",&NoeudMaille);
+      NoeudsMaille[j]=NoeudMaille;
+    }
+    
+    // Analyse des cas de cellules
+    switch (ValElement) {
+    case 102 : ;
+    case 103 : ;
+      {
+      ValDegre=3;
+      nbNoeuds=2;
+      ok = myMesh->AddEdgeWithID(NoeudsMaille[0],NoeudsMaille[1],intNumMaille);
+      break;
+    }
+    case 204 : ;
+    case 208 : ;
+      {
+      ValDegre=9;
+      nbNoeuds=4;
+      ok = myMesh->AddFaceWithID(NoeudsMaille[0],NoeudsMaille[1],NoeudsMaille[2],NoeudsMaille[3],intNumMaille);
+      break;
+    }
+    case 203 : ;
+    case 206 : ;
+      {
+      ValDegre=5;
+      nbNoeuds=3;
+      ok = myMesh->AddFaceWithID(NoeudsMaille[0],NoeudsMaille[1],NoeudsMaille[2],intNumMaille);
+      break;
+    }
+    case 308 : ;
+    case 320 : ;
+      {
+      ValDegre=12;
+      nbNoeuds=8;
+      if (ValElement==320) {
+       //A voir, correspondance VTK
+       NoeudsMaille[4]=NoeudsMaille[8];
+       NoeudsMaille[5]=NoeudsMaille[9];
+       NoeudsMaille[6]=NoeudsMaille[10];
+       NoeudsMaille[7]=NoeudsMaille[11];
+      }
+      ok = myMesh->AddVolumeWithID(NoeudsMaille[0],NoeudsMaille[1],NoeudsMaille[2],NoeudsMaille[3],NoeudsMaille[4],NoeudsMaille[5],NoeudsMaille[6],NoeudsMaille[7],intNumMaille);
+     break;
+    }
+    case 304 : ;
+    case 310 : ;
+      {
+      ValDegre=10;
+      nbNoeuds=4;
+      if (ValElement==310)
+       NoeudsMaille[3]=NoeudsMaille[6];
+      ok = myMesh->AddVolumeWithID(NoeudsMaille[0],NoeudsMaille[1],NoeudsMaille[2],NoeudsMaille[3],intNumMaille);
+     break;
+    }
+    case 306 : ;
+    case 315 : ;
+      {
+      ValDegre=12;
+      nbNoeuds=8;
+      if (ValElement==315) {
+       NoeudsMaille[3]=NoeudsMaille[6];
+       NoeudsMaille[4]=NoeudsMaille[7];
+       NoeudsMaille[5]=NoeudsMaille[8];
+      }
+      NoeudsMaille[7]=NoeudsMaille[5];
+      NoeudsMaille[6]=NoeudsMaille[5];
+      NoeudsMaille[5]=NoeudsMaille[4];
+      NoeudsMaille[4]=NoeudsMaille[3];
+      NoeudsMaille[3]=NoeudsMaille[2];
+      ok = myMesh->AddVolumeWithID(NoeudsMaille[0],NoeudsMaille[1],NoeudsMaille[2],NoeudsMaille[3],NoeudsMaille[4],NoeudsMaille[5],intNumMaille);
+      break;
+    }
+
+    }
+  }
+
+  /****************************************************************************
+  *                      FERMETURE DU FICHIER                      *
+  ****************************************************************************/
+  fclose(myFileId);
+
+}
diff --git a/src/DriverDAT/DriverDAT_R_SMDS_Mesh.h b/src/DriverDAT/DriverDAT_R_SMDS_Mesh.h

new file mode 100644 (file)

index 0000000..0a8a86a
--- /dev/null
+++ b/src/DriverDAT/DriverDAT_R_SMDS_Mesh.h
@@ -0,0 +1,30 @@
+#ifndef _INCLUDE_DRIVERDAT_R_SMDS_MESH
+#define _INCLUDE_DRIVERDAT_R_SMDS_MESH
+
+#include <stdio.h>
+
+#include "SMDS_Mesh.hxx"
+#include "Mesh_Reader.h"
+
+class DriverDAT_R_SMDS_Mesh : public Mesh_Reader {
+
+  public :
+    DriverDAT_R_SMDS_Mesh();
+    ~DriverDAT_R_SMDS_Mesh();
+
+    void Add();
+    void Read();
+    void SetMesh(Handle(SMDS_Mesh)& aMesh);
+    void SetFile(string);
+
+    void SetFileId(FILE*);
+    void SetMeshId(int);
+
+  private :
+    Handle_SMDS_Mesh myMesh;
+    string myFile; 
+    FILE* myFileId;
+    int myMeshId;
+
+};
+#endif
diff --git a/src/DriverDAT/DriverDAT_R_SMESHDS_Document.cxx b/src/DriverDAT/DriverDAT_R_SMESHDS_Document.cxx

new file mode 100644 (file)

index 0000000..a079dbd
--- /dev/null
+++ b/src/DriverDAT/DriverDAT_R_SMESHDS_Document.cxx
@@ -0,0 +1,86 @@
+using namespace std;
+#include "DriverDAT_R_SMESHDS_Document.h"
+#include "DriverDAT_R_SMESHDS_Mesh.h"
+
+#include "utilities.h"
+
+int getOne() {
+  printf("in getOne");
+  return (1);
+}
+
+extern "C" {
+  //  Document_Reader* maker() {
+  DriverDAT_R_SMESHDS_Document* maker() {
+    fprintf(stdout,"here in maker\n");
+    return new DriverDAT_R_SMESHDS_Document;
+  }
+}
+
+DriverDAT_R_SMESHDS_Document::DriverDAT_R_SMESHDS_Document() {
+  myFile = string("");
+}
+
+DriverDAT_R_SMESHDS_Document::~DriverDAT_R_SMESHDS_Document() {
+;
+}
+
+//void DriverDAT_R_SMESHDS_Document::SetFile(string aFile) {
+//myFile = aFile;
+//}
+
+//void DriverDAT_R_SMESHDS_Document::SetDocument(Handle(SMESHDS_Document)& aDoc) {
+//myDocument = aDoc;
+//}
+
+void DriverDAT_R_SMESHDS_Document::Read() {
+
+  int myMeshId;
+  MESSAGE("in read");
+  SCRUTE(myFile);
+  //Handle(SMESHDS_Document) myDocument = new SMESHDS_Document(1);
+
+  /****************************************************************************
+  *                      OUVERTURE DU FICHIER EN LECTURE                      *
+  ****************************************************************************/
+  char* file2Read = (char*)myFile.c_str();
+  FILE* fid = fopen(file2Read,"r");
+  if (fid < 0)
+    {
+      fprintf(stderr,">> ERREUR : ouverture du fichier %s \n",file2Read);
+      exit(EXIT_FAILURE);
+    }
+
+  /****************************************************************************
+  *                      COMBIEN DE MAILLAGES ?                               *
+  ****************************************************************************/
+  int nmaa = 1;
+
+  /****************************************************************************
+  *                      FERMETURE DU FICHIER                                 *
+  ****************************************************************************/
+  fclose(fid);
+  
+  printf("Nombre de maillages = %d\n",nmaa);
+
+  string myClass = string("SMESHDS_Mesh");
+  string myExtension = string("DAT");
+
+  for (int meshIt=1;meshIt<=nmaa;meshIt++) {
+    myMeshId = myDocument->NewMesh();
+
+    Handle(SMDS_Mesh) myMesh = myDocument->GetMesh(myMeshId);
+
+    DriverDAT_R_SMESHDS_Mesh* myReader = new DriverDAT_R_SMESHDS_Mesh;
+
+
+    myReader->SetMesh(myMesh);
+    myReader->SetFile(myFile);
+    //myReader->SetFileId(fid);
+
+    myReader->Read();
+
+  }
+
+
+}
diff --git a/src/DriverDAT/DriverDAT_R_SMESHDS_Document.h b/src/DriverDAT/DriverDAT_R_SMESHDS_Document.h

new file mode 100644 (file)

index 0000000..83b67b5
--- /dev/null
+++ b/src/DriverDAT/DriverDAT_R_SMESHDS_Document.h
@@ -0,0 +1,24 @@
+#ifndef _INCLUDE_DRIVERDAT_R_SMESHDS_DOCUMENT
+#define _INCLUDE_DRIVERDAT_R_SMESHDS_DOCUMENT
+
+#include <stdio.h>
+
+#include "SMESHDS_Document.hxx"
+#include "Document_Reader.h"
+
+class DriverDAT_R_SMESHDS_Document : public Document_Reader {
+  
+public :
+  DriverDAT_R_SMESHDS_Document();
+  ~DriverDAT_R_SMESHDS_Document();
+ 
+  void Read();
+  //void SetFile(string);
+  //void SetDocument(Handle_SMESHDS_Document&);
+
+private :
+  //Handle_SMESHDS_Document myDocument;
+  //string myFile; 
+
+};
+#endif
diff --git a/src/DriverDAT/DriverDAT_R_SMESHDS_Mesh.cxx b/src/DriverDAT/DriverDAT_R_SMESHDS_Mesh.cxx

new file mode 100644 (file)

index 0000000..56fe869
--- /dev/null
+++ b/src/DriverDAT/DriverDAT_R_SMESHDS_Mesh.cxx
@@ -0,0 +1,52 @@
+using namespace std;
+#include "DriverDAT_R_SMESHDS_Mesh.h"
+#include "DriverDAT_R_SMDS_Mesh.h"
+
+#include "utilities.h"
+
+DriverDAT_R_SMESHDS_Mesh::DriverDAT_R_SMESHDS_Mesh() {
+;
+}
+
+DriverDAT_R_SMESHDS_Mesh::~DriverDAT_R_SMESHDS_Mesh() {
+;
+}
+
+void DriverDAT_R_SMESHDS_Mesh::SetMesh(Handle(SMDS_Mesh)& aMesh) {
+  //myMesh = Handle(SMESHDS_Mesh)::DownCast(aMesh);
+  myMesh = aMesh;
+}
+
+void DriverDAT_R_SMESHDS_Mesh::SetFile(string aFile) {
+  myFile = aFile;
+}
+
+void DriverDAT_R_SMESHDS_Mesh::SetFileId(FILE* aFileId) {
+  myFileId = aFileId;
+}
+
+void DriverDAT_R_SMESHDS_Mesh::SetMeshId(int aMeshId) {
+  myMeshId = aMeshId;
+}
+
+void DriverDAT_R_SMESHDS_Mesh::Add() {
+  ;
+}
+
+void DriverDAT_R_SMESHDS_Mesh::Read() {
+  string myClass = string("SMDS_Mesh");
+  string myExtension = string("DAT");
+
+  MESSAGE("in DriverDAT_R_SMESHDS_Mesh::Read() 1");
+  DriverDAT_R_SMDS_Mesh* myReader = new DriverDAT_R_SMDS_Mesh;
+
+  MESSAGE("in DriverDAT_R_SMESHDS_Mesh::Read() 2");
+  myReader->SetMesh(myMesh);
+  MESSAGE("in DriverDAT_R_SMESHDS_Mesh::Read() 3");
+  myReader->SetFile(myFile);
+  //myReader->SetFileId(myFileId);
+
+  MESSAGE("in DriverDAT_R_SMESHDS_Mesh::Read() 4");
+  myReader->Read();
+
+}
diff --git a/src/DriverDAT/DriverDAT_R_SMESHDS_Mesh.h b/src/DriverDAT/DriverDAT_R_SMESHDS_Mesh.h

new file mode 100644 (file)

index 0000000..bbcbd00
--- /dev/null
+++ b/src/DriverDAT/DriverDAT_R_SMESHDS_Mesh.h
@@ -0,0 +1,30 @@
+#ifndef _INCLUDE_DRIVERDAT_R_SMESHDS_MESH
+#define _INCLUDE_DRIVERDAT_R_SMESHDS_MESH
+
+#include <stdio.h>
+
+#include "SMESHDS_Mesh.hxx"
+#include "Mesh_Reader.h"
+
+class DriverDAT_R_SMESHDS_Mesh : public Mesh_Reader {
+
+  public :
+    DriverDAT_R_SMESHDS_Mesh();
+    ~DriverDAT_R_SMESHDS_Mesh();
+
+    void Add();
+    void Read();
+    void SetMesh(Handle(SMDS_Mesh)& aMesh);
+    void SetFile(string);
+
+    void SetFileId(FILE*);
+    void SetMeshId(int);
+
+private :
+    Handle_SMDS_Mesh myMesh;
+    string myFile; 
+    FILE* myFileId;
+    int myMeshId;
+
+};
+#endif
diff --git a/src/DriverDAT/DriverDAT_W_SMDS_Mesh.cxx b/src/DriverDAT/DriverDAT_W_SMDS_Mesh.cxx

new file mode 100644 (file)

index 0000000..0cc429c
--- /dev/null
+++ b/src/DriverDAT/DriverDAT_W_SMDS_Mesh.cxx
@@ -0,0 +1,164 @@
+using namespace std;
+#include "DriverDAT_W_SMDS_Mesh.h"
+
+#include "SMDS_MeshElement.hxx"
+#include "SMDS_MeshNode.hxx"
+#include "SMDS_MeshEdgesIterator.hxx"
+#include "SMDS_MeshFacesIterator.hxx"
+#include "SMDS_MeshNodesIterator.hxx"
+#include "SMDS_MeshVolumesIterator.hxx"
+
+#include "utilities.h"
+
+DriverDAT_W_SMDS_Mesh::DriverDAT_W_SMDS_Mesh() {
+;
+}
+
+DriverDAT_W_SMDS_Mesh::~DriverDAT_W_SMDS_Mesh() {
+;
+}
+
+void DriverDAT_W_SMDS_Mesh::SetMesh(Handle(SMDS_Mesh)& aMesh) {
+  myMesh = aMesh;
+}
+
+void DriverDAT_W_SMDS_Mesh::SetFile(string aFile) {
+  myFile = aFile;
+}
+
+void DriverDAT_W_SMDS_Mesh::SetFileId(FILE* aFileId) {
+  myFileId = aFileId;
+}
+
+void DriverDAT_W_SMDS_Mesh::SetMeshId(int aMeshId) {
+  myMeshId = aMeshId;
+}
+
+void DriverDAT_W_SMDS_Mesh::Add() {
+  ;
+}
+
+void DriverDAT_W_SMDS_Mesh::Write() {
+
+  int nbNodes,nbCells;
+  int i;
+
+  char* file2Read = (char*)myFile.c_str();
+  myFileId = fopen(file2Read,"w+");
+  if (myFileId < 0)
+    {
+      fprintf(stderr,">> ERREUR : ouverture du fichier %s \n",file2Read);
+      exit(EXIT_FAILURE);
+    }
+  SCRUTE(myMesh);
+  /****************************************************************************
+  *                       NOMBRES D'OBJETS                                    *
+  ****************************************************************************/
+  fprintf(stdout,"\n(****************************)\n");
+  fprintf(stdout,"(* INFORMATIONS GENERALES : *)\n");
+  fprintf(stdout,"(****************************)\n");
+
+  /* Combien de noeuds ? */
+  nbNodes = myMesh->NbNodes();
+
+  /* Combien de mailles, faces ou aretes ? */
+  Standard_Integer nb_of_nodes, nb_of_edges,nb_of_faces, nb_of_volumes;
+  nb_of_edges = myMesh->NbEdges();
+  nb_of_faces = myMesh->NbFaces();
+  nb_of_volumes = myMesh->NbVolumes();
+  nbCells = nb_of_edges + nb_of_faces + nb_of_volumes;
+  SCRUTE(nb_of_edges);
+  SCRUTE(nb_of_faces);
+  SCRUTE(nb_of_volumes);
+
+  fprintf(stdout,"%d %d\n",nbNodes,nbCells);
+  fprintf(myFileId,"%d %d\n",nbNodes,nbCells);
+
+  /****************************************************************************
+  *                       ECRITURE DES NOEUDS                                 *
+  ****************************************************************************/
+  fprintf(stdout,"\n(************************)\n");
+  fprintf(stdout,"(* NOEUDS DU MAILLAGE : *)\n");
+  fprintf(stdout,"(************************)\n");
+
+  SMDS_MeshNodesIterator itNodes(myMesh);
+  for (;itNodes.More();itNodes.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itNodes.Value();
+    const Handle(SMDS_MeshNode)& node = myMesh->GetNode(1,elem);
+ 
+    fprintf(myFileId,"%d %e %e %e\n",node->GetID(),node->X(),node->Y(),node->Z());
+  }
+
+  /****************************************************************************
+  *                       ECRITURE DES ELEMENTS                                *
+  ****************************************************************************/
+  fprintf(stdout,"\n(**************************)\n");
+  fprintf(stdout,"(* ELEMENTS DU MAILLAGE : *)\n");
+  fprintf(stdout,"(**************************)");
+  /* Ecriture des connectivites, noms, numeros des mailles */
+
+  SMDS_MeshEdgesIterator itEdges(myMesh);
+  for (;itEdges.More();itEdges.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itEdges.Value();
+
+    switch (elem->NbNodes()) {
+    case 2 : {
+      fprintf(myFileId,"%d %d ",elem->GetID(),102);
+      break;
+    } 
+    case 3 : {
+      fprintf(myFileId,"%d %d ",elem->GetID(),103);
+      break;
+    }
+    }
+
+    for (i=0;i<elem->NbNodes();i++)
+      fprintf(myFileId,"%d ",elem->GetConnection(i+1));
+    
+    fprintf(myFileId,"\n");
+  }
+
+  SMDS_MeshFacesIterator itFaces(myMesh);
+  for (;itFaces.More();itFaces.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itFaces.Value();
+
+    switch (elem->NbNodes()) {
+    case 3 : {
+      fprintf(myFileId,"%d %d ",elem->GetID(),203);
+      break;
+    }
+    case 4 : {
+      fprintf(myFileId,"%d %d ",elem->GetID(),204);
+      break;
+    }
+    case 6 : {
+      fprintf(myFileId,"%d %d ",elem->GetID(),206);
+      break;
+    }
+    }
+
+    for (i=0;i<elem->NbNodes();i++)
+      fprintf(myFileId,"%d ",elem->GetConnection(i+1));
+    
+    fprintf(myFileId,"\n");
+  }
+
+  SMDS_MeshVolumesIterator itVolumes(myMesh);
+  for (;itVolumes.More();itVolumes.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itVolumes.Value();
+
+    switch (elem->NbNodes()) {
+    case 8 : {
+      fprintf(myFileId,"%d %d ",elem->GetID(),308);
+      break;
+    }
+    }
+
+    for (i=0;i<elem->NbNodes();i++)
+      fprintf(myFileId,"%d ",elem->GetConnection(i+1));
+    
+    fprintf(myFileId,"\n");
+  }
+
+  fclose (myFileId);
+}
diff --git a/src/DriverDAT/DriverDAT_W_SMDS_Mesh.h b/src/DriverDAT/DriverDAT_W_SMDS_Mesh.h

new file mode 100644 (file)

index 0000000..06056db
--- /dev/null
+++ b/src/DriverDAT/DriverDAT_W_SMDS_Mesh.h
@@ -0,0 +1,31 @@
+#ifndef _INCLUDE_DRIVERDAT_W_SMDS_MESH
+#define _INCLUDE_DRIVERDAT_W_SMDS_MESH
+
+#include <stdio.h>
+#include <string>
+
+#include "SMDS_Mesh.hxx"
+#include "Mesh_Writer.h"
+
+class DriverDAT_W_SMDS_Mesh : public Mesh_Writer {
+
+  public :
+    DriverDAT_W_SMDS_Mesh();
+    ~DriverDAT_W_SMDS_Mesh();
+
+    void Add();
+    void Write();
+    void SetMesh(Handle(SMDS_Mesh)& aMesh);
+    void SetFile(string);
+
+    void SetFileId(FILE*);
+    void SetMeshId(int);
+
+private :
+    Handle_SMDS_Mesh myMesh;
+    string myFile; 
+    FILE* myFileId;
+    int myMeshId;
+
+};
+#endif
diff --git a/src/DriverDAT/DriverDAT_W_SMESHDS_Document.cxx b/src/DriverDAT/DriverDAT_W_SMESHDS_Document.cxx

new file mode 100644 (file)

index 0000000..1cdb67b
--- /dev/null
+++ b/src/DriverDAT/DriverDAT_W_SMESHDS_Document.cxx
@@ -0,0 +1,78 @@
+using namespace std;
+#include "DriverDAT_W_SMESHDS_Document.h"
+#include "DriverDAT_W_SMESHDS_Mesh.h"
+
+#include "utilities.h"
+
+extern "C" 
+{
+  Document_Writer* Wmaker() {
+    return new DriverDAT_W_SMESHDS_Document;
+  }
+}
+
+DriverDAT_W_SMESHDS_Document::DriverDAT_W_SMESHDS_Document() {
+;
+}
+
+DriverDAT_W_SMESHDS_Document::~DriverDAT_W_SMESHDS_Document() {
+;
+}
+
+//void DriverDAT_W_SMESHDS_Document::SetFile(string aFile) {
+//myFile = aFile;
+//}
+
+//void DriverDAT_W_SMESHDS_Document::SetDocument(Handle(SMESHDS_Document)& aDocument) {
+//myDocument = aDocument;
+//}
+
+void DriverDAT_W_SMESHDS_Document::Write() {
+
+  Handle(SMESHDS_Mesh) myMesh;
+
+  /****************************************************************************
+  *                      OUVERTURE DU FICHIER EN ECRITURE                      *
+  ****************************************************************************/
+  char* file2Write = (char*)myFile.c_str();
+  FILE* fid = fopen(file2Write,"w+");
+  if (fid < 0)
+    {
+      fprintf(stderr,">> ERREUR : ouverture du fichier %s \n",file2Write);
+      exit(EXIT_FAILURE);
+    }
+
+  /****************************************************************************
+  *                      FERMETURE DU FICHIER                                 *
+  ****************************************************************************/
+
+  fclose(fid);
+
+  /******** Nombre de maillages ********/
+  int nb_of_meshes = myDocument->NbMeshes(); //voir avec Yves
+  //nb_of_meshes = 1;
+  int numero = 0;
+
+  string myClass = string("SMESHDS_Mesh");
+  string myExtension = string("DAT");
+
+  //while (numero<nb_of_meshes) {
+  //numero++;
+  //myMesh = myDocument->GetMesh(numero);
+  myDocument->InitMeshesIterator();
+  for (;myDocument->MoreMesh();myDocument->NextMesh()) {
+    numero++;
+    myMesh = myDocument->CurrentMesh();
+
+    DriverDAT_W_SMESHDS_Mesh* myWriter = new DriverDAT_W_SMESHDS_Mesh;
+    //Mesh_Writer* myWriter = Driver::GetMeshWriter(myExtension, myClass);
+
+    myWriter->SetMesh(myMesh);
+    myWriter->SetFile(myFile);
+  SCRUTE(myMesh);
+    //myWriter->SetFileId(fid);
+    myWriter->SetMeshId(numero);
+    myWriter->Write();
+  }
+
+}
diff --git a/src/DriverDAT/DriverDAT_W_SMESHDS_Document.h b/src/DriverDAT/DriverDAT_W_SMESHDS_Document.h

new file mode 100644 (file)

index 0000000..58df6cd
--- /dev/null
+++ b/src/DriverDAT/DriverDAT_W_SMESHDS_Document.h
@@ -0,0 +1,25 @@
+#ifndef _INCLUDE_DRIVERDAT_W_SMESHDS_DOCUMENT
+#define _INCLUDE_DRIVERDAT_W_SMESHDS_DOCUMENT
+
+#include <stdio.h>
+#include <string>
+
+#include "SMESHDS_Document.hxx"
+#include "Document_Writer.h"
+
+class DriverDAT_W_SMESHDS_Document : public Document_Writer {
+  
+public :
+  DriverDAT_W_SMESHDS_Document();
+  ~DriverDAT_W_SMESHDS_Document();
+ 
+  void Write();
+  //void SetFile(string);
+  //void SetDocument(Handle(SMESHDS_Document)&);
+
+private :
+  //Handle_SMESHDS_Document myDocument;
+  //string myFile; 
+
+};
+#endif
diff --git a/src/DriverDAT/DriverDAT_W_SMESHDS_Mesh.cxx b/src/DriverDAT/DriverDAT_W_SMESHDS_Mesh.cxx

new file mode 100644 (file)

index 0000000..9f20237
--- /dev/null
+++ b/src/DriverDAT/DriverDAT_W_SMESHDS_Mesh.cxx
@@ -0,0 +1,175 @@
+using namespace std;
+#include "DriverDAT_W_SMESHDS_Mesh.h"
+#include "DriverDAT_W_SMDS_Mesh.h"
+
+#include "SMDS_MeshElement.hxx"
+#include "SMDS_MeshNode.hxx"
+#include "SMDS_MeshEdgesIterator.hxx"
+#include "SMDS_MeshFacesIterator.hxx"
+#include "SMDS_MeshNodesIterator.hxx"
+#include "SMDS_MeshVolumesIterator.hxx"
+
+#include "utilities.h"
+
+DriverDAT_W_SMESHDS_Mesh::DriverDAT_W_SMESHDS_Mesh() {
+;
+}
+
+DriverDAT_W_SMESHDS_Mesh::~DriverDAT_W_SMESHDS_Mesh() {
+;
+}
+
+void DriverDAT_W_SMESHDS_Mesh::SetMesh(Handle(SMDS_Mesh)& aMesh) {
+  //myMesh = Handle(SMESHDS_Mesh)::DownCast(aMesh);
+  myMesh = aMesh;
+}
+
+void DriverDAT_W_SMESHDS_Mesh::SetFile(string aFile) {
+  myFile = aFile;
+}
+
+void DriverDAT_W_SMESHDS_Mesh::SetFileId(FILE* aFileId) {
+  myFileId = aFileId;
+}
+
+void DriverDAT_W_SMESHDS_Mesh::SetMeshId(int aMeshId) {
+  myMeshId = aMeshId;
+}
+
+void DriverDAT_W_SMESHDS_Mesh::Write() {
+
+  string myClass = string("SMDS_Mesh");
+  string myExtension = string("DAT");
+
+  DriverDAT_W_SMDS_Mesh* myWriter = new DriverDAT_W_SMDS_Mesh;
+
+  myWriter->SetMesh(myMesh);
+  myWriter->SetFile(myFile);
+  myWriter->SetMeshId(myMeshId);
+  //myWriter->SetFileId(myFileId);
+
+  myWriter->Write();
+
+
+}
+void DriverDAT_W_SMESHDS_Mesh::Add() {
+  int nbNodes,nbCells;
+  int i;
+
+  char* file2Read = (char*)myFile.c_str();
+  myFileId = fopen(file2Read,"w+");
+  if (myFileId < 0)
+    {
+      fprintf(stderr,">> ERREUR : ouverture du fichier %s \n",file2Read);
+      exit(EXIT_FAILURE);
+    }
+  
+  /****************************************************************************
+  *                       NOMBRES D'OBJETS                                    *
+  ****************************************************************************/
+  fprintf(stdout,"\n(****************************)\n");
+  fprintf(stdout,"(* INFORMATIONS GENERALES : *)\n");
+  fprintf(stdout,"(****************************)\n");
+
+  /* Combien de noeuds ? */
+  nbNodes = myMesh->NbNodes();
+
+  /* Combien de mailles, faces ou aretes ? */
+  Standard_Integer nb_of_nodes, nb_of_edges,nb_of_faces, nb_of_volumes;
+  nb_of_edges = myMesh->NbEdges();
+  nb_of_faces = myMesh->NbFaces();
+  nb_of_volumes = myMesh->NbVolumes();
+  nbCells = nb_of_edges + nb_of_faces + nb_of_volumes;
+
+  fprintf(stdout,"%d %d\n",nbNodes,nbCells);
+  fprintf(myFileId,"%d %d\n",nbNodes,nbCells);
+
+  /****************************************************************************
+  *                       ECRITURE DES NOEUDS                                 *
+  ****************************************************************************/
+  fprintf(stdout,"\n(************************)\n");
+  fprintf(stdout,"(* NOEUDS DU MAILLAGE : *)\n");
+  fprintf(stdout,"(************************)\n");
+
+  SMDS_MeshNodesIterator itNodes(myMesh);
+  for (;itNodes.More();itNodes.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itNodes.Value();
+    const Handle(SMDS_MeshNode)& node = myMesh->GetNode(1,elem);
+ 
+    fprintf(myFileId,"%d %e %e %e\n",node->GetID(),node->X(),node->Y(),node->Z());
+  }
+
+  /****************************************************************************
+  *                       ECRITURE DES ELEMENTS                                *
+  ****************************************************************************/
+  fprintf(stdout,"\n(**************************)\n");
+  fprintf(stdout,"(* ELEMENTS DU MAILLAGE : *)\n");
+  fprintf(stdout,"(**************************)");
+  /* Ecriture des connectivites, noms, numeros des mailles */
+
+  SMDS_MeshEdgesIterator itEdges(myMesh);
+  for (;itEdges.More();itEdges.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itEdges.Value();
+
+    switch (elem->NbNodes()) {
+    case 2 : {
+      fprintf(myFileId,"%d %d ",elem->GetID(),102);
+      break;
+    } 
+    case 3 : {
+      fprintf(myFileId,"%d %d ",elem->GetID(),103);
+      break;
+    }
+    }
+
+    for (i=0;i<elem->NbNodes();i++)
+      fprintf(myFileId,"%d ",elem->GetConnection(i+1));
+    
+    fprintf(myFileId,"\n");
+  }
+
+  SMDS_MeshFacesIterator itFaces(myMesh);
+  for (;itFaces.More();itFaces.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itFaces.Value();
+
+    switch (elem->NbNodes()) {
+    case 3 : {
+      fprintf(myFileId,"%d %d ",elem->GetID(),203);
+      break;
+    }
+    case 4 : {
+      fprintf(myFileId,"%d %d ",elem->GetID(),204);
+      break;
+    }
+    case 6 : {
+      fprintf(myFileId,"%d %d ",elem->GetID(),206);
+      break;
+    }
+    }
+
+    for (i=0;i<elem->NbNodes();i++)
+      fprintf(myFileId,"%d ",elem->GetConnection(i+1));
+    
+    fprintf(myFileId,"\n");
+  }
+
+  SMDS_MeshVolumesIterator itVolumes(myMesh);
+  for (;itVolumes.More();itVolumes.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itVolumes.Value();
+
+    switch (elem->NbNodes()) {
+    case 8 : {
+      fprintf(myFileId,"%d %d ",elem->GetID(),308);
+      break;
+    }
+    }
+
+    for (i=0;i<elem->NbNodes();i++)
+      fprintf(myFileId,"%d ",elem->GetConnection(i+1));
+    
+    fprintf(myFileId,"\n");
+  }
+
+  fclose (myFileId);
+}
+
diff --git a/src/DriverDAT/DriverDAT_W_SMESHDS_Mesh.h b/src/DriverDAT/DriverDAT_W_SMESHDS_Mesh.h

new file mode 100644 (file)

index 0000000..fe7e02d
--- /dev/null
+++ b/src/DriverDAT/DriverDAT_W_SMESHDS_Mesh.h
@@ -0,0 +1,31 @@
+#ifndef _INCLUDE_DRIVERDAT_W_SMESHDS_MESH
+#define _INCLUDE_DRIVERDAT_W_SMESHDS_MESH
+
+#include <stdio.h>
+#include <string>
+
+#include "SMESHDS_Mesh.hxx"
+#include "Mesh_Writer.h"
+
+class DriverDAT_W_SMESHDS_Mesh : public Mesh_Writer {
+
+  public :
+    DriverDAT_W_SMESHDS_Mesh();
+    ~DriverDAT_W_SMESHDS_Mesh();
+
+    void Add();
+    void Write();
+    void SetMesh(Handle(SMDS_Mesh)& aMesh);
+    void SetFile(string);
+
+    void SetFileId(FILE*);
+    void SetMeshId(int);
+
+private :
+    Handle_SMDS_Mesh myMesh;
+    string myFile; 
+    FILE* myFileId;
+    int myMeshId;
+
+};
+#endif
diff --git a/src/DriverDAT/Makefile.in b/src/DriverDAT/Makefile.in

new file mode 100644 (file)

index 0000000..a4c404d
--- /dev/null
+++ b/src/DriverDAT/Makefile.in
@@ -0,0 +1,39 @@
+# -* Makefile *- 
+#
+# Author : Marc Tajchman (CEA)
+# Date : 5/07/2001
+# $Header$
+#
+
+# source path
+top_srcdir=@top_srcdir@
+top_builddir=../..
+srcdir=@srcdir@
+VPATH=.:@srcdir@
+
+
+@COMMENCE@
+
+# header files 
+EXPORT_HEADERS= DriverDAT_R_SMDS_Mesh.h DriverDAT_R_SMESHDS_Mesh.h DriverDAT_R_SMESHDS_Document.h DriverDAT_W_SMDS_Mesh.h DriverDAT_W_SMESHDS_Mesh.h DriverDAT_W_SMESHDS_Document.h 
+
+# Libraries targets
+LIB = libMeshDriverDAT.la
+LIB_SRC = DriverDAT_R_SMDS_Mesh.cxx DriverDAT_R_SMESHDS_Mesh.cxx DriverDAT_R_SMESHDS_Document.cxx DriverDAT_W_SMDS_Mesh.cxx DriverDAT_W_SMESHDS_Mesh.cxx DriverDAT_W_SMESHDS_Document.cxx 
+
+LIB_CLIENT_IDL =
+
+LIB_SERVER_IDL =
+
+# additionnal information to compil and link file
+CPPFLAGS += $(OCC_INCLUDES)
+CXXFLAGS += $(OCC_CXXFLAGS) $(MED2_INCLUDES)
+LDFLAGS  += $(OCC_LIBS) $(MED2_LIBS) -lMeshDriver
+
+%_moc.cxx: %.h
+       $(MOC) $< -o $@
+
+@CONCLUDE@
+
+
+
diff --git a/src/DriverUNV/DriverUNV_R_SMDS_Mesh.cxx b/src/DriverUNV/DriverUNV_R_SMDS_Mesh.cxx

new file mode 100644 (file)

index 0000000..b9de018
--- /dev/null
+++ b/src/DriverUNV/DriverUNV_R_SMDS_Mesh.cxx
@@ -0,0 +1,153 @@
+using namespace std;
+#include "DriverUNV_R_SMDS_Mesh.h"
+
+#include "utilities.h"
+
+DriverUNV_R_SMDS_Mesh::DriverUNV_R_SMDS_Mesh() {
+;
+}
+
+DriverUNV_R_SMDS_Mesh::~DriverUNV_R_SMDS_Mesh() {
+;
+}
+
+void DriverUNV_R_SMDS_Mesh::SetMesh(Handle(SMDS_Mesh)& aMesh) {
+  myMesh = aMesh;
+}
+
+void DriverUNV_R_SMDS_Mesh::SetFile(string aFile) {
+  myFile = aFile;
+}
+
+void DriverUNV_R_SMDS_Mesh::SetFileId(FILE* aFileId) {
+  myFileId = aFileId;
+}
+
+void DriverUNV_R_SMDS_Mesh::SetMeshId(int aMeshId) {
+  myMeshId = aMeshId;
+}
+
+void DriverUNV_R_SMDS_Mesh::Add() {
+  ;
+}
+
+void DriverUNV_R_SMDS_Mesh::Read() {
+
+  int cell=0,node=0,n1,n2,n3,n4,n_nodes,nodes[6],blockId,i;
+  char *s1,*s2,*s3;
+  string str1,str2,str3;
+  int i1=0;
+  bool ok, found_block2411, found_block2412;
+
+  /****************************************************************************
+  *                      OUVERTURE DU FICHIER EN LECTURE                      *
+  ****************************************************************************/
+  char* file2Read = (char*)myFile.c_str();
+  myFileId = fopen(file2Read,"r");
+  if (myFileId < 0)
+    {
+      fprintf(stderr,">> ERREUR : ouverture du fichier %s \n",file2Read);
+      exit(EXIT_FAILURE);
+    }
+
+  s1 = (char*) malloc(sizeof(char)*100);
+  s2 = (char*) malloc(sizeof(char)*100);
+  s3 = (char*) malloc(sizeof(char)*100);
+
+  found_block2411 = false;
+  found_block2412 = false;
+
+  do {
+
+    while (i1==-1) {
+      fscanf(myFileId,"%d\n",&blockId);
+      switch (blockId) {
+      case 2411 : {
+      MESSAGE("BlockId "<<blockId);
+       fscanf(myFileId,"%d",&node);
+       //MESSAGE("Node "<<node);
+       while (node!=-1) {
+         fscanf(myFileId,"%d %d %d\n",&n1,&n2,&n3);
+         fscanf(myFileId,"%s %s %s\n",s1,s2,s3);
+         str1=string(s1);
+         str2=string(s2);
+         str3=string(s3);
+         if (str1.find("D")!=string::npos) str1.replace(str1.find("D"),1,"E");
+         if (str2.find("D")!=string::npos) str2.replace(str2.find("D"),1,"E");
+         if (str3.find("D")!=string::npos) str3.replace(str3.find("D"),1,"E");
+         ok = myMesh->AddNodeWithID(atof(str1.c_str()),atof(str2.c_str()),atof(str3.c_str()),node);
+         fscanf(myFileId,"%d",&node);
+       }
+       i1=0;
+       found_block2411 = true;
+       break;
+      }        
+      case 2412 : {
+      MESSAGE("BlockId "<<blockId);
+       fscanf(myFileId,"%d",&cell);
+       //MESSAGE("Cell "<<cell);
+       while (cell!=-1) {
+         fscanf(myFileId,"%d %d %d %d %d\n",&n1,&n2,&n3,&n4,&n_nodes);
+
+         if ((n1==71)||(n1==72)||(n1==74)||(n1==91)||(n1==92)) {//203
+           if (n_nodes==3) {
+             for (i=1;i<=n_nodes;i++)
+               fscanf(myFileId,"%d",&nodes[i-1]);
+             ok = myMesh->AddFaceWithID(nodes[0],nodes[1],nodes[2],cell);
+           }
+           else if (n_nodes==6) {//206
+             for (i=1;i<=n_nodes;i++)
+               fscanf(myFileId,"%d",&nodes[i-1]);
+             ok = myMesh->AddFaceWithID(nodes[0],nodes[2],nodes[4],cell);
+           }
+         }
+
+         else if ((n1==11)||(n1==21)||(n1==24)||(n1==25)) {//103
+           fgets(s2,100,myFileId);
+           if (n_nodes==3) {
+             for (i=1;i<=n_nodes;i++)
+               fscanf(myFileId,"%d",&nodes[i-1]);
+             ok = myMesh->AddEdgeWithID(nodes[0],nodes[1],cell);
+             //MESSAGE("in 103 "<<cell);
+           }
+           else if (n_nodes==2) {//102
+             for (i=1;i<=n_nodes;i++)
+               fscanf(myFileId,"%d",&nodes[i-1]);
+             ok = myMesh->AddEdgeWithID(nodes[0],nodes[1],cell);
+             //MESSAGE("in 102 "<<cell);
+           }
+         }
+
+         fscanf(myFileId,"\n");
+         fscanf(myFileId,"%d",&cell);
+       }
+
+       i1=0;
+       found_block2412 = true;
+       break;
+      }
+      case -1 : {
+       break;
+      }
+      default:
+       MESSAGE("BlockId "<<blockId);
+       i1=0;
+       break;
+      }
+    }
+  
+    fscanf(myFileId,"%s\n",s1);
+    i1 = atoi(s1);
+
+  }
+  while ((!feof(myFileId))&&((!found_block2411)||(!found_block2412)));
+
+  /****************************************************************************
+  *                      FERMETURE DU FICHIER                      *
+  ****************************************************************************/
+  free(s1);
+  free(s2);
+  free(s3);
+  fclose(myFileId);
+
+}
diff --git a/src/DriverUNV/DriverUNV_R_SMDS_Mesh.h b/src/DriverUNV/DriverUNV_R_SMDS_Mesh.h

new file mode 100644 (file)

index 0000000..c612677
--- /dev/null
+++ b/src/DriverUNV/DriverUNV_R_SMDS_Mesh.h
@@ -0,0 +1,30 @@
+#ifndef _INCLUDE_DRIVERUNV_R_SMDS_MESH
+#define _INCLUDE_DRIVERUNV_R_SMDS_MESH
+
+#include <stdio.h>
+
+#include "SMDS_Mesh.hxx"
+#include "Mesh_Reader.h"
+
+class DriverUNV_R_SMDS_Mesh : public Mesh_Reader {
+
+  public :
+    DriverUNV_R_SMDS_Mesh();
+    ~DriverUNV_R_SMDS_Mesh();
+
+    void Add();
+    void Read();
+    void SetMesh(Handle(SMDS_Mesh)& aMesh);
+    void SetFile(string);
+
+    void SetFileId(FILE*);
+    void SetMeshId(int);
+
+  private :
+    Handle_SMDS_Mesh myMesh;
+    string myFile; 
+    FILE* myFileId;
+    int myMeshId;
+
+};
+#endif
diff --git a/src/DriverUNV/DriverUNV_R_SMESHDS_Document.cxx b/src/DriverUNV/DriverUNV_R_SMESHDS_Document.cxx

new file mode 100644 (file)

index 0000000..8debd49
--- /dev/null
+++ b/src/DriverUNV/DriverUNV_R_SMESHDS_Document.cxx
@@ -0,0 +1,85 @@
+using namespace std;
+#include "DriverUNV_R_SMESHDS_Document.h"
+#include "DriverUNV_R_SMESHDS_Mesh.h"
+
+#include "utilities.h"
+
+int getOne() {
+  printf("in getOne");
+  return (1);
+}
+
+extern "C" {
+  //  Document_Reader* maker() {
+  DriverUNV_R_SMESHDS_Document* maker() {
+    fprintf(stdout,"here in maker\n");
+    return new DriverUNV_R_SMESHDS_Document;
+  }
+}
+
+DriverUNV_R_SMESHDS_Document::DriverUNV_R_SMESHDS_Document() {
+  myFile = string("");
+}
+
+DriverUNV_R_SMESHDS_Document::~DriverUNV_R_SMESHDS_Document() {
+;
+}
+
+//void DriverUNV_R_SMESHDS_Document::SetFile(string aFile) {
+//myFile = aFile;
+//}
+
+//void DriverUNV_R_SMESHDS_Document::SetDocument(Handle(SMESHDS_Document)& aDoc) {
+//myDocument = aDoc;
+//}
+
+void DriverUNV_R_SMESHDS_Document::Read() {
+
+  int myMeshId;
+  SCRUTE(myFile);
+  //Handle(SMESHDS_Document) myDocument = new SMESHDS_Document(1);
+
+  /****************************************************************************
+  *                      OUVERTURE DU FICHIER EN LECTURE                      *
+  ****************************************************************************/
+  char* file2Read = (char*)myFile.c_str();
+  FILE* fid = fopen(file2Read,"r");
+  if (fid < 0)
+    {
+      fprintf(stderr,">> ERREUR : ouverture du fichier %s \n",file2Read);
+      exit(EXIT_FAILURE);
+    }
+
+  /****************************************************************************
+  *                      COMBIEN DE MAILLAGES ?                               *
+  ****************************************************************************/
+  int nmaa = 1;
+
+  /****************************************************************************
+  *                      FERMETURE DU FICHIER                                 *
+  ****************************************************************************/
+  fclose(fid);
+  
+  printf("Nombre de maillages = %d\n",nmaa);
+
+  string myClass = string("SMESHDS_Mesh");
+  string myExtension = string("UNV");
+
+  for (int meshIt=1;meshIt<=nmaa;meshIt++) {
+    myMeshId = myDocument->NewMesh();
+
+    Handle(SMDS_Mesh) myMesh = myDocument->GetMesh(myMeshId);
+
+    DriverUNV_R_SMESHDS_Mesh* myReader = new DriverUNV_R_SMESHDS_Mesh;
+
+
+    myReader->SetMesh(myMesh);
+    myReader->SetFile(myFile);
+    //myReader->SetFileId(fid);
+
+    myReader->Read();
+
+  }
+
+
+}
diff --git a/src/DriverUNV/DriverUNV_R_SMESHDS_Document.h b/src/DriverUNV/DriverUNV_R_SMESHDS_Document.h

new file mode 100644 (file)

index 0000000..9fbd04f
--- /dev/null
+++ b/src/DriverUNV/DriverUNV_R_SMESHDS_Document.h
@@ -0,0 +1,24 @@
+#ifndef _INCLUDE_DRIVERUNV_R_SMESHDS_DOCUMENT
+#define _INCLUDE_DRIVERUNV_R_SMESHDS_DOCUMENT
+
+#include <stdio.h>
+
+#include "SMESHDS_Document.hxx"
+#include "Document_Reader.h"
+
+class DriverUNV_R_SMESHDS_Document : public Document_Reader {
+  
+public :
+  DriverUNV_R_SMESHDS_Document();
+  ~DriverUNV_R_SMESHDS_Document();
+ 
+  void Read();
+  //void SetFile(string);
+  //void SetDocument(Handle_SMESHDS_Document&);
+
+private :
+  //Handle_SMESHDS_Document myDocument;
+  //string myFile; 
+
+};
+#endif
diff --git a/src/DriverUNV/DriverUNV_R_SMESHDS_Mesh.cxx b/src/DriverUNV/DriverUNV_R_SMESHDS_Mesh.cxx

new file mode 100644 (file)

index 0000000..849a523
--- /dev/null
+++ b/src/DriverUNV/DriverUNV_R_SMESHDS_Mesh.cxx
@@ -0,0 +1,48 @@
+using namespace std;
+#include "DriverUNV_R_SMESHDS_Mesh.h"
+#include "DriverUNV_R_SMDS_Mesh.h"
+
+#include "utilities.h"
+
+DriverUNV_R_SMESHDS_Mesh::DriverUNV_R_SMESHDS_Mesh() {
+;
+}
+
+DriverUNV_R_SMESHDS_Mesh::~DriverUNV_R_SMESHDS_Mesh() {
+;
+}
+
+void DriverUNV_R_SMESHDS_Mesh::SetMesh(Handle(SMDS_Mesh)& aMesh) {
+  //myMesh = Handle(SMESHDS_Mesh)::DownCast(aMesh);
+  myMesh = aMesh;
+}
+
+void DriverUNV_R_SMESHDS_Mesh::SetFile(string aFile) {
+  myFile = aFile;
+}
+
+void DriverUNV_R_SMESHDS_Mesh::SetFileId(FILE* aFileId) {
+  myFileId = aFileId;
+}
+
+void DriverUNV_R_SMESHDS_Mesh::SetMeshId(int aMeshId) {
+  myMeshId = aMeshId;
+}
+
+void DriverUNV_R_SMESHDS_Mesh::Add() {
+  ;
+}
+
+void DriverUNV_R_SMESHDS_Mesh::Read() {
+  string myClass = string("SMDS_Mesh");
+  string myExtension = string("UNV");
+
+  DriverUNV_R_SMDS_Mesh* myReader = new DriverUNV_R_SMDS_Mesh;
+
+  myReader->SetMesh(myMesh);
+  myReader->SetFile(myFile);
+  //myReader->SetFileId(myFileId);
+
+  myReader->Read();
+
+}
diff --git a/src/DriverUNV/DriverUNV_R_SMESHDS_Mesh.h b/src/DriverUNV/DriverUNV_R_SMESHDS_Mesh.h

new file mode 100644 (file)

index 0000000..0c3a7e6
--- /dev/null
+++ b/src/DriverUNV/DriverUNV_R_SMESHDS_Mesh.h
@@ -0,0 +1,30 @@
+#ifndef _INCLUDE_DRIVERUNV_R_SMESHDS_MESH
+#define _INCLUDE_DRIVERUNV_R_SMESHDS_MESH
+
+#include <stdio.h>
+
+#include "SMESHDS_Mesh.hxx"
+#include "Mesh_Reader.h"
+
+class DriverUNV_R_SMESHDS_Mesh : public Mesh_Reader {
+
+  public :
+    DriverUNV_R_SMESHDS_Mesh();
+    ~DriverUNV_R_SMESHDS_Mesh();
+
+    void Add();
+    void Read();
+    void SetMesh(Handle(SMDS_Mesh)& aMesh);
+    void SetFile(string);
+
+    void SetFileId(FILE*);
+    void SetMeshId(int);
+
+private :
+    Handle_SMDS_Mesh myMesh;
+    string myFile; 
+    FILE* myFileId;
+    int myMeshId;
+
+};
+#endif
diff --git a/src/DriverUNV/DriverUNV_W_SMDS_Mesh.cxx b/src/DriverUNV/DriverUNV_W_SMDS_Mesh.cxx

new file mode 100644 (file)

index 0000000..6bf42e6
--- /dev/null
+++ b/src/DriverUNV/DriverUNV_W_SMDS_Mesh.cxx
@@ -0,0 +1,196 @@
+using namespace std;
+#include "DriverUNV_W_SMDS_Mesh.h"
+
+#include "SMDS_MeshElement.hxx"
+#include "SMDS_MeshNode.hxx"
+#include "SMDS_MeshEdgesIterator.hxx"
+#include "SMDS_MeshFacesIterator.hxx"
+#include "SMDS_MeshNodesIterator.hxx"
+#include "SMDS_MeshVolumesIterator.hxx"
+
+#include <utilities.h>
+
+#define sNODE_UNV_ID "  2411"
+#define sELT_UNV_ID  "  2412"
+#define sUNV_SEPARATOR "    -1"
+#define sNODE_UNV_DESCR "%10d         1         1        11\n"
+#define sELT_SURF_DESC  "%10d        %2d         1         1        11         %1d\n"
+#define sELT_VOLU_DESC  "%10d        %2d         1         1         9         %1d\n"
+#define sELT_BEAM_DESC1 "%10d        %2d         1         1         7         %1d\n"
+#define sELT_BEAM_DESC2 "         0         1         1\n"
+
+DriverUNV_W_SMDS_Mesh::DriverUNV_W_SMDS_Mesh() {
+;
+}
+
+DriverUNV_W_SMDS_Mesh::~DriverUNV_W_SMDS_Mesh() {
+;
+}
+
+void DriverUNV_W_SMDS_Mesh::SetMesh(Handle(SMDS_Mesh)& aMesh) {
+  myMesh = aMesh;
+}
+
+void DriverUNV_W_SMDS_Mesh::SetFile(string aFile) {
+  myFile = aFile;
+}
+
+void DriverUNV_W_SMDS_Mesh::SetFileId(FILE* aFileId) {
+  myFileId = aFileId;
+}
+
+void DriverUNV_W_SMDS_Mesh::SetMeshId(int aMeshId) {
+  myMeshId = aMeshId;
+}
+
+void DriverUNV_W_SMDS_Mesh::Add() {
+  ;
+}
+
+void DriverUNV_W_SMDS_Mesh::Write() {
+
+  int nbNodes,nbCells;
+  int i;
+
+  char* file2Read = (char*)myFile.c_str();
+  myFileId = fopen(file2Read,"w+");
+  if (myFileId < 0)
+    {
+      fprintf(stderr,">> ERREUR : ouverture du fichier %s \n",file2Read);
+      exit(EXIT_FAILURE);
+    }
+  SCRUTE(myMesh);
+  /****************************************************************************
+  *                       NOMBRES D'OBJETS                                    *
+  ****************************************************************************/
+  fprintf(stdout,"\n(****************************)\n");
+  fprintf(stdout,"(* INFORMATIONS GENERALES : *)\n");
+  fprintf(stdout,"(****************************)\n");
+
+  /* Combien de noeuds ? */
+  nbNodes = myMesh->NbNodes();
+
+  /* Combien de mailles, faces ou aretes ? */
+  Standard_Integer nb_of_nodes, nb_of_edges,nb_of_faces, nb_of_volumes;
+  nb_of_edges = myMesh->NbEdges();
+  nb_of_faces = myMesh->NbFaces();
+  nb_of_volumes = myMesh->NbVolumes();
+  nbCells = nb_of_edges + nb_of_faces + nb_of_volumes;
+  SCRUTE(nb_of_edges);
+  SCRUTE(nb_of_faces);
+  SCRUTE(nb_of_volumes);
+
+  fprintf(stdout,"%d %d\n",nbNodes,nbCells);
+  fprintf(myFileId,"%d %d\n",nbNodes,nbCells);
+
+  /****************************************************************************
+  *                       ECRITURE DES NOEUDS                                 *
+  ****************************************************************************/
+  fprintf(stdout,"\n(************************)\n");
+  fprintf(stdout,"(* NOEUDS DU MAILLAGE : *)\n");
+  fprintf(stdout,"(************************)\n");
+
+  SMDS_MeshNodesIterator itNodes(myMesh);
+
+  fprintf(myFileId,"%s\n", sUNV_SEPARATOR);
+  fprintf(myFileId,"%s\n", sNODE_UNV_ID  );
+
+  for (;itNodes.More();itNodes.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itNodes.Value();
+    const Handle(SMDS_MeshNode   )& node = myMesh->GetNode(1, elem);
+
+    fprintf(myFileId, sNODE_UNV_DESCR, node->GetID());
+    fprintf(myFileId, "%25.16E%25.16E%25.16E\n", node->X(), node->Y(), node->Z());
+  }
+  fprintf(myFileId,"%s\n", sUNV_SEPARATOR);
+
+  /****************************************************************************
+  *                       ECRITURE DES ELEMENTS                                *
+  ****************************************************************************/
+  fprintf(stdout,"\n(**************************)\n");
+  fprintf(stdout,"(* ELEMENTS DU MAILLAGE : *)\n");
+  fprintf(stdout,"(**************************)");
+  /* Ecriture des connectivites, noms, numeros des mailles */
+
+  fprintf(myFileId,"%s\n", sUNV_SEPARATOR);
+  fprintf(myFileId,"%s\n", sELT_UNV_ID   );
+
+  SMDS_MeshEdgesIterator itEdges(myMesh);
+  for (;itEdges.More();itEdges.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itEdges.Value();
+
+    switch (elem->NbNodes()) {
+      case 2 : {
+        fprintf(myFileId, sELT_BEAM_DESC1, elem->GetID(), 21, elem->NbNodes());
+        fprintf(myFileId, sELT_BEAM_DESC2);
+        fprintf(myFileId, "%10d%10d\n", elem->GetConnection(1), elem->GetConnection(2));
+        break;
+      } 
+      case 3 : {
+        fprintf(myFileId, sELT_BEAM_DESC1, elem->GetID(), 24, elem->NbNodes());
+        fprintf(myFileId, sELT_BEAM_DESC2);
+        fprintf(myFileId, "%10d%10d%10d\n",elem->GetConnection(1), elem->GetConnection(2), elem->GetConnection(3));
+        break;
+      }
+    }
+  }
+
+  SMDS_MeshFacesIterator itFaces(myMesh);
+  for (;itFaces.More();itFaces.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itFaces.Value();
+
+    switch (elem->NbNodes()) {
+      case 3 :
+        // linear triangle
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 74, elem->NbNodes());
+        break;
+      case 4 :
+        // linear quadrilateral
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 71, elem->NbNodes());
+        break;
+      case 6 :
+        // parabolic triangle
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 72, elem->NbNodes());
+        break;
+      case 8 :
+        // parabolic quadrilateral
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 75, elem->NbNodes());
+        break;
+      default:
+        fprintf(myFileId, "element not registered\n");
+    }
+
+    for (i=0;i<elem->NbNodes();i++)
+      fprintf(myFileId,"%10d",elem->GetConnection(i+1));
+    
+    fprintf(myFileId,"\n");
+  }
+
+  SMDS_MeshVolumesIterator itVolumes(myMesh);
+  for (;itVolumes.More();itVolumes.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itVolumes.Value();
+
+    switch (elem->NbNodes()) {
+      case 4 : 
+        // linear tetrahedron
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 111, elem->NbNodes());
+        break;
+      case 6 : 
+        // linear tetrahedron
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 112, elem->NbNodes());
+        break;
+      case 8 : 
+        // linear brick
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 115, elem->NbNodes());
+        break;
+    }
+
+    for (i=0;i<elem->NbNodes();i++)
+      fprintf(myFileId,"%10d",elem->GetConnection(i+1));
+    
+    fprintf(myFileId,"\n");
+  }
+  fprintf(myFileId,"%s\n", sUNV_SEPARATOR);
+
+  fclose (myFileId);
+}
diff --git a/src/DriverUNV/DriverUNV_W_SMDS_Mesh.h b/src/DriverUNV/DriverUNV_W_SMDS_Mesh.h

new file mode 100644 (file)

index 0000000..141ab21
--- /dev/null
+++ b/src/DriverUNV/DriverUNV_W_SMDS_Mesh.h
@@ -0,0 +1,31 @@
+#ifndef _INCLUDE_DRIVERUNV_W_SMDS_MESH
+#define _INCLUDE_DRIVERUNV_W_SMDS_MESH
+
+#include <stdio.h>
+#include <string>
+
+#include "SMDS_Mesh.hxx"
+#include "Mesh_Writer.h"
+
+class DriverUNV_W_SMDS_Mesh : public Mesh_Writer {
+
+  public :
+    DriverUNV_W_SMDS_Mesh();
+    ~DriverUNV_W_SMDS_Mesh();
+
+    void Add();
+    void Write();
+    void SetMesh(Handle(SMDS_Mesh)& aMesh);
+    void SetFile(string);
+
+    void SetFileId(FILE*);
+    void SetMeshId(int);
+
+private :
+    Handle_SMDS_Mesh myMesh;
+    string myFile; 
+    FILE* myFileId;
+    int myMeshId;
+
+};
+#endif
diff --git a/src/DriverUNV/DriverUNV_W_SMESHDS_Document.cxx b/src/DriverUNV/DriverUNV_W_SMESHDS_Document.cxx

new file mode 100644 (file)

index 0000000..7fb740c
--- /dev/null
+++ b/src/DriverUNV/DriverUNV_W_SMESHDS_Document.cxx
@@ -0,0 +1,78 @@
+using namespace std;
+#include "DriverUNV_W_SMESHDS_Document.h"
+#include "DriverUNV_W_SMESHDS_Mesh.h"
+
+#include "utilities.h"
+
+extern "C" 
+{
+  Document_Writer* Wmaker() {
+    return new DriverUNV_W_SMESHDS_Document;
+  }
+}
+
+DriverUNV_W_SMESHDS_Document::DriverUNV_W_SMESHDS_Document() {
+;
+}
+
+DriverUNV_W_SMESHDS_Document::~DriverUNV_W_SMESHDS_Document() {
+;
+}
+
+//void DriverUNV_W_SMESHDS_Document::SetFile(string aFile) {
+//myFile = aFile;
+//}
+
+//void DriverUNV_W_SMESHDS_Document::SetDocument(Handle(SMESHDS_Document)& aDocument) {
+//myDocument = aDocument;
+//}
+
+void DriverUNV_W_SMESHDS_Document::Write() {
+
+  Handle(SMESHDS_Mesh) myMesh;
+
+  /****************************************************************************
+  *                      OUVERTURE DU FICHIER EN ECRITURE                      *
+  ****************************************************************************/
+  char* file2Write = (char*)myFile.c_str();
+  FILE* fid = fopen(file2Write,"w+");
+  if (fid < 0)
+    {
+      fprintf(stderr,">> ERREUR : ouverture du fichier %s \n",file2Write);
+      exit(EXIT_FAILURE);
+    }
+
+  /****************************************************************************
+  *                      FERMETURE DU FICHIER                                 *
+  ****************************************************************************/
+
+  fclose(fid);
+
+  /******** Nombre de maillages ********/
+  int nb_of_meshes = myDocument->NbMeshes(); //voir avec Yves
+  //nb_of_meshes = 1;
+  int numero = 0;
+
+  string myClass = string("SMESHDS_Mesh");
+  string myExtension = string("UNV");
+
+  //while (numero<nb_of_meshes) {
+  //numero++;
+  //myMesh = myDocument->GetMesh(numero);
+  myDocument->InitMeshesIterator();
+  for (;myDocument->MoreMesh();myDocument->NextMesh()) {
+    numero++;
+    myMesh = myDocument->CurrentMesh();
+
+    DriverUNV_W_SMESHDS_Mesh* myWriter = new DriverUNV_W_SMESHDS_Mesh;
+    //Mesh_Writer* myWriter = Driver::GetMeshWriter(myExtension, myClass);
+
+    myWriter->SetMesh(myMesh);
+    myWriter->SetFile(myFile);
+  SCRUTE(myMesh);
+    //myWriter->SetFileId(fid);
+    myWriter->SetMeshId(numero);
+    myWriter->Write();
+  }
+
+}
diff --git a/src/DriverUNV/DriverUNV_W_SMESHDS_Document.h b/src/DriverUNV/DriverUNV_W_SMESHDS_Document.h

new file mode 100644 (file)

index 0000000..7b71287
--- /dev/null
+++ b/src/DriverUNV/DriverUNV_W_SMESHDS_Document.h
@@ -0,0 +1,25 @@
+#ifndef _INCLUDE_DRIVERUNV_W_SMESHDS_DOCUMENT
+#define _INCLUDE_DRIVERUNV_W_SMESHDS_DOCUMENT
+
+#include <stdio.h>
+#include <string>
+
+#include "SMESHDS_Document.hxx"
+#include "Document_Writer.h"
+
+class DriverUNV_W_SMESHDS_Document : public Document_Writer {
+  
+public :
+  DriverUNV_W_SMESHDS_Document();
+  ~DriverUNV_W_SMESHDS_Document();
+ 
+  void Write();
+  //void SetFile(string);
+  //void SetDocument(Handle(SMESHDS_Document)&);
+
+private :
+  //Handle_SMESHDS_Document myDocument;
+  //string myFile; 
+
+};
+#endif
diff --git a/src/DriverUNV/DriverUNV_W_SMESHDS_Mesh.cxx b/src/DriverUNV/DriverUNV_W_SMESHDS_Mesh.cxx

new file mode 100644 (file)

index 0000000..231bba3
--- /dev/null
+++ b/src/DriverUNV/DriverUNV_W_SMESHDS_Mesh.cxx
@@ -0,0 +1,207 @@
+using namespace std;
+#include "DriverUNV_W_SMESHDS_Mesh.h"
+#include "DriverUNV_W_SMDS_Mesh.h"
+
+#include "SMDS_MeshElement.hxx"
+#include "SMDS_MeshNode.hxx"
+#include "SMDS_MeshEdgesIterator.hxx"
+#include "SMDS_MeshFacesIterator.hxx"
+#include "SMDS_MeshNodesIterator.hxx"
+#include "SMDS_MeshVolumesIterator.hxx"
+
+#include "utilities.h"
+
+#define sNODE_UNV_ID "  2411"
+#define sELT_UNV_ID  "  2412"
+#define sUNV_SEPARATOR "    -1"
+#define sNODE_UNV_DESCR "%10d         1         1        11\n"
+#define sELT_SURF_DESC  "%10d        %2d         1         1        11         %1d\n"
+#define sELT_VOLU_DESC  "%10d        %2d         1         1         9         %1d\n"
+#define sELT_BEAM_DESC1 "%10d        %2d         1         1         7         %1d\n"
+#define sELT_BEAM_DESC2 "         0         1         1\n"
+
+DriverUNV_W_SMESHDS_Mesh::DriverUNV_W_SMESHDS_Mesh() {
+;
+}
+
+DriverUNV_W_SMESHDS_Mesh::~DriverUNV_W_SMESHDS_Mesh() {
+;
+}
+
+void DriverUNV_W_SMESHDS_Mesh::SetMesh(Handle(SMDS_Mesh)& aMesh) {
+  //myMesh = Handle(SMESHDS_Mesh)::DownCast(aMesh);
+  myMesh = aMesh;
+}
+
+void DriverUNV_W_SMESHDS_Mesh::SetFile(string aFile) {
+  myFile = aFile;
+}
+
+void DriverUNV_W_SMESHDS_Mesh::SetFileId(FILE* aFileId) {
+  myFileId = aFileId;
+}
+
+void DriverUNV_W_SMESHDS_Mesh::SetMeshId(int aMeshId) {
+  myMeshId = aMeshId;
+}
+
+void DriverUNV_W_SMESHDS_Mesh::Write() {
+
+  string myClass = string("SMDS_Mesh");
+  string myExtension = string("UNV");
+
+  DriverUNV_W_SMDS_Mesh* myWriter = new DriverUNV_W_SMDS_Mesh;
+
+  myWriter->SetMesh(myMesh);
+  myWriter->SetFile(myFile);
+  myWriter->SetMeshId(myMeshId);
+  //myWriter->SetFileId(myFileId);
+
+  myWriter->Write();
+
+
+}
+void DriverUNV_W_SMESHDS_Mesh::Add() {
+  int nbNodes,nbCells;
+  int i;
+
+  char* file2Read = (char*)myFile.c_str();
+  myFileId = fopen(file2Read,"w+");
+  if (myFileId < 0)
+    {
+      fprintf(stderr,">> ERREUR : ouverture du fichier %s \n",file2Read);
+      exit(EXIT_FAILURE);
+    }
+  
+  /****************************************************************************
+  *                       NOMBRES D'OBJETS                                    *
+  ****************************************************************************/
+  fprintf(stdout,"\n(****************************)\n");
+  fprintf(stdout,"(* INFORMATIONS GENERALES : *)\n");
+  fprintf(stdout,"(****************************)\n");
+
+  /* Combien de noeuds ? */
+  nbNodes = myMesh->NbNodes();
+
+  /* Combien de mailles, faces ou aretes ? */
+  Standard_Integer nb_of_nodes, nb_of_edges,nb_of_faces, nb_of_volumes;
+  nb_of_edges = myMesh->NbEdges();
+  nb_of_faces = myMesh->NbFaces();
+  nb_of_volumes = myMesh->NbVolumes();
+  nbCells = nb_of_edges + nb_of_faces + nb_of_volumes;
+
+  fprintf(stdout,"%d %d\n",nbNodes,nbCells);
+//fprintf(myFileId,"%d %d\n",nbNodes,nbCells);
+
+  /****************************************************************************
+  *                       ECRITURE DES NOEUDS                                 *
+  ****************************************************************************/
+  fprintf(stdout,"\n(************************)\n");
+  fprintf(stdout,"(* NOEUDS DU MAILLAGE : *)\n");
+  fprintf(stdout,"(************************)\n");
+
+  SMDS_MeshNodesIterator itNodes(myMesh);
+
+  fprintf(myFileId,"%s\n", sUNV_SEPARATOR);
+  fprintf(myFileId,"%s\n", sNODE_UNV_ID  );
+
+  for (;itNodes.More();itNodes.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itNodes.Value();
+    const Handle(SMDS_MeshNode   )& node = myMesh->GetNode(1, elem);
+
+    fprintf(myFileId, sNODE_UNV_DESCR, node->GetID());
+    fprintf(myFileId, "%25.16E%25.16E%25.16E\n", node->X(), node->Y(), node->Z());
+  }
+  fprintf(myFileId,"%s\n", sUNV_SEPARATOR);
+
+  /****************************************************************************
+  *                       ECRITURE DES ELEMENTS                                *
+  ****************************************************************************/
+  fprintf(stdout,"\n(**************************)\n");
+  fprintf(stdout,"(* ELEMENTS DU MAILLAGE : *)\n");
+  fprintf(stdout,"(**************************)");
+  /* Ecriture des connectivites, noms, numeros des mailles */
+
+  fprintf(myFileId,"%s\n", sUNV_SEPARATOR);
+  fprintf(myFileId,"%s\n", sELT_UNV_ID   );
+
+  SMDS_MeshEdgesIterator itEdges(myMesh);
+  for (;itEdges.More();itEdges.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itEdges.Value();
+
+    switch (elem->NbNodes()) {
+      case 2 : {
+        fprintf(myFileId, sELT_BEAM_DESC1, elem->GetID(), 21, elem->NbNodes());
+        fprintf(myFileId, sELT_BEAM_DESC2);
+        fprintf(myFileId, "%10d%10d\n", elem->GetConnection(1), elem->GetConnection(2));
+        break;
+      }
+      case 3 : {
+        fprintf(myFileId, sELT_BEAM_DESC1, elem->GetID(), 24, elem->NbNodes());
+        fprintf(myFileId, sELT_BEAM_DESC2);
+        fprintf(myFileId, "%10d%10d%10d\n",elem->GetConnection(1), elem->GetConnection(2), elem->GetConnection(3));
+        break;
+      }
+    }
+  }
+
+  SMDS_MeshFacesIterator itFaces(myMesh);
+  for (;itFaces.More();itFaces.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itFaces.Value();
+
+    switch (elem->NbNodes()) {
+      case 3 :
+        // linear triangle
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 74, elem->NbNodes());
+        break;
+      case 4 :
+        // linear quadrilateral
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 71, elem->NbNodes());
+        break;
+      case 6 :
+        // parabolic triangle
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 72, elem->NbNodes());
+        break;
+      case 8 :
+        // parabolic quadrilateral
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 75, elem->NbNodes());
+        break;
+      default:
+        fprintf(myFileId, "element not registered\n");
+    }
+
+    for (i=0;i<elem->NbNodes();i++)
+      fprintf(myFileId,"%10d",elem->GetConnection(i+1));
+
+    fprintf(myFileId,"\n");
+  }
+
+  SMDS_MeshVolumesIterator itVolumes(myMesh);
+  for (;itVolumes.More();itVolumes.Next()) {
+    const Handle(SMDS_MeshElement)& elem = itVolumes.Value();
+
+    switch (elem->NbNodes()) {
+      case 4 :
+        // linear tetrahedron
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 111, elem->NbNodes());
+        break;
+      case 6 :
+        // linear tetrahedron
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 112, elem->NbNodes());
+        break;
+      case 8 :
+        // linear brick
+        fprintf(myFileId, sELT_SURF_DESC, elem->GetID(), 115, elem->NbNodes());
+        break;
+    }
+
+    for (i=0;i<elem->NbNodes();i++)
+      fprintf(myFileId,"%10d",elem->GetConnection(i+1));
+
+    fprintf(myFileId,"\n");
+  }
+  fprintf(myFileId,"%s\n", sUNV_SEPARATOR);
+
+  fclose (myFileId);
+}
+
diff --git a/src/DriverUNV/DriverUNV_W_SMESHDS_Mesh.h b/src/DriverUNV/DriverUNV_W_SMESHDS_Mesh.h

new file mode 100644 (file)

index 0000000..4cd508b
--- /dev/null
+++ b/src/DriverUNV/DriverUNV_W_SMESHDS_Mesh.h
@@ -0,0 +1,31 @@
+#ifndef _INCLUDE_DRIVERUNV_W_SMESHDS_MESH
+#define _INCLUDE_DRIVERUNV_W_SMESHDS_MESH
+
+#include <stdio.h>
+#include <string>
+
+#include "SMESHDS_Mesh.hxx"
+#include "Mesh_Writer.h"
+
+class DriverUNV_W_SMESHDS_Mesh : public Mesh_Writer {
+
+  public :
+    DriverUNV_W_SMESHDS_Mesh();
+    ~DriverUNV_W_SMESHDS_Mesh();
+
+    void Add();
+    void Write();
+    void SetMesh(Handle(SMDS_Mesh)& aMesh);
+    void SetFile(string);
+
+    void SetFileId(FILE*);
+    void SetMeshId(int);
+
+private :
+    Handle_SMDS_Mesh myMesh;
+    string myFile; 
+    FILE* myFileId;
+    int myMeshId;
+
+};
+#endif
diff --git a/src/DriverUNV/Makefile.in b/src/DriverUNV/Makefile.in

new file mode 100644 (file)

index 0000000..c9a4fa8
--- /dev/null
+++ b/src/DriverUNV/Makefile.in
@@ -0,0 +1,41 @@
+# -* Makefile *- 
+#
+# Author : Marc Tajchman (CEA)
+# Date : 5/07/2001
+# $Header$
+#
+
+# source path
+top_srcdir=@top_srcdir@
+top_builddir=../..
+srcdir=@srcdir@
+VPATH=.:@srcdir@
+
+
+@COMMENCE@
+
+# header files 
+EXPORT_HEADERS= DriverUNV_R_SMDS_Mesh.h DriverUNV_R_SMESHDS_Mesh.h DriverUNV_R_SMESHDS_Document.h \
+                DriverUNV_W_SMDS_Mesh.h DriverUNV_W_SMESHDS_Mesh.h DriverUNV_W_SMESHDS_Document.h 
+
+# Libraries targets
+LIB = libMeshDriverUNV.la
+LIB_SRC = DriverUNV_R_SMDS_Mesh.cxx DriverUNV_R_SMESHDS_Mesh.cxx DriverUNV_R_SMESHDS_Document.cxx \
+          DriverUNV_W_SMDS_Mesh.cxx DriverUNV_W_SMESHDS_Mesh.cxx DriverUNV_W_SMESHDS_Document.cxx 
+
+LIB_CLIENT_IDL =
+
+LIB_SERVER_IDL =
+
+# additionnal information to compil and link file
+CPPFLAGS += $(OCC_INCLUDES)
+CXXFLAGS += $(OCC_CXXFLAGS) $(MED2_INCLUDES)
+LDFLAGS  += $(OCC_LIBS) $(MED2_LIBS) -lMeshDriver
+
+%_moc.cxx: %.h
+       $(MOC) $< -o $@
+
+@CONCLUDE@
+
+
+
diff --git a/src/MEFISTO2/Makefile.in b/src/MEFISTO2/Makefile.in

new file mode 100644 (file)

index 0000000..63798ef
--- /dev/null
+++ b/src/MEFISTO2/Makefile.in
@@ -0,0 +1,36 @@
+# -* Makefile *- 
+#
+# Author : 
+# Date : 29/01/2001
+#
+#
+
+# source path
+top_srcdir=@top_srcdir@
+top_builddir=../..
+srcdir=@srcdir@
+VPATH=.:@srcdir@
+
+
+@COMMENCE@
+
+# header files 
+EXPORT_HEADERS = aptrte.h Rn.h
+
+# Libraries targets
+LIB = libMEFISTO2D.la
+LIB_SRC = aptrte.cxx trte.f
+# areteideale.f
+
+LIB_CLIENT_IDL =
+
+LIB_SERVER_IDL = 
+
+# additionnal information to compil and link file
+CPPFLAGS += $(OCC_INCLUDES) -I${KERNEL_ROOT_DIR}/include/salome
+CXXFLAGS += $(OCC_CXXFLAGS) -I${KERNEL_ROOT_DIR}/include/salome
+
+LDFLAGS += $(OCC_LIBS) -lg2c
+
+@CONCLUDE@
+
diff --git a/src/MEFISTO2/Rn.h b/src/MEFISTO2/Rn.h

new file mode 100755 (executable)

index 0000000..a2856ee
--- /dev/null
+++ b/src/MEFISTO2/Rn.h
@@ -0,0 +1,204 @@
+#ifndef Rn__h
+#define Rn__h
+
+#include <gp_Pnt.hxx>      //Dans OpenCascade
+#include <gp_Vec.hxx>      //Dans OpenCascade
+#include <gp_Dir.hxx>      //Dans OpenCascade
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// BUT:   Definir les espaces affines R R2 R3 R4 soit Rn pour n=1,2,3,4
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// AUTEUR : Frederic HECHT      ANALYSE NUMERIQUE UPMC  PARIS   OCTOBRE   2000
+// MODIFS : Alain    PERRONNET  ANALYSE NUMERIQUE UPMC  PARIS   NOVEMBRE  2000
+//...............................................................................
+#include <iostream>
+#include <cmath>
+
+using namespace std;
+
+template<class T> inline T Abs (const T &a){return a <0 ? -a : a;}
+template<class T> inline void Echange (T& a,T& b) {T c=a;a=b;b=c;}
+
+template<class T> inline T Min (const T &a,const T &b)  {return a < b ? a : b;}
+template<class T> inline T Max (const T &a,const T & b) {return a > b ? a : b;}
+
+template<class T> inline T Max (const T &a,const T & b,const T & c){return Max(Max(a,b),c);}
+template<class T> inline T Min (const T &a,const T & b,const T & c){return Min(Min(a,b),c);}
+
+template<class T> inline T Max (const T &a,const T & b,const T & c,const T & d)
+ {return Max(Max(a,b),Max(c,d));}
+template<class T> inline T Min (const T &a,const T & b,const T & c,const T & d)
+ {return Min(Min(a,b),Min(c,d));}
+
+//le type Nom des entites geometriques P L S V O
+//===========
+typedef char Nom[1+24];
+
+//le type N des nombres entiers positifs
+//=========
+typedef unsigned long int N;
+
+//le type Z des nombres entiers relatifs
+//=========
+typedef long int Z;
+
+//le type R des nombres "reels"
+//=========
+typedef double R;
+
+//le type XPoint  des coordonnees d'un pixel dans une fenetre
+//==============
+//typedef struct { short int x,y } XPoint;  //en fait ce type est defini dans X11-Window
+                                            // #include <X11/Xlib.h>
+//la classe R2
+//============
+class R2 
+{
+  friend ostream& operator << (ostream& f, const R2 & P)
+  { f << P.x << ' ' << P.y ; return f; }
+  friend istream& operator >> (istream& f, R2 & P)
+  { f >> P.x >> P.y ; return f; }
+
+  friend ostream& operator << (ostream& f, const R2 * P)
+  { f << P->x << ' ' << P->y ; return f; }
+  friend istream& operator >> (istream& f, R2 * P)
+  { f >> P->x >> P->y ; return f; }
+
+public:
+  R x,y;  //les donnees
+
+  R2 () :x(0),y(0) {}              //les constructeurs
+  R2 (R a,R b)   :x(a),y(b)  {}
+  R2 (R2 A,R2 B) :x(B.x-A.x),y(B.y-A.y)  {} //vecteur defini par 2 points
+
+  R2  operator+(R2 P) const {return R2(x+P.x,y+P.y);}     // Q+P possible
+  R2  operator+=(R2 P)  {x += P.x;y += P.y; return *this;}// Q+=P;
+  R2  operator-(R2 P) const {return R2(x-P.x,y-P.y);}     // Q-P
+  R2  operator-=(R2 P) {x -= P.x;y -= P.y; return *this;} // Q-=P;
+  R2  operator-()const  {return R2(-x,-y);}               // -Q
+  R2  operator+()const  {return *this;}                   // +Q
+  R   operator,(R2 P)const {return x*P.x+y*P.y;} // produit scalaire (Q,P)
+  R   operator^(R2 P)const {return x*P.y-y*P.x;} // produit vectoriel Q^P
+  R2  operator*(R c)const {return R2(x*c,y*c);}  // produit a droite  P*c
+  R2  operator*=(R c)  {x *= c; y *= c; return *this;}
+  R2  operator/(R c)const {return R2(x/c,y/c);}  // division par un reel
+  R2  operator/=(R c)  {x /= c; y /= c; return *this;}
+  R & operator[](int i) {return (&x)[i];}        // la coordonnee i
+  R2  orthogonal() {return R2(-y,x);}    //le vecteur orthogonal dans R2
+  friend R2 operator*(R c,R2 P) {return P*c;}    // produit a gauche c*P
+};
+
+
+//la classe R3
+//============
+class R3
+{
+  friend ostream& operator << (ostream& f, const R3 & P)
+  { f << P.x << ' ' << P.y << ' ' << P.z ; return f; }
+  friend istream& operator >> (istream& f, R3 & P)
+  { f >> P.x >> P.y >> P.z ; return f; }
+
+  friend ostream& operator << (ostream& f, const R3 * P)
+  { f << P->x << ' ' << P->y << ' ' << P->z ; return f; }
+  friend istream& operator >> (istream& f, R3 * P)
+  { f >> P->x >> P->y >> P->z ; return f; }
+
+public:  
+  R  x,y,z;  //les 3 coordonnees
+ 
+  R3 () :x(0),y(0),z(0) {}  //les constructeurs
+  R3 (R a,R b,R c):x(a),y(b),z(c)  {}                  //Point ou Vecteur (a,b,c)
+  R3 (R3 A,R3 B):x(B.x-A.x),y(B.y-A.y),z(B.z-A.z)  {}  //Vecteur AB
+
+  R3 (gp_Pnt P) : x(P.X()), y(P.Y()), z(P.Z()) {}      //Point     d'OpenCascade
+  R3 (gp_Vec V) : x(V.X()), y(V.Y()), z(V.Z()) {}      //Vecteur   d'OpenCascade
+  R3 (gp_Dir P) : x(P.X()), y(P.Y()), z(P.Z()) {}      //Direction d'OpenCascade
+
+  R3   operator+(R3 P)const  {return R3(x+P.x,y+P.y,z+P.z);}
+  R3   operator+=(R3 P)  {x += P.x; y += P.y; z += P.z; return *this;}
+  R3   operator-(R3 P)const  {return R3(x-P.x,y-P.y,z-P.z);}
+  R3   operator-=(R3 P)  {x -= P.x; y -= P.y; z -= P.z; return *this;}
+  R3   operator-()const  {return R3(-x,-y,-z);}
+  R3   operator+()const  {return *this;}
+  R    operator,(R3 P)const {return  x*P.x+y*P.y+z*P.z;} // produit scalaire
+  R3   operator^(R3 P)const {return R3(y*P.z-z*P.y ,P.x*z-x*P.z, x*P.y-y*P.x);} // produit vectoriel
+  R3   operator*(R c)const {return R3(x*c,y*c,z*c);}
+  R3   operator*=(R c)  {x *= c; y *= c; z *= c; return *this;}
+  R3   operator/(R c)const {return R3(x/c,y/c,z/c);}
+  R3   operator/=(R c)  {x /= c; y /= c; z /= c; return *this;}
+  R  & operator[](int i) {return (&x)[i];}
+  friend R3 operator*(R c,R3 P) {return P*c;}
+
+  R3   operator=(gp_Pnt P) {return R3(P.X(),P.Y(),P.Z());}
+  R3   operator=(gp_Dir P) {return R3(P.X(),P.Y(),P.Z());}
+
+  friend gp_Pnt gp_pnt(R3 xyz) { return gp_Pnt(xyz.x,xyz.y,xyz.z); }
+  //friend gp_Pnt operator=() { return gp_Pnt(x,y,z); }
+  friend gp_Dir gp_dir(R3 xyz) { return gp_Dir(xyz.x,xyz.y,xyz.z); }
+
+  bool  DansPave( R3 & xyzMin, R3 & xyzMax )
+    { return xyzMin.x<=x && x<=xyzMax.x &&
+            xyzMin.y<=y && y<=xyzMax.y &&
+            xyzMin.z<=z && z<=xyzMax.z; }
+};
+
+//la classe R4
+//============
+class R4: public R3
+{
+  friend ostream& operator <<(ostream& f, const R4 & P )
+  { f << P.x << ' ' << P.y << ' ' << P.z << ' ' << P.omega; return f; }
+  friend istream& operator >>(istream& f,  R4 & P)
+  { f >> P.x >>  P.y >>  P.z >> P.omega ; return f; }
+
+  friend ostream& operator <<(ostream& f, const R4 * P )
+  { f << P->x << ' ' << P->y << ' ' << P->z << ' ' << P->omega; return f; }
+  friend istream& operator >>(istream& f,  R4 * P)
+  { f >> P->x >>  P->y >>  P->z >> P->omega ; return f; }
+
+public:  
+  R  omega;  //la donnee du poids supplementaire
+ 
+  R4 () :omega(1.0) {}  //les constructeurs
+  R4 (R a,R b,R c,R d):R3(a,b,c),omega(d) {}
+  R4 (R4 A,R4 B) :R3(B.x-A.x,B.y-A.y,B.z-A.z),omega(B.omega-A.omega) {}
+
+  R4   operator+(R4 P)const  {return R4(x+P.x,y+P.y,z+P.z,omega+P.omega);}
+  R4   operator+=(R4 P)  {x += P.x;y += P.y;z += P.z;omega += P.omega;return *this;}
+  R4   operator-(R4 P)const  {return R4(x-P.x,y-P.y,z-P.z,omega-P.omega);}
+  R4   operator-=(R4 P) {x -= P.x;y -= P.y;z -= P.z;omega -= P.omega;return *this;}
+  R4   operator-()const  {return R4(-x,-y,-z,-omega);}
+  R4   operator+()const  {return *this;}
+  R    operator,(R4 P)const {return  x*P.x+y*P.y+z*P.z+omega*P.omega;} // produit scalaire
+  R4   operator*(R c)const {return R4(x*c,y*c,z*c,omega*c);}
+  R4   operator*=(R c)  {x *= c; y *= c; z *= c; omega *= c; return *this;}
+  R4   operator/(R c)const {return R4(x/c,y/c,z/c,omega/c);}
+  R4   operator/=(R c)  {x /= c; y /= c; z /= c; omega /= c; return *this;}
+  R  & operator[](int i) {return (&x)[i];}
+  friend R4 operator*(R c,R4 P) {return P*c;}
+};
+
+//quelques fonctions supplementaires sur ces classes
+//==================================================
+inline R Aire2d(const R2 A,const R2 B,const R2 C){return (B-A)^(C-A);} 
+inline R Angle2d(R2 P){ return atan2(P.y,P.x);}
+
+inline R Norme2_2(const R2 & A){ return (A,A);}
+inline R Norme2(const R2 & A){ return sqrt((A,A));}
+inline R NormeInfinie(const R2 & A){return Max(Abs(A.x),Abs(A.y));}
+
+inline R Norme2_2(const R3 & A){ return (A,A);}
+inline R Norme2(const R3 & A){ return sqrt((A,A));}
+inline R NormeInfinie(const R3 & A){return Max(Abs(A.x),Abs(A.y),Abs(A.z));}
+
+inline R Norme2_2(const R4 & A){ return (A,A);}
+inline R Norme2(const R4 & A){ return sqrt((A,A));}
+inline R NormeInfinie(const R4 & A){return Max(Abs(A.x),Abs(A.y),Abs(A.z),Abs(A.omega));}
+
+inline R2 XY(R3 P) {return R2(P.x, P.y);}  //restriction a R2 d'un R3 par perte de z
+inline R3 Min(R3 P, R3 Q) 
+{return R3(P.x<Q.x ? P.x : Q.x, P.y<Q.y ? P.y : Q.y, P.z<Q.z ? P.z : Q.z);} //Pt de xyz Min
+inline R3 Max(R3 P, R3 Q) 
+{return R3(P.x>Q.x ? P.x : Q.x, P.y>Q.y ? P.y : Q.y, P.z>Q.z ? P.z : Q.z);} //Pt de xyz Max
+
+#endif
diff --git a/src/MEFISTO2/aptrte.cxx b/src/MEFISTO2/aptrte.cxx

new file mode 100755 (executable)

index 0000000..58da981
--- /dev/null
+++ b/src/MEFISTO2/aptrte.cxx
@@ -0,0 +1,760 @@
+using namespace std;
+#include "Rn.h"
+#include "aptrte.h"
+#include "utilities.h"
+
+extern "C"
+{
+  R aretemaxface_;
+  R areteideale_( R3 xyz, R3 direction )
+  {
+    return aretemaxface_;
+  }
+}
+//calcul de la longueur ideale de l'arete au sommet xyz (z ici inactif)
+//dans la direction donnee
+//a ajuster pour chaque surface plane et selon l'entier notysu (voir plus bas)
+
+
+static double cpunew, cpuold=0;
+
+void tempscpu_( double & tempsec )
+//Retourne le temps CPU utilise en secondes
+{  
+  tempsec = ( (double) clock() ) / CLOCKS_PER_SEC;
+  //MESSAGE( "temps cpu=" << tempsec );
+}
+
+
+void deltacpu_( R & dtcpu )
+//Retourne le temps CPU utilise en secondes depuis le precedent appel
+{
+  tempscpu_( cpunew );
+  dtcpu  = R( cpunew - cpuold );
+  cpuold = cpunew;
+  //MESSAGE( "delta temps cpu=" << dtcpu );
+  return;
+}
+
+
+void  aptrte( Z nutysu, R aretmx,
+             Z nblf,   Z * nudslf, R2 * uvslf,
+             Z nbpti,  R2 *uvpti,
+             Z & nbst, R2 * & uvst, Z & nbt, Z * & nust,
+             Z & ierr )
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// but : appel de la triangulation par un arbre-4 recouvrant
+// ----- de triangles equilateraux
+//       le contour du domaine plan est defini par des lignes fermees
+//       la premiere ligne etant l'enveloppe de toutes les autres
+//       la fonction areteideale(s,d) donne la taille d'arete
+//       au point s dans la direction (actuellement inactive) d
+//       des lors toute arete issue d'un sommet s devrait avoir une longueur
+//       comprise entre 0.65 areteideale_(s,d) et 1.3 areteideale_(s,d)
+//
+//Attention:
+//  Les tableaux uvslf et uvpti sont supposes ne pas avoir de sommets identiques!
+//  De meme, un sommet d'une ligne fermee ne peut appartenir a une autre ligne fermee
+//
+// entrees:
+// --------
+// nutysu : numero de traitement de areteideale_(s,d) selon le type de surface
+//          0 pas d'emploi de la fonction areteideale_() et aretmx est active
+//          1 il existe une fonction areteideale_(s,d)
+//            dont seules les 2 premieres composantes de uv sont actives
+//          ... autres options a definir ...
+// aretmx : longueur maximale des aretes de la future triangulation
+// nblf   : nombre de lignes fermees de la surface
+// nudslf : numero du dernier sommet de chacune des nblf lignes fermees
+//          nudslf(0)=0 pour permettre la difference sans test
+//          Attention le dernier sommet de chaque ligne est raccorde au premier
+//          tous les sommets et les points internes ont des coordonnees
+//          UV differentes <=> Pas de point double!
+// uvslf  : uv des nudslf(nblf) sommets des lignes fermees
+// nbpti  : nombre de points internes futurs sommets de la triangulation
+// uvpti  : uv des points internes futurs sommets de la triangulation
+//
+// sorties:
+// --------
+// nbst   : nombre de sommets de la triangulation finale
+// uvst   : coordonnees uv des nbst sommets de la triangulation
+// nbt    : nombre de triangles de la triangulation finale
+// nust   : 4 numeros dans uvst des sommets des nbt triangles
+//          s1, s2, s3, 0: no dans uvst des 3 sommets et 0 car quadrangle!
+// ierr   : 0 si pas d'erreur
+//        > 0 sinon
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// auteur : Alain Perronnet  Analyse Numerique Paris UPMC   decembre 2001
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+{
+  R  d, tcpu=0;
+  R3 direction=R3(0,0,0);  //direction pour areteideale() inactive ici!
+  Z  nbarfr=nudslf[nblf];  //nombre total d'aretes des lignes fermees
+  Z  mxtrou = Max( 1024, nblf );  //nombre maximal de trous dans la surface
+
+  R3 *mnpxyd=NULL;
+  Z  *mnsoar=NULL, mosoar=7, mxsoar, n1soar; //le hachage des aretes
+  Z  *mnartr=NULL, moartr=3, mxartr, n1artr; //le no des 3 aretes des triangles
+  Z  *mntree=NULL, motree=9, mxtree;   //L'arbre 4 de TE et nombre d'entiers par TE
+  Z  *mnqueu=NULL, mxqueu;
+  Z  *mn1arcf=NULL;
+  Z  *mnarcf=NULL, mxarcf;
+  Z  *mnarcf1=NULL;
+  Z  *mnarcf2=NULL;
+  Z  *mnarcf3=NULL;
+  Z  *mntrsu=NULL;
+  Z  *mndalf=NULL;
+  Z  *mnslig=NULL;
+  Z  *mnarst=NULL;
+  Z  *mnlftr=NULL;
+
+  R3 comxmi[2];            //coordonnees UV Min et Maximales
+  R  aremin, aremax;       //longueur minimale et maximale des aretes
+  R  quamoy, quamin;
+
+  Z  noar0, noar, na;
+  Z  i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt;
+  Z  mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn;
+  Z  moins1=-1;
+
+  aretemaxface_ = aretmx;
+
+  // initialisation du temps cpu
+  deltacpu_( d );
+  ierr = 0;
+
+  // quelques reservations de tableaux pour faire les calculs
+  // ========================================================
+  // le tableau pointeur sur la premiere arete de chaque ligne fermee
+  if( mndalf!=NULL ) delete [] mndalf;
+  mndalf = new Z[1+nblf];
+  if( mndalf==NULL ) goto ERREUR;
+  mndalf[0]=0;
+
+  // declaration du tableau des coordonnees des sommets de la frontiere
+  // puis des sommets internes ajoutes
+  // majoration empirique du nombre de sommets de la triangulation
+  i =  4*nbarfr/10;
+  mxsomm = Max( 20000, 64*nbpti+i*i );
+  MESSAGE( "APTRTE: Depart de la triangulation avec " );
+  MESSAGE( "nutysu=" << nutysu << "  aretmx=" << aretmx << "  mxsomm=" << mxsomm );
+
+ NEWDEPART:
+  //mnpxyd( 3, mxsomm ) les coordonnees UV des sommets et la taille d'arete aux sommets
+  if( mnpxyd!=NULL ) delete [] mnpxyd;
+  mnpxyd = new R3[mxsomm];
+  if( mnpxyd==NULL ) goto ERREUR;
+
+  // le tableau mnsoar des aretes des triangles
+  // 1: sommet 1 dans pxyd,
+  // 2: sommet 2 dans pxyd,
+  // 3: numero de 1 a nblf de la ligne qui supporte l'arete
+  // 4: numero dans mnartr du triangle 1 partageant cette arete,
+  // 5: numero dans mnartr du triangle 2 partageant cette arete,
+  // 6: chainage des aretes frontalieres ou internes ou
+  //    des aretes simples des etoiles de triangles,
+  // 7: chainage du hachage des aretes
+  // nombre d'aretes = 3 ( nombre de sommets - 1 + nombre de trous )
+  // pour le hachage des aretes mxsoar doit etre > 3*mxsomm!
+  // h(ns1,ns2) = min( ns1, ns2 )
+  if( mnsoar!=NULL ) delete [] mnsoar;
+  mxsoar = 3 * ( mxsomm + mxtrou );
+  mnsoar = new Z[mosoar*mxsoar];
+  if( mnsoar==NULL ) goto ERREUR;
+  //initialiser le tableau mnsoar pour le hachage des aretes
+  insoar_( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
+
+  // mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets
+  if( mnarst!=NULL ) delete [] mnarst;
+  mnarst = new Z[1+mxsomm];
+  if( mnarst==NULL ) goto ERREUR;
+  n = 1+mxsomm;
+  azeroi_( n, mnarst );
+
+  // mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier
+  //               ou no du point si interne forc'e par l'utilisateur
+  //               ou  0 si interne cree par le module
+  if( mnslig!=NULL ) delete [] mnslig;
+  mnslig = new Z[mxsomm];
+  if( mnslig==NULL ) goto ERREUR;
+  azeroi_( mxsomm, mnslig );
+
+  // initialisation des aretes frontalieres de la triangulation future
+  // renumerotation des sommets des aretes des lignes pour la triangulation
+  // mise a l'echelle des coordonnees des sommets pour obtenir une
+  // meilleure precision lors des calculs + quelques verifications
+  // boucle sur les lignes fermees qui forment la frontiere
+  // ======================================================================
+  noar = 0;
+  aremin = 1e100;
+  aremax = 0;
+
+  for (n=1; n<=nblf; n++)
+  {
+    //l'initialisation de la premiere arete de la ligne n dans la triangulation
+    //-------------------------------------------------------------------------
+    //le sommet ns0 est le numero de l'origine de la ligne
+    ns0 = nudslf[n-1];
+    mnpxyd[ns0].x = uvslf[ns0].x;
+    mnpxyd[ns0].y = uvslf[ns0].y;
+    mnpxyd[ns0].z = areteideale_( mnpxyd[ns0], direction );
+//     cout << "Sommet " << ns0 << ": " << mnpxyd[ns0].x
+//      << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z << endl;
+
+    //carre de la longueur de l'arete 1 de la ligne fermee n
+    d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 ) 
+      + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ;
+    aremin = Min( aremin, d );
+    aremax = Max( aremax, d );
+
+    //le numero des 2 sommets (ns1,ns2) de la premiere arete de la ligne
+    //initialisation de la 1-ere arete ns1-ns1+1 de cette ligne fermee n
+    //le numero des 2 sommets ns1 ns2 de la 1-ere arete
+    //Attention: les numeros ns debutent a 1 (ils ont >0)
+    //           les tableaux c++ demarrent a zero!
+    //           les tableaux fortran demarrent ou l'on veut!
+    ns0++;
+    ns1 = ns0;
+    ns2 = ns1+1;
+
+     //le numero n de la ligne du sommet et son numero ns1 dans la ligne
+    mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1];
+    fasoar_( ns1, ns2, moins1, moins1, n,
+            mosoar, mxsoar, n1soar, mnsoar, mnarst,
+            noar0,  ierr );
+    //pas de test sur ierr car pas de saturation possible a ce niveau
+
+    //le pointeur dans le hachage sur la premiere arete de la ligne fermee n
+    mndalf[n] = noar0;
+
+    //la nouvelle arete est la suivante de l'arete definie juste avant
+    if( noar > 0 )
+      mnsoar[mosoar * noar - mosoar + 5] = noar0;
+
+    //l'initialisation des aretes suivantes de la ligne dans la triangulation
+    //-----------------------------------------------------------------------
+    nbarli = nudslf[n] - nudslf[n-1];  //nombre d'aretes=sommets de la ligne n
+    for (i=2; i<=nbarli; i++)
+    {
+      ns1 = ns2; //le numero de l'arete et le numero du premier sommet de l'arete
+      if( i < nbarli )
+       //nbs+1 est le 2-eme sommet de l'arete i de la ligne fermee n
+       ns2 = ns1+1;
+      else
+       //le 2-eme sommet de la derniere arete est le premier sommet de la ligne
+       ns2 = ns0;
+
+      //l'arete precedente est dotee de sa suivante:celle cree ensuite
+      //les 2 coordonnees du sommet ns2 de la ligne
+      ns = ns1 - 1;
+      mnpxyd[ns].x = uvslf[ns].x;
+      mnpxyd[ns].y = uvslf[ns].y;
+      mnpxyd[ns].z = areteideale_( mnpxyd[ns], direction );
+//       cout << "Sommet " << ns << ": " << mnpxyd[ns].x
+//        << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z << endl;
+
+      //carre de la longueur de l'arete
+      d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2) 
+        + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2);
+      aremin = Min( aremin, d );
+      aremax = Max( aremax, d );
+
+      //le numero n de la ligne du sommet et son numero ns1 dans la ligne
+      mnslig[ns] = 1000000 * n + ns1-nudslf[n-1];
+
+      //ajout de l'arete dans la liste
+      fasoar_( ns1, ns2, moins1, moins1, n,
+              mosoar, mxsoar, n1soar, mnsoar,
+              mnarst, noar, ierr );
+      //pas de test sur ierr car pas de saturation possible a ce niveau
+
+      //chainage des aretes frontalieres en position 6 du tableau mnsoar
+      //la nouvelle arete est la suivante de l'arete definie juste avant
+      mnsoar[ mosoar * noar0 - mosoar + 5 ] = noar;
+      noar0 = noar;
+   }
+    //attention: la derniere arete de la ligne fermee enveloppe
+    //           devient en fait la premiere arete de cette ligne
+    //           dans le chainage des aretes de la frontiere!
+  }
+  if( ierr != 0 ) goto ERREUR;
+
+  aremin = sqrt( aremin );  //longueur minimale d'une arete des lignes fermees
+  aremax = sqrt( aremax );  //longueur maximale d'une arete
+
+  aretmx = Min( aretmx, aremax );  //pour homogeneiser
+  cout << "nutysu=" << nutysu << "  aretmx=" << aretmx 
+       << "  arete min=" << aremin << "  arete max=" << aremax << endl;
+
+  //chainage des aretes frontalieres : la derniere arete frontaliere
+  mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
+
+  //tous les sommets et aretes frontaliers sont numerotes de 1 a nbarfr
+  //reservation du tableau des numeros des 3 aretes de chaque triangle
+  //mnartr( moartr, mxartr )
+  //En nombre: Triangles = Aretes Internes + Aretes Frontalieres - Sommets + 1-Trous
+  //          3Triangles = 2 Aretes internes + Aretes frontalieres
+  //       d'ou 3T/2 < AI + AF => T < 3T/2  - Sommets + 1-Trous
+  //nombre de triangles < 2 ( nombre de sommets - 1 + nombre de trous )
+  if( mnartr!=NULL ) delete [] mnartr;
+  mxartr = 2 * ( mxsomm + mxtrou );
+  mnartr = new Z[moartr*mxartr];
+  if( mnartr==NULL ) goto ERREUR;
+
+  //Ajout des points internes
+  ns1 = nudslf[ nblf ];
+  for (i=0; i<nbpti; i++)
+  {
+    //les 2 coordonnees du point i de sommet nbs
+    mnpxyd[ns1].x = uvpti[i].x;
+    mnpxyd[ns1].y = uvpti[i].y;
+    mnpxyd[ns1].z = areteideale_( mnpxyd[ns1], direction );
+    //le numero i du point interne
+    mnslig[ns1] = i+1;
+    ns1++;
+  }
+
+  //nombre de sommets de la frontiere et internes
+  nbarpi = ns1;
+
+  // creation de l'arbre-4 des te (tableau letree)
+  // ajout dans les te des sommets des lignes et des points internes imposes
+  // =======================================================================
+  // premiere estimation de mxtree
+  mxtree = 2 * mxsomm;
+
+ NEWTREE:  //en cas de saturation de l'un des tableaux, on boucle
+  MESSAGE( "Debut triangulation avec mxsomm=" << mxsomm );
+  if( mntree != NULL ) delete [] mntree;
+  nbsomm = nbarpi;
+  mntree = new Z[motree*(1+mxtree)];
+  if( mntree==NULL ) goto ERREUR;
+
+  //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
+  teajte_( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
+  comxmi[0].z=0;
+  comxmi[1].z=0;
+
+  if( ierr == 51 )
+  {
+    //saturation de letree => sa taille est augmentee et relance
+    mxtree = mxtree * 2;
+    ierr   = 0;
+    MESSAGE( "Nouvelle valeur de mxtree=" << mxtree );
+    goto NEWTREE;
+  }
+
+  deltacpu_( d );
+  tcpu += d;
+  MESSAGE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" );
+  if( ierr != 0 ) goto ERREUR;
+  //ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes
+
+  // homogeneisation de l'arbre des te a un saut de taille au plus
+  // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
+  // ===========================================================================
+  // reservation de la queue pour parcourir les te de l'arbre
+  if( mnqueu != NULL ) delete [] mnqueu;
+  mxqueu = mxtree;
+  mnqueu = new Z[mxqueu];
+  if( mnqueu==NULL) goto ERREUR;
+
+  tehote_( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
+          comxmi, aretmx,
+          mntree, mxqueu, mnqueu,
+          ierr );
+
+  deltacpu_( d );
+  tcpu += d;
+  cout << "Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
+       << d << " secondes" << endl;
+  if( ierr != 0 )
+  {
+    //destruction du tableau auxiliaire et de l'arbre
+    if( ierr == 51 )
+    {
+      //letree sature
+      mxtree = mxtree * 2;
+      MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree );
+      ierr = 0;
+      goto NEWTREE;
+    }
+    else
+      goto ERREUR;
+  }
+
+  // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
+  // et des points de la frontiere, des points internes imposes interieurs
+  // ==========================================================================
+  tetrte_( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
+          mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
+          moartr, mxartr, n1artr, mnartr, mnarst,
+          ierr );
+
+  // destruction de la queue et de l'arbre devenus inutiles
+  delete [] mnqueu;  mnqueu=NULL;
+  delete [] mntree;  mntree=NULL;
+
+  //Temps calcul
+  deltacpu_( d );
+  tcpu += d;
+  MESSAGE( "Temps de la triangulation des TE=" << d << " secondes" );
+
+  // ierr =0 si pas d'erreur
+  //      =1 si le tableau mnsoar est sature
+  //      =2 si le tableau mnartr est sature
+  //      =3 si aucun des triangles ne contient l'un des points internes
+  //      =5 si saturation de la queue de parcours de l'arbre des te
+  if( ierr != 0 ) goto ERREUR;
+
+  //qualites de la triangulation actuelle
+  qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
+              nbt, quamoy, quamin );
+
+  // boucle sur les aretes internes (non sur une ligne de la frontiere)
+  // avec echange des 2 diagonales afin de rendre la triangulation delaunay
+  // ======================================================================
+  // formation du chainage 6 des aretes internes a echanger eventuellement
+  aisoar_( mosoar, mxsoar, mnsoar, na );
+  tedela_( mnpxyd, mnarst,
+          mosoar, mxsoar, n1soar, mnsoar, na,
+          moartr, mxartr, n1artr, mnartr, n );
+
+  MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
+  deltacpu_( d );
+  tcpu += d;
+  cout << "Temps de la triangulation Delaunay par echange des diagonales="
+       << d << " secondes" << endl;
+
+  //qualites de la triangulation actuelle
+  qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
+              nbt, quamoy, quamin );
+
+  // detection des aretes frontalieres initiales perdues
+  // triangulation frontale pour les restaurer
+  // ===================================================
+  mxarcf = mxsomm/5;
+  if( mn1arcf != NULL ) delete [] mn1arcf;
+  if( mnarcf  != NULL ) delete [] mnarcf;
+  if( mnarcf1 != NULL ) delete [] mnarcf1;
+  if( mnarcf2 != NULL ) delete [] mnarcf2;
+  mn1arcf = new Z[1+mxarcf];
+  if( mn1arcf == NULL ) goto ERREUR;
+  mnarcf  = new Z[3*mxarcf];
+  if( mnarcf == NULL ) goto ERREUR;
+  mnarcf1 = new Z[mxarcf];
+  if( mnarcf1 == NULL ) goto ERREUR;
+  mnarcf2 = new Z[mxarcf];
+  if( mnarcf2 == NULL ) goto ERREUR;
+
+  terefr_( nbarpi, mnpxyd,
+          mosoar, mxsoar, n1soar, mnsoar,
+          moartr, n1artr, mnartr, mnarst,
+          mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
+          n, ierr );
+
+  MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere" );
+  deltacpu_( d );
+  tcpu += d;
+  cout << "Temps de la recuperation des aretes perdues de la frontiere="
+       << d << " secondes" << endl;
+
+  if( ierr != 0 ) goto ERREUR;
+
+  //qualites de la triangulation actuelle
+  qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
+              nbt, quamoy, quamin );
+
+  // fin de la triangulation avec respect des aretes initiales frontalieres
+
+  // suppression des triangles externes a la surface
+  // ===============================================
+  // recherche du dernier triangle utilise
+  mn = mxartr * moartr;
+  for ( ndtri0=mxartr; ndtri0<=1; ndtri0-- )
+  {
+    mn -= moartr;
+    if( mnartr[mn] != 0 ) break;
+  }
+
+  if( mntrsu != NULL ) delete [] mntrsu;
+  mntrsu = new Z[ndtri0];
+  if( mntrsu == NULL ) goto ERREUR;
+
+  if( mnlftr != NULL ) delete [] mnlftr;
+  mnlftr = new Z[nblf];
+  if( mnlftr == NULL ) goto ERREUR;
+
+  for (n=0; n<nblf; n++)  //numero de la ligne fermee de 1 a nblf
+    mnlftr[n] = n+1;
+
+  tesuex_( nblf,   mnlftr,
+          ndtri0, nbsomm, mnpxyd, mnslig,
+          mosoar, mxsoar, mnsoar,
+          moartr, mxartr, n1artr, mnartr, mnarst,
+          nbt, mntrsu, ierr );
+
+  delete [] mnlftr; mnlftr=NULL;
+  delete [] mntrsu; mntrsu=NULL;
+
+  deltacpu_( d );
+  tcpu += d;
+  MESSAGE( "Temps de la suppression des triangles externes=" << d );
+  if( ierr != 0 ) goto ERREUR;
+
+  //qualites de la triangulation actuelle
+  qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
+              nbt, quamoy, quamin );
+
+  // amelioration de la qualite de la triangulation par
+  // barycentrage des sommets internes a la triangulation
+  // suppression des aretes trop longues ou trop courtes
+  // modification de la topologie des groupes de triangles
+  // mise en delaunay de la triangulation
+  // =====================================================
+  mnarcf3 = new Z[mxarcf];
+  if( mnarcf3 == NULL ) goto ERREUR;
+
+  teamqt_( nutysu,
+          mnarst, mosoar, mxsoar, n1soar, mnsoar,
+          moartr, mxartr, n1artr, mnartr,
+          mxarcf, mnarcf2, mnarcf3,
+          mn1arcf, mnarcf, mnarcf1,
+          comxmi, nbarpi, nbsomm, mxsomm, mnpxyd, mnslig,
+          ierr );
+  if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
+  if( mnarcf  != NULL ) {delete [] mnarcf;  mnarcf =NULL;}
+  if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
+  if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
+  if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
+
+  deltacpu_( d );
+  tcpu += d;
+  MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
+  if( ierr != 0 ) goto ERREUR;
+
+  //qualites de la triangulation finale
+  qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
+              nbt, quamoy, quamin );
+
+  // renumerotation des sommets internes: mnarst(i)=numero final du sommet
+  // ===================================
+  for (i=0; i<=nbsomm; i++)
+    mnarst[i] = 0;
+
+  for (nt=1; nt<=mxartr; nt++)
+  {
+    if( mnartr[nt*moartr-moartr] != 0 )
+    {
+      //le numero des 3 sommets du triangle nt
+      nusotr_( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
+      //les 3 sommets du triangle sont actifs
+      mnarst[ nosotr[0] ] = 1;
+      mnarst[ nosotr[1] ] = 1;
+      mnarst[ nosotr[2] ] = 1;
+    }
+  }
+  nbst = 0;
+  for (i=1; i<=nbsomm; i++)
+  {
+    if( mnarst[i] >0 )
+      mnarst[i] = ++nbst;
+  }
+
+  // generation du tableau uvst de la surface triangulee
+  // ---------------------------------------------------
+  if( uvst != NULL ) delete [] uvst;
+  uvst = new R2[nbst];
+  if( uvst == NULL ) goto ERREUR;
+
+  nbst=-1;
+  for (i=0; i<nbsomm; i++ )
+  {
+    if( mnarst[i+1]>0 )
+    {
+      nbst++;
+      uvst[nbst].x = mnpxyd[i].x;
+      uvst[nbst].y = mnpxyd[i].y;
+
+      //si le sommet est un point ou appartient a une ligne
+      //ses coordonnees initiales sont restaurees
+      n = mnslig[i];
+      if( n > 0 )
+      {
+       if( n >= 1000000 )
+       {
+         //sommet d'une ligne
+         //retour aux coordonnees initiales dans uvslf
+         l = n / 1000000;
+         n = n - 1000000 * l + nudslf[l-1] - 1;
+         uvst[nbst].x = uvslf[n].x;
+         uvst[nbst].y = uvslf[n].y;
+       }
+       else
+       {
+         //point utilisateur n interne impose
+         //retour aux coordonnees initiales dans uvpti
+         uvst[nbst].x = uvpti[n-1].x;
+         uvst[nbst].y = uvpti[n-1].y;
+       }
+      }
+    }
+  }
+  nbst++;
+
+  // generation du tableau 'nsef' de la surface triangulee
+  // -----------------------------------------------------
+  // boucle sur les triangles occupes (internes et externes)
+  if( nust != NULL ) delete [] nust;
+  nust = new Z[4*nbt];
+  if( nust == NULL ) goto ERREUR;
+  nbt = 0;
+  for (i=1; i<=mxartr; i++)
+  {
+    //le triangle i de mnartr
+    if( mnartr[i*moartr-moartr] != 0 )
+    {
+      //le triangle i est interne => nosotr numero de ses 3 sommets
+      nusotr_( i, mosoar, mnsoar, moartr, mnartr,  nosotr );
+      nust[nbt++] = mnarst[ nosotr[0] ];
+      nust[nbt++] = mnarst[ nosotr[1] ];
+      nust[nbt++] = mnarst[ nosotr[2] ];
+      nust[nbt++] = 0;
+    }
+  }
+  nbt /= 4;  //le nombre final de triangles de la surface
+  cout << "Nombre de sommets=" << nbst
+       << "  Nombre de triangles=" << nbt << endl;
+
+  deltacpu_( d );
+  tcpu += d;
+  MESSAGE( "Temps total de la triangulation=" << tcpu << " secondes" );
+
+  // destruction des tableaux auxiliaires
+  // ------------------------------------
+ NETTOYAGE:
+  if( mnarst != NULL ) delete [] mnarst;
+  if( mnartr != NULL ) delete [] mnartr;
+  if( mnslig != NULL ) delete [] mnslig;
+  if( mnsoar != NULL ) delete [] mnsoar;
+  if( mnpxyd != NULL ) delete [] mnpxyd;
+  if( mndalf != NULL ) delete [] mndalf;
+  if( mntree != NULL ) delete [] mntree;
+  if( mnqueu != NULL ) delete [] mnqueu;
+  if( mntrsu != NULL ) delete [] mntrsu;
+  if( mnlftr != NULL ) delete [] mnlftr;
+  if( mn1arcf != NULL ) delete [] mn1arcf;
+  if( mnarcf  != NULL ) delete [] mnarcf;
+  if( mnarcf1 != NULL ) delete [] mnarcf1;
+  if( mnarcf2 != NULL ) delete [] mnarcf2;
+  if( mnarcf3 != NULL ) delete [] mnarcf3;
+  return;
+
+ ERREUR:
+  if( ierr == 51 || ierr == 52 )
+  {
+    //saturation des sommets => redepart avec 2 fois plus de sommets
+    mxsomm = 2 * mxsomm;
+    ierr   = 0;
+    goto NEWDEPART;
+  }
+  else
+  {
+    MESSAGE( "Triangulation non realisee " << ierr );
+    if( ierr == 0 ) ierr=1;
+    goto NETTOYAGE;
+  }
+}
+
+
+void qualitetrte( R3 *mnpxyd,
+                 Z & mosoar, Z & mxsoar, Z *mnsoar,
+                 Z & moartr, Z & mxartr, Z *mnartr,
+                 Z & nbtria, R & quamoy, R & quamin )
+// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// but :    calculer la qualite moyenne et minimale de la triangulation
+// -----    actuelle definie par les tableaux mnsoar et mnartr
+// entrees:
+// --------
+// mnpxyd : tableau des coordonnees 2d des points
+//          par point : x  y  distance_souhaitee
+// mosoar : nombre maximal d'entiers par arete et
+//          indice dans mnsoar de l'arete suivante dans le hachage
+// mxsoar : nombre maximal d'aretes stockables dans le tableau mnsoar
+//          attention: mxsoar>3*mxsomm obligatoire!
+// mnsoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
+//          chainage des aretes frontalieres, chainage du hachage des aretes
+//          hachage des aretes = mnsoar(1)+mnsoar(2)*2
+//          avec mxsoar>=3*mxsomm
+//          une arete i de mnsoar est vide <=> mnsoar(1,i)=0 et
+//          mnsoar(2,arete vide)=l'arete vide qui precede
+//          mnsoar(3,arete vide)=l'arete vide qui suit
+// moartr : nombre maximal d'entiers par arete du tableau mnartr
+// mxartr : nombre maximal de triangles declarables
+// mnartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
+//          arete1 = 0 si triangle vide => arete2 = triangle vide suivant
+// sorties:
+// --------
+// nbtria : nombre de triangles internes au domaine
+// quamoy : qualite moyenne  des triangles actuels
+// quamin : qualite minimale des triangles actuels
+// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+{
+  R  d, aire, qualite;
+  Z  nosotr[3], mn, nbtrianeg, nt;
+
+  aire   = 0;
+  quamoy = 0;
+  quamin = 2.0;
+  nbtria = 0;
+  nbtrianeg = 0;
+
+  mn = -moartr;
+  for ( nt=1; nt<=mxartr; nt++ )
+  {
+    mn += moartr;
+    if( mnartr[mn]!=0 )
+    {
+      //un triangle occupe de plus
+      nbtria++;
+
+      //le numero des 3 sommets du triangle nt
+      nusotr_( nt, mosoar, mnsoar, moartr, mnartr,  nosotr );
+
+      //la qualite du triangle ns1 ns2 ns3
+      qutr2d_( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
+              qualite );
+
+      //la qualite moyenne
+      quamoy += qualite;
+
+      //la qualite minimale
+      quamin = Min( quamin, qualite );
+
+      //aire signee du triangle nt
+      d = surtd2_( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
+      if( d<0 )
+      {
+       //un triangle d'aire negative de plus
+       nbtrianeg++;
+       cout << "ATTENTION: le triangle " << nt << " de sommets:"
+            << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
+            << " a une aire " << d <<"<=0" << endl;
+      }
+
+      //aire des triangles actuels
+      aire += Abs(d);
+    }
+  }
+
+  //les affichages
+  quamoy /= nbtria;
+  cout << "Qualite moyenne=" << quamoy
+       << "  Qualite minimale=" << quamin
+       << " des " << nbtria << " triangles de surface totale="
+       << aire << endl;
+
+  if( nbtrianeg>0 )
+    MESSAGE( "ATTENTION: nombre de triangles d'aire negative=" << nbtrianeg );
+  return;
+}
diff --git a/src/MEFISTO2/aptrte.h b/src/MEFISTO2/aptrte.h

new file mode 100755 (executable)

index 0000000..9e57131
--- /dev/null
+++ b/src/MEFISTO2/aptrte.h
@@ -0,0 +1,229 @@
+#ifndef aptrte__h
+#define aptrte__h
+
+#include <limits.h>   // limites min max int long real ...
+#include <unistd.h>   // gethostname, ...
+#include <stdio.h>
+#include <iostream.h> // pour cout cin ...
+#include <iomanip.h>  // pour le format des io setw, stx, setfill, ...
+#include <string.h>   // pour les fonctions sur les chaines de caracteres
+#include <ctype.h>
+#include <stdlib.h>
+#include <math.h>     // pour les fonctions mathematiques
+#include <time.h>
+
+#include <sys/types.h>
+#include <sys/time.h>
+
+void qualitetrte( R3 *mnpxyd,
+                 Z & mosoar, Z & mxsoar, Z *mnsoar,
+                 Z & moartr, Z & mxartr, Z *mnartr,
+                 Z & nbtria, R & quamoy, R & quamin );
+// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// but :    calculer la qualite moyenne et minimale de la triangulation
+// -----    actuelle definie par les tableaux nosoar et noartr
+// entrees:
+// --------
+// mnpxyd : tableau des coordonnees 2d des points
+//          par point : x  y  distance_souhaitee
+// mosoar : nombre maximal d'entiers par arete et
+//          indice dans nosoar de l'arete suivante dans le hachage
+// mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
+//          attention: mxsoar>3*mxsomm obligatoire!
+// nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
+//          chainage des aretes frontalieres, chainage du hachage des aretes
+//          hachage des aretes = nosoar(1)+nosoar(2)*2
+//          avec mxsoar>=3*mxsomm
+//          une arete i de nosoar est vide <=> nosoar(1,i)=0 et
+//          nosoar(2,arete vide)=l'arete vide qui precede
+//          nosoar(3,arete vide)=l'arete vide qui suit
+// moartr : nombre maximal d'entiers par arete du tableau noartr
+// mxartr : nombre maximal de triangles declarables
+// noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
+//          arete1 = 0 si triangle vide => arete2 = triangle vide suivant
+// sorties:
+// --------
+// nbtria : nombre de triangles internes au domaine
+// quamoy : qualite moyenne  des triangles actuels
+// quamin : qualite minimale des triangles actuels
+// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+void  aptrte( Z nutysu, R aretmx,
+             Z nblf,   Z *nudslf, R2 *uvslf,
+             Z nbpti,  R2 *uvpti,
+             Z & nbst, R2 * & uvst, Z & nbt, Z * & nust,
+             Z & ierr );
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// but : appel de la triangulation par un arbre-4 recouvrant
+// ----- de triangles equilateraux
+//       le contour du domaine plan est defini par des lignes fermees
+//       la premiere ligne etant l'enveloppe de toutes les autres
+//       la fonction areteideale_(s,d) donne la taille d'arete
+//       au point s dans la direction d (direction inactive pour l'instant)
+//       des lors toute arete issue d'un sommet s devrait avoir une longueur
+//       comprise entre 0.65 areteideale_(s,d) et 1.3 areteideale_(s,d)
+//
+//Attention:
+//  Les tableaux uvslf et uvpti sont supposes ne pas avoir de sommets identiques!
+//  De meme, un sommet d'une ligne fermee ne peut appartenir a une autre ligne fermee
+//
+// entrees:
+// --------
+// nutysu : numero de traitement de areteideale_() selon le type de surface
+//          0 pas d'emploi de la fonction areteideale_() et aretmx est active
+//          1 il existe une fonction areteideale_(s,d)
+//            dont seules les 2 premieres composantes de uv sont actives
+//          ... autres options a definir ...
+// aretmx : longueur maximale des aretes de la future triangulation
+// nblf   : nombre de lignes fermees de la surface
+// nudslf : numero du dernier sommet de chacune des nblf lignes fermees
+//          nudslf(0)=0 pour permettre la difference sans test
+//          Attention le dernier sommet de chaque ligne est raccorde au premier
+//          tous les sommets et les points internes ont des coordonnees
+//          UV differentes <=> Pas de point double!
+// uvslf  : uv des nudslf(nblf) sommets des lignes fermees
+// nbpti  : nombre de points internes futurs sommets de la triangulation
+// uvpti  : uv des points internes futurs sommets de la triangulation
+//
+// sorties:
+// --------
+// nbst   : nombre de sommets de la triangulation finale
+// uvst   : coordonnees uv des nbst sommets de la triangulation
+// nbt    : nombre de triangles de la triangulation finale
+// nust   : 3 numeros dans uvst des sommets des nbt triangles
+// ierr   : 0 si pas d'erreur
+//        > 0 sinon
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// auteur : Alain Perronnet  Analyse Numerique Paris UPMC   decembre 2001
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+extern "C" {  void tempscpu_( double & tempsec );  }
+//Retourne le temps CPU utilise en secondes
+
+extern "C" { void deltacpu_( R & dtcpu ); }
+//Retourne le temps CPU utilise en secondes depuis le precedent appel
+
+//initialiser le tableau mnsoar pour le hachage des aretes
+extern "C" {void insoar_( Z & mxsomm, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar );}
+
+//mettre a zero les nb entiers de tab
+extern "C" {void azeroi_( Z & nb, Z * tab );}
+
+extern "C" {void fasoar_( Z & ns1, Z & ns2, Z & nt1, Z & nt2, Z & nolign,
+                         Z & mosoar,  Z & mxsoar,  Z & n1soar,  Z * mnsoar,  Z * mnarst,
+                         Z & noar, Z & ierr );}
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// but :    former l'arete de sommet ns1-ns2 dans le hachage du tableau
+// -----    nosoar des aretes de la triangulation
+// entrees:
+// --------
+// ns1 ns2: numero pxyd des 2 sommets de l'arete
+// nt1    : numero du triangle auquel appartient l'arete
+//          nt1=-1 si numero inconnu
+// nt2    : numero de l'eventuel second triangle de l'arete si connu
+//          nt2=-1 si numero inconnu
+// nolign : numero de la ligne fermee de l'arete
+//          =0 si l'arete n'est une arete de ligne
+//          ce numero est ajoute seulement si l'arete est creee
+// mosoar : nombre maximal d'entiers par arete du tableau nosoar
+// mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
+// modifies:
+// ---------
+// n1soar : numero de la premiere arete vide dans le tableau nosoar
+//          une arete i de nosoar est vide  <=>  nosoar(1,i)=0
+//          chainage des aretes vides amont et aval
+//          l'arete vide qui precede=nosoar(4,i)
+//          l'arete vide qui suit   =nosoar(5,i)
+// nosoar : numero des 2 sommets, no ligne, 2 triangles de l'arete,
+//          chainage momentan'e d'aretes, chainage du hachage des aretes
+//          hachage des aretes = min( nosoar(1), nosoar(2) )
+// noarst : noarst(np) numero d'une arete du sommet np
+
+// ierr   : si < 0  en entree pas d'affichage en cas d'erreur du type
+//         "arete appartenant a plus de 2 triangles et a creer!"
+//          si >=0  en entree       affichage de ce type d'erreur
+// sorties:
+// --------
+// noar   : >0 numero de l'arete retrouvee ou ajoutee
+// ierr   : =0 si pas d'erreur
+//          =1 si le tableau nosoar est sature
+//          =2 si arete a creer et appartenant a 2 triangles distincts
+//             des triangles nt1 et nt2
+//          =3 si arete appartenant a 2 triangles distincts
+//             differents des triangles nt1 et nt2
+//          =4 si arete appartenant a 2 triangles distincts
+//             dont le second n'est pas le triangle nt2
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
+extern "C" {void teajte_( Z & mxsomm, Z &  nbsomm, R3 * mnpxyd,  R3 * comxmi,
+                         R & aretmx,  Z & mxtree, Z * letree,
+                         Z & ierr );}
+
+
+extern "C" {void tehote_( Z & nutysu, Z & nbarpi, Z &  mxsomm, Z &  nbsomm, R3 * mnpxyd,
+                         R3 * comxmi, R & aretmx,
+                         Z * letree, Z & mxqueu, Z * mnqueu,
+                         Z & ierr );}
+// homogeneisation de l'arbre des te a un saut de taille au plus
+// prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
+
+extern "C" {void tetrte_(  R3 * comxmi, R & aretmx, Z & nbarpi, Z & mxsomm, R3 * mnpxyd,
+                          Z & mxqueu,  Z * mnqueu,  Z * mntree,
+                          Z & mosoar,  Z & mxsoar,  Z & n1soar, Z * mnsoar,
+                          Z & moartr, Z &  mxartr,  Z & n1artr,  Z * mnartr,  Z * mnarst,
+                          Z & ierr );}
+// trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
+// et des points de la frontiere, des points internes imposes interieurs
+
+extern "C" {void aisoar_( Z & mosoar, Z & mxsoar, Z * mnsoar, Z & na );}
+  // formation du chainage 6 des aretes internes a echanger eventuellement
+
+extern "C" {void tedela_( R3 * mnpxyd, Z * mnarst,
+                         Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar, Z & na,
+                         Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z & n );}
+  // boucle sur les aretes internes (non sur une ligne de la frontiere)
+  // avec echange des 2 diagonales afin de rendre la triangulation delaunay
+ 
+extern "C" {void terefr_( Z & nbarpi, R3 * mnpxyd,
+                         Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
+                         Z & moartr, Z & n1artr, Z * mnartr, Z * mnarst,
+                         Z & mxarcf, Z * mnarc1, Z * mnarc2,
+                         Z * mnarc3, Z * mnarc4,
+                         Z & n, Z & ierr );}
+// detection des aretes frontalieres initiales perdues
+// triangulation frontale pour les restaurer
+
+extern "C" {void tesuex_( Z & nblf, Z * nulftr,
+                         Z & ndtri0, Z & nbsomm, R3 * mnpxyd, Z * mnslig,
+                         Z & mosoar, Z & mxsoar, Z * mnsoar,
+                         Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr, Z * mnarst,
+                         Z & nbtria, Z * mntrsu, Z & ierr );}
+// suppression des triangles externes a la surface
+
+extern "C" {void teamqt_( Z & nutysu,
+                         Z * mnarst, Z & mosoar, Z & mxsoar, Z & n1soar, Z * mnsoar,
+                         Z & moartr, Z & mxartr, Z & n1artr, Z * mnartr,
+                         Z & mxarcf, Z * mntrcf, Z * mnstbo,
+                         Z * n1arcf, Z * mnarcf, Z * mnarc1,
+                         R3 * comxmi, Z & nbarpi, Z & nbsomm, Z & mxsomm,
+                         R3 * mnpxyd, Z * mnslig,
+                         Z & ierr );}
+// amelioration de la qualite de la triangulation par
+// barycentrage des sommets internes a la triangulation
+// suppression des aretes trop longues ou trop courtes
+// modification de la topologie des groupes de triangles
+// mise en delaunay de la triangulation
+ 
+extern "C" {void nusotr_( Z & nt, Z & mosoar, Z * mnsoar, Z & moartr, Z * mnartr,
+                         Z * nosotr );}
+//retrouver les numero des 3 sommets du triangle nt
+
+extern "C" {void qutr2d_( R3 & p1, R3 & p2, R3 & p3, R & qualite );}
+//calculer la qualite d'un triangle de R2 de sommets p1, p2, p3
+
+extern "C" { R surtd2_( R3 & p1, R3 & p2, R3 & p3 ); }
+//calcul de la surface d'un triangle defini par 3 points de r**2
+
+#endif
diff --git a/src/MEFISTO2/areteideale.f b/src/MEFISTO2/areteideale.f

new file mode 100755 (executable)

index 0000000..cabc8d4
--- /dev/null
+++ b/src/MEFISTO2/areteideale.f
@@ -0,0 +1,5 @@
+      double precision function areteideale( xyz, direction )
+      double precision xyz(3), direction(3)
+      areteideale = 10
+      return
+      end
diff --git a/src/MEFISTO2/trte.f b/src/MEFISTO2/trte.f

new file mode 100755 (executable)

index 0000000..f66104e
--- /dev/null
+++ b/src/MEFISTO2/trte.f
@@ -0,0 +1,8319 @@
+      subroutine qutr2d( p1, p2, p3, qualite )
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but :     calculer la qualite d'un triangle de r**2
+c -----     2 coordonnees des 3 sommets en double precision
+c
+c entrees :
+c ---------
+c p1,p2,p3 : les 3 coordonnees des 3 sommets du triangle
+c            sens direct pour une surface et qualite >0
+c sorties :
+c ---------
+c qualite: valeur de la qualite du triangle entre 0 et 1 (equilateral)
+c          1 etant la qualite optimale
+c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : alain perronnet analyse numerique upmc paris     janvier 1995
+c2345x7..............................................................012
+      parameter  ( d2uxr3 = 3.4641016151377544d0 )
+c                  d2uxr3 = 2 * sqrt(3)
+      double precision  p1(2), p2(2), p3(2), qualite, a, b, c, p
+c
+c     la longueur des 3 cotes
+      a = sqrt( (p2(1)-p1(1))**2 + (p2(2)-p1(2))**2 )
+      b = sqrt( (p3(1)-p2(1))**2 + (p3(2)-p2(2))**2 )
+      c = sqrt( (p1(1)-p3(1))**2 + (p1(2)-p3(2))**2 )
+c
+c     demi perimetre
+      p = (a+b+c) * 0.5d0
+c
+      if ( (a*b*c) .ne. 0d0 ) then
+c        critere : 2 racine(3) * rayon_inscrit / plus longue arete
+         qualite = d2uxr3 * sqrt( abs( (p-a) / p * (p-b) * (p-c) ) )
+     %          / max(a,b,c)
+      else
+         qualite = 0d0
+      endif
+c
+c
+c     autres criteres possibles:
+c     critere : 2 * rayon_inscrit / rayon_circonscrit
+c     qualite = 8d0 * (p-a) * (p-b) * (p-c) / (a * b * c)
+c
+c     critere : 3*sqrt(3.) * ray_inscrit / demi perimetre
+c     qualite = 3*sqrt(3.) * sqrt ((p-a)*(p-b)*(p-c) / p**3)
+c
+c     critere : 2*sqrt(3.) * ray_inscrit / max( des aretes )
+c     qualite = 2*sqrt(3.) * sqrt( (p-a)*(p-b)*(p-c) / p ) / max(a,b,c)
+      end
+
+
+      double precision function surtd2( p1 , p2 , p3 )
+c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but : calcul de la surface d'un triangle defini par 3 points de R**2
+c -----
+c parametres d entree :
+c ---------------------
+c p1 p2 p3 : les 3 fois 2 coordonnees des sommets du triangle
+c
+c parametre resultat :
+c --------------------
+c surtd2 : surface du triangle
+c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : alain perronnet analyse numerique upmc paris     fevrier 1992
+c2345x7..............................................................012
+      double precision  p1(2), p2(2), p3(2)
+c
+c     la surface du triangle
+      surtd2 = ( ( p2(1)-p1(1) ) * ( p3(2)-p1(2) )
+     %         - ( p2(2)-p1(2) ) * ( p3(1)-p1(1) ) ) * 0.5d0
+      end
+
+      integer function nopre3( i )
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but :   numero precedent i dans le sens circulaire  1 2 3 1 ...
+c -----
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
+c2345x7..............................................................012
+      if( i .eq. 1 ) then
+         nopre3 = 3
+      else
+         nopre3 = i - 1
+      endif
+      end
+
+      integer function nosui3( i )
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but :   numero suivant i dans le sens circulaire  1 2 3 1 ...
+c -----
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
+c2345x7..............................................................012
+      if( i .eq. 3 ) then
+         nosui3 = 1
+      else
+         nosui3 = i + 1
+      endif
+      end
+
+      subroutine provec( v1 , v2 , v3 )
+c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but :    v3 vecteur = produit vectoriel de 2 vecteurs de r ** 3
+c -----
+c entrees:
+c --------
+c v1, v2 : les 2 vecteurs de 3 composantes
+c
+c sortie :
+c --------
+c v3     : vecteur = v1  produit vectoriel v2
+cc++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : perronnet alain upmc analyse numerique paris        mars 1987
+c2345x7..............................................................012
+      double precision    v1(3), v2(3), v3(3)
+c
+      v3( 1 ) = v1( 2 ) * v2( 3 ) - v1( 3 ) * v2( 2 )
+      v3( 2 ) = v1( 3 ) * v2( 1 ) - v1( 1 ) * v2( 3 )
+      v3( 3 ) = v1( 1 ) * v2( 2 ) - v1( 2 ) * v2( 1 )
+c
+      return
+      end
+
+      subroutine norme1( n, v, ierr )
+c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but :   normalisation euclidienne a 1 d un vecteur v de n composantes
+c -----
+c entrees :
+c ---------
+c n       : nombre de composantes du vecteur
+c
+c modifie :
+c ---------
+c v       : le vecteur a normaliser a 1
+c
+c sortie  :
+c ---------
+c ierr    : 1 si la norme de v est egale a 0
+c           0 si pas d'erreur
+c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : alain perronnet analyse numerique paris             mars 1987
+c ......................................................................
+      double precision  v( n ), s, sqrt
+c
+      s = 0.0d0
+      do 10 i=1,n
+         s = s + v( i ) * v( i )
+   10 continue
+c
+c     test de nullite de la norme du vecteur
+c     --------------------------------------
+      if( s .le. 0.0d0 ) then
+c        norme nulle du vecteur non normalisable a 1
+         ierr = 1
+         return
+      endif
+c
+      s = 1.0d0 / sqrt( s )
+      do 20 i=1,n
+         v( i ) = v ( i ) * s
+   20 continue
+c
+      ierr = 0
+      end
+
+
+      subroutine insoar( mxsomm, mosoar, mxsoar, n1soar, nosoar )
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but :    initialiser le tableau nosoar pour le hachage des aretes
+c -----
+c
+c entrees:
+c --------
+c mxsomm : plus grand numero de sommet d'une arete au cours du calcul
+c mosoar : nombre maximal d'entiers par arete du tableau nosoar
+c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
+c          avec mxsoar>=3*mxsomm
+c
+c sorties:
+c --------
+c n1soar : numero de la premiere arete vide dans le tableau nosoar
+c          une arete i de nosoar est vide  <=>  nosoar(1,i)=0
+c          chainage des aretes vides amont et aval
+c          l'arete vide qui precede=nosoar(4,i)
+c          l'arete vide qui suit   =nosoar(5,i)
+c nosoar : numero des 2 sommets, no ligne, 2 triangles de l'arete,
+c          chainage momentan'e d'aretes, chainage du hachage des aretes
+c          hachage des aretes = min( nosoar(1), nosoar(2) )
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
+c2345x7..............................................................012
+      integer   nosoar(mosoar,mxsoar)
+c
+c     initialisation des aretes 1 a mxsomm
+      do 10 i=1,mxsomm
+c
+c        sommet 1 = 0 <=> temoin d'arete vide pour le hachage
+         nosoar( 1, i ) = 0
+c
+c        arete sur aucune ligne
+         nosoar( 3, i ) = 0
+c
+c        la position de l'arete interne ou frontaliere est inconnue
+         nosoar( 6, i ) = -2
+c
+c        fin de chainage du hachage pas d'arete suivante
+         nosoar( mosoar, i ) = 0
+c
+ 10   continue
+c
+c     la premiere arete vide chainee est la mxsomm+1 du tableau
+c     car ces aretes ne sont pas atteignables par le hachage direct
+      n1soar = mxsomm + 1
+c
+c     initialisation des aretes vides et des chainages
+      do 20 i = n1soar, mxsoar
+c
+c        sommet 1 = 0 <=> temoin d'arete vide pour le hachage
+         nosoar( 1, i ) = 0
+c
+c        arete sur aucune ligne
+         nosoar( 3, i ) = 0
+c
+c        chainage sur l'arete vide qui precede
+c        (si arete occupee cela deviendra le no du triangle 1 de l'arete)
+         nosoar( 4, i ) = i-1
+c
+c        chainage sur l'arete vide qui suit
+c        (si arete occupee cela deviendra le no du triangle 2 de l'arete)
+         nosoar( 5, i ) = i+1
+c
+c        chainages des aretes frontalieres ou internes ou ...
+         nosoar( 6, i ) = -2
+c
+c        fin de chainage du hachage
+         nosoar( mosoar, i ) = 0
+c
+ 20   continue
+c
+c     la premiere arete vide n'a pas de precedent
+      nosoar( 4, n1soar ) = 0
+c
+c     la derniere arete vide est mxsoar sans arete vide suivante
+      nosoar( 5, mxsoar ) = 0
+      end
+
+
+      subroutine azeroi ( l , ntab )
+c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but : initialisation a zero d un tableau ntab de l variables entieres
+c -----
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : alain perronnet analyse numerique upmc paris septembre 1988
+c23456---------------------------------------------------------------012
+      integer ntab(l)
+      do 1 i = 1 , l
+         ntab( i ) = 0
+    1 continue
+      end
+
+
+      subroutine fasoar( ns1,    ns2,    nt1,    nt2,    nolign,
+     %                   mosoar, mxsoar, n1soar, nosoar, noarst,
+     %                   noar,   ierr )
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but :    former l'arete de sommet ns1-ns2 dans le hachage du tableau
+c -----    nosoar des aretes de la triangulation
+c
+c entrees:
+c --------
+c ns1 ns2: numero pxyd des 2 sommets de l'arete
+c nt1    : numero du triangle auquel appartient l'arete
+c          nt1=-1 si numero inconnu
+c nt2    : numero de l'eventuel second triangle de l'arete si connu
+c          nt2=-1 si numero inconnu
+c nolign : numero de la ligne de l'arete dans ladefi(wulftr-1+nolign)
+c          =0 si l'arete n'est une arete de ligne
+c          ce numero est ajoute seulement si l'arete est creee
+c mosoar : nombre maximal d'entiers par arete du tableau nosoar
+c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
+c
+c modifies:
+c ---------
+c n1soar : numero de la premiere arete vide dans le tableau nosoar
+c          une arete i de nosoar est vide  <=>  nosoar(1,i)=0
+c          chainage des aretes vides amont et aval
+c          l'arete vide qui precede=nosoar(4,i)
+c          l'arete vide qui suit   =nosoar(5,i)
+c nosoar : numero des 2 sommets, no ligne, 2 triangles de l'arete,
+c          chainage momentan'e d'aretes, chainage du hachage des aretes
+c          hachage des aretes = min( nosoar(1), nosoar(2) )
+c noarst : noarst(np) numero d'une arete du sommet np
+c
+c ierr   : si < 0  en entree pas d'affichage en cas d'erreur du type
+c         "arete appartenant a plus de 2 triangles et a creer!"
+c          si >=0  en entree       affichage de ce type d'erreur
+c
+c sorties:
+c --------
+c noar   : >0 numero de l'arete retrouvee ou ajoutee
+c ierr   : =0 si pas d'erreur
+c          =1 si le tableau nosoar est sature
+c          =2 si arete a creer et appartenant a 2 triangles distincts
+c             des triangles nt1 et nt2
+c          =3 si arete appartenant a 2 triangles distincts
+c             differents des triangles nt1 et nt2
+c          =4 si arete appartenant a 2 triangles distincts
+c             dont le second n'est pas le triangle nt2
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
+c2345x7..............................................................012
+      common / unites / lecteu, imprim, nunite(30)
+      integer           nosoar(mosoar,mxsoar), noarst(*)
+      integer           nu2sar(2)
+c
+c     ajout eventuel de l'arete s1 s2 dans nosoar
+      nu2sar(1) = ns1
+      nu2sar(2) = ns2
+c
+c     hachage de l'arete de sommets nu2sar
+      call hasoar( mosoar, mxsoar, n1soar, nosoar, nu2sar, noar )
+c     en sortie: noar>0 => no arete retrouvee
+c                    <0 => no arete ajoutee
+c                    =0 => saturation du tableau nosoar
+c
+      if( noar .eq. 0 ) then
+c
+c        saturation du tableau nosoar
+         write(imprim,*) 'fasoar: tableau nosoar sature'
+         ierr = 1
+         return
+c
+      else if( noar .lt. 0 ) then
+c
+c        l'arete a ete ajoutee. initialisation des autres informations
+         noar = -noar
+c        le numero de la ligne de l'arete
+         nosoar(3,noar) = nolign
+c        le triangle 1 de l'arete => le triangle nt1
+         nosoar(4,noar) = nt1
+c        le triangle 2 de l'arete => le triangle nt2
+         nosoar(5,noar) = nt2
+c
+c        le sommet appartient a l'arete noar
+         noarst( nu2sar(1) ) = noar
+         noarst( nu2sar(2) ) = noar
+c
+      else
+c
+c        l'arete a ete retrouvee.
+c        si elle appartient a 2 triangles differents de nt1 et nt2
+c        alors il y a une erreur
+         if( nosoar(4,noar) .gt. 0 .and.
+     %       nosoar(5,noar) .gt. 0 ) then
+             if( nosoar(4,noar) .ne. nt1 .and.
+     %           nosoar(4,noar) .ne. nt2 .or.
+     %           nosoar(5,noar) .ne. nt1 .and.
+     %           nosoar(5,noar) .ne. nt2 ) then
+c                arete appartenant a plus de 2 triangles => erreur
+                 if( ierr .ge. 0 ) then
+                    write(imprim,*) 'erreur fasoar: arete ',noar,
+     %              ' dans 2 triangles et a creer!'
+                 endif
+                 ierr = 2
+                 return
+             endif
+         endif
+c
+c        mise a jour du numero des triangles de l'arete noar
+c        le triangle 2 de l'arete => le triangle nt1
+         if( nosoar(4,noar) .lt. 0 ) then
+c            pas de triangle connu pour cette arete
+             n = 4
+         else
+c            deja un triangle connu. ce nouveau est le second
+             if( nosoar(5,noar) .gt. 0  .and.  nt1 .gt. 0 .and.
+     %          nosoar(5,noar) .ne. nt1 ) then
+c               arete appartenant a plus de 2 triangles => erreur
+                write(imprim,*) 'erreur fasoar: arete ',noar,
+     %          ' dans plus de 2 triangles'
+                ierr = 3
+                return
+             endif
+             n = 5
+         endif
+         nosoar(n,noar) = nt1
+c
+c        cas de l'arete frontaliere retrouvee comme diagonale d'un quadrangle
+         if( nt2 .gt. 0 ) then
+c           l'arete appartient a 2 triangles
+            if( nosoar(5,noar) .gt. 0  .and.
+     %          nosoar(5,noar) .ne. nt2 ) then
+c               arete appartenant a plus de 2 triangles => erreur
+                write(imprim,*) 'erreur fasoar: arete ',noar,
+     %         ' dans plus de 2 triangles'
+                ierr = 4
+                return
+            endif
+            nosoar(5,noar) = nt2
+         endif
+c
+      endif
+c
+c     pas d'erreur
+      ierr = 0
+      end
+
+      subroutine fq1inv( x, y, s, xc, yc, ierr )
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but :   calcul des 2 coordonnees (xc,yc) dans le carre (0,1)
+c -----   image par f:carre unite-->quadrangle appartenant a q1**2
+c         par une resolution directe due a nicolas thenault
+c
+c entrees:
+c --------
+c x,y   : coordonnees du point image dans le quadrangle de sommets s
+c s     : les 2 coordonnees des 4 sommets du quadrangle
+c
+c sorties:
+c --------
+c xc,yc : coordonnees dans le carre dont l'image par f vaut (x,y)
+c ierr  : 0 si calcul sans erreur, 1 si quadrangle degenere
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteurs: thenault tulenew  analyse numerique paris        janvier 1998
+c modifs : perronnet alain   analyse numerique paris        janvier 1998
+c234567..............................................................012
+      real             s(1:2,1:4), dist(2)
+      double precision a,b,c,d,alpha,beta,gamma,delta,x0,y0,t(2),u,v,w
+c
+      a = s(1,1)
+      b = s(1,2) - s(1,1)
+      c = s(1,4) - s(1,1)
+      d = s(1,1) - s(1,2) + s(1,3) - s(1,4)
+c
+      alpha = s(2,1)
+      beta  = s(2,2) - s(2,1)
+      gamma = s(2,4) - s(2,1)
+      delta = s(2,1) - s(2,2) + s(2,3) - s(2,4)
+c
+      u = beta  * c - b * gamma
+      if( u .eq. 0 ) then
+c        quadrangle degenere
+         ierr = 1
+         return
+      endif
+      v = delta * c - d * gamma
+      w = b * delta - beta * d
+c
+      x0 = c * (y-alpha) - gamma * (x-a)
+      y0 = b * (y-alpha) - beta  * (x-a)
+c
+      a = v  * w
+      b = u  * u - w * x0 - v * y0
+      c = x0 * y0
+c
+      if( a .ne. 0 ) then
+c
+         delta = sqrt( b*b-4*a*c )
+         if( b .ge. 0.0 ) then
+            t(2) = -b - delta
+         else
+            t(2) = -b + delta
+         endif
+c        la racine de plus grande valeur absolue
+c       (elle donne le plus souvent le point exterieur au carre unite
+c        donc a tester en second pour reduire les calculs)
+         t(2) = t(2) / ( 2 * a )
+c        calcul de la seconde racine a partir de la somme => plus stable
+         t(1) = - b/a - t(2)
+c
+         do 10 i=1,2
+c
+c           la solution i donne t elle un point interne au carre unite?
+            xc = ( x0 - v * t(i) ) / u
+            yc = ( w * t(i) - y0 ) / u
+            if( 0.0 .le. xc .and. xc .le. 1.0 ) then
+               if( 0.0 .le. yc .and. yc .le. 1.0 ) goto 9000
+            endif
+c
+c           le point (xc,yc) n'est pas dans le carre unite
+c           cela peut etre du aux erreurs d'arrondi
+c           => choix par le minimum de la distance aux bords du carre
+            dist(i) = max( 0.0, -xc, xc-1.0, -yc, yc-1.0 )
+c
+ 10      continue
+c
+         if( dist(1) .gt. dist(2) ) then
+c           f(xc,yc) pour la racine 2 est plus proche de x,y
+c           xc yc sont deja calcules
+            goto 9000
+         endif
+c
+      else if ( b .ne. 0 ) then
+         t(1) = - c / b
+      else
+         t(1) = 0
+      endif
+c
+c     les 2 coordonnees du point dans le carre unite
+      xc = ( x0 - v * t(1) ) / u
+      yc = ( w * t(1) - y0 ) / u
+c
+ 9000 ierr = 0
+      return
+      end
+
+
+      subroutine ptdatr( point, pxyd, nosotr, nsigne )
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but :    le point est il dans le triangle de sommets nosotr
+c -----
+c
+c entrees:
+c --------
+c point  : les 2 coordonnees du point
+c pxyd   : les 2 coordonnees et distance souhaitee des points du maillage
+c nosotr : le numero des 3 sommets du triangle
+c
+c sorties:
+c --------
+c nsigne : >0 si le point est dans le triangle ou sur une des 3 aretes
+c          =0 si le triangle est degenere ou indirect ou ne contient pas le poin
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
+c....................................................................012
+      integer           nosotr(3)
+      double precision  point(2), pxyd(3,*)
+      double precision  xp,yp, x1,x2,x3, y1,y2,y3, d,dd, cb1,cb2,cb3
+c
+      xp = point( 1 )
+      yp = point( 2 )
+c
+      n1 = nosotr( 1 )
+      x1 = pxyd( 1 , n1 )
+      y1 = pxyd( 2 , n1 )
+c
+      n2 = nosotr( 2 )
+      x2 = pxyd( 1 , n2 )
+      y2 = pxyd( 2 , n2 )
+c
+      n3 = nosotr( 3 )
+      x3 = pxyd( 1 , n3 )
+      y3 = pxyd( 2 , n3 )
+c
+c     2 fois la surface du triangle = determinant de la matrice
+c     de calcul des coordonnees barycentriques du point p
+      d  = ( x2 - x1 ) * ( y3 - y1 ) - ( x3 - x1 ) * ( y2 - y1 )
+c
+      if( d .gt. 0 ) then
+c
+c        triangle non degenere
+c        =====================
+c        calcul des 3 coordonnees barycentriques du
+c        point xp yp dans le triangle
+         cb1 = ( ( x2-xp ) * ( y3-yp ) - ( x3-xp ) * ( y2-yp ) ) / d
+         cb2 = ( ( x3-xp ) * ( y1-yp ) - ( x1-xp ) * ( y3-yp ) ) / d
+         cb3 = 1d0 - cb1 -cb2
+ccc         cb3 = ( ( x1-xp ) * ( y2-yp ) - ( x2-xp ) * ( y1-yp ) ) / d
+c
+ccc         if( cb1 .ge. -0.00005d0 .and. cb1 .le. 1.00005d0 .and.
+         if( cb1 .ge. 0d0 .and. cb1 .le. 1d0 .and.
+     %       cb2 .ge. 0d0 .and. cb2 .le. 1d0 .and.
+     %       cb3 .ge. 0d0 .and. cb3 .le. 1d0 ) then
+c
+c           le triangle nosotr contient le point
+            nsigne = 1
+         else
+            nsigne = 0
+         endif
+c
+      else
+c
+c        triangle degenere
+c        =================
+c        le point est il du meme cote que le sommet oppose de chaque arete?
+         nsigne = 0
+         do 10 i=1,3
+c           le sinus de l'angle p1 p2-p1 point
+            x1  = pxyd(1,n1)
+            y1  = pxyd(2,n1)
+            d   = ( pxyd(1,n2) - x1 ) * ( point(2) - y1 )
+     %          - ( pxyd(2,n2) - y1 ) * ( point(1) - x1 )
+            dd  = ( pxyd(1,n2) - x1 ) * ( pxyd(2,n3) - y1 )
+     %          - ( pxyd(2,n2) - y1 ) * ( pxyd(1,n3) - x1 )
+            cb1 = ( pxyd(1,n2) - x1 ) ** 2
+     %          + ( pxyd(2,n2) - y1 ) ** 2
+            cb2 = ( point(1) - x1 ) ** 2
+     %          + ( point(2) - y1 ) ** 2
+            cb3 = ( pxyd(1,n3) - x1 ) ** 2
+     %          + ( pxyd(2,n3) - y1 ) ** 2
+            if( abs( dd ) .le. 1e-4 * sqrt( cb1 * cb3 ) ) then
+c              le point 3 est sur l'arete 1-2
+c              le point doit y etre aussi
+               if( abs( d ) .le. 1e-4 * sqrt( cb1 * cb2 ) ) then
+c                 point sur l'arete
+                  nsigne = nsigne + 1
+               endif
+            else
+c              le point 3 n'est pas sur l'arete . test des signes
+               if( d * dd .ge. 0 ) then
+                  nsigne = nsigne + 1
+               endif
+            endif
+c           permutation circulaire des 3 sommets et aretes
+            n  = n1
+            n1 = n2
+            n2 = n3
+            n3 = n
+ 10      continue
+         if( nsigne .ne. 3 ) nsigne = 0
+      endif
+      end
+
+      integer function nosstr( p, pxyd, nt, letree )
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but :    calculer le numero 0 a 3 du sous-triangle te contenant
+c -----    le point p
+c
+c entrees:
+c --------
+c p      : point de r**2 contenu dans le te nt de letree
+c pxyd   : x y distance des points
+c nt     : numero letree du te de te voisin a calculer
+c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
+c      letree(0,0)  no du 1-er te vide dans letree
+c      letree(0,1) : maximum du 1-er indice de letree (ici 8)
+c      letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
+c      letree(0:8,1) : racine de l'arbre  (triangle sans sur triangle)
+c      si letree(0,.)>0 alors
+c         letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
+c      sinon
+c         letree(0:3,j) :-no pxyd des 1 \85a 4 points internes au triangle j
+c                         0  si pas de point
+c                       ( j est alors une feuille de l'arbre )
+c      letree(4,j) : no letree du sur-triangle du triangle j
+c      letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
+c      letree(6:8,j) : no pxyd des 3 sommets du triangle j
+c
+c sorties :
+c ---------
+c nosstr : 0 si le sous-triangle central contient p
+c          i =1,2,3 numero du sous-triangle contenant p
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
+c2345x7..............................................................012
+      integer           letree(0:8,0:*)
+      double precision  pxyd(3,*), p(2),
+     %                  x1, y1, x21, y21, x31, y31, d, xe, ye
+c
+c     le numero des 3 sommets du triangle
+      ns1 = letree( 6, nt )
+      ns2 = letree( 7, nt )
+      ns3 = letree( 8, nt )
+c
+c     les coordonnees entre 0 et 1 du point p
+      x1  = pxyd(1,ns1)
+      y1  = pxyd(2,ns1)
+c
+      x21 = pxyd(1,ns2) - x1
+      y21 = pxyd(2,ns2) - y1
+c
+      x31 = pxyd(1,ns3) - x1
+      y31 = pxyd(2,ns3) - y1
+c
+      d   = 1.0 / ( x21 * y31 - x31 * y21 )
+c
+      xe  = ( ( p(1) - x1 ) * y31 - ( p(2) - y1 ) * x31 ) * d
+      ye  = ( ( p(2) - y1 ) * x21 - ( p(1) - x1 ) * y21 ) * d
+c
+      if( xe .gt. 0.5d0 ) then
+c        sous-triangle droit
+         nosstr = 2
+      else if( ye .gt. 0.5d0 ) then
+c        sous-triangle haut
+         nosstr = 3
+      else if( xe+ye .lt. 0.5d0 ) then
+c        sous-triangle gauche
+         nosstr = 1
+      else
+c        sous-triangle central
+         nosstr = 0
+      endif
+      end
+
+
+      integer function notrpt( p, pxyd, notrde, letree )
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but :    calculer le numero letree du sous-triangle feuille contenant
+c -----    le point p a partir du te notrde de letree
+c
+c entrees:
+c --------
+c p      : point de r**2 contenu dans le te nt de letree
+c pxyd   : x y distance des points
+c notrde : numero letree du triangle depart de recherche (1=>racine)
+c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
+c      letree(0,0)  no du 1-er te vide dans letree
+c      letree(0,1) : maximum du 1-er indice de letree (ici 8)
+c      letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
+c      letree(0:8,1) : racine de l'arbre  (triangle sans sur triangle)
+c      si letree(0,.)>0 alors
+c         letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
+c      sinon
+c         letree(0:3,j) :-no pxyd des 1 \85 4 points internes au triangle j
+c                         0  si pas de point
+c                        ( j est alors une feuille de l'arbre )
+c      letree(4,j) : no letree du sur-triangle du triangle j
+c      letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
+c      letree(6:8,j) : no pxyd des 3 sommets du triangle j
+c
+c sorties :
+c ---------
+c notrpt : numero letree du triangle contenant le point p
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
+c2345x7..............................................................012
+      integer           letree(0:8,0:*)
+      double precision  pxyd(1:3,*), p(2)
+c
+c     la racine depart de la recherche
+      notrpt = notrde
+c
+c     tant que la feuille n'est pas atteinte descendre l'arbre
+ 10   if( letree(0,notrpt) .gt. 0 ) then
+c
+c        recherche du sous-triangle contenant p
+         nsot = nosstr( p, pxyd, notrpt, letree )
+c
+c        le numero letree du sous-triangle
+         notrpt = letree( nsot, notrpt )
+         goto 10
+c
+      endif
+      end
+
+
+      subroutine teajpt( ns,   nbsomm, mxsomm, pxyd, letree,
+     &                   ntrp, ierr )
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but :    ajout du point ns de pxyd dans letree
+c -----
+c
+c entrees:
+c --------
+c ns     : numero du point a ajouter dans letree
+c mxsomm : nombre maximal de points declarables dans pxyd
+c pxyd   : tableau des coordonnees des points
+c          par point : x  y  distance_souhaitee
+c
+c modifies :
+c ----------
+c nbsomm : nombre actuel de points dans pxyd
+c
+c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
+c      letree(0,0) : no du 1-er te vide dans letree
+c      letree(0,1) : maximum du 1-er indice de letree (ici 8)
+c      letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
+c      letree(0:8,1) : racine de l'arbre  (triangle sans sur triangle)
+c      si letree(0,.)>0 alors
+c         letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
+c      sinon
+c         letree(0:3,j) :-no pxyd des 1 \85a 4 points internes au triangle j
+c                         0  si pas de point
+c                        ( j est alors une feuille de l'arbre )
+c      letree(4,j) : no letree du sur-triangle du triangle j
+c      letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
+c      letree(6:8,j) : no pxyd des 3 sommets du triangle j
+c
+c sorties :
+c ---------
+c ntrp    : numero letree du triangle te ou a ete ajoute le point
+c ierr    : 0 si pas d'erreur,  51 saturation letree, 52 saturation pxyd
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
+c2345x7..............................................................012
+      integer           letree(0:8,0:*)
+      double precision  pxyd(3,mxsomm)
+c
+c     depart de la racine
+      ntrp = 1
+c
+c     recherche du triangle contenant le point pxyd(ns)
+ 1    ntrp = notrpt( pxyd(1,ns), pxyd, ntrp, letree )
+c
+c     existe t il un point libre
+      do 10 i=0,3
+         if( letree(i,ntrp) .eq. 0 ) then
+c           la place i est libre
+            letree(i,ntrp) = -ns
+            return
+         endif
+ 10   continue
+c
+c     pas de place libre => 4 sous-triangles sont crees
+c                           a partir des 3 milieux des aretes
+      call te4ste( nbsomm, mxsomm, pxyd, ntrp, letree, ierr )
+      if( ierr .ne. 0 ) return
+c
+c     ajout du point ns
+      goto 1
+      end
+
+      subroutine n1trva( nt, lar, letree, notrva, lhpile )
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but :    calculer le numero letree du triangle voisin du te nt
+c -----    par l'arete lar (1 a 3 ) de nt
+c          attention : notrva n'est pas forcement minimal
+c
+c entrees:
+c --------
+c nt     : numero letree du te de te voisin a calculer
+c lar    : numero 1 a 3 de l'arete du triangle nt
+c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
+c   letree(0,0)  no du 1-er te vide dans letree
+c   letree(0,1) : maximum du 1-er indice de letree (ici 8)
+c   letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
+c   letree(0:8,1) : racine de l'arbre  (triangle sans sur-triangle)
+c   si letree(0,.)>0 alors
+c      letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
+c   sinon
+c      letree(0:3,j) :-no pxyd des 1 a 4 points internes au triangle j
+c                      0  si pas de point
+c                     ( j est alors une feuille de l'arbre )
+c   letree(4,j) : no letree du sur-triangle du triangle j
+c   letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
+c   letree(6:8,j) : no pxyd des 3 sommets du triangle j
+c
+c sorties :
+c ---------
+c notrva  : >0 numero letree du te voisin par l'arete lar
+c           =0 si pas de te voisin (racine , ... )
+c lhpile  : =0 si nt et notrva ont meme taille
+c           >0 nt est 4**lhpile fois plus petit que notrva
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
+c2345x7..............................................................012
+      integer   letree(0:8,0:*)
+      integer   lapile(1:64)
+c
+c     initialisation de la pile
+c     le triangle est empile
+      lapile(1) = nt
+      lhpile = 1
+c
+c     tant qu'il existe un sur-triangle
+ 10   ntr  = lapile( lhpile )
+      if( ntr .eq. 1 ) then
+c        racine atteinte => pas de triangle voisin
+         notrva = 0
+         lhpile = lhpile - 1
+         return
+      endif
+c
+c     le type du triangle ntr
+      nty  = letree( 5, ntr )
+c     l'eventuel sur-triangle
+      nsut = letree( 4, ntr )
+c
+      if( nty .eq. 0 ) then
+c
+c        triangle de type 0 => triangle voisin de type precedent(lar)
+c                              dans le sur-triangle de ntr
+c                              ce triangle remplace ntr dans lapile
+         lapile( lhpile ) = letree( nopre3(lar), nsut )
+         goto 20
+      endif
+c
+c     triangle ntr de type nty>0
+      if( nosui3(nty) .eq. lar ) then
+c
+c        le triangle voisin par lar est le triangle 0
+         lapile( lhpile ) = letree( 0, nsut )
+         goto 20
+      endif
+c
+c     triangle sans voisin direct => passage par le sur-triangle
+      if( nsut .eq. 0 ) then
+c
+c        ntr est la racine => pas de triangle voisin par cette arete
+         notrva = 0
+         return
+      else
+c
+c        le sur-triangle est empile
+         lhpile = lhpile + 1
+         lapile(lhpile) = nsut
+         goto 10
+      endif
+c
+c     descente aux sous-triangles selon la meme arete
+ 20   notrva = lapile( lhpile )
+c
+ 30   lhpile = lhpile - 1
+      if( letree(0,notrva) .le. 0 ) then
+c        le triangle est une feuille de l'arbre 0 sous-triangle
+c        lhpile = nombre de differences de niveaux dans l'arbre
+         return
+      else
+c        le triangle a 4 sous-triangles
+         if( lhpile .gt. 0 ) then
+c
+c           bas de pile non atteint
+            nty  = letree( 5, lapile(lhpile) )
+            if( nty .eq. lar ) then
+c              l'oppose est suivant(nty) de notrva
+               notrva = letree( nosui3(nty) , notrva )
+            else
+c              l'oppose est precedent(nty) de notrva
+               notrva = letree( nopre3(nty) , notrva )
+            endif
+            goto 30
+         endif
+      endif
+c
+c     meme niveau dans l'arbre lhpile = 0
+      end
+
+
+      subroutine cenced( xy1, xy2, xy3, cetria, ierr )
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c but : calcul des coordonnees du centre du cercle circonscrit
+c ----- du triangle defini par ses 3 sommets de coordonnees
+c       xy1 xy2 xy3 ainsi que le carre du rayon de ce cercle
+c
+c entrees :
+c ---------
+c xy1 xy2 xy3 : les 2 coordonnees des 3 sommets du triangle
+c ierr   : <0  => pas d'affichage si triangle degenere
+c          >=0 =>       affichage si triangle degenere
+c
+c sortie :
+c --------
+c cetria : cetria(1)=abcisse  du centre
+c          cetria(2)=ordonnee du centre
+c          cetria(3)=carre du rayon   1d28 si triangle degenere
+c ierr   : 0 si triangle non degenere
+c          1 si triangle degenere
+c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c auteur : perronnet alain upmc analyse numerique paris        juin 1995
+c2345x7..............................................................012
+      parameter        (epsurf=1d-7)
+      common / unites / lecteu,imprim,nunite(30)
+      double precision  x1,y1,x21,y21,x31,y31,
+     %                  aire2,xc,yc,rot,
+     %                  xy1(2),xy2(2),xy3(2),cetria(3)
+c
+c     le calcul de 2 fois l'aire du triangle
+c     attention l'ordre des 3 sommets est direct ou non
+      x1  = xy1(1)
+      x21 = xy2(1) - x1
+      x31 = xy3(1) - x1
+c
+      y1  = xy1(2)
+      y21 = xy2(2) - y1
+      y31 = xy3(2) - y1
+c
+      aire2  = x21 * y31 - x31 * y21
+c
+c     recherche d'un test relatif peu couteux
+c     pour reperer la degenerescence du triangle
+      if( abs(aire2) .le.
+     %    epsurf*(abs(x21)+abs(x31))*(abs(y21)+abs(y31)) ) then
+c        triangle de qualite trop faible
+         if( ierr .ge. 0 ) then
+c            nblgrc(nrerr) = 1
+c            kerr(1) = 'erreur cenced: triangle degenere'
+c            call lereur
+            write(imprim,*) 'erreur cenced: triangle degenere'
+            write(imprim,10000)  xy1,xy2,xy3,aire2
+         endif
+10000 format( 3(' x=',g24.16,' y=',g24.16/),' aire*2=',g24.16)
+         cetria(1) = 0d0
+         cetria(2) = 0d0
+         cetria(3) = 1d28
+         ierr = 1
+         return
+      endif
+c
+c     les 2 coordonnees du centre intersection des 2 mediatrices
+c     x = (x1+x2)/2 + lambda * (y2-y1)
+c     y = (y1+y2)/2 - lambda * (x2-x1)
+c     x = (x1+x3)/2 + rot    * (y3-y1)
+c
author	admin <salome-admin@opencascade.com>
	Mon, 19 May 2003 14:07:01 +0000 (14:07 +0000)
committer	admin <salome-admin@opencascade.com>
	Mon, 19 May 2003 14:07:01 +0000 (14:07 +0000)