Blub

AMDiS/CMakeLists.txt

@@ -36,6 +36,7 @@ if(CMAKE_CXX_COMPILER MATCHES ".*icpc")
 	set(CMAKE_CXX_FLAGS "-diag-disable 654 -diag-disable 858") 
 endif()
 
+
 SET(ENABLE_PARALLEL_DOMAIN "OFF" CACHE STRING "use parallel domain decomposition. please set to one of: PMTL, PETSC, OFF" )
 option(USE_PETSC_DEV false)
 option(ENABLE_ZOLTAN false)
@@ -43,6 +44,12 @@ option(ENABLE_UMFPACK "Use of UMFPACK solver" false)
 option(ENABLE_PNG "use png reader/writer" false)
 option(ENABLE_BDDCML "Use of BDDCML library" false)
 option(ENABLE_EXTENSIONS "Use extensions" false)
+option(ENABLE_OPENMP "Use OpenMP" false)
+option(ENABLE_OUTPUT "AMDiS output priniting, disable only for debugging!" true)
+
+
+mark_as_advanced(ENABLE_OUTPUT)
+
 
 find_package(Boost 1.42 REQUIRED)
 if(Boost_FOUND)
@@ -307,6 +314,12 @@ if(ENABLE_PNG)
 endif(ENABLE_PNG)
 
 
+if(NOT ENABLE_OUTPUT)
+       message(WARNING "AMDiS cout output disabled!")
+       list(APPEND COMPILEFLAGS "-DSUPPRESS_OUTPUT")
+endif(NOT ENABLE_OUTPUT)
+
+
 if(ENABLE_BDDCML)
 	SET(BDDCML_LINK_LIST "" CACHE STRING "Further libraries to be linked with BDDCML")
 		
@@ -440,6 +453,10 @@ if(ENABLE_EXTENSIONS)
 endif(ENABLE_EXTENSIONS)
 
 
+if(ENABLE_OPENMP)
+	list(APPEND COMPILEFLAGS "-fopenmp")
+endif(ENABLE_OPENMP)
+
 SET(COMPOSITE_SOURCE_DIR ${SOURCE_DIR}/compositeFEM)
 SET(COMPOSITE_FEM_SRC ${COMPOSITE_SOURCE_DIR}/CFE_Integration.cc 
 		      ${COMPOSITE_SOURCE_DIR}/CFE_NormAndErrorFcts.cc 

AMDiS/src/AMDiS_fwd.h

@@ -25,6 +25,7 @@
 #define AMDIS_AMDIS_FWD_INCLUDE
 
 #include <boost/numeric/mtl/mtl.hpp>
+#include "OpenMP.h"
 
 namespace AMDiS {
 

AMDiS/src/Assembler.h

@@ -252,16 +252,12 @@ namespace AMDiS {
     ///
     ElementMatrix tmpMat;
   
-    /** \brief
-     * Used to check whether \ref initElement() must be called, because
-     * a new Element is visited.
-     */
+    /// Used to check whether \ref initElement() must be called, because
+    /// a new Element is visited.
     Element* lastMatEl;
 
-    /** \brief
-     * Used to check whether \ref initElement() must be called, because
-     * a new Element is visited.
-     */
+    /// Used to check whether \ref initElement() must be called, because
+    /// a new Element is visited.
     Element* lastVecEl;
 
     /// Used to check for new traverse.

AMDiS/src/BasisFunction.h

@@ -111,23 +111,6 @@ namespace AMDiS {
     virtual int* orderOfPositionIndices(const Element* el, GeoIndex position, 
 					int positionIndex) const = 0;
 
-    /** \brief
-     * Pointer to a function which connects the set of local basis functions
-     * with its global DOFs.
-     * getDOFIndices(el, admin, dof) returns a pointer to a const vector of 
-     * length \ref nBasFcts where the i-th entry is the index of the DOF 
-     * associated to the i-th basis function; arguments are the actual element 
-     * el and the DOF admin admin of the corresponding finite element space 
-     * (these indices depend on all defined DOF admins and thus on the 
-     * corresponding admin); if the last argument dof is NULL, getDOFndices 
-     * has to provide memory for storing this vector, which is overwritten on the
-     * next call of getDOFIndices; if dof is not NULL, dof is a pointer to a 
-     * vector which has to be filled;   
-     */
-    virtual const DegreeOfFreedom* getDOFIndices(const Element*,
-						 const DOFAdmin&, 
-						 DegreeOfFreedom*) const = 0;
-
     /** \brief
      * The second argument 'bound' has to be a pointer to a vector which has 
      * to be filled. Its length is \ref nBasFcts (the number of basis functions
@@ -190,8 +173,7 @@ namespace AMDiS {
      * , indices[n-1] are the local indices of the basis functions where the
      * coefficients have to be calculated, and the i-th entry in the return 
      * vector is then the coefficient of the indices[i]-th basis function; coeff 
-     * may be a pointer to a vector which has to be filled 
-     * (compare the dof argument of \ref getDOFIndices());
+     * may be a pointer to a vector which has to be filled.
      * such a function usually needs vertex coordinate information; thus, all 
      * routines using this function on the elements need the FILL COORDS flag 
      * during mesh traversal.
@@ -285,35 +267,20 @@ namespace AMDiS {
     {}
 
     /// Returns local dof indices of the element for the given fe space.
-    virtual const DegreeOfFreedom *getLocalIndices(const Element *el,
-						   const DOFAdmin *admin,
-						   DegreeOfFreedom *dofPtr) const
-    {
-      return NULL;
-    }
-
-    inline void getLocalIndices(const Element *el,
-				const DOFAdmin *admin,
-				vector<DegreeOfFreedom> &indices) const
-    {
-      FUNCNAME("BasisFunction::getLocalIndices()");
-      
-      TEST_EXIT_DBG(static_cast<int>(indices.size()) >= nBasFcts)
-	("Index vector is too small!\n");
-
-      getLocalIndices(el, admin, &(indices[0]));
-    }
+    virtual void getLocalIndices(const Element *el,
+				 const DOFAdmin *admin,
+				 vector<DegreeOfFreedom> &indices) const
+    {}
 
+    ///
     virtual void getLocalDofPtrVec(const Element *el, 
 				   const DOFAdmin *admin,
 				   vector<const DegreeOfFreedom*>& vec) const
     {}
 
 
-    /** \brief
-     * Evaluates elements value at barycentric coordinates lambda with local 
-     * coefficient vector uh.
-     */
+    /// Evaluates elements value at barycentric coordinates lambda with local 
+    /// coefficient vector uh.
     template<typename T>
     T evalUh(const DimVec<double>& lambda, const mtl::dense_vector<T>& uh) const;
 

AMDiS/src/DOFIndexed.h

@@ -60,28 +60,20 @@ namespace AMDiS {
     virtual void compressDOFIndexed(int first, int last, 
 				    std::vector<DegreeOfFreedom> &newDOF) = 0;
 
-    /** \brief
-     * Performs needed action when a DOF index is freed. Can be overriden in
-     * sub classes. The default behavior is to do nothing.
-     */
+    /// Performs needed action when a DOF index is freed. Can be overriden in
+    /// sub classes. The default behavior is to do nothing.
     virtual void freeDOFContent(int) {}
 
-    /** \brief
-     * Interpolation after refinement. Can be overriden in subclasses.
-     * The default behavior is to do nothing.
-     */
+    /// Interpolation after refinement. Can be overriden in subclasses.
+    /// The default behavior is to do nothing.
     virtual void refineInterpol(RCNeighbourList&, int) {}
 
-    /** \brief
-     * Restriction after coarsening. Can be overriden in subclasses.
-     * The default behavior is to do nothing.
-     */
+    /// Restriction after coarsening. Can be overriden in subclasses.
+    /// The default behavior is to do nothing.
     virtual void coarseRestrict(RCNeighbourList&, int) {}
 
-    /** \brief
-     * Returns the finite element space of this DOFIndexed object. Must be
-     * overriden in sub classes. 
-     */
+    /// Returns the finite element space of this DOFIndexed object. Must be
+    /// overriden in sub classes. 
     virtual const FiniteElemSpace* getFeSpace() const = 0;
   };
 

AMDiS/src/DOFMatrix.cc

@@ -231,16 +231,8 @@ namespace AMDiS {
 	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
 
       if (condition && condition->isDirichlet()) {	
-	if (condition->applyBoundaryCondition()) {
-
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-	  if (dofMap->isRankDof(rowIndices[i])) {
-	    dirichletDofs.insert(row);
-	  }
-#else
+	if (condition->applyBoundaryCondition())
 	  dirichletDofs.insert(row);
-#endif
-	}
       } else {
 	for (int j = 0; j < nCol; j++) {
 	  DegreeOfFreedom col = colIndices[j];
@@ -497,7 +489,10 @@ namespace AMDiS {
     inserter_type &ins = *inserter;  
     for (std::set<int>::iterator it = dirichletDofs.begin(); 
 	 it != dirichletDofs.end(); ++it)
-      ins[*it][*it] = 1.0;    
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+      if (dofMap->isRankDof(*it))
+#endif
+	ins[*it][*it] = 1.0;    
   }
 
 

AMDiS/src/DOFMatrix.h

@@ -313,6 +313,11 @@ namespace AMDiS {
       return rowFeSpace->getAdmin()->getUsedSize();
     }
 
+    std::set<DegreeOfFreedom>& getDirichletRows()
+    {
+      return dirichletDofs;
+    }
+
     /// Returns \ref name
     inline string getName() const 
     { 

AMDiS/src/DOFVector.cc

@@ -79,7 +79,7 @@ namespace AMDiS {
     int dim = mesh->getDim();
     int nBasFcts = basFcts->getNumber();
   
-    DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
+    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
     DimVec<double> lambda(dim, NO_INIT);
   
     ElInfo *elInfo = mesh->createNewElInfo();
@@ -104,7 +104,6 @@ namespace AMDiS {
     } else
       throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));
 
-    delete [] localIndices;
     if (oldElInfo == NULL)
       delete elInfo;
 
@@ -124,11 +123,9 @@ namespace AMDiS {
     int dim = mesh->getDim();
     int nBasFcts = basFcts->getNumber();
   
-    DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
-    DimVec<double> lambda(dim, NO_INIT);
-  
+    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
+    DimVec<double> lambda(dim, NO_INIT);  
     ElInfo *elInfo = mesh->createNewElInfo();
-
     WorldVector<double> value(DEFAULT_VALUE, 0.0);
     bool inside = false;
   
@@ -150,7 +147,6 @@ namespace AMDiS {
     } else
       throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));
 
-    delete [] localIndices;
     if (oldElInfo == NULL)
       delete elInfo;
 

AMDiS/src/DOFVector.h

@@ -237,6 +237,18 @@ namespace AMDiS {
       boundaryManager = bm;
     }
 
+    inline void setDirichletDofValue(DegreeOfFreedom dof,
+				     T value)
+    {
+      dirichletDofValues[dof] = value;
+    }
+
+    map<DegreeOfFreedom, T>& getDirichletValues()
+    {
+      return dirichletDofValues;
+    }
+
+
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
     void setDofMap(FeSpaceDofMap& m)
     {
@@ -249,17 +261,6 @@ namespace AMDiS {
 
       return dofMap->isRankDof(dof);
     }
-
-    inline void setDirichletDofValue(DegreeOfFreedom dof,
-				     T value)
-    {
-      dirichletDofValues[dof] = value;
-    }
-
-    map<DegreeOfFreedom, T>& getDirichletValues()
-    {
-      return dirichletDofValues;
-    }
 #endif
 
   protected:
@@ -290,10 +291,10 @@ namespace AMDiS {
     /// Dimension of the mesh this DOFVectorBase belongs to
     int dim;
 
+    map<DegreeOfFreedom, T> dirichletDofValues;
+
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
     FeSpaceDofMap *dofMap;
-
-    map<DegreeOfFreedom, T> dirichletDofValues;
 #endif
   };
 

AMDiS/src/DOFVector.hh

@@ -1374,9 +1374,9 @@ namespace AMDiS {
     TEST_EXIT_DBG(feSpace->getMesh() == el->getMesh())
       ("Element is defined on a different mesh than the DOF vector!\n");
 
+    std::vector<int> localIndices(nBasFcts);
     const DOFAdmin* admin = feSpace->getAdmin();
-    const DegreeOfFreedom *localIndices = 
-      feSpace->getBasisFcts()->getLocalIndices(el, admin, NULL);   
+    feSpace->getBasisFcts()->getLocalIndices(el, admin, localIndices);
 
     for (int i = 0; i < nBasFcts; i++)
       d[i] = (*this)[localIndices[i]];

AMDiS/src/DirichletBC.cc

@@ -65,7 +65,6 @@ namespace AMDiS {
 
     for (int i = 0; i < nBasFcts; i++) {
 
-#if 1
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
       if (vector->isRankDof(dofIndices[i]))
 #endif
@@ -80,29 +79,10 @@ namespace AMDiS {
 	    else
 	      ERROR_EXIT("There is something wrong!\n");
 	  }
-
+	  vector->setDirichletDofValue(dofIndices[i], value);
 	  (*vector)[dofIndices[i]] = value;
 	}
-#else
-      if (localBound[i] == boundaryType) {
-	double value = 0.0;
-	if (f) {
-          elInfo->coordToWorld(*(basFcts->getCoords(i)), worldCoords);
-	  value = (*f)(worldCoords);
-	}
-	if (dofVec)
-	  value = (*dofVec)[dofIndices[i]];
-
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-	vector->setDirichletDofValue(dofIndices[i], value);
-#else
-	(*vector)[dofIndices[i]] = value;
-#endif
-      }
-#endif
-
     }
-
   }
 
 }

AMDiS/src/FirstOrderAssembler.cc

@@ -22,11 +22,15 @@
 
 namespace AMDiS {
 
-  std::vector<SubAssembler*> FirstOrderAssembler::optimizedSubAssemblersGrdPhi;
-  std::vector<SubAssembler*> FirstOrderAssembler::optimizedSubAssemblersGrdPsi;
+  ThreadPrivate<vector<SubAssembler*> > 
+  FirstOrderAssembler::optimizedSubAssemblersGrdPhi;
+  ThreadPrivate<vector<SubAssembler*> >
+  FirstOrderAssembler::optimizedSubAssemblersGrdPsi;
 
-  std::vector<SubAssembler*> FirstOrderAssembler::standardSubAssemblersGrdPhi;
-  std::vector<SubAssembler*> FirstOrderAssembler::standardSubAssemblersGrdPsi;
+  ThreadPrivate<vector<SubAssembler*> >
+  FirstOrderAssembler::standardSubAssemblersGrdPhi;
+  ThreadPrivate<vector<SubAssembler*> >
+  FirstOrderAssembler::standardSubAssemblersGrdPsi;
 
 
   FirstOrderAssembler::FirstOrderAssembler(Operator *op,
@@ -51,16 +55,16 @@ namespace AMDiS {
 							    FirstOrderType type,
 							    bool optimized)
   {
-    std::vector<SubAssembler*> *subAssemblers =
+    vector<SubAssembler*> &subAssemblers =
       optimized ?
       (type == GRD_PSI ? 
-       &optimizedSubAssemblersGrdPsi : 
-       &optimizedSubAssemblersGrdPhi) :
+       optimizedSubAssemblersGrdPsi.get() : 
+       optimizedSubAssemblersGrdPhi.get()) :
       (type == GRD_PSI ? 
-       &standardSubAssemblersGrdPsi :
-       &standardSubAssemblersGrdPhi);
+       standardSubAssemblersGrdPsi.get() :
+       standardSubAssemblersGrdPhi.get());
     
-    std::vector<OperatorTerm*> opTerms = 
+    vector<OperatorTerm*> opTerms = 
       (type == GRD_PSI) ? op->firstOrderGrdPsi : op->firstOrderGrdPhi;
 
     // check if a assembler is needed at all
@@ -72,13 +76,13 @@ namespace AMDiS {
     FirstOrderAssembler *newAssembler;
 
     // check if a new assembler is needed
-    for (unsigned int i = 0; i < subAssemblers->size(); i++) {
-      std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
+    for (unsigned int i = 0; i < subAssemblers.size(); i++) {
+      vector<OperatorTerm*> assTerms = *(subAssemblers[i]->getTerms());
     
       sort(assTerms.begin(), assTerms.end());   
 
-      if (opTerms == assTerms && (*subAssemblers)[i]->getQuadrature() == quad)
-	return dynamic_cast<FirstOrderAssembler*>((*subAssemblers)[i]);
+      if (opTerms == assTerms && subAssemblers[i]->getQuadrature() == quad)
+	return dynamic_cast<FirstOrderAssembler*>(subAssemblers[i]);
     }
 
     // check if all terms are pw_const
@@ -110,7 +114,7 @@ namespace AMDiS {
        }
     }
 
-    subAssemblers->push_back(newAssembler);
+    subAssemblers.push_back(newAssembler);
     return newAssembler;
   }
 
@@ -130,7 +134,7 @@ namespace AMDiS {
     mtl::dense_vector<double> grdPsi(dim + 1, 0.0);
     int nPoints = quadrature->getNumPoints();
     Lb.resize(nPoints);
-    std::vector<double> phival(nCol);
+    vector<double> phival(nCol);
 
     for (int iq = 0; iq < nPoints; iq++) {
       Lb[iq].change_dim(dim + 1);

AMDiS/src/FirstOrderAssembler.h

@@ -23,11 +23,14 @@
 #ifndef AMDIS_FIRST_ORDER_ASSEMBLER_H
 #define AMDIS_FIRST_ORDER_ASSEMBLER_H
 
+#include <vector>
 #include "AMDiS_fwd.h"
 #include "SubAssembler.h"
 
 namespace AMDiS {
 
+  using namespace std;
+
   /**
    * \ingroup Assembler
    * 
@@ -62,19 +65,19 @@ namespace AMDiS {
 
   protected:
     /// Vector of DimMats for calculation in function calculateElementMatrix().
-    std::vector<mtl::dense_vector<double> > Lb;
+    vector<mtl::dense_vector<double> > Lb;
 
     /// List of all yet created optimized zero order assemblers for grdPsi.
-    static std::vector<SubAssembler*> optimizedSubAssemblersGrdPsi;
+    static ThreadPrivate<vector<SubAssembler*> > optimizedSubAssemblersGrdPsi;
 
     /// List of all yet created standard zero order assemblers for grdPsi.
-    static std::vector<SubAssembler*> standardSubAssemblersGrdPsi;
+    static ThreadPrivate<vector<SubAssembler*> > standardSubAssemblersGrdPsi;
 
     /// List of all yet created optimized zero order assemblers for grdPhi.
-    static std::vector<SubAssembler*> optimizedSubAssemblersGrdPhi;
+    static ThreadPrivate<vector<SubAssembler*> > optimizedSubAssemblersGrdPhi;
 
     /// List of all yet created standard zero order assemblers for grdPhi.
-    static std::vector<SubAssembler*> standardSubAssemblersGrdPhi;
+    static ThreadPrivate<vector<SubAssembler*> > standardSubAssemblersGrdPhi;
   };
 
 

AMDiS/src/Global.cc

@@ -29,7 +29,7 @@ namespace AMDiS {
   bool Msg::outputMainRank = false;
 #endif
 
-  const char *Msg::oldFuncName = NULL;
+  ThreadPrivate<const char *> Msg::oldFuncName(NULL);
   std::ostream* Msg::out = NULL;
   std::ostream* Msg::error = NULL;
   int Global::dimOfWorld = 0;
@@ -125,14 +125,14 @@ namespace AMDiS {
     if (!out) 
       out = &std::cout;
 
-    if (funcName &&  oldFuncName != funcName) {
+    if (funcName &&  oldFuncName.get() != funcName) {
       PRINT_LINE((*out), funcName << ":" << std::endl);
     } else if (!funcName) {
       PRINT_LINE((*out), "*unknown function*" << std::endl);
     }
     PRINT_LINE((*out), "               ");
 
-    oldFuncName = funcName;
+    oldFuncName.set(funcName);
   }
 
 
@@ -145,16 +145,16 @@ namespace AMDiS {
 
     std::stringstream oss;
 
-    if (funcName && oldFuncName != funcName) {
+    if (funcName && oldFuncName.get() != funcName) {
       oss << funcName << ": ";
     } else if (!funcName) {
       if (line-old_line > 5) 
 	oss << "*unknown function*";
     }
 
-    if (oldFuncName != funcName) {
+    if (oldFuncName.get() != funcName) {
       oss << "ERROR in " << file << ", line " << line << std::endl;;
-      oldFuncName = funcName;
+      oldFuncName.set(funcName);
     } else if (line - old_line > 5)
       oss << "ERROR in " << file << ", line " << line << "\n" << std::endl;
 
@@ -207,15 +207,15 @@ namespace AMDiS {
 
     std::stringstream oss;
 
-    if (funcName  &&  oldFuncName != funcName) {
+    if (funcName  &&  oldFuncName.get() != funcName) {
       oss << funcName << ": ";
     } else if (!funcName) {
       oss << "*unknown function*";
     }
 
-    if (oldFuncName != funcName) {
+    if (oldFuncName.get() != funcName) {
       oss << "WARNING in " << file << ", line " << line << std::endl;
-      oldFuncName = funcName;
+      oldFuncName.set(funcName);
     } else if (line - old_line > 5) {
       oss << "WARNING in " << file << ", line " << line << std::endl;
     }
@@ -244,6 +244,7 @@ namespace AMDiS {
 
   void Msg::print(const char *format, ...)
   {
+#ifndef SUPPRESS_OUTPUT
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
     if (outputMainRank && MPI::COMM_WORLD.Get_rank() != 0)
       return;
@@ -259,6 +260,7 @@ namespace AMDiS {
     vsprintf(buff, format, arg);
     (*out) << buff;
     va_end(arg);
+#endif
   }
 
 

AMDiS/src/Global.h

@@ -60,6 +60,7 @@
 #include <boost/algorithm/string/trim.hpp>
 #include "boost/tuple/tuple.hpp"
 #include "AMDiS_fwd.h"
+#include "OpenMP.h"
 
 namespace AMDiS {
 
@@ -100,10 +101,15 @@ namespace AMDiS {
   /// Returns the dimension of GeoIndex ind for dimension dim
 #define DIM_OF_INDEX(ind, dim) ((static_cast<int>(ind) == 0) ? dim : static_cast<int>(ind) - 1)
 
+
+#if SUPPRESS_OUTPUT
+#define PRINT_LINE(stream, line)
+#else
 #if HAVE_PARALLEL_DOMAIN_AMDIS
 #define PRINT_LINE(stream, line) stream << "[" << MPI::COMM_WORLD.Get_rank() << "] " << line
 #else
 #define PRINT_LINE(stream, line) stream << line
+#endif
 #endif
 
   /// Calculates factorial of i
@@ -298,7 +304,7 @@ namespace AMDiS {
 
     /// Remember funcName to avoid multiple output of funcName within the same
     /// function call
-    static const char *oldFuncName;
+    static ThreadPrivate<const char*> oldFuncName;
 
     /// Global info level
     static int msgInfo;

AMDiS/src/Lagrange.h

@@ -73,11 +73,6 @@ namespace AMDiS {
     /// Returns the barycentric coordinates of the i-th basis function.
     DimVec<double> *getCoords(int i) const;
 
-    /// Implements BasisFunction::getDOFIndices
-    const DegreeOfFreedom* getDOFIndices(const Element*,
-					 const DOFAdmin&, 
-					 DegreeOfFreedom*) const;