diff --git a/AMDiS/other/include/Makefile_AMDiS.mk b/AMDiS/other/include/Makefile_AMDiS.mk
index dc5315c15e0d5ec8b54d4a740a3f83dc633b316a..5623d5193fc87b96b9f8871cee5629da299cc91a 100644
--- a/AMDiS/other/include/Makefile_AMDiS.mk
+++ b/AMDiS/other/include/Makefile_AMDiS.mk
@@ -97,7 +97,7 @@ endif
# ============================================================================
ifeq ($(strip $(DEBUG)), 0)
- CPPFLAGS += -O2
+ CPPFLAGS += -O3
else
CPPFLAGS += -g -O0
endif
diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 611fd67a5e12b6ff1f3010edb72748f6e886f978..13e5893c8b600057419c8d16ee807a0b3bc5293a 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -252,12 +252,11 @@ namespace AMDiS {
{
FUNCNAME("Assembler::matVecAssemble()");
- Element *el = elInfo->getElement();
double *uhOldLoc = new double[nRow];
- operat->uhOld->getLocalVector(el, uhOldLoc);
+ operat->uhOld->getLocalVector(elInfo, uhOldLoc);
- if (el != lastMatEl) {
+ if (elInfo->getElement() != lastMatEl) {
set_to_zero(elementMatrix);
calculateElementMatrix(elInfo, elementMatrix);
}
diff --git a/AMDiS/src/BasisFunction.h b/AMDiS/src/BasisFunction.h
index 6203bb6c6a4e54db78f7a97916f359e908f39c8b..ce045c73ac2511673cb22ce949b2e41f169b9010 100644
--- a/AMDiS/src/BasisFunction.h
+++ b/AMDiS/src/BasisFunction.h
@@ -275,7 +275,7 @@ namespace AMDiS {
virtual void coarseRestr(DOFVector<WorldVector<double> >*, RCNeighbourList*, int)
{}
- /// Returns local dof indices of the element for the given fe space.
+ /// Returns local dof indices of the element for the given DOF admin.
virtual const DegreeOfFreedom *getLocalIndices(const Element *el,
const DOFAdmin *admin,
DegreeOfFreedom *dofPtr) const
@@ -283,17 +283,11 @@ namespace AMDiS {
return NULL;
}
- inline void getLocalIndices(const Element *el,
- const DOFAdmin *admin,
- std::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]));
- }
+ /// Calculates local dof indices of the element for the given DOF admin.
+ virtual void getLocalIndices(const Element *el,
+ const DOFAdmin *admin,
+ std::vector<DegreeOfFreedom> &indices) const
+ {}
virtual void getLocalDofPtrVec(const Element *el,
const DOFAdmin *admin,
diff --git a/AMDiS/src/BoundaryManager.cc b/AMDiS/src/BoundaryManager.cc
index 6fd1a42631da8d1755fbf40d9caba266624aadd3..57a464a850ee89364397eac4b90ab434834d6818 100644
--- a/AMDiS/src/BoundaryManager.cc
+++ b/AMDiS/src/BoundaryManager.cc
@@ -15,7 +15,6 @@ namespace AMDiS {
BoundaryManager::BoundaryManager(const FiniteElemSpace *feSpace)
{
localBounds.resize(omp_get_overall_max_threads());
- dofIndices.resize(omp_get_overall_max_threads());
allocatedMemoryLocalBounds = feSpace->getBasisFcts()->getNumber();
for (int i = 0; i < static_cast<int>(localBounds.size()); i++)
localBounds[i] = new BoundaryType[allocatedMemoryLocalBounds];
@@ -26,7 +25,6 @@ namespace AMDiS {
localBCs = bm.localBCs;
allocatedMemoryLocalBounds = bm.allocatedMemoryLocalBounds;
localBounds.resize(bm.localBounds.size());
- dofIndices.resize(bm.localBounds.size());
for (int i = 0; i < static_cast<int>(localBounds.size()); i++)
localBounds[i] = new BoundaryType[allocatedMemoryLocalBounds];
}
@@ -49,54 +47,61 @@ namespace AMDiS {
return result;
}
+
void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo,
DOFVectorBase<double> *vec)
{
- if (localBCs.size() > 0) {
- const FiniteElemSpace *feSpace = vec->getFESpace();
- std::vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
- const BasisFunction *basisFcts = feSpace->getBasisFcts();
- int nBasFcts = basisFcts->getNumber();
- dofVec.resize(nBasFcts);
-
- // get boundaries of all DOFs
- BoundaryType *localBound = localBounds[omp_get_thread_num()];
- basisFcts->getBound(elInfo, localBound);
-
- // get dof indices
- basisFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
-
- // apply non dirichlet boundary conditions
- for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second && !(*it).second->isDirichlet())
- (*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0],
- localBound, nBasFcts);
-
- // apply dirichlet boundary conditions
- for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second && (*it).second->isDirichlet())
- (*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0],
- localBound, nBasFcts);
- }
+ FUNCNAME("BoundaryManager::fillBoundaryConditions()");
+
+ if (localBCs.size() <= 0)
+ return;
+
+ const FiniteElemSpace *feSpace = vec->getFESpace();
+ const BasisFunction *basisFcts = feSpace->getBasisFcts();
+ int nBasFcts = basisFcts->getNumber();
+
+ // get boundaries of all DOFs
+ BoundaryType *localBound = localBounds[omp_get_thread_num()];
+ basisFcts->getBound(elInfo, localBound);
+
+ // get dof indices
+ std::vector<DegreeOfFreedom> &dofVec = elInfo->getLocalIndices(feSpace);
+ TEST_EXIT_DBG(static_cast<int>(dofVec.size()) == nBasFcts)
+ ("Local index vector is too small!\n");
+
+ // apply non dirichlet boundary conditions
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && !(*it).second->isDirichlet())
+ (*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0],
+ localBound, nBasFcts);
+
+ // apply dirichlet boundary conditions
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && (*it).second->isDirichlet())
+ (*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0],
+ localBound, nBasFcts);
}
+
void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo, DOFMatrix *mat)
{
+ FUNCNAME("BoundaryManager::fillBoundaryConditions()");
+
if (localBCs.size() <= 0)
return;
const FiniteElemSpace *feSpace = mat->getRowFESpace();
- std::vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
const BasisFunction *basisFcts = feSpace->getBasisFcts();
int nBasFcts = basisFcts->getNumber();
- dofVec.resize(nBasFcts);
// get boundaries of all DOFs
BoundaryType *localBound = localBounds[omp_get_thread_num()];
basisFcts->getBound(elInfo, localBound);
// get dof indices
- basisFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
+ std::vector<DegreeOfFreedom> &dofVec = elInfo->getLocalIndices(feSpace);
+ TEST_EXIT_DBG(static_cast<int>(dofVec.size()) == nBasFcts)
+ ("Local index vector is too small!\n");
// apply non dirichlet boundary conditions
for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
diff --git a/AMDiS/src/BoundaryManager.h b/AMDiS/src/BoundaryManager.h
index c9f391a0c5daebae407c74ae8466920008e24533..98615d230dc885e07d217e48ada0fac9b1086d84 100644
--- a/AMDiS/src/BoundaryManager.h
+++ b/AMDiS/src/BoundaryManager.h
@@ -125,9 +125,6 @@ namespace AMDiS {
/// Temporary thread-safe variable for functions fillBoundaryconditions.
std::vector<BoundaryType*> localBounds;
- /// Temporary thread-safe variable for functions fillBoundaryconditions.
- std::vector<std::vector<DegreeOfFreedom> > dofIndices;
-
/** \brief
* Stores the number of byte that were allocated in the constructor for
* each localBounds value. Is used to free the memory in the destructor.
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index fcc63e7684cc8abba0e20e577a23a93c9c53f84f..163f7aa6e365220ed3f1545a2cb420430b0f840c 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -168,9 +168,12 @@ namespace AMDiS {
// === Get indices mapping from local to global matrix indices. ===
- rowFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
- rowFESpace->getAdmin(),
- rowIndices);
+ std::vector<DegreeOfFreedom> &rowIndices2 = rowElInfo->getLocalIndices(rowFESpace);
+ std::vector<DegreeOfFreedom> &colIndices2 = rowIndices2;
+
+// rowFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
+// rowFESpace->getAdmin(),
+// rowIndices);
if (rowFESpace == colFESpace) {
colIndices = rowIndices;
} else {
@@ -188,21 +191,8 @@ namespace AMDiS {
using namespace mtl;
-#if 0
- std::cout << "----- PRINT MAT --------" << std::endl;
- std::cout << elMat << std::endl;
- std::cout << "rows: ";
- for (int i = 0; i < rowIndices.size(); i++)
- std::cout << rowIndices[i] << " ";
- std::cout << std::endl;
- std::cout << "cols: ";
- for (int i = 0; i < colIndices.size(); i++)
- std::cout << colIndices[i] << " ";
- std::cout << std::endl;
-#endif
-
for (int i = 0; i < nRow; i++) {
- DegreeOfFreedom row = rowIndices[i];
+ DegreeOfFreedom row = rowIndices2[i];
BoundaryCondition *condition =
bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
@@ -210,20 +200,15 @@ namespace AMDiS {
if (condition && condition->isDirichlet()) {
if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- if ((*rankDofs)[rowIndices[i]])
+ if ((*rankDofs)[row])
applyDBCs.insert(static_cast<int>(row));
#else
applyDBCs.insert(static_cast<int>(row));
#endif
}
} else {
- for (int j = 0; j < nCol; j++) {
- DegreeOfFreedom col = colIndices[j];
- double entry = elMat[i][j];
-
- // std::cout << "ADD at " << row << " " << col << std::endl;
- ins[row][col] += entry;
- }
+ for (int j = 0; j < nCol; j++)
+ ins[row][colIndices2[j]] += elMat[i][j];
}
}
}
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index cb978648331bbb4e5ce4ef25a05d365ba584be72..61389406bb220184f44bebbd910f2cad562cd5fa 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -108,13 +108,16 @@ namespace AMDiS {
// traverse mesh
std::vector<bool> visited(getUsedSize(), false);
TraverseStack stack;
- Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
+ stack.addFeSpace(feSpace);
+ Flag fillFlag =
+ Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA |
+ Mesh::FILL_COORDS | Mesh::FILL_LOCAL_INDICES;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
while (elInfo) {
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
- getLocalVector(elInfo->getElement(), localUh);
+ getLocalVector(elInfo, localUh);
int localDOFNr = 0;
for (int i = 0; i < nNodes; i++) { // for all nodes
@@ -177,9 +180,11 @@ namespace AMDiS {
// traverse mesh
Mesh *mesh = feSpace->getMesh();
TraverseStack stack;
+ stack.addFeSpace(feSpace);
ElInfo *elInfo = stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
+ Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS |
+ Mesh::FILL_LOCAL_INDICES);
double *localUh = new double[basFcts->getNumber()];
@@ -187,7 +192,7 @@ namespace AMDiS {
double det = elInfo->getDet();
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
- getLocalVector(elInfo->getElement(), localUh);
+ getLocalVector(elInfo, localUh);
basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
for (int i = 0; i < dim + 1; i++) {
@@ -249,7 +254,7 @@ namespace AMDiS {
}
double *localVec = localVectors[myRank];
- getLocalVector(elInfo->getElement(), localVec);
+ getLocalVector(elInfo, localVec);
DimVec<double> &grd1 = *grdTmp[myRank];
int parts = Global::getGeo(PARTS, dim);
@@ -334,7 +339,7 @@ namespace AMDiS {
}
double *localVec = localVectors[myRank];
- getLocalVector(largeElInfo->getElement(), localVec);
+ getLocalVector(largeElInfo, localVec);
const BasisFunction *basFcts = feSpace->getBasisFcts();
mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
@@ -391,8 +396,6 @@ namespace AMDiS {
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
- Element *el = elInfo->getElement();
-
int myRank = omp_get_thread_num();
int dow = Global::getGeo(WORLD);
int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
@@ -412,7 +415,7 @@ namespace AMDiS {
}
double *localVec = localVectors[myRank];
- getLocalVector(el, localVec);
+ getLocalVector(elInfo, localVec);
DimMat<double> D2Tmp(dim, DEFAULT_VALUE, 0.0);
int parts = Global::getGeo(PARTS, dim);
@@ -776,7 +779,7 @@ namespace AMDiS {
}
double *localVec = localVectors[omp_get_thread_num()];
- getLocalVector(largeElInfo->getElement(), localVec);
+ getLocalVector(largeElInfo, localVec);
mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 2d42e1286eaa00db0b0df5605fea9f5a6bb5333e..7e88f4b95e22f18773f4a23d6e99c1c4c80c8ee8 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -67,7 +67,9 @@ namespace AMDiS {
* For the given element, this function returns an array of all DOFs of this
* DOFVector that are defined on this element.
*/
- virtual const T *getLocalVector(const Element *el, T* localVec) const;
+ const T *getLocalVector(const Element *el, T* localVec) const;
+
+ const T *getLocalVector(const ElInfo *elInfo, T* localVec) const;
/** \brief
* Evaluates the DOF vector at a set of quadrature points defined on the
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index d97dede0e2fb5138877ea9e3ea87db160548121c..b8435cca8dc3fe37a335fef4a3d8d7167713b765 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -147,9 +147,9 @@ namespace AMDiS {
{
FUNCNAME("DOFVector::addElementVector()");
- std::vector<DegreeOfFreedom> indices(nBasFcts);
- feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
- indices);
+ std::vector<DegreeOfFreedom> &indices = elInfo->getLocalIndices(feSpace);
+ TEST_EXIT_DBG(static_cast<int>(indices.size()) == nBasFcts)
+ ("Local index vector is too small!\n");
for (DegreeOfFreedom i = 0; i < nBasFcts; i++) {
BoundaryCondition *condition =
@@ -454,9 +454,11 @@ namespace AMDiS {
int nPoints = quadFast->getNumPoints();
std::vector<T> uh_vec(nPoints);
TraverseStack stack;
+ stack.addFeSpace(this->feSpace);
ElInfo *elInfo =
stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
+ Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
+ Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
while (elInfo) {
double det = elInfo->getDet();
double normT = 0.0;
@@ -490,9 +492,11 @@ namespace AMDiS {
int nPoints = quadFast->getNumPoints();
std::vector<T> uh_vec(nPoints);
TraverseStack stack;
+ stack.addFeSpace(this->feSpace);
ElInfo *elInfo =
stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
+ Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
+ Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
while (elInfo) {
double det = elInfo->getDet();
double normT = 0.0;
@@ -526,9 +530,11 @@ namespace AMDiS {
int nPoints = quadFast->getNumPoints();
std::vector<T> uh_vec(nPoints);
TraverseStack stack;
+ stack.addFeSpace(this->feSpace);
ElInfo *elInfo =
stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
+ Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
+ Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
while (elInfo) {
double det = elInfo->getDet();
double normT = 0.0;
@@ -986,6 +992,49 @@ namespace AMDiS {
}
+ template<typename T>
+ const T *DOFVectorBase<T>::getLocalVector(const ElInfo *elInfo, T *d) const
+ {
+ FUNCNAME("DOFVectorBase<T>::getLocalVector()");
+
+ TEST_EXIT_DBG(feSpace->getMesh() == elInfo->getElement()->getMesh())
+ ("Element is defined on a different mesh than the DOF vector!\n");
+
+ static T* localVec = NULL;
+ static int localVecSize = 0;
+ T *result;
+
+ if (d) {
+ result = d;
+ } else {
+#ifdef _OPENMP
+ ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
+#endif
+ if (localVec && nBasFcts > localVecSize) {
+ delete [] localVec;
+ localVec = new T[nBasFcts];
+ }
+
+ if (!localVec)
+ localVec = new T[nBasFcts];
+
+ localVecSize = nBasFcts;
+ result = localVec;
+ }
+
+ std::vector<DegreeOfFreedom> &localIndices =
+ const_cast<ElInfo*>(elInfo)->getLocalIndices(feSpace);
+
+ TEST_EXIT_DBG(static_cast<int>(localIndices.size()) == nBasFcts)
+ ("Local indices vector is too small!\n");
+
+ for (int i = 0; i < nBasFcts; i++)
+ result[i] = (*this)[localIndices[i]];
+
+ return result;
+ }
+
+
template<typename T>
const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo,
const Quadrature *quad,
@@ -1004,7 +1053,6 @@ namespace AMDiS {
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
- Element *el = elInfo->getElement();
const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
const BasisFunction *basFcts = feSpace->getBasisFcts();
int nPoints = quadrature->getNumPoints();
@@ -1028,7 +1076,7 @@ namespace AMDiS {
}
T *localVec = localVectors[omp_get_thread_num()];
- getLocalVector(el, localVec);
+ getLocalVector(elInfo, localVec);
for (int i = 0; i < nPoints; i++) {
result[i] = 0.0;
@@ -1093,9 +1141,11 @@ namespace AMDiS {
int nPoints = quadFast->getNumPoints();
std::vector<T> uh_vec(nPoints);
TraverseStack stack;
+ stack.addFeSpace(this->feSpace);
ElInfo *elInfo =
stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
+ Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
+ Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
while (elInfo) {
double det = elInfo->getDet();
double normT = 0.0;
diff --git a/AMDiS/src/DualTraverse.h b/AMDiS/src/DualTraverse.h
index 8226414014cc2828246ae0517b84eea2e475b3ca..d5edd66b71d45ad0cd22afef4d95c93570063b23 100644
--- a/AMDiS/src/DualTraverse.h
+++ b/AMDiS/src/DualTraverse.h
@@ -125,6 +125,19 @@ namespace AMDiS {
* \param[in] largeFace A local edge/face number on the large element.
*/
static int getFace(DualElInfo *dualElInfo, int largeFace);
+
+ void addFeSpace(int stack, const FiniteElemSpace *feSpace)
+ {
+ FUNCNAME("DualTraverse::addFeSpace()");
+
+ TEST_EXIT_DBG(stack == 0 || stack == 1)("Wrong stack number!\n");
+
+ if (stack == 0)
+ stack1.addFeSpace(feSpace);
+
+ if (stack == 1)
+ stack2.addFeSpace(feSpace);
+ }
protected:
/** \brief
diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h
index bd0f395530d0fc945357db966860a02568d4afc3..6e0cc6f61e37f965ee25ef9fe97110e0511ff662 100644
--- a/AMDiS/src/ElInfo.h
+++ b/AMDiS/src/ElInfo.h
@@ -227,6 +227,17 @@ namespace AMDiS {
return parametric;
}
+ /// Returns the local indices, \ref localIndices, for given FE space.
+ inline std::vector<DegreeOfFreedom>& getLocalIndices(const FiniteElemSpace* feSpace)
+ {
+ FUNCNAME("ElInfo::getLocalIndices()");
+
+ TEST_EXIT_DBG(element->isLeaf() == true)
+ ("Local indices are computed only for leaf elements!\n");
+
+ return localIndices[feSpace];
+ }
+
/// Returns \ref refinementPath
inline unsigned long getRefinementPath() const
{
@@ -393,7 +404,7 @@ namespace AMDiS {
* parent data parentInfo.
* pure virtual => must be overriden in sub-class.
*/
- virtual void fillElInfo(int ichild, const ElInfo *parentInfo) = 0;
+ virtual void fillElInfo(int iChild, const ElInfo *parentInfo) = 0;
/** \brief
* calculates the Jacobian of the barycentric coordinates on \element and stores
@@ -525,6 +536,13 @@ namespace AMDiS {
/// Gradient of lambda.
DimVec<WorldVector<double> > grdLambda;
+ /** \brief
+ * If the local indices of elements should be set during mesh traverse, they are
+ * stored in this variable for all differnt FE spaces that are defined on the
+ * mesh. Node that local indices are computed only for leaf elements of the mesh.
+ */
+ std::map<const FiniteElemSpace*, std::vector<DegreeOfFreedom> > localIndices;
+
/// True, if this elInfo stores parametrized information. False, otherwise.
bool parametric;
diff --git a/AMDiS/src/Error.hh b/AMDiS/src/Error.hh
index 9dd255b15c7e8851acc1a8e3f0ea0d4cffbcc42c..11df84a7f4dd5e06ab9ba74fe4c9150d2fe6daec 100644
--- a/AMDiS/src/Error.hh
+++ b/AMDiS/src/Error.hh
@@ -28,29 +28,33 @@ namespace AMDiS {
{
FUNCNAME("Error<T>::maxErrAtQp()");
- const FiniteElemSpace *fe_space;
+ const FiniteElemSpace *feSpace = uh->getFESpace();
+
if (!(pU = &u)) {
ERROR("no function u specified; doing nothing\n");
return(-1.0);
}
- if (!(errUh = &uh) || !(fe_space = uh->getFESpace())) {
- ERROR("no discrete function or no fe_space for it; doing nothing\n");
+ if (!(errUh = &uh) || !feSpace) {
+ ERROR("no discrete function or no feSpace for it; doing nothing\n");
return(-1.0);
}
- if (!(basFct = fe_space->getBasisFcts())) {
+ if (!(basFct = feSpace->getBasisFcts())) {
ERROR("no basis functions at discrete solution ; doing nothing\n");
return(-1.0);
}
if (!q)
- q = Quadrature::provideQuadrature(fe_space->getMesh()->getDim(),
- 2 * fe_space->getBasisFcts()->getDegree() - 2);
+ q = Quadrature::provideQuadrature(feSpace->getMesh()->getDim(),
+ 2 * feSpace->getBasisFcts()->getDegree() - 2);
quadFast = FastQuadrature::provideFastQuadrature(basFct, *q, INIT_PHI);
double maxErr = 0.0;
TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
- Mesh::FILL_COORDS | Mesh::CALL_LEAF_EL);
+ stack.addFeSpace(feSpace);
+ ElInfo *elInfo =
+ stack.traverseFirst(feSpace->getMesh(), -1,
+ Mesh::FILL_COORDS | Mesh::FILL_LOCAL_INDICES |
+ Mesh::CALL_LEAF_EL);
while (elInfo) {
elinfo = elInfo;
double err = 0.0;
@@ -81,25 +85,25 @@ namespace AMDiS {
{
FUNCNAME("Error<T>::H1Err()");
- const FiniteElemSpace *fe_space;
+ const FiniteElemSpace *feSpace;
writeInLeafData = writeLeafData;
component = comp;
Quadrature *q = NULL;
pGrdU = &grdU;
errUh = &uh;
- if (!(fe_space = uh.getFESpace())) {
- ERROR("no fe_space for uh; doing nothing\n");
+ if (!(feSpace = uh.getFESpace())) {
+ ERROR("no feSpace for uh; doing nothing\n");
return(0.0);
}
- if (!(basFct = fe_space->getBasisFcts())) {
+ if (!(basFct = feSpace->getBasisFcts())) {
ERROR("no basis functions at discrete solution ; doing nothing\n");
return(0.0);
}
- int dim = fe_space->getMesh()->getDim();
+ int dim = feSpace->getMesh()->getDim();
int deg = grdU.getDegree();
- int degree = deg ? deg : 2 * fe_space->getBasisFcts()->getDegree() - 2;
+ int degree = deg ? deg : 2 * feSpace->getBasisFcts()->getDegree() - 2;
q = Quadrature::provideQuadrature(dim, degree);
quadFast = FastQuadrature::provideFastQuadrature(basFct,
@@ -111,7 +115,7 @@ namespace AMDiS {
TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
+ ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
Mesh::FILL_COORDS |
Mesh::CALL_LEAF_EL |
Mesh::FILL_DET |
@@ -159,7 +163,7 @@ namespace AMDiS {
if (relative) {
double relNorm2 = h1Norm2 + 1.e-15;
- elInfo = stack.traverseFirst(fe_space->getMesh(), -1, Mesh::CALL_LEAF_EL);
+ elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
double exact = elInfo->getElement()->getEstimation(component) / relNorm2;
if (writeInLeafData)
@@ -187,7 +191,7 @@ namespace AMDiS {
{
FUNCNAME("Error<T>::L2Err()");
- const FiniteElemSpace *fe_space;
+ const FiniteElemSpace *feSpace;
Quadrature *q = NULL;
writeInLeafData = writeLeafData;
component = comp;
@@ -197,19 +201,19 @@ namespace AMDiS {
return(0.0);
}
- if (!(errUh = &uh) || !(fe_space = uh.getFESpace())) {
- ERROR("no discrete function or no fe_space for it; doing nothing\n");
+ if (!(errUh = &uh) || !(feSpace = uh.getFESpace())) {
+ ERROR("no discrete function or no feSpace for it; doing nothing\n");
return(0.0);
}
- if (!(basFct = fe_space->getBasisFcts())) {
+ if (!(basFct = feSpace->getBasisFcts())) {
ERROR("no basis functions at discrete solution ; doing nothing\n");
return(0.0);
}
- int dim = fe_space->getMesh()->getDim();
+ int dim = feSpace->getMesh()->getDim();
int deg = u.getDegree();
- int degree = deg ? deg : 2 * fe_space->getBasisFcts()->getDegree() - 2;
+ int degree = deg ? deg : 2 * feSpace->getBasisFcts()->getDegree() - 2;
q = Quadrature::provideQuadrature(dim, degree);
quadFast = FastQuadrature::provideFastQuadrature(basFct, *q, INIT_PHI);
@@ -219,11 +223,12 @@ namespace AMDiS {
double *u_vec = new double[quadFast->getQuadrature()->getNumPoints()];
TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
+ ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
Mesh::FILL_COORDS |
Mesh::CALL_LEAF_EL |
Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA);
+ Mesh::FILL_GRD_LAMBDA |
+ Mesh::FILL_LOCAL_INDICES);
while (elInfo) {
elinfo = elInfo;
const double *uh_vec;
@@ -261,7 +266,7 @@ namespace AMDiS {
if (relative) {
double relNorm2 = l2Norm2 + 1.e-15;
- elInfo = stack.traverseFirst(fe_space->getMesh(), -1, Mesh::CALL_LEAF_EL);
+ elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
double exact = elInfo->getElement()->getEstimation(component) / relNorm2;
if (writeInLeafData)
diff --git a/AMDiS/src/Estimator.cc b/AMDiS/src/Estimator.cc
index ba4ca4bca896ee53816c8a632c6d2453376687ed..c35f2f44c4a551b2aaf548c2e6aa64d3d0ad6b5c 100644
--- a/AMDiS/src/Estimator.cc
+++ b/AMDiS/src/Estimator.cc
@@ -76,6 +76,11 @@ namespace AMDiS {
FUNCNAME("Estimator::singleMeshTraverse()");
TraverseStack stack;
+
+ if (traverseFlag.isSet(Mesh::FILL_LOCAL_INDICES))
+ for (unsigned int i = 0; i < matrix.size(); i++)
+ stack.addFeSpace(matrix[i]->getFESpace());
+
ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
while (elInfo) {
estimateElement(elInfo);
@@ -91,6 +96,16 @@ namespace AMDiS {
DualTraverse dualTraverse;
DualElInfo dualElInfo;
+ if (traverseFlag.isSet(Mesh::FILL_LOCAL_INDICES)) {
+ for (unsigned int i = 0; i < matrix.size(); i++) {
+ if (matrix[i]->getFESpace()->getMesh() == mesh)
+ dualTraverse.addFeSpace(0, matrix[i]->getFESpace());
+
+ if (matrix[i]->getFESpace()->getMesh() == auxMesh)
+ dualTraverse.addFeSpace(1, matrix[i]->getFESpace());
+ }
+ }
+
bool cont = dualTraverse.traverseFirst(mesh, auxMesh, -1, -1,
traverseFlag, traverseFlag,
dualElInfo);
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index 9d10f8cc07aa738d39d08b49d196cfef15ca49aa..8c50657a96f4e14f6223884ad211cbd087ec11c2 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -692,6 +692,7 @@ namespace AMDiS {
return getLocalIndices(el, &admin, idof);
}
+
int* Lagrange::orderOfPositionIndices(const Element* el,
GeoIndex position,
int positionIndex) const
@@ -717,11 +718,10 @@ namespace AMDiS {
int vertex[3];
int** dof = const_cast<int**>(el->getDOF());
- int verticesOfPosition = dimOfPosition + 1;
+ const Element *refElement = Global::getReferenceElement(dim);
- for (int i = 0; i < verticesOfPosition; i++)
- vertex[i] = Global::getReferenceElement(dim)->
- getVertexOfPosition(position, positionIndex, i);
+ for (int i = 0; i <= dimOfPosition; i++)
+ vertex[i] = refElement->getVertexOfPosition(position, positionIndex, i);
if (dimOfPosition == 1) {
if (degree == 3) {
@@ -769,6 +769,7 @@ namespace AMDiS {
return NULL;
}
+
void Lagrange::getBound(const ElInfo* el_info, BoundaryType* bound) const
{
el_info->testFlag(Mesh::FILL_BOUND);
@@ -953,6 +954,38 @@ namespace AMDiS {
return result;
}
+
+ void Lagrange::getLocalIndices(const Element* el,
+ const DOFAdmin *admin,
+ std::vector<DegreeOfFreedom> &indices) const
+ {
+ indices.resize(nBasFcts);
+
+ const DegreeOfFreedom **dof = el->getDOF();
+
+ int nrDOFs, n0, node0, num = 0;
+ GeoIndex posIndex;
+
+ for (int pos = 0, j = 0; pos <= dim; pos++) {
+ posIndex = INDEX_OF_DIM(pos, dim);
+ nrDOFs = admin->getNumberOfDOFs(posIndex);
+
+ if (nrDOFs) {
+ n0 = admin->getNumberOfPreDOFs(posIndex);
+ node0 = admin->getMesh()->getNode(posIndex);
+ num = Global::getGeo(posIndex, dim);
+
+ for (int i = 0; i < num; node0++, i++) {
+ const int *indi = orderOfPositionIndices(el, posIndex, i);
+
+ for (int k = 0; k < nrDOFs; k++)
+ indices[j++] = dof[node0][n0 + indi[k]];
+ }
+ }
+ }
+ }
+
+
void Lagrange::getLocalDofPtrVec(const Element *el,
const DOFAdmin *admin,
std::vector<const DegreeOfFreedom*>& vec) const
diff --git a/AMDiS/src/Lagrange.h b/AMDiS/src/Lagrange.h
index f4a7bddbbfab202a9e03f2d756f1b040b7ceced4..4cb390c4962b05c611b2da65f91d836824c047c8 100644
--- a/AMDiS/src/Lagrange.h
+++ b/AMDiS/src/Lagrange.h
@@ -123,6 +123,11 @@ namespace AMDiS {
const DOFAdmin *admin,
DegreeOfFreedom *dofs) const;
+ /// Implements BasisFunction::getLocalIndices().
+ void getLocalIndices(const Element *el,
+ const DOFAdmin *admin,
+ std::vector<DegreeOfFreedom> &indices) const;
+
void getLocalDofPtrVec(const Element *el,
const DOFAdmin *admin,
std::vector<const DegreeOfFreedom*>& vec) const;
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index 5b03933e94ba60e7067fc0d6dda242e59588a24e..0bae77b69f307bffca55420d2f9bfa8ec552c2dd 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -33,20 +33,21 @@ namespace AMDiS {
// flags, which information should be present in the elInfo structure
//**************************************************************************
- const Flag Mesh::FILL_NOTHING = 0X00L;
- const Flag Mesh::FILL_COORDS = 0X01L;
- const Flag Mesh::FILL_BOUND = 0X02L;
- const Flag Mesh::FILL_NEIGH = 0X04L;
- const Flag Mesh::FILL_OPP_COORDS = 0X08L;
- const Flag Mesh::FILL_ORIENTATION= 0X10L;
- const Flag Mesh::FILL_DET = 0X20L;
- const Flag Mesh::FILL_GRD_LAMBDA = 0X40L;
- const Flag Mesh::FILL_ADD_ALL = 0X80L;
-
-
- const Flag Mesh::FILL_ANY_1D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L);
- const Flag Mesh::FILL_ANY_2D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L);
- const Flag Mesh::FILL_ANY_3D = (0X01L|0X02L|0X04L|0X08L|0X10L|0x20L|0X40L|0X80L);
+ const Flag Mesh::FILL_NOTHING = 0X000L;
+ const Flag Mesh::FILL_COORDS = 0X001L;
+ const Flag Mesh::FILL_BOUND = 0X002L;
+ const Flag Mesh::FILL_NEIGH = 0X004L;
+ const Flag Mesh::FILL_OPP_COORDS = 0X008L;
+ const Flag Mesh::FILL_ORIENTATION = 0X010L;
+ const Flag Mesh::FILL_DET = 0X020L;
+ const Flag Mesh::FILL_GRD_LAMBDA = 0X040L;
+ const Flag Mesh::FILL_LOCAL_INDICES = 0X080L;
+ const Flag Mesh::FILL_ADD_ALL = 0X100L;
+
+
+ const Flag Mesh::FILL_ANY_1D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X100L);
+ const Flag Mesh::FILL_ANY_2D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X100L);
+ const Flag Mesh::FILL_ANY_3D = (0X01L|0X02L|0X04L|0X08L|0X10L|0x20L|0X40L|0X100L);
//**************************************************************************
// flags for Mesh traversal
@@ -60,9 +61,6 @@ namespace AMDiS {
const Flag Mesh::CALL_EL_LEVEL = 0X2000L;
const Flag Mesh::CALL_MG_LEVEL = 0X4000L ; // used in mg methods
-
- // const Flag Mesh::USE_PARAMETRIC = 0X8000L ; // used in mg methods
-
std::vector<DegreeOfFreedom> Mesh::dof_used;
const int Mesh::MAX_DOF = 100;
std::map<std::pair<DegreeOfFreedom, int>, DegreeOfFreedom*> Mesh::serializedDOFs;
diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h
index 57369844df89a6bfc464aa3eecc2354079251768..87b62a7cd217a92b94e64b5beb523eeae47bccd9 100644
--- a/AMDiS/src/Mesh.h
+++ b/AMDiS/src/Mesh.h
@@ -616,6 +616,9 @@ namespace AMDiS {
///
static const Flag FILL_GRD_LAMBDA;
+ ///
+ static const Flag FILL_LOCAL_INDICES;
+
//**************************************************************************
// flags for Mesh traversal
//**************************************************************************
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 827d6c5c93eb3820f80aa40ee86db339eb8209d4..68e7af285b8a6a36b660e8aada86ebc491de205b 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -564,7 +564,7 @@ namespace AMDiS {
// here is reached already because of time adaption
allowFirstRefinement();
- TEST_EXIT_DBG(static_cast<unsigned int>(nComponents == marker.size()))
+ TEST_EXIT_DBG(static_cast<unsigned int>(nComponents) == marker.size())
("Wrong number of markers!\n");
Flag markFlag = 0;
@@ -655,10 +655,11 @@ namespace AMDiS {
flag |
(*systemMatrix)[0][0]->getAssembleFlag() |
rhs->getDOFVector(0)->getAssembleFlag() |
- Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA |
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET |
+ Mesh::FILL_LOCAL_INDICES |
+ Mesh::FILL_GRD_LAMBDA |
Mesh::FILL_NEIGH;
if (useGetBound)
@@ -807,10 +808,11 @@ namespace AMDiS {
flag |
(*systemMatrix)[0][0]->getAssembleFlag() |
rhs->getDOFVector(0)->getAssembleFlag() |
- Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA |
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET |
+ Mesh::FILL_GRD_LAMBDA |
+ Mesh::FILL_LOCAL_INDICES |
Mesh::FILL_NEIGH;
if (useGetBound)
@@ -894,6 +896,8 @@ namespace AMDiS {
useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;
DualTraverse dualTraverse;
+ dualTraverse.addFeSpace(0, feSpaces[0]);
+ dualTraverse.addFeSpace(1, feSpaces[1]);
DualElInfo dualElInfo;
int oldElIndex0 = -1;
int oldElIndex1 = -1;
@@ -1264,9 +1268,12 @@ namespace AMDiS {
#ifdef _OPENMP
TraverseParallelStack stack(0, 1);
#else
- TraverseStack stack;
+ TraverseStack stack;
#endif
+ for (unsigned int i = 0; i < feSpaces.size(); i++)
+ stack.addFeSpace(feSpaces[i]);
+
// == Initialize matrix and vector. If we have to assemble in parallel, ===
// == temporary thread owned matrix and vector must be created. ===
@@ -1423,6 +1430,8 @@ namespace AMDiS {
// === Parallel boundary assemblage ===
TraverseParallelStack stack;
+ for (unsigned int i = 0; i < feSpaces.size(); i++)
+ stack.addFeSpace(feSpaces[i]);
#pragma omp parallel
{
@@ -1469,6 +1478,9 @@ namespace AMDiS {
// === Sequentiell boundary assemblage ===
TraverseStack stack;
+ for (unsigned int i = 0; i < feSpaces.size(); i++)
+ stack.addFeSpace(feSpaces[i]);
+
ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
while (elInfo) {
if (rhs->getBoundaryManager())
diff --git a/AMDiS/src/QPInfo.h b/AMDiS/src/QPInfo.h
index dac5de3600b3d18014f70a3c62912ea1651c8bf2..a718294553d1e67f18cbabb62ec88f1b124f5f04 100644
--- a/AMDiS/src/QPInfo.h
+++ b/AMDiS/src/QPInfo.h
@@ -34,14 +34,10 @@ namespace AMDiS {
class QPInfo
{
public:
- /** \brief
- * Sets \ref currentElInfo_ to elInfo and all valid flags to false.
- */
+ /// Sets \ref currentElInfo_ to elInfo and all valid flags to false.
void initElement(const ElInfo *elInfo);
- /** \brief
- * Returns coordinates at quadrature points.
- */
+ /// Returns coordinates at quadrature points.
WorldVector<double> *getCoordsAtQPs(int numPoints);
/** \brief
@@ -68,53 +64,33 @@ namespace AMDiS {
int numPoints,
const FastQuadrature *quadFast = NULL);
-
- /** \brief
- * Returns element normals at quadrature points.
- */
+ /// Returns element normals at quadrature points.
WorldVector<double> **getElementNormalAtQPs(int numPoints);
- /** \brief
- *
- */
+ ///
DimVec<WorldVector<double> > **getGrdLambdaAtQPs(int numPoints);
- /** \brief
- * Returns a QPInfo instance for the given quadrature.
- */
+ /// Returns a QPInfo instance for the given quadrature.
static QPInfo *provideQPInfo(const Quadrature*, const FastQuadrature*);
- /** \brief
- * Deletes the QPInfo instance for the given quadrature.
- */
+ /// Deletes the QPInfo instance for the given quadrature.
static void clearQPInfo(const Quadrature*, const FastQuadrature*);
- /** \brief
- * Deletes all QPInfo instances.
- */
+ /// Deletes all QPInfo instances.
static void clearAllQPInfos();
protected:
- /** \brief
- * Constructor. Called by \ref provideQPInfo().
- */
+ /// Constructor. Called by \ref provideQPInfo().
QPInfo(const Quadrature*);
- /** \brief
- * Destructor. Called by \ref clearQPInfo() and \ref clearAllQPInfos().
- */
+ /// Destructor. Called by \ref clearQPInfo() and \ref clearAllQPInfos().
~QPInfo();
protected:
- /** \brief
- * Structure which stores infos about one DOFVector.
- */
+ /// Structure which stores infos about one DOFVector.
class VecQPInfo
{
public:
- /** \brief
- * Constructor.
- */
VecQPInfo()
: valAtQPs_(NULL),
grdAtQPs_(NULL),
@@ -124,100 +100,62 @@ namespace AMDiS {
D2NumPointsValid_(0)
{}
- /** \brief
- * Values at quadrature points
- */
+ /// Values at quadrature points
double *valAtQPs_;
- /** \brief
- * Gradients at quadrature points
- */
+ /// Gradients at quadrature points
WorldVector<double> *grdAtQPs_;
- /** \brief
- * D2 at quadrature points
- */
+ /// D2 at quadrature points
WorldMatrix<double> *D2AtQPs_;
- /** \brief
- * valid flag for values
- */
+ /// Valid flag for values
int valNumPointsValid_;
- /** \brief
- * valid flag for gradients
- */
+ /// Valid flag for gradients
int grdNumPointsValid_;
- /** \brief
- * valid flag for D2
- */
+ /// Valid flag for D2
bool D2NumPointsValid_;
};
- /** \brief
- * Quadrature of this QPInfo
- */
+ /// Quadrature of this QPInfo
const Quadrature *quadrature_;
- /** \brief
- * Set to \ref quadrature_->getNumPoints().
- */
+ /// Set to \ref quadrature_->getNumPoints().
int numPoints_;
- /** \brief
- * ElInfo of the current element
- */
+ /// ElInfo of the current element
const ElInfo *currentElInfo_;
- /** \brief
- * Coords at quadrature points
- */
+ /// Coords at quadrature points
WorldVector<double> *coordsAtQPs_;
- /** \brief
- * Valid flag for coords
- */
+ /// Valid flag for coords
int coordsNumPointsValid_;
- /** \brief
- * Map of all vector infos
- */
+ /// Map of all vector infos
std::map<const DOFVector<double>*, VecQPInfo*> vecQPInfos_;
- /** \brief
- * element normal at quadrature points (array of pointers)
- */
+ /// element normal at quadrature points (array of pointers)
WorldVector<double> **elementNormalAtQPs_;
- /** \brief
- * for constant values at all QPs (all entries point to same memory)
- */
+ /// for constant values at all QPs (all entries point to same memory)
WorldVector<double> **elementNormalConst_;
- /** \brief
- * valid flag for element normals
- */
+ /// valid flag for element normals
int elementNormalNumPointsValid_;
- /** \brief
- * gradient of barycentric coordinates at QPs (array of pointers)
- */
+ /// gradient of barycentric coordinates at QPs (array of pointers)
DimVec<WorldVector<double> > **grdLambdaAtQPs_;
- /** \brief
- * for constant values at all QPs (all entries point to same memory)
- */
+ /// for constant values at all QPs (all entries point to same memory)
DimVec<WorldVector<double> > **grdLambdaConst_;
- /** \brief
- * number of valid points of grdLambdaAtQPs_
- */
+ /// number of valid points of grdLambdaAtQPs_
int grdLambdaNumPointsValid_;
- /** \brief
- * Static map of all QPInfos. Used by \ref provideQPInfo().
- */
+ /// Static map of all QPInfos. Used by \ref provideQPInfo().
static std::map<const Quadrature*, QPInfo*> qpInfos_;
};
diff --git a/AMDiS/src/Recovery.cc b/AMDiS/src/Recovery.cc
index f9a669033456438b51e0a4c0d4cb87dca2c3062c..8e05184cf65e20dec3b07ef95db01c3fb2b95265 100644
--- a/AMDiS/src/Recovery.cc
+++ b/AMDiS/src/Recovery.cc
@@ -501,21 +501,23 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
Mesh::FILL_COORDS |
Mesh::FILL_BOUND |
Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA;
+ Mesh::FILL_GRD_LAMBDA |
+ Mesh::FILL_LOCAL_INDICES;
if (degree == 2 && dim > 1)
fill_flag |= Mesh::FILL_NEIGH | Mesh::FILL_OPP_COORDS;
TraverseStack stack;
- ElInfo *el_info = stack.traverseFirst(mesh, -1, fill_flag);
+ stack.addFeSpace(feSpace);
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, fill_flag);
- while (el_info) { // traversing the mesh.
- const DegreeOfFreedom **dof = el_info->getElement()->getDOF();
+ while (elInfo) { // traversing the mesh.
+ const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
int n_neighbors = 0; // counting interior vertices of element
for (int i = 0; i < n_vertices; i++) {
DegreeOfFreedom k = dof[i][preDOFs[VERTEX]];
- if (el_info->getBoundary(VERTEX, i) == INTERIOR)
+ if (elInfo->getBoundary(VERTEX, i) == INTERIOR)
interior_vertices[n_neighbors++] = k;
}
@@ -528,10 +530,10 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
// Setting world coordinates of node.
if (!(*struct_vec)[k].coords) {
(*struct_vec)[k].coords = new WorldVector<double>;
- *(*struct_vec)[k].coords = el_info->getCoord(i);
+ *(*struct_vec)[k].coords = elInfo->getCoord(i);
}
- if (el_info->getBoundary(VERTEX, i) == INTERIOR) {
+ if (elInfo->getBoundary(VERTEX, i) == INTERIOR) {
// Allocating memory for matrix and right hand side.
if (!(*struct_vec)[k].A) {
(*struct_vec)[k].A = new Matrix<double>(n_monomials, n_monomials);
@@ -549,20 +551,20 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
// Computing the integrals.
if (method)
- compute_integrals(uh, el_info, &(*struct_vec)[k],
+ compute_integrals(uh, elInfo, &(*struct_vec)[k],
f_vec, f_scal, aux_vec);
else if (gradient)
- compute_interior_sums(uh, el_info, &(*struct_vec)[k], quad,
+ compute_interior_sums(uh, elInfo, &(*struct_vec)[k], quad,
f_vec, f_scal, aux_vec);
else {
if (!(*struct_vec)[k].neighbors)
(*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
if (degree == 2 && dim > 1)
- compute_sums_linear(uh, el_info, &(*struct_vec)[k],
+ compute_sums_linear(uh, elInfo, &(*struct_vec)[k],
i, pre_dofs, n_vertices);
else
- compute_node_sums(uh, el_info, &(*struct_vec)[k],
+ compute_node_sums(uh, elInfo, &(*struct_vec)[k],
pre_dofs, n_vertices, n_edges, n_faces);
}
} else { // Setting list of adjacent interior vertices.
@@ -583,7 +585,7 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
if (!(*struct_vec)[k].coords) {
// Setting world coordinates of node.
- el_info->coordToWorld(*basis_fcts->getCoords(n_count),
+ elInfo->coordToWorld(*basis_fcts->getCoords(n_count),
coordinates);
(*struct_vec)[k].coords = new WorldVector<double>;
*(*struct_vec)[k].coords = coordinates;
@@ -591,10 +593,10 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
// Setting list of adjacent interior vertices.
(*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
- if (el_info->getBoundary(EDGE, i) == INTERIOR)
+ if (elInfo->getBoundary(EDGE, i) == INTERIOR)
for (int m = 0; m < 2; m++) {
int l = Global::getReferenceElement(dim)->getVertexOfEdge(i, m);
- if (el_info->getBoundary(VERTEX, l) == INTERIOR)
+ if (elInfo->getBoundary(VERTEX, l) == INTERIOR)
(*struct_vec)[k].neighbors->insert(dof[l][preDOFs[VERTEX]]);
} else
for (int m = 0; m < n_neighbors; m++)
@@ -611,7 +613,7 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
if (!(*struct_vec)[k].coords) {
// Setting world coordinates of node.
- el_info->coordToWorld(*basis_fcts->getCoords(n_count),
+ elInfo->coordToWorld(*basis_fcts->getCoords(n_count),
coordinates);
(*struct_vec)[k].coords = new WorldVector<double>;
*(*struct_vec)[k].coords = coordinates;
@@ -619,10 +621,10 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
// Setting list of adjacent interior vertices.
(*struct_vec)[k].neighbors = new std::set<DegreeOfFreedom>;
- if (el_info->getBoundary(FACE, i) == INTERIOR)
+ if (elInfo->getBoundary(FACE, i) == INTERIOR)
for (int m = 0; m < 3; m++) {
int l = Global::getReferenceElement(dim)->getVertexOfPosition(FACE, i, m);
- if (el_info->getBoundary(VERTEX, l) == INTERIOR)
+ if (elInfo->getBoundary(VERTEX, l) == INTERIOR)
(*struct_vec)[k].neighbors->insert(dof[l][preDOFs[VERTEX]]);
}
else
@@ -638,7 +640,7 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
DegreeOfFreedom k = dof[n_vertices+n_edges+n_faces][preDOFs[CENTER] + j];
// Setting world coordinates of node.
- el_info->coordToWorld(*basis_fcts->getCoords(n_count), coordinates);
+ elInfo->coordToWorld(*basis_fcts->getCoords(n_count), coordinates);
(*struct_vec)[k].coords = new WorldVector<double>;
*(*struct_vec)[k].coords = coordinates;
@@ -651,7 +653,7 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
n_count++;
}
- el_info = stack.traverseNext(el_info);
+ elInfo = stack.traverseNext(elInfo);
}
}
@@ -884,8 +886,11 @@ Recovery::recovery(DOFVector<double> *uh,
// traverse mesh
TraverseStack stack;
- Flag fillFlag =
- Mesh::CALL_LEAF_EL | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
+ stack.addFeSpace(feSpace);
+ Flag fillFlag =
+ Mesh::CALL_LEAF_EL | Mesh::FILL_DET |
+ Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS |
+ Mesh::FILL_LOCAL_INDICES;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
while (elInfo) {
diff --git a/AMDiS/src/RecoveryEstimator.cc b/AMDiS/src/RecoveryEstimator.cc
index 7eac009fda4e268d96f9ce58cd122f7499899a2e..4c28a5873164c5c47ef67f8047c41a311fad8ed7 100644
--- a/AMDiS/src/RecoveryEstimator.cc
+++ b/AMDiS/src/RecoveryEstimator.cc
@@ -94,6 +94,7 @@ namespace AMDiS {
// clear error indicators
TraverseStack stack;
+ stack.addFeSpace(uh_->getFESpace());
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
elInfo->getElement()->setEstimation(0.0, row);
@@ -105,7 +106,8 @@ namespace AMDiS {
elInfo = stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
- Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
+ Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA |
+ Mesh::FILL_LOCAL_INDICES);
while (elInfo) {
Element *el = elInfo->getElement();
double det = elInfo->getDet();
diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/ResidualEstimator.cc
index 271a6ae980475956ab84be329035aa4fea36b73a..387ee37bc0d73d66c0faecccff6b8ff868581871 100644
--- a/AMDiS/src/ResidualEstimator.cc
+++ b/AMDiS/src/ResidualEstimator.cc
@@ -102,14 +102,15 @@ namespace AMDiS {
est_max = 0.0;
est_t_sum = 0.0;
est_t_max = 0.0;
-
+
traverseFlag =
- Mesh::FILL_NEIGH |
- Mesh::FILL_COORDS |
- Mesh::FILL_OPP_COORDS |
- Mesh::FILL_BOUND |
- Mesh::FILL_GRD_LAMBDA |
- Mesh::FILL_DET |
+ Mesh::FILL_NEIGH |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_OPP_COORDS |
+ Mesh::FILL_BOUND |
+ Mesh::FILL_GRD_LAMBDA |
+ Mesh::FILL_DET |
+ Mesh::FILL_LOCAL_INDICES |
Mesh::CALL_LEAF_EL;
neighInfo = mesh->createNewElInfo();
@@ -286,7 +287,7 @@ namespace AMDiS {
if (timestep && uhOld[system]) {
TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
- uhOld[system]->getLocalVector(elInfo->getElement(), uhOldEl[system]);
+ uhOld[system]->getLocalVector(elInfo, uhOldEl[system]);
// === Compute time error. ===
@@ -445,7 +446,7 @@ namespace AMDiS {
if (matrix[system] == NULL || secondOrderTerms[system] == false)
continue;
- uh[system]->getLocalVector(el, uhEl[system]);
+ uh[system]->getLocalVector(elInfo, uhEl[system]);
uh[system]->getLocalVector(neigh, uhNeigh[system]);
for (int iq = 0; iq < nPointsSurface; iq++) {
diff --git a/AMDiS/src/RobinBC.cc b/AMDiS/src/RobinBC.cc
index 9df83ac1838735235480daf5b74b8970979b2f93..6a4ec41edfea7a7b7c946a7e0f9060ea9be9e567 100644
--- a/AMDiS/src/RobinBC.cc
+++ b/AMDiS/src/RobinBC.cc
@@ -22,10 +22,10 @@ namespace AMDiS {
// create barycentric coords for each vertex of each side
const Element *refElement = Global::getReferenceElement(dim);
- coords = new VectorOfFixVecs<DimVec<double> >*[dim+1];
+ coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
// for all element sides
- for (int i = 0; i < dim+1; i++) {
+ for (int i = 0; i < dim + 1; i++) {
coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
DimVec<double>(dim,
DEFAULT_VALUE,
@@ -74,7 +74,7 @@ namespace AMDiS {
// create barycentric coords for each vertex of each side
const Element *refElement = Global::getReferenceElement(dim);
- coords = new VectorOfFixVecs<DimVec<double> >*[dim+1];
+ coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
// for all element sides
for (int i = 0; i < dim + 1; i++) {
@@ -95,7 +95,7 @@ namespace AMDiS {
jOp->addZeroOrderTerm(new VecAtQP_ZOT(j, NULL));
neumannOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
- for (int i=0; i < dim+1; i++)
+ for (int i=0; i < dim + 1; i++)
(*neumannOperators)[i] = new SurfaceOperator(jOp, *coords[i]);
delete jOp;
@@ -126,7 +126,7 @@ namespace AMDiS {
// create barycentric coords for each vertex of each side
const Element *refElement = Global::getReferenceElement(dim);
- coords = new VectorOfFixVecs<DimVec<double> >*[dim+1];
+ coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
// for all element sides
for (int i = 0; i < dim + 1; i++) {
@@ -248,10 +248,7 @@ namespace AMDiS {
(*neumannOperators)[face]->getAssembler(omp_get_thread_num())->getZeroOrderAssembler()->
initElement(elInfo);
- const double *uhAtQp = dv->getVecAtQPs(elInfo,
- quadrature,
- NULL,
- NULL);
+ const double *uhAtQp = dv->getVecAtQPs(elInfo, quadrature, NULL, NULL);
std::vector<double> f(numPoints, 0.0);
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index 8c9cb62f78c84fabeaf38075dabff0252fe86d86..9b8482276c274886be4b6df033a4c4344c08454e 100644
--- a/AMDiS/src/Traverse.cc
+++ b/AMDiS/src/Traverse.cc
@@ -187,7 +187,7 @@ namespace AMDiS {
int i = info_stack[stack_used];
el = const_cast<Element*>(((i == 0) ? el->getFirstChild() : el->getSecondChild()));
info_stack[stack_used]++;
- elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
+ fillElInfo(i);
stack_used++;
TEST_EXIT_DBG(stack_used < stack_size)
@@ -279,7 +279,7 @@ namespace AMDiS {
int i = info_stack[stack_used];
info_stack[stack_used]++;
- elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
+ fillElInfo(i);
stack_used++;
TEST_EXIT_DBG(stack_used < stack_size)
@@ -347,7 +347,7 @@ namespace AMDiS {
int i = info_stack[stack_used];
info_stack[stack_used]++;
- elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
+ fillElInfo(i);
stack_used++;
TEST_EXIT_DBG(stack_used < stack_size)
@@ -412,12 +412,12 @@ namespace AMDiS {
while (elinfo_stack[stack_used]->getElement()->getFirstChild() &&
info_stack[stack_used] < 2) {
- if (stack_used >= stack_size-1)
+ if (stack_used >= stack_size - 1)
enlargeTraverseStack();
int i = info_stack[stack_used];
info_stack[stack_used]++;
- elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
+ fillElInfo(i);
stack_used++;
info_stack[stack_used] = 0;
}
@@ -482,7 +482,7 @@ namespace AMDiS {
if (stack_used >= stack_size - 1)
enlargeTraverseStack();
i = 1 - neighbour;
- elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
+ fillElInfo(i);
info_stack[stack_used] = i + 1;
stack_used++;
info_stack[stack_used] = 0;
@@ -540,11 +540,11 @@ namespace AMDiS {
elinfo2 = save_elinfo_stack[stack2_used];
el2 = elinfo2->getElement();
- if (stack_used >= stack_size-1)
+ if (stack_used >= stack_size - 1)
enlargeTraverseStack();
i = 2 - info_stack[stack_used];
- info_stack[stack_used] = i+1;
- elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
+ info_stack[stack_used] = i + 1;
+ fillElInfo(i);
stack_used++;
info_stack[stack_used] = 0;
nb = 0;
@@ -558,8 +558,8 @@ namespace AMDiS {
if (nb < 2) { /* go down one level in hierarchy */
if (stack_used >= stack_size-1)
enlargeTraverseStack();
- info_stack[stack_used] = 2-nb;
- elinfo_stack[stack_used+1]->fillElInfo(1 - nb, elinfo_stack[stack_used]);
+ info_stack[stack_used] = 2 - nb;
+ fillElInfo(1 - nb);
stack_used++;
info_stack[stack_used] = 0;
elinfo = elinfo_stack[stack_used];
@@ -594,10 +594,10 @@ namespace AMDiS {
nb = 2;
}
- info_stack[stack_used] = i+1;
- if (stack_used >= stack_size-1)
+ info_stack[stack_used] = i + 1;
+ if (stack_used >= stack_size - 1)
enlargeTraverseStack();
- elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
+ fillElInfo(i);
stack_used++;
info_stack[stack_used] = 0;
@@ -630,11 +630,11 @@ namespace AMDiS {
if (elinfo->getNeighbour(opp_vertex) != old_elinfo->getElement()) {
MSG(" looking for neighbour %d of element %d at %p\n",
- neighbour, old_elinfo->getElement()->getIndex(), reinterpret_cast<void*>( old_elinfo->getElement()));
+ neighbour, old_elinfo->getElement()->getIndex(), reinterpret_cast<void*>(old_elinfo->getElement()));
MSG(" originally: neighbour %d of element %d at %p\n",
- sav_neighbour, sav_index, reinterpret_cast<void*>( old_elinfo->getElement()));
+ sav_neighbour, sav_index, reinterpret_cast<void*>(old_elinfo->getElement()));
MSG(" got element %d at %p with opp_vertex %d neigh %d\n",
- elinfo->getElement()->getIndex(), reinterpret_cast<void*>( elinfo->getElement()),
+ elinfo->getElement()->getIndex(), reinterpret_cast<void*>(elinfo->getElement()),
opp_vertex, (elinfo->getNeighbour(opp_vertex))?(elinfo->getNeighbour(opp_vertex))->getIndex():-1);
TEST_EXIT_DBG(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
("didn't succeed !?!?!?\n");
@@ -679,15 +679,16 @@ namespace AMDiS {
TEST_EXIT_DBG(elinfo_old == elinfo_stack[stack_used])("invalid old elinfo");
elinfo_stack[stack_used]->testFlag(Mesh::FILL_NEIGH);
- el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo_stack[stack_used]->getElement()));
+ el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo_stack[stack_used]->getElement()));
sav_index = el->getIndex();
sav_el = el;
/* first, goto to leaf level, if necessary... */
if (!(el->isLeaf()) && (neighbour < 2)) {
- if (stack_used >= static_cast<int>( elinfo_stack.size())-1) enlargeTraverseStack();
+ if (stack_used >= static_cast<int>(elinfo_stack.size()) - 1)
+ enlargeTraverseStack();
i = 1 - neighbour;
- elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
+ fillElInfo(i);
info_stack[stack_used] = i + 1;
stack_used++;
neighbour = 2;
@@ -736,7 +737,7 @@ namespace AMDiS {
stack2_used = 1;
}
elinfo2 = save_elinfo_stack[stack2_used];
- el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
+ el2 = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo2->getElement()));
i = traverse_mel->getOppVertex(nb);
traverse_mel = traverse_mel->getNeighbour(nb);
@@ -745,7 +746,7 @@ namespace AMDiS {
nb = i;
stack_used = 1;
- elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>( traverse_mel));
+ elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>(traverse_mel));
info_stack[stack_used] = 0;
} else { /* goto other child */
@@ -755,17 +756,17 @@ namespace AMDiS {
stack2_used++; /* go down one level in OLD hierarchy */
}
elinfo2 = save_elinfo_stack[stack2_used];
- el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
+ el2 = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo2->getElement()));
if (stack_used >= stack_size-1) {
enlargeTraverseStack();
}
i = 2 - info_stack[stack_used];
- info_stack[stack_used] = i+1;
- elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
+ info_stack[stack_used] = i + 1;
+ fillElInfo(i);
stack_used++;
- nb = 1-i;
+ nb = 1 - i;
}
/****************************************************************************/
@@ -775,15 +776,15 @@ namespace AMDiS {
/****************************************************************************/
elinfo = elinfo_stack[stack_used];
- el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
+ el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo->getElement()));
while (el->getFirstChild()) {
if (nb < 2) { /* go down one level in hierarchy */
- if (stack_used >= stack_size-1)
+ if (stack_used >= stack_size - 1)
enlargeTraverseStack();
- elinfo_stack[stack_used+1]->fillElInfo(1-nb, elinfo_stack[stack_used]);
- info_stack[stack_used] = 2-nb;
+ fillElInfo(1 - nb);
+ info_stack[stack_used] = 2 - nb;
stack_used++;
nb = 2;
}
@@ -794,25 +795,23 @@ namespace AMDiS {
/* use child i, neighbour of el2->child[nb-1] */
i = 2 - save_info_stack[stack2_used];
TEST_EXIT_DBG(i < 2)("invalid OLD refinement?");
- info_stack[stack_used] = i+1;
- elinfo_stack[stack_used+1]->fillElInfo(i, elinfo_stack[stack_used]);
+ info_stack[stack_used] = i + 1;
+ fillElInfo(i);
stack_used++;
nb = i;
elinfo = elinfo_stack[stack_used];
- el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
+ el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo->getElement()));
stack2_used++;
- if (save_stack_used > stack2_used) {
- stack2_used++; /* go down one level in OLD hierarchy */
- }
+ if (save_stack_used > stack2_used)
+ stack2_used++; /* go down one level in OLD hierarchy */
elinfo2 = save_elinfo_stack[stack2_used];
- el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
+ el2 = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo2->getElement()));
- }
- else { /* now we're done... */
+ } else { /* now we're done... */
elinfo = elinfo_stack[stack_used];
- el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
+ el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo->getElement()));
}
}
@@ -874,5 +873,4 @@ namespace AMDiS {
elInfo.setRefinementPath(rPath);
elInfo.setRefinementPathLength(levelDif);
}
-
}
diff --git a/AMDiS/src/Traverse.h b/AMDiS/src/Traverse.h
index b4ded5bd6966441337abfe71a9bc2db36619a01e..19479de134f5d208a8575a0b387d6ab6f65b4706 100644
--- a/AMDiS/src/Traverse.h
+++ b/AMDiS/src/Traverse.h
@@ -38,6 +38,8 @@
#include "ElInfo.h"
#include "ElInfoStack.h"
#include "OpenMP.h"
+#include "Mesh.h"
+#include "BasisFunction.h"
#include "AMDiS_fwd.h"
namespace AMDiS {
@@ -136,6 +138,42 @@ namespace AMDiS {
return stack_used;
}
+ /** \brief
+ * Adds new FE space to \ref feSpaces. Is used, when local indices of leaf
+ * elements should be computed during mesh traverse.
+ */
+ void addFeSpace(const FiniteElemSpace *feSpace)
+ {
+ feSpaces.insert(feSpace);
+ }
+
+ private:
+ /** \brief
+ * Fills the element info data of a new ElInfo object on the top of the stack.
+ *
+ * \param[in] iChild Defines if the new element on the stack is the left or the
+ * right children of its parent element.
+ */
+ inline void fillElInfo(int iChild)
+ {
+ FUNCNAME("Traverse::fillElInfo()");
+
+ ElInfo *elInfo = elinfo_stack[stack_used + 1];
+ elInfo->fillElInfo(iChild, elinfo_stack[stack_used]);
+
+ if (traverse_fill_flag.isSet(Mesh::FILL_LOCAL_INDICES) &&
+ elInfo->getElement()->isLeaf()) {
+
+ TEST_EXIT_DBG(feSpaces.size() > 0)("There are no FE space defined!\n");
+
+ for (std::set<const FiniteElemSpace*>::iterator it = feSpaces.begin();
+ it != feSpaces.end(); ++it)
+ (*it)->getBasisFcts()->getLocalIndices(elInfo->getElement(),
+ (*it)->getAdmin(),
+ elInfo->getLocalIndices(*it));
+ }
+ }
+
private:
/// Enlargement of the stack
void enlargeTraverseStack();
@@ -224,6 +262,12 @@ namespace AMDiS {
///
int id;
+ /** \brief
+ * If local indices should be field for leaf elements, here all used FE spaces
+ * are stored, for which the local indices should be computed.
+ */
+ std::set<const FiniteElemSpace*> feSpaces;
+
/** \brief
* Is used for parallel mesh traverse. The thread with the id
* myThreadId is only allowed to access coarse elements, which id
diff --git a/AMDiS/src/TraverseParallel.h b/AMDiS/src/TraverseParallel.h
index 44046eaaa0179047bc78cdffbb98f334db0f7499..aa99479746d2703df2cf5c6f01e2629d2bead268 100644
--- a/AMDiS/src/TraverseParallel.h
+++ b/AMDiS/src/TraverseParallel.h
@@ -44,6 +44,12 @@ namespace AMDiS {
ElInfo* traverseNext(ElInfo* elInfoOld);
+ void addFeSpace(const FiniteElemSpace* feSpace)
+ {
+ for (unsigned int i = 0; i < stacks.size(); i++)
+ stacks[i].addFeSpace(feSpace);
+ }
+
private:
/// Number of threads using the stack in parallel.
int nThreads;