pmtl usage

implicit problem read arh file instead of dat files
memory error warning in dofvector

AMDiS/AMDISUse.cmake

@@ -14,11 +14,17 @@ if(AMDIS_HAS_PARALLEL_DOMAIN)
 			list(APPEND AMDIS_LIBRARIES ${PETSC_LIBRARY_SYS} ${PETSC_LIBRARIES})
 			list(APPEND AMDIS_INCLUDE_DIRS ${PETSC_INCLUDES})
 		else()
-			message(FATAL_ERROR "Could nit find PETSc!")
+			message(FATAL_ERROR "Could not find PETSc!")
 		endif(PETSC_FOUND)
 	elseif(AMDIS_HAS_PARALLEL_DOMAIN STREQUAL "PMTL")
 		find_package(MTL REQUIRED)
 		list(APPEND AMDIS_LIBRARIES ${MTL_LIBRARIES})
+		find_library(PARMETIS_LIB parmetis)
+		if(PARMETIS_LIB)
+			list(APPEND AMDIS_LIBRARIES ${PARMETIS_LIB})
+		else(PARMETIS_LIB)
+			message(ERROR "could not find the parmetis libraries needed by amdis")
+		endif(PARMETIS_LIB)
 	endif()
 endif(AMDIS_HAS_PARALLEL_DOMAIN)
 

AMDiS/src/AdaptInstationary.cc

@@ -19,8 +19,10 @@
 
 #if HAVE_PARALLEL_DOMAIN_AMDIS
 #include "parallel/MeshDistributor.h"
+#ifndef HAVE_PARALLEL_MTL4
 #include <petsc.h>
 #endif
+#endif
 
 namespace AMDiS {
 

AMDiS/src/AdaptStationary.cc

@@ -18,8 +18,10 @@
 
 #if HAVE_PARALLEL_DOMAIN_AMDIS
 #include "parallel/MeshDistributor.h"
+#ifndef HAVE_PARALLEL_MTL4
 #include <petsc.h>
 #endif
+#endif
 
 namespace AMDiS {
 

AMDiS/src/DOFVector.cc

@@ -68,7 +68,7 @@ namespace AMDiS {
 
 
   template<>
-  const double& DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* values) const
+  const double DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* values) const
   {
     FUNCNAME("DOFVector<double>::evalAtCoords()");
   
@@ -114,7 +114,7 @@ namespace AMDiS {
 
 
   template<>
-  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* values) const
+  const WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* values) const
   {
     FUNCNAME("DOFVector<double>::evalAtCoords()");
   

AMDiS/src/DOFVector.h

@@ -542,7 +542,7 @@ namespace AMDiS {
 
     /// eval DOFVector at given point p. If oldElInfo != NULL the search for the element, where p is inside,
     /// starts from oldElInfo.
-    const T& evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo = NULL, T* value = NULL) const;
+    const T evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo = NULL, T* value = NULL) const;
 
     /// Writes the data of the DOFVector to an output stream.
     void serialize(std::ostream &out) 
@@ -582,8 +582,8 @@ namespace AMDiS {
   }; 
 
 
-  template<>
-  const double& DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
+  /*template<>
+  const double DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
 
   template<>
   const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* value) const;
@@ -592,7 +592,7 @@ namespace AMDiS {
   void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 
 
   template<>
-  void DOFVector<double>::coarseRestrict(RCNeighbourList&, int); 
+  void DOFVector<double>::coarseRestrict(RCNeighbourList&, int); */
 
   inline double min(const DOFVector<double>& v) 
   {

AMDiS/src/MTL4Solver.h

@@ -55,7 +55,7 @@ namespace AMDiS {
 
       error = worker.solve(matrix ,mtl_x, mtl_b);
 
-      mtl::dense_vector<typename MTLMatrix::value_type> r(mtl_b); 
+      MTLVector r(mtl_b); 
       r -= matrix * mtl_x; 
       double residual = two_norm(r);
       MSG("MTL4Solver: ||b-Ax||= %e\n", residual);

AMDiS/src/ProblemImplicit.cc

@@ -12,14 +12,15 @@
 
 #include "ProblemImplicit.h"
 #include "MathFunctions.h"
-
+#include "io/ArhReader.h"
 using namespace std;
 
 namespace AMDiS {
 
-  void ProblemImplicitBase::readDofVec(std::istream& in, DOFVector<double>* vec, 
+  void ProblemImplicit::readDofVec(const std::string& filename, DOFVector<double>* vec, 
                                        Mesh* mesh) 
   {
+#if 0
     long size = mesh->getNumberOfVertices();
 
     TEST_EXIT(vec->getSize()>=size)("dofvector smaller than vertex data vector\n");
@@ -34,13 +35,15 @@ namespace AMDiS {
       in >> buf;
       (*vec)[i] = atof(buf.c_str());
     }
+#endif
+    ArhReader::read(filename, vec->getFeSpace()->getMesh(), vec);
   }
 
 
-  void ProblemImplicitBase::readR(std::istream& inStream, double eps, 
+  void ProblemImplicit::readR(const std::string& inStream, double eps, 
 	                          Mesh* mesh, int implMesh , int comp) 
   {
-    FUNCNAME("ProblemImplicitBase::readR()");
+    FUNCNAME("ProblemImplicit::readR()");
 
     DOFVector<double>* r = getSignedDistance(implMesh, comp);
     DOFVector<double>* phi1 = getPhi1(implMesh, comp);
@@ -67,7 +70,7 @@ namespace AMDiS {
   }
 
 
-  void ProblemImplicitBase::readPhi1(istream& inStream, double eps, 
+  void ProblemImplicit::readPhi1(const std::string& inStream, double eps, 
                                      Mesh* mesh, int implMesh, int comp)
   {
     DOFVector<double>* r = getSignedDistance(implMesh, comp);
@@ -95,7 +98,7 @@ namespace AMDiS {
   }
 
 
-  void ProblemImplicitBase::readPhi2(istream& inStream, double eps, 
+  void ProblemImplicit::readPhi2(const std::string& inStream, double eps, 
 	                             Mesh* mesh, int implMesh, int comp)
   { 
     DOFVector<double>* r = getSignedDistance(implMesh, comp);
@@ -122,21 +125,27 @@ namespace AMDiS {
   }
   
 
-  std::string ProblemImplicitBase::getDofFilename(std::string path, int implMesh)
+  std::string ProblemImplicit::getDofFilename(std::string path, int implMesh)
   {
+    FUNCNAME("ProblemImplicit::getDofFilename()");
     std::string suffix = "[" + 
       boost::lexical_cast< std::string >(implMesh) +
       "]";
 
     std::string dofFilename("");
     Parameters::get(path + "dof file" + suffix, dofFilename);
-    if (implMesh == 0 && dofFilename.length() == 0)
+    TEST_EXIT(dofFilename.length() == 0)("dof file was changed to \"arh file\" and reads an arh file now");
+    Parameters::get(path + "arh file" + suffix, dofFilename);
+    if (implMesh == 0 && dofFilename.length() == 0) {
       Parameters::get(path + "dof file", dofFilename);
+      TEST_EXIT(dofFilename.length() == 0)("dof file was changed to \"arh file\" and reads an arh file now");
+      Parameters::get(path + "arh file", dofFilename);
+    }
     return dofFilename;    
   }
 
 
-  double ProblemImplicitBase::getEpsilon(std::string path, int implMesh)
+  double ProblemImplicit::getEpsilon(std::string path, int implMesh)
   {
     std::string suffix = "[" + 
       boost::lexical_cast< std::string >(implMesh) +
@@ -150,7 +159,7 @@ namespace AMDiS {
   }
 
 
-  int ProblemImplicitBase::getType(std::string path, int implMesh)
+  int ProblemImplicit::getType(std::string path, int implMesh)
   {
     std::string suffix = "[" + 
       boost::lexical_cast< std::string >(implMesh) +
@@ -215,17 +224,16 @@ namespace AMDiS {
   {
     FUNCNAME("ProblemImplicit::createImplicitMesh()");
 
-    if (!isImplicitMesh(p))
-      return false;
     std::string path = name + "->implicit mesh[" 
       + boost::lexical_cast< std::string >(p) + "]->";
     int nImplMeshes(0);
     Parameters::get(path + "nr meshes", nImplMeshes);
-
-    r[p].resize(nImplMeshes);
-    phi1[p].resize(nImplMeshes);
-    phi2[p].resize(nImplMeshes);
-    levelSet[p].resize(nImplMeshes);
+    if (nImplMeshes == 0)
+      return false;
+    r[p].resize(nImplMeshes, NULL);
+    phi1[p].resize(nImplMeshes, NULL);
+    phi2[p].resize(nImplMeshes, NULL);
+    levelSet[p].resize(nImplMeshes, NULL);
 
     for ( int i = 0; i < nImplMeshes ; ++i ) {
       (r[p])[i] = new DOFVector< double >(getFeSpace(p), "r");
@@ -256,21 +264,19 @@ namespace AMDiS {
 
     TEST_EXIT(meshes[comp] != NULL)("the mesh was not created\n");
     
-    std::ifstream inStream(dofFilename.c_str());
     switch (serType) {
     case 0:
-      readR(inStream, eps, meshes[comp], implMesh, comp);
+      readR(dofFilename, eps, meshes[comp], implMesh, comp);
       break;
     case 1:
-      readPhi1(inStream, eps, meshes[comp], implMesh, comp);
+      readPhi1(dofFilename, eps, meshes[comp], implMesh, comp);
       break;
     case 2:
-      readPhi2(inStream, eps, meshes[comp], implMesh, comp);
+      readPhi2(dofFilename, eps, meshes[comp], implMesh, comp);
       break;
     default:
       WARNING("unknown implicit mesh type\n");
     }
-    inStream.close();
     return true;
   }
 
@@ -279,8 +285,8 @@ namespace AMDiS {
   {
     FUNCNAME("ProblemImplicit::createMesh()");
 
-    ProblemStatSeq::createMesh();
-
+    ProblemStat::createMesh();
+#if 0
     implMesh.resize(meshes.size());
     for (unsigned int i = 0 ; i < meshes.size() ; ++i ) {
       std::string path = name + "->implicit mesh[" +
@@ -299,6 +305,7 @@ namespace AMDiS {
       } else
 	implMesh[i] = false;
     }
+#endif    
   }
 
 
@@ -308,7 +315,27 @@ namespace AMDiS {
   {
     FUNCNAME("ProblemImplicit::initialize()");
 
-    ProblemStatSeq::initialize(initFlag);    
+    ProblemStat::initialize(initFlag);
+    if ( initFlag.isSet(INIT_IMPLICIT_MESH) )
       createImplicitMesh();
   }
+
+  ProblemImplicit::~ProblemImplicit()
+  {
+    for ( unsigned int p(0); p < meshes.size(); ++p ) {
+      for ( unsigned int i = 0; i < r[p].size() ; ++i )
+	  if ( r[p][i] != NULL)
+	    delete r[p][i];
+       for ( unsigned int i(0); i < phi1[p].size(); ++i )
+	 if ( phi1[p][i] != NULL)
+	    delete phi1[p][i];
+       for ( unsigned int i(0); i < phi2[p].size(); ++i )
+	 if ( phi2[p][i] != NULL)
+	    delete phi2[p][i];
+       for ( unsigned int i(0); i < levelSet[p].size(); ++i )
+	 if ( levelSet[p][i] != NULL)
+	    delete levelSet[p][i];
+    }
+
+  }
 }

AMDiS/src/ProblemImplicit.h

@@ -23,81 +23,69 @@
 #ifndef AMDIS_PROBLEM_IMPLICIT_H
 #define AMDIS_PROBLEM_IMPLICIT_H
 
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+#ifdef HAVE_PARALLEL_MTL4
+#include "parallel/Mtl4Solver.h"
+#else
+#include "parallel/PetscProblemStat.h"
+#endif
+#else
 #include "ProblemStat.h"
-
+#endif
 namespace AMDiS {
   
   using namespace std;
+  const Flag INIT_IMPLICIT_MESH = 0X2000L;
 
-  class ProblemImplicitBase 
-  {
-  public:
-    virtual ~ProblemImplicitBase() {}
-
-    virtual bool createImplicitMesh() = 0;
-
-    virtual DOFVector<double>* getSignedDistance(unsigned int im = 0, 
-						 unsigned int m = 0) = 0;
-
-    virtual DOFVector<double>* getPhi1(unsigned int im = 0, 
-				       unsigned int m = 0) = 0;
-
-    virtual DOFVector<double>* getPhi2(unsigned int im = 0, 
-				       unsigned int m = 0) = 0;
-
-    virtual DOFVector<double>* getLevelset(unsigned int im = 0, 
-					   unsigned int m = 0) = 0;
-
-  protected:
-    void readDofVec(istream& , DOFVector<double>* , Mesh*);
-
-    void readR(istream&, double, Mesh*, int implMesh = 0, int comp = 0);
-
-    void readPhi1(istream&, double, Mesh*, int implMesh = 0, int comp = 0);
-
-    void readPhi2(istream&, double, Mesh*, int implMesh = 0, int comp = 0);
-
-    string getDofFilename(string path, int implMesh = 0);
-
-    double getEpsilon(string path, int implMesh = 0);
-
-    int getType(string path, int implMesh = 0);
-
-    virtual bool isImplicitMesh(int comp = 0) = 0;
-  };
- 
-
-  class ProblemImplicit : public ProblemStatSeq, 
-			  public ProblemImplicitBase
+  class ProblemImplicit : public ProblemStat
   {
   public:
     ProblemImplicit(string name, 
 		    ProblemIterationInterface* problem = NULL)
-      : ProblemStatSeq(name, problem),
+      : ProblemStat(name, problem),
 	r(0),
 	phi1(0),
 	phi2(0),
-	levelSet(0),
-        implMesh(0)
+	levelSet(0)
     {}
 
-    virtual void createMesh();
+    virtual ~ProblemImplicit();
 
-    virtual bool createImplicitMesh();
+    virtual void createMesh();
 
     virtual void initialize(Flag initFlag, 
 			    ProblemStatSeq* adoptProblem = NULL, 
 		            Flag adoptFlag = INIT_NOTHING);
 
-    DOFVector<double>* getSignedDistance(unsigned int im = 0, unsigned int m = 0);
+    bool createImplicitMesh();
 
-    DOFVector<double>* getPhi1(unsigned int im = 0, unsigned int m = 0); 
+    DOFVector<double>* getSignedDistance(unsigned int im = 0, 
+						 unsigned int m = 0) ;
 
-    DOFVector<double>* getPhi2(unsigned int im = 0, unsigned int m = 0);
+    DOFVector<double>* getPhi1(unsigned int im = 0, 
+				       unsigned int m = 0) ;
 
-    DOFVector<double>* getLevelset(unsigned int im = 0, unsigned int m = 0);
+    DOFVector<double>* getPhi2(unsigned int im = 0, 
+				       unsigned int m = 0) ;
+
+    DOFVector<double>* getLevelset(unsigned int im = 0, 
+					   unsigned int m = 0) ;
 
   protected:
+    void readDofVec(const std::string& , DOFVector<double>* , Mesh*);
+
+    void readR(const std::string&, double, Mesh*, int implMesh = 0, int comp = 0);
+
+    void readPhi1(const std::string&, double, Mesh*, int implMesh = 0, int comp = 0);
+
+    void readPhi2(const std::string&, double, Mesh*, int implMesh = 0, int comp = 0);
+
+    string getDofFilename(string path, int implMesh = 0);
+
+    double getEpsilon(string path, int implMesh = 0);
+
+    int getType(string path, int implMesh = 0);
+
     bool createImplicitMesh(int p);
 
     bool createImplicitMesh(string, int, int); 
@@ -115,12 +103,8 @@ namespace AMDiS {
     /// (levelSet(x): x \in \Omega: 1, x \not \in Omega: -1, x \in \Gamma: 0)
     vector< vector<DOFVector<double>*> >  levelSet;
 
-    vector<bool> implMesh;
-    
-    bool isImplicitMesh(int comp = 0) 
-    { 
-      return implMesh[comp]; 
-    }
   };
+ 
+
 }
 #endif

AMDiS/src/nonlin/NonLinSolver.h

@@ -33,8 +33,12 @@
 #include "ProblemStat.h"
 
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+#ifdef HAVE_PARALLEL_MTL4
+#include "parallel/Mtl4Solver.h"
+#else
 #include "parallel/PetscProblemStat.h"
 #endif
+#endif
 
 namespace AMDiS {
 

AMDiS/src/nonlin/ProblemNonLin.h

@@ -30,8 +30,12 @@
 #include "MatrixVector.h"
 
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+#ifdef HAVE_PARALLEL_MTL4
+#include "parallel/Mtl4Solver.h"
+#else
 #include "parallel/PetscProblemStat.h"
 #endif
+#endif
 
 namespace AMDiS {