diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt
index a6f4ce880dc67d6fbfe8196ad6e8459cc6ade700..af9f197a14da69471080dfa7d09b9c36c861729f 100644
--- a/AMDiS/CMakeLists.txt
+++ b/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
diff --git a/AMDiS/src/AMDiS_fwd.h b/AMDiS/src/AMDiS_fwd.h
index 9078dfb18a14a4f7eb43393cf1e526b715c64a4f..7afcc8a0efa2c3cd02154a092a63c24909b657a6 100644
--- a/AMDiS/src/AMDiS_fwd.h
+++ b/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 {
diff --git a/AMDiS/src/Assembler.h b/AMDiS/src/Assembler.h
index fabbb2e00d86df08741fbf61502f0fddf8f2e6e6..367cdfa98033f226adb95b9e9c99735bca176c2e 100644
--- a/AMDiS/src/Assembler.h
+++ b/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.
diff --git a/AMDiS/src/BasisFunction.h b/AMDiS/src/BasisFunction.h
index 55f681632e7980e80d21342ee6ed3860a3f0e35c..553fd8960d64fd3d4dd62aa8b28e09f3a060f6ac 100644
--- a/AMDiS/src/BasisFunction.h
+++ b/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;
diff --git a/AMDiS/src/DOFIndexed.h b/AMDiS/src/DOFIndexed.h
index 494cc41fc74b9f2e11fdc5d17c24e05f8abd0567..bfc0fc885f880367d8a1a2d3269a1244133c8250 100644
--- a/AMDiS/src/DOFIndexed.h
+++ b/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;
};
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 9850d8f2cd8a8824cb9b0c03d4e7b19b27861b8b..20f90cf4cda1f8d72fedf7b6ce07d8f4ee45a0aa 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/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;
}
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index 2f3b082095271be758a45476525800205b9e9eb6..d000661864cdc5328aa2bf5c811b2e807de0d331 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/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
{
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 57212a3dcda8ba01296b7360aa52d952f98fa97a..9023f6ec429f94ce94028e482637f16804371128 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/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;
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index a8335726f0fcf59105a90ad9e83be9803016c9e9..e7fc40b349cb5116bb33f31748aa4723b9fcdbe9 100644
--- a/AMDiS/src/DOFVector.h
+++ b/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
};
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 39ecb551407ed2488f7969748c491d35ebc57ef2..e0cc38656d999a782bd9ca37c2ca5dc28774e228 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/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]];
diff --git a/AMDiS/src/DirichletBC.cc b/AMDiS/src/DirichletBC.cc
index face0ddd3ece92c35d86f11a294aeeb2b101cd39..38952e597980ab43e29916470a1b187f9eb2c5f8 100644
--- a/AMDiS/src/DirichletBC.cc
+++ b/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
-
}
-
}
}
diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc
index 4359539fe4a1275a59ce0760255a682f8e05b7f4..1c5760543ddcc25b96a4ee1050a92bc5b6591f39 100644
--- a/AMDiS/src/FirstOrderAssembler.cc
+++ b/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);
diff --git a/AMDiS/src/FirstOrderAssembler.h b/AMDiS/src/FirstOrderAssembler.h
index 298522537ae005ca712722e9fef513dce202fad8..61f2aaae476647c04dd7915bbcd973f52ce7a996 100644
--- a/AMDiS/src/FirstOrderAssembler.h
+++ b/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;
};
diff --git a/AMDiS/src/Global.cc b/AMDiS/src/Global.cc
index 1b8499423fcbe3860e8efef0fbbf42e40a05ede1..55ad68b6e69b1cd0ed6ce9aea3f0c023476c055a 100644
--- a/AMDiS/src/Global.cc
+++ b/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
}
diff --git a/AMDiS/src/Global.h b/AMDiS/src/Global.h
index abfb3fbddc6d7593417582d5bb7b46aceadbd1d1..d11c33c9fa8a215f265a40cee768957b4b769a72 100644
--- a/AMDiS/src/Global.h
+++ b/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;
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index 21d597780e3f4f1f6049849ae69e1f67006beaf5..064cae81d76bec559e81b0851a3a5dde4c67fedb 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -391,7 +391,9 @@ namespace AMDiS {
}
}
- Lagrange::Phi::~Phi() {
+
+ Lagrange::Phi::~Phi()
+ {
delete [] vertices;
}
@@ -701,14 +703,6 @@ namespace AMDiS {
}
- const DegreeOfFreedom *Lagrange::getDOFIndices(const Element *el,
- const DOFAdmin& admin,
- DegreeOfFreedom *idof) const
- {
- WARNING("depricated: use getLocalIndices() instead\n");
- return getLocalIndices(el, &admin, idof);
- }
-
int* Lagrange::orderOfPositionIndices(const Element* el,
GeoIndex position,
int positionIndex) const
@@ -887,32 +881,16 @@ namespace AMDiS {
}
- const DegreeOfFreedom *Lagrange::getLocalIndices(const Element* el,
- const DOFAdmin *admin,
- DegreeOfFreedom *indices) const
+ void Lagrange::getLocalIndices(const Element *el,
+ const DOFAdmin *admin,
+ std::vector<DegreeOfFreedom>& dofs) const
{
- static DegreeOfFreedom *localVec = NULL;
- static int localVecSize = 0;
+ if (static_cast<int>(dofs.size()) < nBasFcts)
+ dofs.resize(nBasFcts);
const DegreeOfFreedom **dof = el->getDof();
-
- int nrDOFs, n0, node0, num = 0;
GeoIndex posIndex;
- DegreeOfFreedom* result;
-
- if (indices) {
- result = indices;
- } else {
- if (localVec && nBasFcts > localVecSize) {
- delete [] localVec;
- localVec = new DegreeOfFreedom[nBasFcts];
- }
- if (!localVec)
- localVec = new DegreeOfFreedom[nBasFcts];
-
- localVecSize = nBasFcts;
- result = localVec;
- }
+ int nrDOFs, n0, node0, num = 0;
for (int pos = 0, j = 0; pos <= dim; pos++) {
posIndex = INDEX_OF_DIM(pos, dim);
@@ -927,20 +905,19 @@ namespace AMDiS {
const int *indi = orderOfPositionIndices(el, posIndex, i);
for (int k = 0; k < nrDOFs; k++)
- result[j++] = dof[node0][n0 + indi[k]];
+ dofs[j++] = dof[node0][n0 + indi[k]];
}
}
}
-
- return result;
}
+
void Lagrange::getLocalDofPtrVec(const Element *el,
const DOFAdmin *admin,
- std::vector<const DegreeOfFreedom*>& vec) const
+ std::vector<const DegreeOfFreedom*>& dofs) const
{
- if (static_cast<int>(vec.size()) < nBasFcts)
- vec.resize(nBasFcts);
+ if (static_cast<int>(dofs.size()) < nBasFcts)
+ dofs.resize(nBasFcts);
const DegreeOfFreedom **dof = el->getDof();
GeoIndex posIndex;
@@ -959,12 +936,13 @@ namespace AMDiS {
const int *indi = orderOfPositionIndices(el, posIndex, i);
for (int k = 0; k < nrDOFs; k++)
- vec[j++] = &(dof[node0][n0 + indi[k]]);
+ dofs[j++] = &(dof[node0][n0 + indi[k]]);
}
}
}
}
+
void Lagrange::l2ScpFctBas(Quadrature *q,
AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
DOFVector<WorldVector<double> >* fh)
@@ -972,6 +950,7 @@ namespace AMDiS {
ERROR_EXIT("not yet\n");
}
+
void Lagrange::l2ScpFctBas(Quadrature *quad,
AbstractFunction<double, WorldVector<double> >* f,
DOFVector<double>* fh)
@@ -991,16 +970,17 @@ namespace AMDiS {
const FastQuadrature *quad_fast =
FastQuadrature::provideFastQuadrature(this, *quad, INIT_PHI);
- double *wdetf_qp = new double[quad->getNumPoints()];
+ vector<double> wdetf_qp(quad->getNumPoints());
int nPoints = quad->getNumPoints();
DOFAdmin *admin = fh->getFeSpace()->getAdmin();
WorldVector<double> x;
+ vector<DegreeOfFreedom> dof;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(fh->getFeSpace()->getMesh(), -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
while (elInfo) {
- const DegreeOfFreedom *dof = getLocalIndices(elInfo->getElement(), admin, NULL);
+ getLocalIndices(elInfo->getElement(), admin, dof);
double det = elInfo->getDet();
for (int iq = 0; iq < nPoints; iq++) {
@@ -1018,8 +998,6 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
-
- delete [] wdetf_qp;
}
@@ -1083,7 +1061,8 @@ namespace AMDiS {
Element *el = list->getElement(0);
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
- const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
+ vector<DegreeOfFreedom> pdof;
+ basFct->getLocalIndices(el, admin, pdof);
int node = drv->getFeSpace()->getMesh()->getNode(VERTEX);
int n0 = admin->getNumberOfPreDofs(VERTEX);
@@ -1127,7 +1106,8 @@ namespace AMDiS {
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
Element *el = list->getElement(0);
- const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
+ vector<DegreeOfFreedom> pdof;
+ basFct->getLocalIndices(el, admin, pdof);
int node = drv->getFeSpace()->getMesh()->getNode(VERTEX);
int n0 = admin->getNumberOfPreDofs(VERTEX);
@@ -1173,7 +1153,7 @@ namespace AMDiS {
/* adjust the value at the midpoint of the common edge of neigh's children */
/****************************************************************************/
el = list->getElement(1);
- pdof = basFct->getLocalIndices(el, admin, NULL);
+ basFct->getLocalIndices(el, admin, pdof);
cdof = el->getChild(0)->getDof(node + 1, n0);
(*drv)[cdof] =
@@ -1193,8 +1173,7 @@ namespace AMDiS {
Element *el = list->getElement(0);
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
-
- DegreeOfFreedom pdof[10];
+ vector<DegreeOfFreedom> pdof(10);
basFct->getLocalIndices(el, admin, pdof);
int node0 = drv->getFeSpace()->getMesh()->getNode(EDGE);
@@ -1204,7 +1183,8 @@ namespace AMDiS {
/* values on child[0] */
/****************************************************************************/
- const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
+ vector<DegreeOfFreedom> cdof(10);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
(*drv)[cdof[3]] = ((*drv)[pdof[4]]);
(*drv)[cdof[6]] =
@@ -1275,14 +1255,13 @@ namespace AMDiS {
Element *el = list->getElement(0);
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
-
- DegreeOfFreedom *pdof = new DegreeOfFreedom[basFct->getNumber()];
+ vector<DegreeOfFreedom> pdof, cdof;
basFct->getLocalIndices(el, admin, pdof);
-
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
+
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
(*drv)[cdof[2]] =
(-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]])
@@ -1306,7 +1285,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[cdof[5]] = (*drv)[pdof[8]];
(*drv)[cdof[6]] =
@@ -1317,10 +1296,9 @@ namespace AMDiS {
- 0.125 * (*drv)[pdof[6]] + 0.1875 * (-(*drv)[pdof[7]] + (*drv)[pdof[8]])
+ 0.75 * (*drv)[pdof[9]]);
- if (n <= 1) {
- delete [] pdof;
+ if (n <= 1)
return;
- }
+
/****************************************************************************/
/* adjust the value on the neihgbour */
/****************************************************************************/
@@ -1331,7 +1309,7 @@ namespace AMDiS {
/****************************************************************************/
/* values on neigh's child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
(*drv)[cdof[5]] = (*drv)[pdof[9]];
(*drv)[cdof[6]] =
@@ -1355,8 +1333,6 @@ namespace AMDiS {
(0.0625 * ((*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.375 * (*drv)[pdof[3]]
- 0.125 * (*drv)[pdof[6]] + 0.1875 * (-(*drv)[pdof[7]] + (*drv)[pdof[8]])
+ 0.75 *(*drv)[pdof[9]]);
-
- delete [] pdof;
}
void Lagrange::refineInter3_2d(DOFIndexed<double> *drv,
@@ -1372,13 +1348,13 @@ namespace AMDiS {
Element *el = list->getElement(0);
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
- DegreeOfFreedom pdof[10];
+ vector<DegreeOfFreedom> pdof, cdof;
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
(*drv)[cdof[2]] =
(-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]])
@@ -1402,7 +1378,7 @@ namespace AMDiS {
/* values on child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[cdof[5]] = (*drv)[pdof[8]];
(*drv)[cdof[6]] =
@@ -1422,11 +1398,11 @@ namespace AMDiS {
el = list->getElement(1);
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on neigh's child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
(*drv)[cdof[5]] = (*drv)[pdof[9]];
(*drv)[cdof[6]] =
@@ -1461,25 +1437,19 @@ namespace AMDiS {
if (n < 1)
return;
- Element *el;
- const DegreeOfFreedom *cd;
- DegreeOfFreedom pd[20], cdi;
- int i, typ, lr_set, node0, n0;
- const DOFAdmin *admin;
-
- el = list->getElement(0);
- typ = list->getType(0);
-
- admin = drv->getFeSpace()->getAdmin();
-
+ DegreeOfFreedom cdi;
+ int i, lr_set, node0, n0;
+ Element *el = list->getElement(0);
+ int typ = list->getType(0);
+ const DOFAdmin * admin = drv->getFeSpace()->getAdmin();
+ vector<DegreeOfFreedom> pd(20), cd(20);
basFct->getLocalIndices(el, admin, pd);
+ basFct->getLocalIndices(el->getChild(0), admin, cd);
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[cd[3]] =
(0.0625*(-(*drv)[pd[0]] - (*drv)[pd[1]])
+ 0.5625*((*drv)[pd[4]] + (*drv)[pd[5]]));
@@ -1515,7 +1485,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cd);
if (typ == 0) {
/****************************************************************************/
@@ -1578,7 +1548,7 @@ namespace AMDiS {
/* values on child[0] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
+ basFct->getLocalIndices(el->getChild(0), admin, cd);
switch (lr_set) {
case 1:
@@ -1626,7 +1596,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cd);
if (typ == 0) {
switch (lr_set) {
@@ -1679,19 +1649,17 @@ namespace AMDiS {
if (n < 1)
return;
- DegreeOfFreedom *pdof = new DegreeOfFreedom[basFct->getNumber()];
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
Element *el = list->getElement(0);
+ vector<DegreeOfFreedom> pdof, cdof;
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- const DegreeOfFreedom *cdof =
- basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[cdof[2]] = (*drv)[pdof[10]];
(*drv)[cdof[3]] =
(0.2734375*(*drv)[pdof[0]] - 0.0390625*(*drv)[pdof[1]]
@@ -1732,7 +1700,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[cdof[6]] =
(0.0234375*(*drv)[pdof[0]] - 0.0390625*(*drv)[pdof[1]]
@@ -1756,24 +1724,21 @@ namespace AMDiS {
- 0.03125*(*drv)[pdof[11]] + 0.75*(*drv)[pdof[14]]);
(*drv)[cdof[14]] = (*drv)[pdof[13]];
- if (n <= 1) {
- delete [] pdof;
+ if (n <= 1)
return;
- }
+
/****************************************************************************/
/* adjust neighbour values */
/****************************************************************************/
el = list->getElement(1);
-
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on neighbour's child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[cdof[6]] =
(0.0234375*((*drv)[pdof[0]] + (*drv)[pdof[1]])
+ 0.0625*(-(*drv)[pdof[3]] - (*drv)[pdof[8]])
@@ -1804,7 +1769,7 @@ namespace AMDiS {
/* values on neighbour's child[1] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[cdof[12]] =
(0.0390625*(-(*drv)[pdof[0]] - (*drv)[pdof[1]])
@@ -1818,8 +1783,6 @@ namespace AMDiS {
+ 0.09375*(*drv)[pdof[9]] - 0.046875*(*drv)[pdof[10]]
- 0.03125*(*drv)[pdof[11]] + 0.75*(*drv)[pdof[14]]);
(*drv)[cdof[14]] = (*drv)[pdof[13]];
-
- delete [] pdof;
}
void Lagrange::refineInter4_2d(DOFIndexed<double> *drv,
@@ -1827,24 +1790,20 @@ namespace AMDiS {
int n, BasisFunction* basFct)
{
FUNCNAME("Lagrange::refineInter4_2d");
- Element *el;
- DegreeOfFreedom pdof[15];
- const DegreeOfFreedom *cdof;
- const DOFAdmin *admin;
if (n < 1)
return;
- el = list->getElement(0);
- admin = drv->getFeSpace()->getAdmin();
+ Element * el = list->getElement(0);
+ const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
+ vector<DegreeOfFreedom> pdof(15), cdof(15);
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[cdof[2]] = (*drv)[pdof[10]];
(*drv)[cdof[3]] =
(0.2734375*(*drv)[pdof[0]] - 0.0390625*(*drv)[pdof[1]]
@@ -1887,7 +1846,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[cdof[6]] =
(0.0234375*(*drv)[pdof[0]] - 0.0390625*(*drv)[pdof[1]]
@@ -1919,15 +1878,13 @@ namespace AMDiS {
/****************************************************************************/
el = list->getElement(1);
-
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on neighbour's child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[cdof[6]] =
(0.0234375*((*drv)[pdof[0]] + (*drv)[pdof[1]])
+ 0.0625*(-(*drv)[pdof[3]] - (*drv)[pdof[8]])
@@ -1958,7 +1915,7 @@ namespace AMDiS {
/* values on neighbour's child[1] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[cdof[12]] =
(0.0390625*(-(*drv)[pdof[0]] - (*drv)[pdof[1]])
@@ -1979,9 +1936,6 @@ namespace AMDiS {
int n, BasisFunction* basFct)
{
FUNCNAME("Lagrange::refineInter4_3d");
- DegreeOfFreedom pd[35];
- const DegreeOfFreedom *cd;
- int i, lr_set;
if (n < 1)
return;
@@ -1989,15 +1943,14 @@ namespace AMDiS {
Element *el = list->getElement(0);
int typ = list->getType(0);
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
-
+ vector<DegreeOfFreedom> pd(35), cd(35);
basFct->getLocalIndices(el, admin, pd);
+ basFct->getLocalIndices(el->getChild(0), admin, cd);
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[cd[3]] = ((*drv)[pd[5]]);
(*drv)[cd[10]] =
(0.2734375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]]
@@ -2090,7 +2043,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cd);
if (typ == 0) {
/****************************************************************************/
/* parent of el_type 0 */
@@ -2201,12 +2154,12 @@ namespace AMDiS {
/* adjust neighbour values */
/****************************************************************************/
- for (i = 1; i < n; i++) {
+ for (int i = 1; i < n; i++) {
el = list->getElement(i);
typ = list->getType(i);
basFct->getLocalIndices(el, admin, pd);
- lr_set = 0;
+ int lr_set = 0;
if (list->getNeighbourElement(i, 0) && list->getNeighbourNr(i, 0) < i)
lr_set = 1;
@@ -2219,7 +2172,7 @@ namespace AMDiS {
/* values on child[0] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
+ basFct->getLocalIndices(el->getChild(0), admin, cd);
switch (lr_set) {
case 1:
@@ -2387,7 +2340,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cd);
if (typ == 0) {
switch (lr_set) {
@@ -2582,7 +2535,8 @@ namespace AMDiS {
Element *el = list->getElement(0);
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
- const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
+ vector<DegreeOfFreedom> pdof;
+ basFct->getLocalIndices(el, admin, pdof);
/****************************************************************************/
/* contributions of dofs located on child[0] */
@@ -2615,7 +2569,7 @@ namespace AMDiS {
if (n > 1) {
el = list->getElement(1);
- pdof = basFct->getLocalIndices(el, admin, NULL);
+ basFct->getLocalIndices(el, admin, pdof);
/****************************************************************************/
/* first set those values not effected by previous element */
@@ -2647,28 +2601,21 @@ namespace AMDiS {
if (n < 1)
return;
- Element *el;
- const DegreeOfFreedom *cdof;
- DegreeOfFreedom pdof[10], cdofi;
- int i, lr_set;
- int node0, n0;
- const DOFAdmin *admin;
-
- el = list->getElement(0);
-
- admin = drv->getFeSpace()->getAdmin();
+ Element *el = list->getElement(0);
+ const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
+ vector<DegreeOfFreedom> pdof(10), cdof(10);
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
- node0 = drv->getFeSpace()->getMesh()->getNode(EDGE);
- n0 = admin->getNumberOfPreDofs(EDGE);
+ int node0 = drv->getFeSpace()->getMesh()->getNode(EDGE);
+ int n0 = admin->getNumberOfPreDofs(EDGE);
+ int i, lr_set;
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[pdof[0]] +=
(0.375*(*drv)[cdof[6]] + 0.125*(-(*drv)[cdof[8]] - (*drv)[cdof[9]]));
(*drv)[pdof[1]] +=
@@ -2689,8 +2636,8 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
- cdofi = el->getChild(1)->getDof(node0 + 2, n0);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
+ DegreeOfFreedom cdofi = el->getChild(1)->getDof(node0 + 2, n0);
(*drv)[pdof[0]] += (-0.125*(*drv)[cdofi]);
(*drv)[pdof[1]] += (0.375*(*drv)[cdofi]);
@@ -2717,7 +2664,7 @@ namespace AMDiS {
/* values on child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
switch (lr_set) {
case 1:
@@ -2750,19 +2697,20 @@ namespace AMDiS {
TEST_EXIT(drv->getFeSpace()->getBasisFcts())("No basis functions in fe_space!\n");
int node, n0;
- DegreeOfFreedom pdof[10], dof9;
+ DegreeOfFreedom dof9;
if (n < 1)
return;
Element *el = list->getElement(0);
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
+ vector<DegreeOfFreedom> pdof(10), cdof(10);
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
(*drv)[pdof[0]] +=
(0.0625*(-(*drv)[cdof[2]] + (*drv)[cdof[6]] - (*drv)[cdof[9]])
@@ -2791,7 +2739,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[pdof[0]] += 0.0625*(*drv)[cdof[6]] + 0.0625*(*drv)[cdof[9]];
(*drv)[pdof[1]] += 0.3125*(*drv)[cdof[6]] - 0.0625*(*drv)[cdof[9]];
@@ -2811,11 +2759,11 @@ namespace AMDiS {
el = list->getElement(1);
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on neigh's child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
(*drv)[pdof[0]] += 0.0625*((*drv)[cdof[6]] - (*drv)[cdof[9]]);
(*drv)[pdof[1]] += 0.0625*((*drv)[cdof[6]] + (*drv)[cdof[9]]);
@@ -2853,23 +2801,22 @@ namespace AMDiS {
TEST_EXIT(drv->getFeSpace())("No fe_space in dof_real_vec!\n");
TEST_EXIT(drv->getFeSpace()->getBasisFcts())("No basis functions in fe_space!\n");
- const Element *el;
- int node, n0;
- const DegreeOfFreedom *cdof;
- DegreeOfFreedom pdof[10], dof9;
- const DOFAdmin *admin;
+ if (n < 1)
+ return;
- if (n < 1) return;
- el = list->getElement(0);
- admin = drv->getFeSpace()->getAdmin();
+ int node, n0;
+ DegreeOfFreedom dof9;
+ const Element *el = list->getElement(0);
+ const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
+ vector<DegreeOfFreedom> pdof(10), cdof(10);
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
(*drv)[pdof[0]] +=
(0.0625*(-(*drv)[cdof[2]] + (*drv)[cdof[6]] - (*drv)[cdof[9]])
@@ -2894,7 +2841,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[pdof[0]] += 0.0625*(*drv)[cdof[6]] + 0.0625*(*drv)[cdof[9]];
(*drv)[pdof[1]] += 0.3125*(*drv)[cdof[6]] - 0.0625*(*drv)[cdof[9]];
@@ -2913,12 +2860,12 @@ namespace AMDiS {
el = list->getElement(1);
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on neigh's child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
+
(*drv)[pdof[0]] += 0.0625*((*drv)[cdof[6]] - (*drv)[cdof[9]]);
(*drv)[pdof[1]] += 0.0625*((*drv)[cdof[6]] + (*drv)[cdof[9]]);
(*drv)[pdof[3]] += -0.25*(*drv)[cdof[6]] - 0.12500*(*drv)[cdof[9]];
@@ -2958,18 +2905,18 @@ namespace AMDiS {
if (n < 1)
return;
- DegreeOfFreedom pd[20], cdi;
+ DegreeOfFreedom cdi;
const Element *el = list->getElement(0);
int typ = list->getType(0);
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
+ vector<DegreeOfFreedom> pd(20), cd(20);
basFct->getLocalIndices(el, admin, pd);
+ basFct->getLocalIndices(el->getChild(0), admin, cd);
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- const DegreeOfFreedom *cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[pd[0]] +=
(0.0625*(-(*drv)[cd[3]] + (*drv)[cd[12]] + (*drv)[cd[14]]
+ (*drv)[cd[16]] - (*drv)[cd[17]] - (*drv)[cd[18]])
@@ -3010,7 +2957,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cd);
if (typ == 0) {
/****************************************************************************/
@@ -3077,7 +3024,7 @@ namespace AMDiS {
/* values on child[0] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
+ basFct->getLocalIndices(el->getChild(0), admin, cd);
cdi = cd[16];
switch (lr_set) {
@@ -3145,7 +3092,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cd);
if (typ == 0) {
switch (lr_set) {
@@ -3212,18 +3159,16 @@ namespace AMDiS {
if (n < 1)
return;
- DegreeOfFreedom pdof[15];
const Element *el = list->getElement(0);
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
-
+ vector<DegreeOfFreedom> pdof(15), cdof(15);
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[pdof[0]] +=
(0.2734375*(*drv)[cdof[3]]
+ 0.0390625*(-(*drv)[cdof[5]] - (*drv)[cdof[8]] - (*drv)[cdof[13]])
@@ -3266,7 +3211,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[pdof[0]] +=
(0.0234375*(*drv)[cdof[6]]
@@ -3299,15 +3244,13 @@ namespace AMDiS {
/****************************************************************************/
el = list->getElement(1);
-
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on neighbour's child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[pdof[0]] +=
(0.0234375*((*drv)[cdof[6]] + (*drv)[cdof[12]])
+ 0.0390625*(-(*drv)[cdof[8]] - (*drv)[cdof[13]]));
@@ -3346,7 +3289,7 @@ namespace AMDiS {
/* values on neighbour's child[1] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[pdof[0]] += 0.0390625*(-(*drv)[cdof[12]] - (*drv)[cdof[13]]);
(*drv)[pdof[1]] += -0.0390625*(*drv)[cdof[12]] + 0.0234375*(*drv)[cdof[13]];
@@ -3371,24 +3314,19 @@ namespace AMDiS {
TEST_EXIT(drv->getFeSpace())("No fe_space in dof_real_vec!\n");
TEST_EXIT(drv->getFeSpace()->getBasisFcts())("No basis functions in fe_space!\n");
- const Element *el;
- DegreeOfFreedom pdof[15];
- const DegreeOfFreedom *cdof;
- const DOFAdmin *admin;
-
- if (n < 1) return;
- el = list->getElement(0);
-
- admin = drv->getFeSpace()->getAdmin();
+ if (n < 1)
+ return;
+ const Element *el = list->getElement(0);
+ const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
+ vector<DegreeOfFreedom> pdof(15), cdof(15);
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[pdof[0]] +=
(0.2734375*(*drv)[cdof[3]]
+ 0.0390625*(-(*drv)[cdof[5]] - (*drv)[cdof[8]] - (*drv)[cdof[13]])
@@ -3431,7 +3369,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[pdof[0]] +=
(0.0234375*(*drv)[cdof[6]]
@@ -3464,15 +3402,13 @@ namespace AMDiS {
/****************************************************************************/
el = list->getElement(1);
-
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
/****************************************************************************/
/* values on neighbour's child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[pdof[0]] +=
(0.0234375*((*drv)[cdof[6]] + (*drv)[cdof[12]])
+ 0.0390625*(-(*drv)[cdof[8]] - (*drv)[cdof[13]]));
@@ -3511,7 +3447,7 @@ namespace AMDiS {
/* values on neighbour's child[1] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[pdof[0]] += 0.0390625*(-(*drv)[cdof[12]] - (*drv)[cdof[13]]);
(*drv)[pdof[1]] += -0.0390625*(*drv)[cdof[12]] + 0.0234375*(*drv)[cdof[13]];
@@ -3539,19 +3475,17 @@ namespace AMDiS {
if (n < 1)
return;
- DegreeOfFreedom pd[35];
+ vector<DegreeOfFreedom> pd(35), cd(35);
const Element *el = list->getElement(0);
int typ = list->getType(0);
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
-
basFct->getLocalIndices(el, admin, pd);
+ basFct->getLocalIndices(el->getChild(0), admin, cd);
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- const DegreeOfFreedom *cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[pd[0]] +=
(0.2734375*(*drv)[cd[10]]
+ 0.0390625*(-(*drv)[cd[12]] - (*drv)[cd[16]] - (*drv)[cd[19]]
@@ -3661,7 +3595,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cd);
if (typ == 0) {
/****************************************************************************/
@@ -3785,6 +3719,7 @@ namespace AMDiS {
el = list->getElement(i);
typ = list->getType(i);
basFct->getLocalIndices(el, admin, pd);
+ basFct->getLocalIndices(el->getChild(0), admin, cd);
int lr_set = 0;
if (list->getNeighbourElement(i,0) && list->getNeighbourNr(i,0) < i)
@@ -3799,8 +3734,6 @@ namespace AMDiS {
/* values on child[0] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
switch (lr_set) {
case 1:
(*drv)[pd[0]] +=
@@ -4013,7 +3946,7 @@ namespace AMDiS {
/* values on child[1] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cd);
if (typ == 0) {
switch (lr_set) {
@@ -4473,7 +4406,8 @@ namespace AMDiS {
const DegreeOfFreedom **pds = el->getDof();
const DegreeOfFreedom **cds = el->getChild(0)->getDof();
- DegreeOfFreedom pd_o[20], cd, pd;
+ DegreeOfFreedom cd, pd;
+ vector<DegreeOfFreedom> pd_o(20);
basFct->getLocalIndices(el, admin, pd_o);
pd = (pds[0][0] < pds[1][0]) ? pds[node_e][n0_e] : pds[node_e][n0_e + 1];
cd = cds[0][0] < cds[3][0] ? cds[node_e + 2][n0_e + 1] : cds[node_e + 2][n0_e];
@@ -4541,13 +4475,12 @@ namespace AMDiS {
Element *el = list->getElement(0);
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
- DegreeOfFreedom pdof[15];
+ vector<DegreeOfFreedom> pdof(15), cdof(15);
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
// values on child[0]
- const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[pdof[9]] = (*drv)[cdof[4]];
(*drv)[pdof[10]] = (*drv)[cdof[2]];
(*drv)[pdof[12]] = (*drv)[cdof[14]];
@@ -4555,7 +4488,7 @@ namespace AMDiS {
// values on child[1]
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[pdof[11]] = (*drv)[cdof[7]];
(*drv)[pdof[13]] = (*drv)[cdof[14]];
@@ -4568,17 +4501,16 @@ namespace AMDiS {
el = list->getElement(1);
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
// values on neighbour's child[0]
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[pdof[12]] = (*drv)[cdof[14]];
(*drv)[pdof[14]] = (*drv)[cdof[7]];
// values on neighbour's child[1]
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[pdof[13]] = (*drv)[cdof[14]];
}
@@ -4596,13 +4528,12 @@ namespace AMDiS {
const Element *el = list->getElement(0);
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
- DegreeOfFreedom pdof[15];
+ vector<DegreeOfFreedom> pdof(15), cdof(15);
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
// values on child[0]
- const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[pdof[9]] = (*drv)[cdof[4]];
(*drv)[pdof[10]] = (*drv)[cdof[2]];
(*drv)[pdof[12]] = (*drv)[cdof[14]];
@@ -4610,7 +4541,7 @@ namespace AMDiS {
// values on child[1]
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[pdof[11]] = (*drv)[cdof[7]];
(*drv)[pdof[13]] = (*drv)[cdof[14]];
@@ -4623,17 +4554,16 @@ namespace AMDiS {
el = list->getElement(1);
basFct->getLocalIndices(el, admin, pdof);
+ basFct->getLocalIndices(el->getChild(0), admin, cdof);
// values on neighbour's child[0]
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[pdof[12]] = (*drv)[cdof[14]];
(*drv)[pdof[14]] = (*drv)[cdof[7]];
// values on neighbour's child[1]
- cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cdof);
(*drv)[pdof[13]] = (*drv)[cdof[14]];
}
@@ -4651,15 +4581,13 @@ namespace AMDiS {
const Element* el = list->getElement(0);
int typ = list->getType(0);
-
const DOFAdmin *admin = drv->getFeSpace()->getAdmin();
- DegreeOfFreedom pd[35];
+ vector<DegreeOfFreedom> pd(35), cd(35);
basFct->getLocalIndices(el, admin, pd);
+ basFct->getLocalIndices(el->getChild(0), admin, cd);
// values on child[0]
- const DegreeOfFreedom *cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
-
(*drv)[pd[4]] = (*drv)[cd[11]];
(*drv)[pd[5]] = (*drv)[cd[3]];
(*drv)[pd[28]] = (*drv)[cd[27]];
@@ -4670,7 +4598,7 @@ namespace AMDiS {
// values on child[1]
- cd = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cd);
if (typ == 0) {
// parent of el_type 0
@@ -4704,7 +4632,7 @@ namespace AMDiS {
// values on child[0]
- cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
+ basFct->getLocalIndices(el->getChild(0), admin, cd);
switch (lr_set) {
case 1:
@@ -4724,7 +4652,7 @@ namespace AMDiS {
// values on child[1]
- cd = basFct->getLocalIndices(el->getChild(1), admin, NULL);
+ basFct->getLocalIndices(el->getChild(1), admin, cd);
if (typ == 0) {
switch (lr_set) {
diff --git a/AMDiS/src/Lagrange.h b/AMDiS/src/Lagrange.h
index 52b56780cd5e778b082f47aa98d072a0f6f8ec13..310f7ab90eca813c95206712e42e04e52442885c 100644
--- a/AMDiS/src/Lagrange.h
+++ b/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;
-
/// Implements BasisFunction::getBound
void getBound(const ElInfo*, BoundaryType *) const;
@@ -116,13 +111,13 @@ namespace AMDiS {
}
/// Implements BasisFunction::getLocalIndices().
- const DegreeOfFreedom *getLocalIndices(const Element *el,
- const DOFAdmin *admin,
- DegreeOfFreedom *dofs) const;
+ void getLocalIndices(const Element *el,
+ const DOFAdmin *admin,
+ std::vector<DegreeOfFreedom> &dofs) const;
void getLocalDofPtrVec(const Element *el,
const DOFAdmin *admin,
- std::vector<const DegreeOfFreedom*>& vec) const;
+ vector<const DegreeOfFreedom*>& vec) const;
/// Implements BasisFunction::l2ScpFctBas
void l2ScpFctBas(Quadrature* q,
@@ -154,7 +149,7 @@ namespace AMDiS {
/// Sets used function pointers
void setFunctionPointer();
- /// Used by \ref getDOFIndices and \ref getVec
+ /// Used by \ref getVec
int* orderOfPositionIndices(const Element* el,
GeoIndex position,
int positionIndex) const;
diff --git a/AMDiS/src/OpenMP.h b/AMDiS/src/OpenMP.h
new file mode 100644
index 0000000000000000000000000000000000000000..c1a512e6ec3c58a6bd24a7966d10d2665927eb4c
--- /dev/null
+++ b/AMDiS/src/OpenMP.h
@@ -0,0 +1,103 @@
+// ============================================================================
+// == ==
+// == AMDiS - Adaptive multidimensional simulations ==
+// == ==
+// == http://www.amdis-fem.org ==
+// == ==
+// ============================================================================
+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+/** \file OpenMP.h */
+
+#ifndef AMDIS_OPENMP_H
+#define AMDIS_OPENMP_H
+
+#ifdef _OPENMP
+#include <omp.h>
+#include <vector>
+#endif
+
+namespace AMDiS {
+
+ using namespace std;
+
+#ifdef _OPENMP
+
+ template<typename T>
+ class ThreadPrivate {
+ public:
+ ThreadPrivate()
+ : data(omp_get_max_threads())
+ {}
+
+ ThreadPrivate(T val)
+ : data(omp_get_max_threads(), val)
+ {}
+
+ inline T& get()
+ {
+#if (DEBUG != 0)
+ if (omp_get_thread_num() >= data.size()) {
+ cout << "Error in ThreadPrivate::get()!\n";
+ exit(0);
+ }
+#endif
+ return data[omp_get_thread_num()];
+ }
+
+ inline void set(T& val)
+ {
+#if (DEBUG != 0)
+ if (omp_get_thread_num() >= data.size()) {
+ cout << "Error in ThreadPrivate::set()!\n";
+ exit(0);
+ }
+#endif
+ data[omp_get_thread_num()] = val;
+ }
+
+ private:
+ vector<T> data;
+ };
+
+#else
+
+ template<typename T>
+ class ThreadPrivate {
+ public:
+ ThreadPrivate() {}
+
+ ThreadPrivate(T val)
+ : data(val)
+ {}
+
+ inline T& get()
+ {
+ return data;
+ }
+
+ inline void set(T& val)
+ {
+ data = val;
+ }
+
+ private:
+ T data;
+ };
+
+#endif
+
+}
+
+#endif
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index f0d5ca7cf749b9acf19ee144aab416c993b326cd..7f192f36160c13e9e3f0032ac69d5fb240c16385 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -19,70 +19,16 @@
namespace AMDiS {
- const Flag Operator::MATRIX_OPERATOR = 1;
- const Flag Operator::VECTOR_OPERATOR = 2;
-
-
- int Operator::getQuadratureDegree(int order, FirstOrderType firstOrderType)
- {
- std::vector<OperatorTerm*>* terms = NULL;
-
- switch(order) {
- case 0:
- terms = &zeroOrder;
- break;
- case 1:
- if (firstOrderType == GRD_PHI)
- terms = &firstOrderGrdPhi;
- else
- terms = &firstOrderGrdPsi;
- break;
- case 2:
- terms = &secondOrder;
- break;
- }
-
- const BasisFunction *psi = rowFeSpace->getBasisFcts();
- const BasisFunction *phi = colFeSpace->getBasisFcts();
-
- int psiDegree = psi->getDegree();
- int phiDegree = phi->getDegree();
- int maxTermDegree = 0;
-
- for (int i = 0; i < static_cast<int>(terms->size()); i++) {
- OperatorTerm *term = (*terms)[i];
- maxTermDegree = std::max(maxTermDegree, term->degree);
- }
-
- return psiDegree + phiDegree - order + maxTermDegree;
- }
-
-
- Operator::Operator(Flag operatorType,
- const FiniteElemSpace *row,
- const FiniteElemSpace *col)
- {
- FUNCNAME("Operator::Operator()");
-
- std::cout << "\n\n\n============== ERROR MESSAGE ===============\n";
- std::cout << "You make use of the obsolete constructor Operator::Operator(Flag, FiniteElemSpace*, FiniteElemSpace*). Please\n";
- std::cout << "make use of the new constructor Operator::Operatpr(FiniteElemSpace*, FiniteElemSpace*). So you have just to\n";
- std::cout << "remove the MATRIX_OPERATOR/VECTOR_OPERATOR values from the constuctor call.\n\n\n";
-
- exit(0);
- }
-
-
Operator::Operator(const FiniteElemSpace *row, const FiniteElemSpace *col)
: rowFeSpace(row),
colFeSpace(col ? col : row),
fillFlag(Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA),
needDualTraverse(false),
+ assembler(NULL),
uhOld(NULL),
optimized(true)
{
- assembler = NULL;
secondOrder.resize(0);
firstOrderGrdPsi.resize(0);
firstOrderGrdPhi.resize(0);
@@ -104,10 +50,10 @@ namespace AMDiS {
ElementMatrix& userMat,
double factor)
{
- if (!assembler)
+ if (!assembler.get())
initAssembler(NULL, NULL, NULL, NULL);
- assembler->calculateElementMatrix(elInfo, userMat, factor);
+ assembler.get()->calculateElementMatrix(elInfo, userMat, factor);
}
@@ -117,13 +63,13 @@ namespace AMDiS {
ElementMatrix& userMat,
double factor)
{
- if (!assembler)
+ if (!assembler.get())
initAssembler(NULL, NULL, NULL, NULL);
- assembler->calculateElementMatrix(rowElInfo, colElInfo,
- smallElInfo, largeElInfo,
- rowColFeSpaceEqual,
- userMat, factor);
+ assembler.get()->calculateElementMatrix(rowElInfo, colElInfo,
+ smallElInfo, largeElInfo,
+ rowColFeSpaceEqual,
+ userMat, factor);
}
@@ -131,10 +77,10 @@ namespace AMDiS {
ElementVector& userVec,
double factor)
{
- if (!assembler)
+ if (!assembler.get())
initAssembler(NULL, NULL, NULL, NULL);
- assembler->calculateElementVector(elInfo, userVec, factor);
+ assembler.get()->calculateElementVector(elInfo, userVec, factor);
}
@@ -143,12 +89,12 @@ namespace AMDiS {
ElementVector& userVec,
double factor)
{
- if (!assembler)
+ if (!assembler.get())
initAssembler(NULL, NULL, NULL, NULL);
- assembler->calculateElementVector(mainElInfo, auxElInfo,
- smallElInfo, largeElInfo,
- userVec, factor);
+ assembler.get()->calculateElementVector(mainElInfo, auxElInfo,
+ smallElInfo, largeElInfo,
+ userVec, factor);
}
@@ -157,34 +103,74 @@ namespace AMDiS {
Quadrature *quad1GrdPhi,
Quadrature *quad0)
{
- if (optimized)
- assembler = new OptimizedAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi, quad0,
- rowFeSpace, colFeSpace);
- else
- assembler = new StandardAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi, quad0,
- rowFeSpace, colFeSpace);
+ if (optimized) {
+ Assembler *aptr =
+ new OptimizedAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi,
+ quad0, rowFeSpace, colFeSpace);
+ assembler.set(aptr);
+ } else {
+ Assembler *aptr =
+ new StandardAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi,
+ quad0, rowFeSpace, colFeSpace);
+ assembler.set(aptr);
+ }
+ }
+
+
+ int Operator::getQuadratureDegree(int order, FirstOrderType firstOrderType)
+ {
+ std::vector<OperatorTerm*>* terms = NULL;
+
+ switch(order) {
+ case 0:
+ terms = &zeroOrder;
+ break;
+ case 1:
+ if (firstOrderType == GRD_PHI)
+ terms = &firstOrderGrdPhi;
+ else
+ terms = &firstOrderGrdPsi;
+ break;
+ case 2:
+ terms = &secondOrder;
+ break;
+ }
+
+ const BasisFunction *psi = rowFeSpace->getBasisFcts();
+ const BasisFunction *phi = colFeSpace->getBasisFcts();
+
+ int psiDegree = psi->getDegree();
+ int phiDegree = phi->getDegree();
+ int maxTermDegree = 0;
+
+ for (int i = 0; i < static_cast<int>(terms->size()); i++) {
+ OperatorTerm *term = (*terms)[i];
+ maxTermDegree = std::max(maxTermDegree, term->degree);
+ }
+
+ return psiDegree + phiDegree - order + maxTermDegree;
}
void Operator::finishAssembling()
{
- if (assembler)
- assembler->finishAssembling();
+ if (assembler.get())
+ assembler.get()->finishAssembling();
}
- Assembler *Operator::getAssembler()
+ Assembler* Operator::getAssembler()
{
- if (!assembler)
+ if (!assembler.get())
initAssembler(NULL, NULL, NULL, NULL);
- return assembler;
+ return assembler.get();
}
void Operator::setAssembler(Assembler *a)
{
- assembler = a;
+ assembler.set(a);
}
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 8af1747e1264bb70ed0208aa6b4e8a32838e6bfe..068265c6bb9351e3bc9bb8dd8bc63eed28f06b50 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -53,12 +53,6 @@ namespace AMDiS {
class Operator
{
public:
- /// Obsolete consructor. Calling this constructor yields an error. Will be removed
- /// in on of the next revisions.
- Operator(Flag operatorType,
- const FiniteElemSpace *rowFeSpace,
- const FiniteElemSpace *colFeSpace = NULL);
-
/// Constructs an empty Operator of type operatorType for the given FiniteElemSpace.
Operator(const FiniteElemSpace *rowFeSpace,
const FiniteElemSpace *colFeSpace = NULL);
@@ -94,8 +88,7 @@ namespace AMDiS {
void addTerm(ZeroOrderTerm *term);
/// Adds a FirstOrderTerm to the Operator
- void addTerm(FirstOrderTerm *term,
- FirstOrderType type);
+ void addTerm(FirstOrderTerm *term, FirstOrderType type);
/// Adds a SecondOrderTerm to the Operator
void addTerm(SecondOrderTerm *term);
@@ -315,13 +308,6 @@ namespace AMDiS {
return zeroOrder.size() != 0;
}
- public:
- /// Constant type flag for matrix operators. Obsolete, will be removed!
- static const Flag MATRIX_OPERATOR;
-
- /// Constant type flag for vector operators. Obsolete, will be removed!
- static const Flag VECTOR_OPERATOR;
-
protected:
/// FiniteElemSpace for matrix rows and element vector
const FiniteElemSpace *rowFeSpace;
@@ -347,7 +333,7 @@ namespace AMDiS {
/// Calculates the element matrix and/or the element vector. It is
/// created especially for this Operator, when \ref getElementMatrix()
/// or \ref getElementVector is called for the first time.
- Assembler* assembler;
+ ThreadPrivate<Assembler*> assembler;
/// List of all second order terms
vector<OperatorTerm*> secondOrder;
diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc
index fc3a905813a369c19bd2f95858078ddeb0e1b2a8..ece5db9b64a2dc3d5422a4e1a52eb9c1d8585dbc 100644
--- a/AMDiS/src/ProblemStat.cc
+++ b/AMDiS/src/ProblemStat.cc
@@ -1247,9 +1247,8 @@ namespace AMDiS {
}
}
- OperatorPos opPos = {i, j, factor, estFactor, Operator::MATRIX_OPERATOR};
+ OperatorPos opPos = {i, j, factor, estFactor};
operators[op].push_back(opPos);
- opFlags[op].setFlag(Operator::MATRIX_OPERATOR);
}
@@ -1284,9 +1283,8 @@ namespace AMDiS {
}
}
- OperatorPos opPos = {i, -1, factor, estFactor, Operator::VECTOR_OPERATOR};
+ OperatorPos opPos = {i, -1, factor, estFactor};
operators[op].push_back(opPos);
- opFlags[op].setFlag(Operator::VECTOR_OPERATOR);
}
diff --git a/AMDiS/src/ProblemStat.h b/AMDiS/src/ProblemStat.h
index 4fd293e238cca5284d6d18df6392a54240f6614d..707e3b87dbefa97c16ca25e9ccc8e8b5d35cb046 100644
--- a/AMDiS/src/ProblemStat.h
+++ b/AMDiS/src/ProblemStat.h
@@ -45,7 +45,6 @@ namespace AMDiS {
{
int row, col;
double *factor, *estFactor;
- Flag operatorType;
};
@@ -656,8 +655,6 @@ namespace AMDiS {
bool writeAsmInfo;
map<Operator*, vector<OperatorPos> > operators;
-
- map<Operator*, Flag> opFlags;
};
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
diff --git a/AMDiS/src/Quadrature.cc b/AMDiS/src/Quadrature.cc
index c8f72b5675875a871d6c191221bf36103cbb3347..53169448ff78c6960a9c99dd9c27c00326c5bc59 100644
--- a/AMDiS/src/Quadrature.cc
+++ b/AMDiS/src/Quadrature.cc
@@ -1463,9 +1463,9 @@ namespace AMDiS {
double result = 0.0;
// calculate weighted sum over all quadrature-points
- for (int i = 0; i < n_points; i++) {
+ for (int i = 0; i < n_points; i++)
result += w[i] * (*f)((*lambda)[i]);
- }
+
return result;
}
@@ -1476,31 +1476,35 @@ namespace AMDiS {
{
FUNCNAME("FastQuadrature::provideFastQuuadrature()");
- FastQuadrature *quad_fast;
-
- list<FastQuadrature*>::iterator fast = fastQuadList.begin();
-
- for (fast = fastQuadList.begin(); fast != fastQuadList.end(); fast++)
- if ((*fast)->basisFunctions == bas_fcts && (*fast)->quadrature == &quad)
- break;
-
- if (fast != fastQuadList.end() && ((*fast)->init_flag & init_flag) == init_flag) {
- quad_fast = *fast;
- } else {
- if (fast == fastQuadList.end()) {
- quad_fast =
- new FastQuadrature(const_cast<BasisFunction*>(bas_fcts),
- const_cast<Quadrature*>(&quad), 0);
-
- fastQuadList.push_front(quad_fast);
-
- max_points = std::max(max_points, quad.getNumPoints());
+ FastQuadrature *quad_fast = NULL;
+
+#pragma omp critical
+ {
+ list<FastQuadrature*>::iterator fast = fastQuadList.begin();
+ for (; fast != fastQuadList.end(); fast++)
+ if ((*fast)->basisFunctions == bas_fcts &&
+ (*fast)->quadrature == &quad)
+ break;
+
+ if (fast != fastQuadList.end() &&
+ ((*fast)->init_flag & init_flag) == init_flag) {
+ quad_fast = *fast;
} else {
- quad_fast = (*fast);
+ if (fast == fastQuadList.end()) {
+ quad_fast =
+ new FastQuadrature(const_cast<BasisFunction*>(bas_fcts),
+ const_cast<Quadrature*>(&quad), 0);
+
+ fastQuadList.push_front(quad_fast);
+
+ max_points = std::max(max_points, quad.getNumPoints());
+ } else {
+ quad_fast = (*fast);
+ }
}
+
+ quad_fast->init(init_flag);
}
-
- quad_fast->init(init_flag);
return quad_fast;
}
diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc
index 9dd029e05bcd89121c7049f90a58113b6884e271..9b446d0ffb9e5d563b419c6099f76133ab3f0905 100644
--- a/AMDiS/src/SecondOrderAssembler.cc
+++ b/AMDiS/src/SecondOrderAssembler.cc
@@ -22,8 +22,12 @@
namespace AMDiS {
- std::vector<SubAssembler*> SecondOrderAssembler::optimizedSubAssemblers;
- std::vector<SubAssembler*> SecondOrderAssembler::standardSubAssemblers;
+ ThreadPrivate<vector<SubAssembler*> >
+ SecondOrderAssembler::optimizedSubAssemblers;
+
+ ThreadPrivate<vector<SubAssembler*> >
+ SecondOrderAssembler::standardSubAssemblers;
+
SecondOrderAssembler::SecondOrderAssembler(Operator *op,
Assembler *assembler,
@@ -44,23 +48,21 @@ namespace AMDiS {
SecondOrderAssembler *newAssembler;
- std::vector<SubAssembler*> *subAssemblers =
- optimized ?
- &optimizedSubAssemblers :
- &standardSubAssemblers;
+ vector<SubAssembler*> &subAssemblers =
+ optimized ? optimizedSubAssemblers.get() : standardSubAssemblers.get();
- std::vector<OperatorTerm*> opTerms = op->zeroOrder;
+ vector<OperatorTerm*> opTerms = op->zeroOrder;
sort(opTerms.begin(), opTerms.end());
// 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<SecondOrderAssembler*>((*subAssemblers)[i]);
+ if (opTerms == assTerms && subAssemblers[i]->getQuadrature() == quad)
+ return dynamic_cast<SecondOrderAssembler*>(subAssemblers[i]);
}
// check if all terms are pw_const
@@ -83,7 +85,7 @@ namespace AMDiS {
}
}
- subAssemblers->push_back(newAssembler);
+ subAssemblers.push_back(newAssembler);
return newAssembler;
}
@@ -286,7 +288,7 @@ namespace AMDiS {
int nPoints = quadrature->getNumPoints();
- std::vector<mtl::dense2D<double> > LALt(nPoints);
+ vector<mtl::dense2D<double> > LALt(nPoints);
for (int iq = 0; iq < nPoints; iq++) {
LALt[iq].change_dim(dim + 1, dim + 1);
LALt[iq] = 0.0;
diff --git a/AMDiS/src/SecondOrderAssembler.h b/AMDiS/src/SecondOrderAssembler.h
index 8fc6bc6ed26c9eb9b43284b914ace191d74eb61e..140deae3ab8410b15fb2bfe6c6889a09a5b12388 100644
--- a/AMDiS/src/SecondOrderAssembler.h
+++ b/AMDiS/src/SecondOrderAssembler.h
@@ -65,10 +65,10 @@ namespace AMDiS {
protected:
/// List of all yet created optimized second order assemblers.
- static vector<SubAssembler*> optimizedSubAssemblers;
+ static ThreadPrivate<vector<SubAssembler*> > optimizedSubAssemblers;
/// List of all yet created standard second order assemblers.
- static vector<SubAssembler*> standardSubAssemblers;
+ static ThreadPrivate<vector<SubAssembler*> > standardSubAssemblers;
};
diff --git a/AMDiS/src/SubAssembler.cc b/AMDiS/src/SubAssembler.cc
index 92ab6acfb7e7a95e5ebb309a503d3b6dc2f1ef8e..4a596bb012eb87ff660db6a5fa153d3ecf117f25 100644
--- a/AMDiS/src/SubAssembler.cc
+++ b/AMDiS/src/SubAssembler.cc
@@ -93,8 +93,11 @@ namespace AMDiS {
if (!quadFast) {
quadFast = FastQuadrature::provideFastQuadrature(psi, *quadrature, updateFlag);
} else {
- if (!quadFast->initialized(updateFlag))
- quadFast->init(updateFlag);
+#pragma omp critical
+ {
+ if (!quadFast->initialized(updateFlag))
+ quadFast->init(updateFlag);
+ }
}
return quadFast;
diff --git a/AMDiS/src/UmfPackSolver.h b/AMDiS/src/UmfPackSolver.h
index 74461e3aa6fd6a35b1753c74281d9034bc43222f..3a6643a635c095893b0df7d535bd682182de78e1 100644
--- a/AMDiS/src/UmfPackSolver.h
+++ b/AMDiS/src/UmfPackSolver.h
@@ -92,7 +92,7 @@ namespace AMDiS {
double residual = two_norm(r);
oem.setResidual(residual);
oem.setErrorCode(code);
- std::cout << "UmfPackSolver: ||b-Ax|| = " << residual << "\n";
+ MSG("UmfPackSolver: ||b-Ax|| = %e\n", residual);
TEST_EXIT(residual <= oem.getTolerance())("Tolerance tol = %e could not be reached!\n Set tolerance by '->solver->tolerance:' \n", oem.getTolerance());
}
return code;
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index a83b09023074249c389fe7c38f6e8ddc3e67ec47..853724edd42a4eef925501ee61b737c06dfbb610 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -22,8 +22,12 @@
namespace AMDiS {
- std::vector<SubAssembler*> ZeroOrderAssembler::optimizedSubAssemblers;
- std::vector<SubAssembler*> ZeroOrderAssembler::standardSubAssemblers;
+ ThreadPrivate<vector<SubAssembler*> >
+ ZeroOrderAssembler::optimizedSubAssemblers;
+
+ ThreadPrivate<vector<SubAssembler*> >
+ ZeroOrderAssembler::standardSubAssemblers;
+
ZeroOrderAssembler::ZeroOrderAssembler(Operator *op,
Assembler *assembler,
@@ -44,22 +48,22 @@ namespace AMDiS {
ZeroOrderAssembler *newAssembler;
- std::vector<SubAssembler*> *subAssemblers =
- optimized ? &optimizedSubAssemblers : &standardSubAssemblers;
+ vector<SubAssembler*> subAssemblers =
+ optimized ? optimizedSubAssemblers.get() : standardSubAssemblers.get();
- std::vector<OperatorTerm*> opTerms = op->zeroOrder;
+ vector<OperatorTerm*> opTerms = op->zeroOrder;
sort(opTerms.begin(), opTerms.end());
// check if a new assembler is needed
if (quad) {
- for (int i = 0; i < static_cast<int>(subAssemblers->size()); i++) {
- std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
+ for (int i = 0; i < static_cast<int>(subAssemblers.size()); i++) {
+ vector<OperatorTerm*> assTerms = *(subAssemblers[i]->getTerms());
sort(assTerms.begin(), assTerms.end());
- if (opTerms == assTerms && (*subAssemblers)[i]->getQuadrature() == quad)
- return dynamic_cast<ZeroOrderAssembler*>((*subAssemblers)[i]);
+ if (opTerms == assTerms && subAssemblers[i]->getQuadrature() == quad)
+ return dynamic_cast<ZeroOrderAssembler*>(subAssemblers[i]);
}
}
@@ -82,7 +86,8 @@ namespace AMDiS {
newAssembler = new FastQuadZOA(op, assembler, quad);
}
- subAssemblers->push_back(newAssembler);
+ subAssemblers.push_back(newAssembler);
+
return newAssembler;
}
@@ -102,9 +107,9 @@ namespace AMDiS {
int nPoints = quadrature->getNumPoints();
ElementVector c(nPoints, 0.0);
- std::vector<double> phival(nCol);
+ vector<double> phival(nCol);
- for (std::vector<OperatorTerm*>::iterator termIt = terms.begin();
+ for (vector<OperatorTerm*>::iterator termIt = terms.begin();
termIt != terms.end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
@@ -151,7 +156,7 @@ namespace AMDiS {
int nPoints = quadrature->getNumPoints();
ElementVector c(nPoints, 0.0);
- std::vector<OperatorTerm*>::iterator termIt;
+ vector<OperatorTerm*>::iterator termIt;
for (termIt = terms.begin(); termIt != terms.end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
@@ -190,7 +195,7 @@ namespace AMDiS {
c.change_dim(nPoints);
c = 0.0;
- for (std::vector<OperatorTerm*>::iterator termIt = terms.begin();
+ for (vector<OperatorTerm*>::iterator termIt = terms.begin();
termIt != terms.end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
@@ -240,7 +245,7 @@ namespace AMDiS {
}
ElementVector c(nPoints, 0.0);
- std::vector<OperatorTerm*>::iterator termIt;
+ vector<OperatorTerm*>::iterator termIt;
for (termIt = terms.begin(); termIt != terms.end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
@@ -313,7 +318,7 @@ namespace AMDiS {
firstCall = false;
}
- std::vector<OperatorTerm*>::iterator termIt;
+ vector<OperatorTerm*>::iterator termIt;
ElementVector c(1, 0.0);
for (termIt = terms.begin(); termIt != terms.end(); ++termIt)
diff --git a/AMDiS/src/ZeroOrderAssembler.h b/AMDiS/src/ZeroOrderAssembler.h
index 703479fadaa88d2e7392615e3de2fb5629d5de74..43970c64b54835acb8dc964051f9bc6a3ef01d23 100644
--- a/AMDiS/src/ZeroOrderAssembler.h
+++ b/AMDiS/src/ZeroOrderAssembler.h
@@ -29,6 +29,8 @@
namespace AMDiS {
+ using namespace std;
+
/**
* \ingroup Assembler
*
@@ -63,10 +65,10 @@ namespace AMDiS {
protected:
/// List of all yet created optimized SubAssembler objects.
- static std::vector<SubAssembler*> optimizedSubAssemblers;
+ static ThreadPrivate<vector<SubAssembler*> > optimizedSubAssemblers;
/// List of all yet created standard SubAssembler objects.
- static std::vector<SubAssembler*> standardSubAssemblers;
+ static ThreadPrivate<vector<SubAssembler*> > standardSubAssemblers;
};
diff --git a/AMDiS/src/compositeFEM/CompositeFEMMethods.cc b/AMDiS/src/compositeFEM/CompositeFEMMethods.cc
index a6df47aa38489034dcc3d6d0df88b1e9a22ee534..5ea0c82adb7ccb1323ce4b72b6cfb8c95c746052 100644
--- a/AMDiS/src/compositeFEM/CompositeFEMMethods.cc
+++ b/AMDiS/src/compositeFEM/CompositeFEMMethods.cc
@@ -41,6 +41,7 @@ void CompositeFEMMethods::setPosLsToVal(DOFVector<double> *dof,
}
}
+
void CompositeFEMMethods::setPosLsToFct(DOFVector<double> *dof,
const AbstractFunction<double, WorldVector<double> > *fct,
const DOFVector<double> *lsFct_dof)
@@ -48,12 +49,12 @@ void CompositeFEMMethods::setPosLsToFct(DOFVector<double> *dof,
const BasisFunction *basisFcts = dof->getFeSpace()->getBasisFcts();
const DOFAdmin *admin = dof->getFeSpace()->getAdmin();
const int dim = dof->getFeSpace()->getMesh()->getDim();
- const DegreeOfFreedom *locInd;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(dof->getFeSpace()->getMesh(), -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
ElementVector locVec(basisFcts->getNumber());
+ std::vector<DegreeOfFreedom> locInd(basisFcts->getNumber());
while (elInfo) {
const Element *el = elInfo->getElement();
@@ -62,7 +63,7 @@ void CompositeFEMMethods::setPosLsToFct(DOFVector<double> *dof,
lsFct_dof->getLocalVector(el, locVec);
// Get dof indices of vertices.
- locInd = basisFcts->getLocalIndices(el, admin, NULL);
+ basisFcts->getLocalIndices(el, admin, locInd);
for (int i = 0; i <= dim; i++) {
// Is vertex in domain with positive level set function values ?
@@ -74,6 +75,7 @@ void CompositeFEMMethods::setPosLsToFct(DOFVector<double> *dof,
}
}
+
void CompositeFEMMethods::printBoundaryElements(const std::string fn_str,
ElementLevelSet *elLS,
FiniteElemSpace *feSpace)
@@ -95,7 +97,7 @@ void CompositeFEMMethods::printBoundaryElements(const std::string fn_str,
WorldVector<double> coord;
const int nBasFcts = feSpace->getBasisFcts()->getNumber();
- DegreeOfFreedom *locInd = new DegreeOfFreedom[nBasFcts];
+ std::vector<DegreeOfFreedom> locInd(nBasFcts);
ElInfo *loc_elInfo = stack.traverseFirst(feSpace->getMesh(),
-1,
@@ -105,10 +107,9 @@ void CompositeFEMMethods::printBoundaryElements(const std::string fn_str,
while (loc_elInfo) {
// Get local indices of vertices.
- feSpace->getBasisFcts()->getLocalIndices(
- const_cast<Element *>(loc_elInfo->getElement()),
- const_cast<DOFAdmin *>(feSpace->getAdmin()),
- locInd);
+ feSpace->getBasisFcts()->getLocalIndices(const_cast<Element *>(loc_elInfo->getElement()),
+ const_cast<DOFAdmin *>(feSpace->getAdmin()),
+ locInd);
// Get element status.
elStatus = elLS->createElementLevelSet(loc_elInfo);
@@ -140,8 +141,6 @@ void CompositeFEMMethods::printBoundaryElements(const std::string fn_str,
} // end of: mesh traverse
- delete [] locInd;
-
boundaryOut << "\nNumber of boundary elements: \t" << boundEl_cntr << "\n";
boundaryOut.close();
diff --git a/AMDiS/src/compositeFEM/CompositeFEMOperator.cc b/AMDiS/src/compositeFEM/CompositeFEMOperator.cc
index cb572ba3d7d335f90e4e01ea68afd7a7e615f715..2c8c32e7a00a289c9ec59638cd551dff5e4ee430 100644
--- a/AMDiS/src/compositeFEM/CompositeFEMOperator.cc
+++ b/AMDiS/src/compositeFEM/CompositeFEMOperator.cc
@@ -119,9 +119,9 @@ void CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
subElementAssembler->getNCol());
set_to_zero(subPolMat2);
- if (!assembler) {
- assembler =
- new StandardAssembler(this, NULL, NULL, NULL, NULL, rowFeSpace, colFeSpace);
+ if (!assembler.get()) {
+ Assembler *aptr = new StandardAssembler(this, NULL, NULL, NULL, NULL, rowFeSpace, colFeSpace);
+ assembler.set(aptr);
}
if (elLS->getLevelSetDomain() ==
@@ -132,7 +132,7 @@ void CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
}
- assembler->calculateElementMatrix(elInfo, elMat, 1.0);
+ assembler.get()->calculateElementMatrix(elInfo, elMat, 1.0);
subElementAssembler->getSubPolytopeMatrix(subPolytope,
subElementAssembler,
elInfo,
@@ -244,17 +244,16 @@ void CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
elInfo,
subPolVec1);
- /**
- * Integration on second subpolytope produced by the intersection.
- */
+ // Integration on second subpolytope produced by the intersection.
ElementVector elVec(subElementAssembler->getNRow());
set_to_zero(elVec);
ElementVector subPolVec2(subElementAssembler->getNRow());
set_to_zero(subPolVec2);
- if (!assembler)
- assembler =
- new StandardAssembler(this, NULL, NULL, NULL, NULL, rowFeSpace, colFeSpace);
+ if (!assembler.get()) {
+ Assembler *aptr = new StandardAssembler(this, NULL, NULL, NULL, NULL, rowFeSpace, colFeSpace);
+ assembler.set(aptr);
+ }
if (elLS->getLevelSetDomain() ==
ElementLevelSet::LEVEL_SET_INTERIOR) {
@@ -263,7 +262,7 @@ void CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
}
- assembler->calculateElementVector(elInfo, elVec, 1.0);
+ assembler.get()->calculateElementVector(elInfo, elVec, 1.0);
subElementAssembler->getSubPolytopeVector(subPolytope,
subElementAssembler,
elInfo,
diff --git a/AMDiS/src/compositeFEM/CompositeFEMOperator.h b/AMDiS/src/compositeFEM/CompositeFEMOperator.h
index b7346873487b10b868bfa43fd303ec4c817f0e37..a1330d29fd0c1c9899833a0c2af584db7f078658 100644
--- a/AMDiS/src/compositeFEM/CompositeFEMOperator.h
+++ b/AMDiS/src/compositeFEM/CompositeFEMOperator.h
@@ -55,11 +55,10 @@ class CompositeFEMOperator : public Operator
{
public:
/// Constructor.
- CompositeFEMOperator(Flag operatorType_,
- ElementLevelSet *elLS_,
+ CompositeFEMOperator(ElementLevelSet *elLS_,
const FiniteElemSpace *rowFeSpace_,
const FiniteElemSpace *colFeSpace_ = NULL)
- : Operator(operatorType_, rowFeSpace_, colFeSpace_),
+ : Operator(rowFeSpace_, colFeSpace_),
elLS(elLS_),
subElementAssembler(NULL),
elStatus(ElementLevelSet::LEVEL_SET_UNDEFINED)
diff --git a/AMDiS/src/compositeFEM/PenaltyOperator.h b/AMDiS/src/compositeFEM/PenaltyOperator.h
index 79b88c4a3abac9830f88b5455ac098915ef523aa..da31a3ee6b44ac998fd562c53bd2eab34ee2f157 100644
--- a/AMDiS/src/compositeFEM/PenaltyOperator.h
+++ b/AMDiS/src/compositeFEM/PenaltyOperator.h
@@ -41,13 +41,12 @@ class PenaltyOperator : public Operator
{
public:
/// Constructor.
- PenaltyOperator(Flag operatorType_,
- ElementLevelSet *elLS_,
+ PenaltyOperator(ElementLevelSet *elLS_,
double factor_,
bool penaltyCoeffFlag_,
FiniteElemSpace *rowFeSpace_,
FiniteElemSpace *colFeSpace_ = NULL)
- : Operator(operatorType_, rowFeSpace_, colFeSpace_),
+ : Operator(rowFeSpace_, colFeSpace_),
elLS(elLS_),
elStatus(ElementLevelSet::LEVEL_SET_UNDEFINED),
factor(factor_),
diff --git a/AMDiS/src/io/DataCollector.h b/AMDiS/src/io/DataCollector.h
index c25fc7982ca87bc11417e5fad2f89822526baf4e..e468b95b1f4cc07a1b069348f743d58a9b37639c 100644
--- a/AMDiS/src/io/DataCollector.h
+++ b/AMDiS/src/io/DataCollector.h
@@ -181,10 +181,8 @@ namespace AMDiS {
/// Stores for each simplex the interpolation points.
std::vector<std::vector<int> > interpPoints;
- /** \brief
- * Stores for each DOF a list of its coordinates. If there are now periodic
- * boundaries than there is also only one coordinate per DOF.
- */
+ /// Stores for each DOF a list of its coordinates. If there are now periodic
+ /// boundaries than there is also only one coordinate per DOF.
DOFVector<std::list<WorldVector<double> > > *interpPointCoords;
/// list of coords for each dof
@@ -199,10 +197,8 @@ namespace AMDiS {
/// Stores a list of vertex infos for each dof.
DOFVector<std::list<VertexInfo> > *vertexInfos;
- /** \brief
- * periodicConnections[i][j] stores whether the connection at side j of
- * the element with output index i has already been written.
- */
+ /// periodicConnections[i][j] stores whether the connection at side j of
+ /// the element with output index i has already been written.
std::vector<DimVec<bool> > periodicConnections;
/// Stores if element data was collected before.
@@ -217,9 +213,6 @@ namespace AMDiS {
/// Pointer to a function which decides whether an element is considered.
bool (*writeElem)(ElInfo*);
- /// Temporary variable used in functions addValueData() and addInterpData().
- DegreeOfFreedom *localDOFs;
-
/// Temporary variable used in functions addValueData() and addInterpData().
BasisFunction *basisFcts;
diff --git a/AMDiS/src/io/DataCollector.hh b/AMDiS/src/io/DataCollector.hh
index f253913c7fc6a9e93528e3e9afee1f6404e0d93a..5dab35394b74b441d78118759f7f7b83c4875672 100644
--- a/AMDiS/src/io/DataCollector.hh
+++ b/AMDiS/src/io/DataCollector.hh
@@ -138,12 +138,10 @@ namespace AMDiS {
basisFcts = const_cast<BasisFunction*>(feSpace->getBasisFcts());
nBasisFcts = basisFcts->getNumber();
- localDOFs = new DegreeOfFreedom[nBasisFcts];
- TraverseStack stack;
-
// Traverse elements to add value information and to mark
// interpolation points.
+ TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, level,
traverseFlag | Mesh::FILL_COORDS);
while (elInfo) {
@@ -171,7 +169,6 @@ namespace AMDiS {
}
}
- delete [] localDOFs;
valueDataCollected = true;
}
@@ -299,9 +296,10 @@ namespace AMDiS {
{
FUNCNAME("DataCollector<T>::addValueData()");
- WorldVector<double> vertexCoords;
+ vector<DegreeOfFreedom> localDOFs(basisFcts->getNumber());
basisFcts->getLocalIndices(elInfo->getElement(), localAdmin, localDOFs);
+ WorldVector<double> vertexCoords;
// First, traverse all DOFs at the vertices of the element, determine
// their coordinates and add them to the corresponding entry in dofCoords.
for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
@@ -358,6 +356,7 @@ namespace AMDiS {
{
FUNCNAME("DataCollector<T>::addInterpData()");
+ vector<DegreeOfFreedom> localDOFs(basisFcts->getNumber());
basisFcts->getLocalIndices(elInfo->getElement(), localAdmin, localDOFs);
std::vector<int> elemInterpPoints(0);
diff --git a/AMDiS/src/io/GNUPlotWriter.cc b/AMDiS/src/io/GNUPlotWriter.cc
index 211c4f4a7e892d97e40981f6ecbf95f743de9fc8..c6f0484d9b77684aef21d0d1936077c616d1f645 100644
--- a/AMDiS/src/io/GNUPlotWriter.cc
+++ b/AMDiS/src/io/GNUPlotWriter.cc
@@ -40,7 +40,7 @@ namespace AMDiS {
int dow = Global::getGeo(WORLD);
const BasisFunction *basFcts = feSpace_->getBasisFcts();
int numFcts = basFcts->getNumber();
- DegreeOfFreedom *localDOFs = new DegreeOfFreedom[numFcts];
+ std::vector<DegreeOfFreedom> localDOFs(numFcts);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1,
@@ -56,8 +56,6 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
- delete [] localDOFs;
-
FILE *file = NULL;
if (!(file = fopen(filename_.c_str(), "w")))
ERROR("could not open file %s for writing\n", filename_.c_str());
diff --git a/AMDiS/src/io/PovrayWriter.cc b/AMDiS/src/io/PovrayWriter.cc
index 244eb26428b6ea0ca065016f76e40cd9ea2617cd..585bab629d70c9548177468b6905b22bae22e47b 100644
--- a/AMDiS/src/io/PovrayWriter.cc
+++ b/AMDiS/src/io/PovrayWriter.cc
@@ -26,9 +26,10 @@ using namespace std;
namespace AMDiS {
/* destructor */
- PovrayWriter::~PovrayWriter(){
+ PovrayWriter::~PovrayWriter()
+ {
// free the memory for the bounding box
- if(bBox==NULL){
+ if (bBox==NULL) {
delete bBox;
bBox = NULL;
}
@@ -46,7 +47,7 @@ namespace AMDiS {
DOFVector<double> *values = dataCollector->getValues();
const FiniteElemSpace *feSpace = dataCollector->getFeSpace();
const BasisFunction *basFcts = feSpace->getBasisFcts();
- DegreeOfFreedom *dofs = new DegreeOfFreedom[basFcts->getNumber()];
+ std::vector<DegreeOfFreedom> dofs(basFcts->getNumber());
/* create an 'iterator' */
TraverseStack stack;
@@ -93,8 +94,6 @@ namespace AMDiS {
out << endl;
}
out << "\t}" << endl; // end of texture_list
-
- delete [] dofs;
}
// provides the bounding box of the mesh
diff --git a/AMDiS/src/parallel/PetscSolver.cc b/AMDiS/src/parallel/PetscSolver.cc
index 58b2be7afd11ca40361c965b5dd641f0d3e4d047..904f8cba3122e40d64ba8b3442d5f48f299dd4a4 100644
--- a/AMDiS/src/parallel/PetscSolver.cc
+++ b/AMDiS/src/parallel/PetscSolver.cc
@@ -88,4 +88,72 @@ namespace AMDiS {
VecScatterDestroy(&scatter);
}
+
+ bool PetscSolver::testMatrixSymmetric(Mat mat, bool advancedTest)
+ {
+ FUNCNAME("PetscSolver::testMatrixSymmetric()");
+
+ Mat matTrans;
+ MatTranspose(mat, MAT_INITIAL_MATRIX, &matTrans);
+
+ if (advancedTest) {
+ int rowStart, rowEnd;
+ MatGetOwnershipRange(mat, &rowStart, &rowEnd);
+
+ PetscInt nCols, nColsTrans;
+ const PetscInt *cols, *colsTrans;
+ const PetscScalar *vals, *valsTrans;
+
+ for (int i = rowStart; i < rowEnd; i++) {
+ bool isSym = true;
+ MatGetRow(mat, i, &nCols, &cols, &vals);
+ MatGetRow(matTrans, i, &nColsTrans, &colsTrans, &valsTrans);
+
+ if (nCols != nColsTrans) {
+ MSG("Symmetric test fails: mat row %d: %d nnz mat col %d: %d nnz\n",
+ i, nCols, i, nColsTrans);
+ isSym = false;
+ } else {
+ for (int j = 0; j < nCols; j++) {
+ if (cols[j] != colsTrans[j]) {
+ if (cols[j] > colsTrans[j]) {
+ MSG("mat[%d][%d] does not exists: mat[%d][%d] = %e\n",
+ i, colsTrans[j],
+ colsTrans[j], i, valsTrans[j]);
+ } else {
+ MSG("mat[%d][%d] does not exists: mat[%d][%d] = %e\n",
+ cols[j], i,
+ i, cols[j], vals[j]);
+ }
+ isSym = false;
+ } else {
+ double n = fabs(vals[j] - valsTrans[j]);
+ if (n > 1e-10) {
+ MSG("value diff: mat[%d][%d] = %e mat[%d][%d] = %e\n",
+ i, cols[j], vals[j], colsTrans[j], i, valsTrans[j]);
+
+ isSym = false;
+ }
+ }
+ }
+ }
+
+ MatRestoreRow(mat, i, &nCols, &cols, &vals);
+ MatRestoreRow(matTrans, i, &nColsTrans, &colsTrans, &valsTrans);
+
+ if (!isSym)
+ return false;
+ }
+ }
+
+ MatAXPY(matTrans, -1.0, mat, DIFFERENT_NONZERO_PATTERN);
+ double norm1, norm2;
+ MatNorm(matTrans, NORM_FROBENIUS, &norm1);
+ MatNorm(matTrans, NORM_INFINITY, &norm2);
+ MatDestroy(&matTrans);
+
+ MSG("Matrix norm test: %e %e\n", norm1, norm2);
+
+ return (norm1 < 1e-10 && norm2 < 1e-10);
+ }
}
diff --git a/AMDiS/src/parallel/PetscSolver.h b/AMDiS/src/parallel/PetscSolver.h
index 1156591395c071d95cd2c443ec475197fbb715b0..2a541b8f0e72fbec2a539154717f8f2bd9c68db9 100644
--- a/AMDiS/src/parallel/PetscSolver.h
+++ b/AMDiS/src/parallel/PetscSolver.h
@@ -141,6 +141,8 @@ namespace AMDiS {
void copyVec(Vec& originVec, Vec& destVec,
vector<int>& originIndex, vector<int>& destIndex);
+ /// Run test, if matrix is symmetric.
+ bool testMatrixSymmetric(Mat mat, bool advancedTest = false);
protected:
/// PETSc solver object
KSP kspInterior;
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index d0d091cdfdfb20fadec83b962258d01ab7e0b7a7..64bdc7e27b2000d1dd258194fb7367bd2b3e57d8 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -25,6 +25,32 @@ namespace AMDiS {
using namespace std;
+ PetscErrorCode KSPMonitorFeti(KSP ksp, PetscInt n, PetscReal rnorm, void *dummy)
+ {
+ Vec Br,v,w;
+ Mat A;
+
+ KSPGetOperators(ksp, &A, PETSC_NULL, PETSC_NULL);
+ MatGetVecs(A, &w, &v);
+ KSPBuildResidual(ksp, v, w, &Br);
+
+ Vec nest0, nest1;
+ VecNestGetSubVec(Br, 0, &nest0);
+ VecNestGetSubVec(Br, 1, &nest1);
+
+ PetscScalar norm0, norm1;
+ VecNorm(nest0, NORM_2, &norm0);
+ VecNorm(nest1, NORM_2, &norm1);
+
+ VecDestroy(&v);
+ VecDestroy(&w);
+
+ PetscPrintf(PETSC_COMM_WORLD, "%3D KSP Component residual norm [%1.12e %1.12e] and relative [%1.12e]\n",
+ norm0, norm1, rnorm);
+
+ return 0;
+ }
+
PetscSolverFeti::PetscSolverFeti()
: PetscSolver(),
schurPrimalSolver(0),
@@ -122,6 +148,27 @@ namespace AMDiS {
}
+ void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
+ {
+ FUNCNAME("PetscSolverFeti::createDirichletData()");
+
+ int nComponents = mat.getSize();
+ for (int i = 0; i < nComponents; i++) {
+ DOFMatrix* dofMat = mat[i][i];
+ if (!dofMat)
+ continue;
+
+ const FiniteElemSpace *feSpace = dofMat->getRowFeSpace();
+ std::set<DegreeOfFreedom>& dRows = dofMat->getDirichletRows();
+ if (dirichletRows.count(feSpace)) {
+ WARNING("Implement test that for all components dirichlet rows are equal!\n");
+ } else {
+ dirichletRows[feSpace] = dRows;
+ }
+ }
+ }
+
+
void PetscSolverFeti::createFetiData()
{
FUNCNAME("PetscSolverFeti::createFetiData()");
@@ -235,9 +282,9 @@ namespace AMDiS {
lagrangeMap[feSpace].nRankDofs,
lagrangeMap[feSpace].nOverallDofs);
- TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() ==
- static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
- ("Should not happen!\n");
+// TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() ==
+// static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
+// ("Should not happen!\n");
}
}
@@ -272,15 +319,11 @@ namespace AMDiS {
DofIndexSet primals;
for (DofContainerSet::iterator it = vertices.begin();
it != vertices.end(); ++it) {
- if (meshLevel == 0) {
- MSG("WARNING: Modified primal detection algorithm!\n");
+ if (dirichletRows[feSpace].count(**it))
+ continue;
- double e = 1e-3;
- WorldVector<double> c;
- feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);
-
- if (fabs(c[0]) > e && fabs(c[1]) > e && fabs(c[0] - 1.0) > e && fabs(c[1] - 1.0) > e)
- primals.insert(**it);
+ if (meshLevel == 0) {
+ primals.insert(**it);
} else {
double e = 1e-8;
WorldVector<double> c;
@@ -318,16 +361,22 @@ namespace AMDiS {
DofContainer allBoundaryDofs;
meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
-
+
for (DofContainer::iterator it = allBoundaryDofs.begin();
- it != allBoundaryDofs.end(); ++it)
- if (!isPrimal(feSpace, **it))
- if (meshLevel == 0) {
+ it != allBoundaryDofs.end(); ++it) {
+ if (dirichletRows[feSpace].count(**it))
+ continue;
+
+ if (isPrimal(feSpace, **it))
+ continue;
+
+ if (meshLevel == 0) {
+ dualDofMap[feSpace].insertRankDof(**it);
+ } else {
+ if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
dualDofMap[feSpace].insertRankDof(**it);
- } else {
- if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
- dualDofMap[feSpace].insertRankDof(**it);
- }
+ }
+ }
}
@@ -342,11 +391,15 @@ namespace AMDiS {
meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
for (DofContainer::iterator it = allBoundaryDofs.begin();
- it != allBoundaryDofs.end(); ++it)
+ it != allBoundaryDofs.end(); ++it) {
+ if (dirichletRows[feSpace].count(**it))
+ continue;
+
if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
interfaceDofMap[feSpace].insertRankDof(**it);
else
- interfaceDofMap[feSpace].insertNonRankDof(**it);
+ interfaceDofMap[feSpace].insertNonRankDof(**it);
+ }
}
@@ -503,41 +556,44 @@ namespace AMDiS {
int nLocalInterior = 0;
for (int i = 0; i < admin->getUsedSize(); i++) {
- if (admin->isDofFree(i) == false &&
- isPrimal(feSpace, i) == false &&
- isDual(feSpace, i) == false &&
- isInterface(feSpace, i) == false) {
- if (meshLevel == 0) {
- localDofMap[feSpace].insertRankDof(i, nLocalInterior);
+ if (admin->isDofFree(i) ||
+ isPrimal(feSpace, i) ||
+ isDual(feSpace, i) ||
+ isInterface(feSpace, i) ||
+ dirichletRows[feSpace].count(i))
+ continue;
- if (fetiPreconditioner == FETI_DIRICHLET)
- interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
-
- nLocalInterior++;
- } else {
- if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
- localDofMap[feSpace].insertRankDof(i);
- else
- localDofMap[feSpace].insertNonRankDof(i);
-
- TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
- ("Not yet implemnted!\n");
- }
+ if (meshLevel == 0) {
+ localDofMap[feSpace].insertRankDof(i, nLocalInterior);
+
+ if (fetiPreconditioner == FETI_DIRICHLET)
+ interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
+
+ nLocalInterior++;
+ } else {
+ if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
+ localDofMap[feSpace].insertRankDof(i);
+ else
+ localDofMap[feSpace].insertNonRankDof(i);
+
+ TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
+ ("Not yet implemnted!\n");
}
}
// === And finally, add the global indicies of all dual nodes. ===
for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
- it != dualDofMap[feSpace].getMap().end(); ++it)
- if (meshLevel == 0)
+ it != dualDofMap[feSpace].getMap().end(); ++it) {
+ if (meshLevel == 0) {
localDofMap[feSpace].insertRankDof(it->first);
- else {
+ } else {
if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it->first))
localDofMap[feSpace].insertRankDof(it->first);
else
localDofMap[feSpace].insertNonRankDof(it->first);
}
+ }
}
@@ -687,7 +743,7 @@ namespace AMDiS {
VecAssemblyBegin(tmpVec);
VecAssemblyEnd(tmpVec);
-
+
subdomain->solve(tmpVec, tmpVec);
PetscScalar *tmpValues;
@@ -810,6 +866,8 @@ namespace AMDiS {
KSPSetFromOptions(ksp_feti);
+
+
// === Create FETI-DP preconditioner object. ===
if (fetiPreconditioner != FETI_NONE || stokesMode) {
@@ -1458,8 +1516,12 @@ namespace AMDiS {
initialize(feSpaces);
+ createDirichletData(*mat);
+
createFetiData();
+ dirichletRows.clear();
+
// === Create matrices for the FETI-DP method. ===
@@ -1470,6 +1532,14 @@ namespace AMDiS {
subdomain->fillPetscMatrix(mat);
+ // MSG("MATRIX IS SYM: %d\n", testMatrixSymmetric(subdomain->getMatInterior(), true));
+
+// if (MPI::COMM_WORLD.Get_rank() == 0)
+// MatView(subdomain->getMatInterior(), PETSC_VIEWER_STDOUT_SELF);
+
+// AMDiS::finalize();
+// exit(0);
+
if (printTimings) {
MPI::COMM_WORLD.Barrier();
MSG("FETI-DP timing 02: %.5f seconds (creation of interior matrices)\n",
@@ -1488,6 +1558,11 @@ namespace AMDiS {
// === Create matrices for FETI-DP preconditioner. ===
createPreconditionerMatrix(mat);
+
+ // MSG("DUALS MATRIX IS SYM: %d\n", testMatrixSymmetric(mat_duals_duals, true));
+
+// AMDiS::finalize();
+// exit(0);
// === Create and fill PETSc matrix for Lagrange constraints. ===
createMatLagrange(feSpaces);
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index 3313bab362b41aa15be5aba3328639e40c874528..affcf1ec141604fa6dd02834eb9399cb8d5c5aa7 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/AMDiS/src/parallel/PetscSolverFeti.h
@@ -72,11 +72,6 @@ namespace AMDiS {
/// spaces of all components.
void initialize(vector<const FiniteElemSpace*> feSpaces);
- /// After mesh changes, or if the solver is called the first time, this
- /// function creates all information about primal nodes, dual nodes and
- /// lagrange constraints.
- void createFetiData();
-
int getNumberOfPrimals()
{
return primalDofMap.getOverallDofs();
@@ -103,6 +98,14 @@ namespace AMDiS {
}
protected:
+ ///
+ void createDirichletData(Matrix<DOFMatrix*> &mat);
+
+ /// After mesh changes, or if the solver is called the first time, this
+ /// function creates all information about primal nodes, dual nodes and
+ /// lagrange constraints.
+ void createFetiData();
+
/// Defines which boundary nodes are primal. Creates global index of
/// the primal variables.
void createPrimals(const FiniteElemSpace *feSpace);
@@ -167,6 +170,7 @@ namespace AMDiS {
Vec &vec_sol_primal,
SystemVector &vec);
+ ///
void recoverInterfaceSolution(Vec& vecInterface,
SystemVector &vec);
@@ -322,6 +326,8 @@ namespace AMDiS {
const FiniteElemSpace* pressureFeSpace;
int pressureComponent;
+
+ map<const FiniteElemSpace*, std::set<DegreeOfFreedom> > dirichletRows;
};
}
diff --git a/AMDiS/src/parallel/PetscSolverFetiOperators.cc b/AMDiS/src/parallel/PetscSolverFetiOperators.cc
index 582053c82bcbf37bd514bfa845e8480d39d8ebda..a3bae20089123975698c5daecc2eed5c8c7ace01 100644
--- a/AMDiS/src/parallel/PetscSolverFetiOperators.cc
+++ b/AMDiS/src/parallel/PetscSolverFetiOperators.cc
@@ -288,6 +288,7 @@ namespace AMDiS {
PetscErrorCode petscApplyFetiInterfaceLumpedPrecon(PC pc, Vec xvec, Vec yvec)
{
+ FUNCNAME("precon");
// Get data for the preconditioner
void *ctx;
PCShellGetContext(pc, &ctx);
@@ -301,16 +302,20 @@ namespace AMDiS {
VecNestGetSubVec(yvec, 1, &y_lagrange);
#if 0
- MatMult(data->subSolver->getMatCoarse(0, 1), x_interface, data->tmp_primal);
- KSPSolve(data->ksp_primal, data->tmp_primal, data->tmp_primal);
- MatMult(data->subSolver->getMatCoarse(1, 0), data->tmp_primal, y_interface);
+// MatMult(data->subSolver->getMatCoarse(0, 1), x_interface, data->tmp_primal);
+// KSPSolve(data->ksp_primal, data->tmp_primal, data->tmp_primal);
+// MatMult(data->subSolver->getMatCoarse(1, 0), data->tmp_primal, y_interface);
MatMult(data->subSolver->getMatInteriorCoarse(1), x_interface, data->tmp_vec_b0);
MatMultTranspose(*(data->mat_lagrange_scaled), x_lagrange, data->tmp_vec_b1);
VecAXPY(data->tmp_vec_b0, 1.0, data->tmp_vec_b1);
+
+ //MatMultTranspose(*(data->mat_lagrange_scaled), x_lagrange, data->tmp_vec_b0);
#else
VecCopy(x_interface, y_interface);
- VecScale(y_interface, 1024.0);
+ double scaleFactor = 1.0;
+ Parameters::get("scal", scaleFactor);
+ VecScale(y_interface, scaleFactor);
// VecPointwiseMult(y_interface, x_interface, data->h2vec);
MatMultTranspose(*(data->mat_lagrange_scaled), x_lagrange, data->tmp_vec_b0);
@@ -358,10 +363,10 @@ namespace AMDiS {
MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b0, y_lagrange);
#if 0
- Vec tmp_interface;
- VecDuplicate(y_interface, &tmp_interface);
- MatMult(data->subSolver->getMatCoarseInterior(1), data->tmp_vec_b0, tmp_interface);
- VecAXPY(y_interface, 1.0, tmp_interface);
+// Vec tmp_interface;
+// VecDuplicate(y_interface, &tmp_interface);
+ MatMult(data->subSolver->getMatCoarseInterior(1), data->tmp_vec_b0, y_interface);
+ // VecAXPY(y_interface, 1.0, tmp_interface);
#endif
return 0;
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 55e96dacf0309364e68eda993f30ea60b5efcca6..24b2da760256dc4949da42397b5f4b77397da489 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -158,6 +158,10 @@ namespace AMDiS {
bool isColCoarse =
isCoarseSpace(colComponent, feSpaces[colComponent], col(*icursor));
+ if (isColCoarse == false)
+ if ((*interiorMap)[dofMat->getColFeSpace()].isSet(col(*icursor)) == false)
+ continue;
+
if (isColCoarse == isRowCoarse) {
cols.push_back(col(*icursor));
values.push_back(value(*icursor));
@@ -183,13 +187,16 @@ namespace AMDiS {
for (unsigned int i = 0; i < colsOther.size(); i++)
colsOther[i] =
interiorMap->getMatIndex(colComponent, colsOther[i]) +
- rStartInterior;
+ rStartInterior;
MatSetValues(getMatCoarseInteriorByComponent(rowComponent),
1, &rowIndex, colsOther.size(),
&(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
}
} else {
+ if ((*interiorMap)[dofMat->getRowFeSpace()].isSet(*cursor) == false)
+ continue;
+
int localRowIndex =
(subdomainLevel == 0 ?
interiorMap->getLocalMatIndex(rowComponent, *cursor) :
@@ -682,6 +689,9 @@ namespace AMDiS {
int index = rowCoarseSpace->getMatIndex(nRowVec, dofIt.getDOFIndex());
VecSetValue(vecCoarse, index, *dofIt, ADD_VALUES);
} else {
+ if ((*interiorMap)[feSpace].isSet(dofIt.getDOFIndex()) == false)
+ continue;
+
// Calculate global row index of the DOF.
DegreeOfFreedom globalRowDof =
(*interiorMap)[feSpace][dofIt.getDOFIndex()].global;
diff --git a/AMDiS/src/reinit/HL_SignedDistTraverse.cc b/AMDiS/src/reinit/HL_SignedDistTraverse.cc
index a0bb681f81a5cef93a8c3abe210dde7afb7cc661..9359890e1c43e720a55bebb445815ac6874e7188 100644
--- a/AMDiS/src/reinit/HL_SignedDistTraverse.cc
+++ b/AMDiS/src/reinit/HL_SignedDistTraverse.cc
@@ -52,10 +52,9 @@ void HL_SignedDistTraverse::initializeBoundary()
((VelocityExtFromVelocityField *)velExt)->setElInfo(elInfo);
// Get local indices of vertices.
- feSpace->getBasisFcts()->getLocalIndices(
- const_cast<Element*>(elInfo->getElement()),
- const_cast<DOFAdmin*>(feSpace->getAdmin()),
- &(locInd[0]));
+ feSpace->getBasisFcts()->getLocalIndices(const_cast<Element*>(elInfo->getElement()),
+ const_cast<DOFAdmin*>(feSpace->getAdmin()),
+ locInd);
// Get element status.
elStatus = elLS->createElementLevelSet(elInfo);
@@ -161,10 +160,9 @@ void HL_SignedDistTraverse::HL_elementUpdate(ElInfo *elInfo)
if (static_cast<int>(locInd.size()) < nBasFcts)
locInd.resize(nBasFcts);
- feSpace->getBasisFcts()->getLocalIndices(
- const_cast<Element *>(elInfo->getElement()),
- const_cast<DOFAdmin *>(feSpace->getAdmin()),
- &(locInd[0]));
+ feSpace->getBasisFcts()->getLocalIndices(const_cast<Element *>(elInfo->getElement()),
+ const_cast<DOFAdmin *>(feSpace->getAdmin()),
+ locInd);
// ===== Hopf-Lax element update for each vertex of element. =====
for (int i = 0; i <= dim; i++) {