diff --git a/AMDiS/compositeFEM/src/CFE_Integration.cc b/AMDiS/compositeFEM/src/CFE_Integration.cc
index 8c4e0b29e6a2bfca8ecfef2ab0f406051c4d2083..9fe952d9b94c881047514eb9d8fbc445bf38f508 100644
--- a/AMDiS/compositeFEM/src/CFE_Integration.cc
+++ b/AMDiS/compositeFEM/src/CFE_Integration.cc
@@ -279,12 +279,12 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> > &surfVert)
{
double surfDet;
- int dim = surfVert[0].getSize()-1;
- FixVec<WorldVector<double>, VERTEX> worldCoords(dim-1, NO_INIT);
+ int dim = surfVert[0].getSize() - 1;
+ FixVec<WorldVector<double>, VERTEX> worldCoords(dim - 1, NO_INIT);
// transform barycentric coords to world coords
- for(int i = 0; i < dim; i++) {
- loc_elInfo->coordToWorld(surfVert[i], &worldCoords[i]);
+ for (int i = 0; i < dim; i++) {
+ loc_elInfo->coordToWorld(surfVert[i], worldCoords[i]);
}
// calculate determinant for surface
diff --git a/AMDiS/compositeFEM/src/CFE_NormAndErrorFcts.cc b/AMDiS/compositeFEM/src/CFE_NormAndErrorFcts.cc
index 6ac83b9a1969055e9b74130609428d4a369dbb08..fd1895ca1519b07f1e3b4c690b6e195c2179b652 100644
--- a/AMDiS/compositeFEM/src/CFE_NormAndErrorFcts.cc
+++ b/AMDiS/compositeFEM/src/CFE_NormAndErrorFcts.cc
@@ -17,11 +17,11 @@ namespace AMDiS {
const double &fac)
{
double val = 0.0;
- const WorldVector<double> *worldCoordsAtQP;
+ WorldVector<double> worldCoordsAtQP;
for (int iq = 0; iq < nQPts; ++iq) {
- worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
- val += q->getWeight(iq) * fabs((*f)(*worldCoordsAtQP));
+ elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
+ val += q->getWeight(iq) * fabs((*f)(worldCoordsAtQP));
}
double nrm = det * val;
@@ -34,11 +34,11 @@ namespace AMDiS {
const double &fac)
{
double val = 0.0;
- const WorldVector<double> *worldCoordsAtQP;
+ WorldVector<double> worldCoordsAtQP;
for (int iq = 0; iq < nQPts; ++iq) {
- worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
- val += q->getWeight(iq) * sqr((*f)(*worldCoordsAtQP));
+ elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
+ val += q->getWeight(iq) * sqr((*f)(worldCoordsAtQP));
}
double nrm = det * val;
@@ -51,15 +51,15 @@ namespace AMDiS {
const double &fac)
{
double val = 0.0;
- double norm_grd2;
- const WorldVector<double> *worldCoordsAtQP;
+ double norm_grd2 = 0.0;
+ WorldVector<double> worldCoordsAtQP;
for (int iq = 0; iq < nQPts; ++iq) {
- worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
+ elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
norm_grd2 = 0.0;
for (int j = 0; j < dim; j++)
- norm_grd2 += sqr(((*grd)(*worldCoordsAtQP))[j]);
+ norm_grd2 += sqr(((*grd)(worldCoordsAtQP))[j]);
val += q->getWeight(iq) * norm_grd2;
}
@@ -142,14 +142,14 @@ namespace AMDiS {
q,
NULL,
NULL);
- const WorldVector<double> *worldCoordsAtQP;
+ WorldVector<double> worldCoordsAtQP;
for (int iq = 0; iq < nQPts; ++iq) {
- worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
- val += q->getWeight(iq) * sqr((*u)(*worldCoordsAtQP) - uhAtQPs[iq]);
+ elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
+ val += q->getWeight(iq) * sqr((*u)(worldCoordsAtQP) - uhAtQPs[iq]);
if (relErr)
- val_nrm += q->getWeight(iq) * sqr((*u)(*worldCoordsAtQP));
+ val_nrm += q->getWeight(iq) * sqr((*u)(worldCoordsAtQP));
}
double err = det * val;
@@ -173,23 +173,22 @@ namespace AMDiS {
q,
NULL,
NULL);
- const WorldVector<double> *worldCoordsAtQP;
+ WorldVector<double> worldCoordsAtQP;
for (int iq = 0; iq < nQPts; ++iq) {
-
- worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
+ elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
norm_err_grd2 = 0.0;
for (int j = 0; j < dim; ++j)
norm_err_grd2 +=
- sqr(((*grdu)(*worldCoordsAtQP))[j] - grdUhAtQPs[iq][j]);
+ sqr(((*grdu)(worldCoordsAtQP))[j] - grdUhAtQPs[iq][j]);
val += q->getWeight(iq) * norm_err_grd2;
if (relErr) {
norm_grd2 = 0.0;
for (int j = 0; j < dim; ++j)
- norm_grd2 += sqr(((*grdu)(*worldCoordsAtQP))[j]);
+ norm_grd2 += sqr(((*grdu)(worldCoordsAtQP))[j]);
val_nrm += q->getWeight(iq) * norm_grd2;
}
diff --git a/AMDiS/compositeFEM/src/ElementLevelSet.cc b/AMDiS/compositeFEM/src/ElementLevelSet.cc
index 0f4176361207cc9702ae5ce4ac4b7097eba49621..5ef6b0763188c8aa2f55c83e9450419e869841f7 100644
--- a/AMDiS/compositeFEM/src/ElementLevelSet.cc
+++ b/AMDiS/compositeFEM/src/ElementLevelSet.cc
@@ -249,8 +249,8 @@ ElementLevelSet::calcIntersecNormal_2d(WorldVector<double> &normalVec)
// ===== Get world coordinates of intersection points. =====
WorldVector<double> sP1;
WorldVector<double> sP2;
- elInfo->coordToWorld((*elIntersecPoints)[0], &sP1);
- elInfo->coordToWorld((*elIntersecPoints)[1], &sP2);
+ elInfo->coordToWorld((*elIntersecPoints)[0], sP1);
+ elInfo->coordToWorld((*elIntersecPoints)[1], sP2);
// ===== Calculate normal vector. =====
double norm2 = 0.0;
@@ -326,9 +326,9 @@ ElementLevelSet::calcIntersecNormal_3d(WorldVector<double> &normalVec)
WorldVector<double> sP1;
WorldVector<double> sP2;
WorldVector<double> sP3;
- elInfo->coordToWorld((*elIntersecPoints)[0], &sP1);
- elInfo->coordToWorld((*elIntersecPoints)[1], &sP2);
- elInfo->coordToWorld((*elIntersecPoints)[2], &sP3);
+ elInfo->coordToWorld((*elIntersecPoints)[0], sP1);
+ elInfo->coordToWorld((*elIntersecPoints)[1], sP2);
+ elInfo->coordToWorld((*elIntersecPoints)[2], sP3);
// ===== Calculate normal vector. =====
WorldVector<double> A = sP2-sP1;
diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 182cb004e37eb806232c68466f437d7b9259b7bb..91c80f46d784ed5cb484486d056e8d8649da82bc 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -406,13 +406,11 @@ namespace AMDiS {
ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);
}
- ElementMatrix *Assembler::initElementMatrix(ElementMatrix *elMat,
- const ElInfo *rowElInfo,
- const ElInfo *colElInfo)
+ void Assembler::initElementMatrix(ElementMatrix *elMat,
+ const ElInfo *rowElInfo,
+ const ElInfo *colElInfo)
{
- if (!elMat) {
- elMat = NEW ElementMatrix(nRow, nCol);
- }
+ TEST_EXIT_DBG(elMat)("No ElementMatrix!\n");
elMat->set(0.0);
@@ -434,25 +432,17 @@ namespace AMDiS {
&(elMat->colIndices));
}
}
-
- return elMat;
}
- ElementVector *Assembler::initElementVector(ElementVector *elVec,
- const ElInfo *elInfo)
+ void Assembler::initElementVector(ElementVector *elVec, const ElInfo *elInfo)
{
- if (!elVec) {
- elVec = NEW ElementVector(nRow);
- }
+ TEST_EXIT_DBG(elVec)("No ElementVector!\n");
elVec->set(0.0);
-
- DOFAdmin *rowAdmin = rowFESpace->getAdmin();
- Element *element = elInfo->getElement();
-
- rowFESpace->getBasisFcts()->getLocalIndicesVec(element, rowAdmin, &(elVec->dofIndices));
-
- return elVec;
+
+ rowFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
+ rowFESpace->getAdmin(),
+ &(elVec->dofIndices));
}
void Assembler::checkQuadratures()
diff --git a/AMDiS/src/Assembler.h b/AMDiS/src/Assembler.h
index dc489aff32f519ce5badbdedb6c0821f9ebcbf2b..cc537d0ddd5649c2915fece542e69ef8b3dfddc0 100644
--- a/AMDiS/src/Assembler.h
+++ b/AMDiS/src/Assembler.h
@@ -69,13 +69,13 @@ namespace AMDiS {
virtual ~Assembler() {}
- ElementMatrix *initElementMatrix(ElementMatrix *elMat,
- const ElInfo *rowElInfo,
- const ElInfo *colElInfo = NULL);
+ void initElementMatrix(ElementMatrix *elMat,
+ const ElInfo *rowElInfo,
+ const ElInfo *colElInfo = NULL);
- ElementVector *initElementVector(ElementVector *elVec,
- const ElInfo *elInfo);
+ void initElementVector(ElementVector *elVec,
+ const ElInfo *elInfo);
/** \brief
* Assembles the element matrix for the given ElInfo
@@ -230,9 +230,7 @@ namespace AMDiS {
const ElInfo *smallElInfo, const ElInfo *largeElInfo,
ElementVector *vec);
- /** \brief
- * Checks whether this is a new travese.
- */
+ /// Checks whether this is a new travese.
inline void checkForNewTraverse() {
if (lastTraverseId != ElInfo::traverseId[omp_get_thread_num()]) {
lastVecEl = lastMatEl = NULL;
diff --git a/AMDiS/src/BoundaryManager.cc b/AMDiS/src/BoundaryManager.cc
index 091c39325ab07de1656b267906adec763c45165f..8b7298b8d7ebaaabc7e30217e726ada026ee7367 100644
--- a/AMDiS/src/BoundaryManager.cc
+++ b/AMDiS/src/BoundaryManager.cc
@@ -12,6 +12,7 @@ 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] = GET_MEMORY(BoundaryType, allocatedMemoryLocalBounds);
@@ -23,6 +24,7 @@ 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] = GET_MEMORY(BoundaryType, allocatedMemoryLocalBounds);
}
@@ -51,29 +53,26 @@ namespace AMDiS {
void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo,
DOFVectorBase<double> *vec)
{
- // ===== fill local conditions ==============================================
- const FiniteElemSpace *feSpace = vec->getFESpace();
- Vector<DegreeOfFreedom> dofIndices;
- const BasisFunction *basisFcts = feSpace->getBasisFcts();
- int nBasFcts = basisFcts->getNumber();
- DOFAdmin *admin = feSpace->getAdmin();
-
- std::map<BoundaryType, BoundaryCondition*>::iterator it;
-
if (localBCs.size() > 0) {
+ const FiniteElemSpace *feSpace = vec->getFESpace();
+ Vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
+ const BasisFunction *basisFcts = feSpace->getBasisFcts();
+ int nBasFcts = basisFcts->getNumber();
+ std::map<BoundaryType, BoundaryCondition*>::iterator it;
// get boundaries of all DOFs
BoundaryType *localBound = localBounds[omp_get_thread_num()];
basisFcts->getBound(elInfo, localBound);
// get dof indices
- basisFcts->getLocalIndicesVec(elInfo->getElement(), admin, &dofIndices);
+ basisFcts->getLocalIndicesVec(elInfo->getElement(),
+ feSpace->getAdmin(), &dofVec);
// apply non dirichlet boundary conditions
for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) {
if (!(*it).second->isDirichlet()) {
- (*it).second->fillBoundaryCondition(vec, elInfo, &dofIndices[0], localBound, nBasFcts);
+ (*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0], localBound, nBasFcts);
}
}
}
@@ -82,7 +81,7 @@ namespace AMDiS {
for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) {
if ((*it).second->isDirichlet()) {
- (*it).second->fillBoundaryCondition(vec, elInfo, &dofIndices[0], localBound, nBasFcts);
+ (*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0], localBound, nBasFcts);
}
}
}
@@ -92,28 +91,26 @@ namespace AMDiS {
void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo,
DOFMatrix *mat)
{
- // ===== fill local conditions ==============================================
- const FiniteElemSpace *feSpace = mat->getRowFESpace();
- Vector<DegreeOfFreedom> dofIndices;
- const BasisFunction *basisFcts = feSpace->getBasisFcts();
- int nBasFcts = basisFcts->getNumber();
- DOFAdmin *admin = feSpace->getAdmin();
-
- std::map<BoundaryType, BoundaryCondition*>::iterator it;
-
if (localBCs.size() > 0) {
+ const FiniteElemSpace *feSpace = mat->getRowFESpace();
+ Vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
+ const BasisFunction *basisFcts = feSpace->getBasisFcts();
+ int nBasFcts = basisFcts->getNumber();
+ std::map<BoundaryType, BoundaryCondition*>::iterator it;
+
// get boundaries of all DOFs
BoundaryType *localBound = localBounds[omp_get_thread_num()];
basisFcts->getBound(elInfo, localBound);
// get dof indices
- basisFcts->getLocalIndicesVec(elInfo->getElement(), admin, &dofIndices);
+ basisFcts->getLocalIndicesVec(elInfo->getElement(),
+ feSpace->getAdmin(), &dofVec);
// apply non dirichlet boundary conditions
for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) {
if (!(*it).second->isDirichlet()) {
- (*it).second->fillBoundaryCondition(mat, elInfo, &dofIndices[0], localBound, nBasFcts);
+ (*it).second->fillBoundaryCondition(mat, elInfo, &dofVec[0], localBound, nBasFcts);
}
}
}
@@ -122,7 +119,7 @@ namespace AMDiS {
for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) {
if ((*it).second->isDirichlet()) {
- (*it).second->fillBoundaryCondition(mat, elInfo, &dofIndices[0], localBound, nBasFcts);
+ (*it).second->fillBoundaryCondition(mat, elInfo, &dofVec[0], localBound, nBasFcts);
}
}
}
diff --git a/AMDiS/src/BoundaryManager.h b/AMDiS/src/BoundaryManager.h
index f90ac30d987a6a651989a9bc80952a446ab84340..a0e3bd0aa9c554a1466b8e3ee29f88f5b6b44a80 100644
--- a/AMDiS/src/BoundaryManager.h
+++ b/AMDiS/src/BoundaryManager.h
@@ -31,6 +31,7 @@ namespace AMDiS {
class DOFMatrix;
class FiniteElemSpace;
+ template<typename T> class Vector;
template<typename T> class DOFVectorBase;
// ============================================================================
@@ -107,16 +108,15 @@ namespace AMDiS {
};
protected:
- /** \brief
- * Map of managed local boundary conditions.
- */
+ /// Map of managed local boundary conditions.
std::map<BoundaryType, BoundaryCondition*> localBCs;
- /** \brief
- * Temporary thread-safe variable for functions fillBoundaryconditions.
- */
+ /// Temporary thread-safe variable for functions fillBoundaryconditions.
std::vector<BoundaryType*> localBounds;
+ /// Temporary thread-safe variable for functions fillBoundaryconditions.
+ std::vector<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 7994e2a45a8806ca061d9d9188bd851c9284c7a3..0ad6285634994e1389563e7341b10336fb9f9aaa 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -249,10 +249,6 @@ namespace AMDiS {
int nRow = elMat.rowIndices.getSize();
int nCol = elMat.colIndices.getSize();
-// elMat.rowIndices.print();
-// elMat.colIndices.print();
-// elMat.print();
-
for (int i = 0; i < nRow; i++) { // for all rows of element matrix
row = elMat.rowIndices[i];
BoundaryCondition *condition =
@@ -433,37 +429,40 @@ namespace AMDiS {
}
- void DOFMatrix::assemble(double factor, ElInfo *elInfo,
- const BoundaryType *bound, Operator *op)
+ void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
{
FUNCNAME("DOFMatrix::assemble()");
- if (!op && operators.size() == 0) {
- return;
- }
+ if (operators.size() == 0)
+ return;
- Operator *operat = op ? op : operators[0];
- operat->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
+ operators[0]->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
- if (op) {
- op->getElementMatrix(elInfo, elementMatrix);
- } else {
- std::vector<Operator*>::iterator it;
- std::vector<double*>::iterator factorIt;
- for (it = operators.begin(), factorIt = operatorFactor.begin();
- it != operators.end();
- ++it, ++factorIt) {
+ std::vector<Operator*>::iterator it = operators.begin();
+ std::vector<double*>::iterator factorIt = operatorFactor.begin();
+ for (; it != operators.end(); ++it, ++factorIt) {
if ((*it)->getNeedDualTraverse() == false) {
- (*it)->getElementMatrix(elInfo,
- elementMatrix,
- *factorIt ? **factorIt : 1.0);
+ (*it)->getElementMatrix(elInfo,
+ elementMatrix,
+ *factorIt ? **factorIt : 1.0);
}
- }
- }
+ }
addElementMatrix(factor, *elementMatrix, bound);
}
+ void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
+ Operator *op)
+ {
+ FUNCNAME("DOFMatrix::assemble()");
+
+ TEST_EXIT_DBG(op)("No operator!\n");
+
+ op->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
+ op->getElementMatrix(elInfo, elementMatrix);
+ addElementMatrix(factor, *elementMatrix, bound);
+ }
+
void DOFMatrix::assemble(double factor,
ElInfo *rowElInfo, ElInfo *colElInfo,
ElInfo *smallElInfo, ElInfo *largeElInfo,
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index 26706cb05cddd0218851bd720955eb17038de9be..00fe35d611f4ae33e3b4c69aa6820239d6cb6138 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -357,9 +357,10 @@ namespace AMDiS {
* matrix should be cleared by calling clear dof matrix().
*/
- void assemble(double factor, ElInfo *elInfo,
- const BoundaryType *bound, Operator *op = NULL);
+ void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound);
+ void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
+ Operator *op);
void assemble(double factor,
ElInfo *rowElInfo, ElInfo *colElInfo,
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 96d96a1b8390e2d5824666fe204acfcb27b846aa..b2da7a69445ce8195ddaa2690eef1034c52b2bf2 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -511,7 +511,6 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
} else {
- DimVec<double> *coords1 = NULL;
DimVec<double> coords2(feSpace->getMesh()->getDim(), NO_INIT);
DualTraverse dualStack;
ElInfo *elInfo1, *elInfo2;
@@ -534,8 +533,7 @@ namespace AMDiS {
for (int i = 0; i < nBasisFcts; i++) {
if (vec[myLocalIndices[i]] == 0.0) {
- coords1 = basisFcts->getCoords(i);
- elInfo1->coordToWorld(*coords1, &worldVec);
+ elInfo1->coordToWorld(*(basisFcts->getCoords(i)), worldVec);
elInfo2->worldToCoord(worldVec, &coords2);
bool isPositive = true;
@@ -821,13 +819,14 @@ namespace AMDiS {
{
FUNCNAME("DOFVector::assemble()");
+ TEST_EXIT_DBG(this->elementVector)("No ElementVector!\n");
+
if (!(op || this->operators.size()))
return;
Operator *operat = op ? op : this->operators[0];
- this->elementVector =
- operat->getAssembler(omp_get_thread_num())->initElementVector(this->elementVector, elInfo);
+ operat->getAssembler(omp_get_thread_num())->initElementVector(this->elementVector, elInfo);
if (op) {
op->getElementVector(elInfo, this->elementVector);
@@ -862,9 +861,8 @@ namespace AMDiS {
Operator *operat = op ? op : this->operators[0];
- this->elementVector =
- operat->getAssembler(omp_get_thread_num())->
- initElementVector(this->elementVector, mainElInfo);
+ operat->getAssembler(omp_get_thread_num())->
+ initElementVector(this->elementVector, mainElInfo);
if (op) {
ERROR_EXIT("TODO");
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index c4f011864e1273415a30d354229ff6ab6324c9b1..9b3378072737279471975cd3ee12c6f07bf968c1 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -560,7 +560,7 @@ namespace AMDiS {
/** \brief
* Assignment operator between two vectors
*/
- DOFVector<T>& operator=(const DOFVector<T>& );
+ DOFVector<T>& operator=(const DOFVector<T>&);
/** \brief
* vec[i] = v.vec[i]
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 533c5da27fb9facba7f79bb5129b48fec850cc1a..9bf1874b1f720c2e827a45e90bdabf406386e5c1 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -26,7 +26,6 @@ namespace AMDiS {
DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
: feSpace(f),
name(n),
- elementVector(NULL),
boundaryManager(NULL)
{
nBasFcts = feSpace->getBasisFcts()->getNumber();
@@ -45,6 +44,8 @@ namespace AMDiS {
grdTmp[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
D2Phis[i] = NEW DimMat<double>(dim, NO_INIT);
}
+
+ elementVector = NEW ElementVector(this->nBasFcts);
}
template<typename T>
@@ -743,8 +744,9 @@ namespace AMDiS {
DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs )
{
feSpace = rhs.feSpace;
+ this->nBasFcts = rhs.nBasFcts;
vec = rhs.vec;
- this->elementVector = NULL;
+ this->elementVector = new ElementVector(this->nBasFcts);
interFct = rhs.interFct;
coarsenOperation = rhs.coarsenOperation;
this->operators = rhs.operators;
@@ -756,7 +758,7 @@ namespace AMDiS {
// boundaryManager->setDOFVector(this);
} else {
- this->boundaryManager=NULL;
+ this->boundaryManager = NULL;
}
return *this;
diff --git a/AMDiS/src/DataCollector.cc b/AMDiS/src/DataCollector.cc
index a80867f3a5c472fa0ca6fd59e0ca1240f827fc17..b5bb4c6e03b3752eb37bc6dd5bcade82f35088ff 100644
--- a/AMDiS/src/DataCollector.cc
+++ b/AMDiS/src/DataCollector.cc
@@ -220,7 +220,7 @@ namespace AMDiS {
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
DegreeOfFreedom vertexDOF;
WorldVector<double> vertexCoords;
-
+
// create ElementInfo
ElementInfo elementInfo(dim_);
@@ -326,7 +326,7 @@ namespace AMDiS {
for (int i = mesh_->getGeo(VERTEX); i < nBasisFcts_; i++) {
DegreeOfFreedom dofi = localDOFs_[i];
- elInfo->coordToWorld(*basisFcts_->getCoords(i), vertexCoords);
+ elInfo->coordToWorld(*basisFcts_->getCoords(i), *vertexCoords);
if ((*interpPointInd_)[dofi] == -1) {
// mark as interpolation point
diff --git a/AMDiS/src/DataCollector.h b/AMDiS/src/DataCollector.h
index 88e44be2c80fd18791e552fa1a02298a6a0c2bc5..d917b49be96400875824f8842127a127a5c634a4 100644
--- a/AMDiS/src/DataCollector.h
+++ b/AMDiS/src/DataCollector.h
@@ -213,44 +213,28 @@ namespace AMDiS {
*/
int nPreDofs_;
- /** \brief
- * Number of vertices.
- */
+ /// Number of vertices.
int nVertices_;
- /** \brief
- * Number of elements.
- */
+ /// Number of elements.
int nElements_;
- /** \brief
- * Total number of interpolation points.
- */
+ /// Total number of interpolation points.
int nInterpPoints_;
- /** \brief
- * Number of connections in periodic problems.
- */
+ /// Number of connections in periodic problems.
int nConnections_;
- /** \brief
- * Dimension of \ref mesh
- */
+ /// Dimension of \ref mesh
int dim_;
- /** \brief
- * Maps internal element indices to global output indices.
- */
+ /// Maps internal element indices to global output indices.
std::map<int, int> outputIndices_;
- /** \brief
- * Global interpolation point indexing
- */
+ /// Global interpolation point indexing
DOFVector<int> *interpPointInd_;
- /** \brief
- * Stores for each simplex the interpolation points.
- */
+ /// Stores for each simplex the interpolation points.
std::vector< std::vector<int> > interpPoints_;
/** \brief
@@ -259,24 +243,16 @@ namespace AMDiS {
*/
DOFVector< std::list<WorldVector<double> > > *interpPointCoords_;
- /** \brief
- * list of coords for each dof
- */
+ /// list of coords for each dof
DOFVector< std::list<WorldVector<double> > > *dofCoords_;
- /** \brief
- * List that stores an ElementInfo for each element.
- */
+ /// List that stores an ElementInfo for each element.
std::list<ElementInfo> elements_;
- /** \brief
- * List stat stores information about all periodics.
- */
+ /// List stat stores information about all periodics.
std::list<PeriodicInfo> periodicInfos_;
- /** \brief
- * Stores a list of vertex infos for each dof.
- */
+ /// Stores a list of vertex infos for each dof.
DOFVector< std::list<VertexInfo> > *vertexInfos_;
/** \brief
@@ -285,44 +261,28 @@ namespace AMDiS {
*/
std::vector<DimVec<bool> > periodicConnections_;
- /** \brief
- * Stores if element data was collected before.
- */
+ /// Stores if element data was collected before.
bool elementDataCollected_;
- /** \brief
- * Stores if value data was collected before.
- */
+ /// Stores if value data was collected before.
bool valueDataCollected_;
- /** \brief
- * Stores if periodic data was collected before.
- */
+ /// Stores if periodic data was collected before.
bool periodicDataCollected_;
- /** \brief
- * Pointer to a function which decides whether an element is considered.
- */
+ /// Pointer to a function which decides whether an element is considered.
bool (*writeElem_)(ElInfo*);
- /** \brief
- * Temporary variable used in functions addValueData() and addInterpData().
- */
+ /// Temporary variable used in functions addValueData() and addInterpData().
DegreeOfFreedom *localDOFs_;
-
- /** \brief
- * Temporary variable used in functions addValueData() and addInterpData().
- */
+
+ /// Temporary variable used in functions addValueData() and addInterpData().
BasisFunction *basisFcts_;
- /** \brief
- * Temporary variable used in functions addValueData() and addInterpData().
- */
+ /// Temporary variable used in functions addValueData() and addInterpData().
int nBasisFcts_;
- /** \brief
- * Temporary variable used in function \ref addValueData.
- */
+ /// Temporary variable used in function \ref addValueData.
WorldVector<double> *vertexCoords;
};
}
diff --git a/AMDiS/src/DirichletBC.cc b/AMDiS/src/DirichletBC.cc
index 35421b312ae89571591c8852856941e23e9c2243..705a882982424e411adb0107c058c21fec3ef914 100644
--- a/AMDiS/src/DirichletBC.cc
+++ b/AMDiS/src/DirichletBC.cc
@@ -55,8 +55,7 @@ namespace AMDiS {
for (int i = 0; i < nBasFcts; i++) {
if (localBound[i] == boundaryType) {
if (f) {
- DimVec<double> *coords = basFcts->getCoords(i);
- elInfo->coordToWorld(*coords, &(worldCoords[myRank]));
+ elInfo->coordToWorld(*(basFcts->getCoords(i)), worldCoords[myRank]);
double fAtCoords = (*f)(worldCoords[myRank]);
(*vector)[dofIndices[i]] = fAtCoords;
}
diff --git a/AMDiS/src/DualTraverse.cc b/AMDiS/src/DualTraverse.cc
index b890b8d0c19bc2ab6438e0c5f12bbb2e36612f94..8e5e7909b023c9901211e84ea4dfc7b1a95274bf 100644
--- a/AMDiS/src/DualTraverse.cc
+++ b/AMDiS/src/DualTraverse.cc
@@ -139,8 +139,6 @@ namespace AMDiS {
rest = 1.0;
inc1 = true;
inc2 = true;
- stack1.stopDisplacementCalculation();
- stack2.stopDisplacementCalculation();
} else {
// increment only small element
inc1 = (*elInfo1 == *elInfoSmall) ? true : false;
diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc
index a7058c5ea05ee3b61abcd921bb68850e2e36ae35..ca7fb7d07cbe5ff407600340f062ef29885f1176 100644
--- a/AMDiS/src/ElInfo.cc
+++ b/AMDiS/src/ElInfo.cc
@@ -24,7 +24,6 @@ namespace AMDiS {
parent_(NULL),
macroElement_(NULL),
level(0),
- displacement(0),
iChild(0),
coord_(mesh_->getDim(), NO_INIT),
boundary_(mesh_->getDim(), DEFAULT_VALUE, INTERIOR),
@@ -53,32 +52,24 @@ namespace AMDiS {
}
- /****************************************************************************/
- /* transform local coordintes l to world coordinates; if w is non NULL */
- /* store them at w otherwise return a pointer to some local static */
- /* area containing world coordintes */
- /****************************************************************************/
- const WorldVector<double> *ElInfo::coordToWorld(const DimVec<double>& l,
- WorldVector<double>* w) const
+ void ElInfo::coordToWorld(const DimVec<double>& l,
+ WorldVector<double>& w) const
{
int dim = l.getSize() - 1;
- static WorldVector<double> world;
- WorldVector<double> *ret = w ? w : &world;
double c = l[0];
for (int j = 0; j < dimOfWorld; j++)
- (*ret)[j] = c * coord_[0][j];
+ w[j] = c * coord_[0][j];
int vertices = Global::getGeo(VERTEX, dim);
for (int i = 1; i < vertices; i++) {
c = l[i];
for (int j = 0; j < dimOfWorld; j++)
- (*ret)[j] += c * coord_[i][j];
+ w[j] += c * coord_[i][j];
}
- return ret;
}
/****************************************************************************/
diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h
index 9f5c9bc2ee24a084d9de13242824ac5582f260bd..aeb16a496246d148a838139d762598a316558815 100644
--- a/AMDiS/src/ElInfo.h
+++ b/AMDiS/src/ElInfo.h
@@ -149,13 +149,6 @@ namespace AMDiS {
return level;
}
- /** \brief
- * Get ElInfo's \ref displacement. Is used only during dual traverse.
- */
- inline int getDisplacement() const {
- return displacement;
- }
-
/** \brief
* Get ElInfo's \ref iChild
*/
@@ -339,13 +332,6 @@ namespace AMDiS {
level = l;
}
- /** \brief
- * Set ElInfo's \ref displacement. Used only during dual traverse.
- */
- inline void setDisplacement(int d) {
- displacement = d;
- }
-
/** \brief
* Set ElInfo's \ref boundary_[i]
*/
@@ -406,14 +392,11 @@ namespace AMDiS {
void testFlag(const Flag& flag) const;
/** \brief
- * Returns a pointer to a vector, which contains the world coordinates
- * of a point in barycentric coordinates lambda with respect to \ref element.
- * If world is not NULL the world coordinates are stored in this vector.
- * Otherwise the function itself provides memory for this vector. In this
- * case the vector is overwritten during the next call of coordToWorld.
+ * Transforms local barycentric coordinates of a point defined on this element
+ * to global world coordinates.
*/
- virtual const WorldVector<double> *coordToWorld(const DimVec<double>& lambda,
- WorldVector<double>* world) const;
+ void coordToWorld(const DimVec<double>& lambda,
+ WorldVector<double>& world) const;
/** \brief
@@ -518,33 +501,6 @@ namespace AMDiS {
*/
unsigned char level;
- /** \brief
- * This value is set only during dual traverse, and only if this is the smaller
- * element of an element pair. Then the value describes the displacement of
- * this element to the first element in the refinement hierachy starting from
- * the bigger element of the element pair.
- *
- * Example, if the level difference of the two elements is 1:
- *
- * a
- * / \
- * b c
- *
- * The element b has displacement 0, the element c has displacement 1.
- *
- * Example, if the level difference of the two elements is 2:
- *
- * a
- * / \
- * / \
- * b c
- * / \ / \
- * d e f g
- *
- * The displacements are as follows: d = 0, e = 1, f = 2, g = 3.
- */
- int displacement;
-
/** \brief
* This ElInfo is the iChild-th child of the parent element.
*/
diff --git a/AMDiS/src/ElInfo1d.cc b/AMDiS/src/ElInfo1d.cc
index 9ae0c0a03074771091a5111bdf3cf955f25ec7b7..07215df766623685484e1891c9a646019df45e40 100644
--- a/AMDiS/src/ElInfo1d.cc
+++ b/AMDiS/src/ElInfo1d.cc
@@ -50,7 +50,7 @@ namespace AMDiS {
while (!(nb->isLeaf())) { // make nb nearest element
if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
- if (nb->getNewCoord(-1)) {
+ if (nb->isNewCoordSet()) {
oppC = *(nb->getNewCoord());
} else {
oppC = (mel->coord[i] + oppC) * 0.5;
@@ -211,7 +211,7 @@ namespace AMDiS {
const FixVec<WorldVector<double>, VERTEX> *old_coord = &(elInfoOld->coord_);
coord_[ichild] = (*old_coord)[ichild];
- if (elem->getNewCoord(-1)) {
+ if (elem->isNewCoordSet()) {
coord_[1 - ichild] = *(elem->getNewCoord());
} else {
coord_[1 - ichild] = ((*old_coord)[0] + (*old_coord)[1]) * 0.5;
@@ -241,7 +241,7 @@ namespace AMDiS {
if (nb) {
while (nb->getChild(0)) { // make nb nearest element
if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
- if (nb->getNewCoord(-1)) {
+ if (nb->isNewCoordSet()) {
oppC = *(nb->getNewCoord());
} else {
oppC = (coord_[i] + oppC) * 0.5;
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index fe0efee1f8670f83bf3cd08a041d4e0e2e408072..65e90afaac8564162caf62fae057a0b4431b51f4 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -73,7 +73,7 @@ namespace AMDiS {
oppVertex_[i] = 2;
if (fill_opp_coords) {
- if (nb->getNewCoord(-1)) {
+ if (nb->isNewCoordSet()) {
oppCoord_[i] = *(nb->getNewCoord());
} else {
oppCoord_[i] = (macroNeighbour->coord[0] + macroNeighbour->coord[1]) * 0.5;
@@ -180,7 +180,7 @@ namespace AMDiS {
fillFlag_.isSet(Mesh::FILL_DET) ||
fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
- if (elem->getNewCoord(-1)) {
+ if (elem->isNewCoordSet()) {
coord_[2] = *(elem->getNewCoord());
} else {
coord_[2].setMidpoint(elInfoOld->coord_[0], elInfoOld->coord_[1]);
@@ -222,7 +222,7 @@ namespace AMDiS {
oppVertex_[1] = 2;
if (fill_opp_coords) {
- if (elem->getSecondChild()->getNewCoord(-1)) {
+ if (elem->getSecondChild()->isNewCoordSet()) {
oppCoord_[1] = *(elem->getSecondChild()->getNewCoord());
} else {
oppCoord_[1].setMidpoint(elInfoOld->coord_[1],
@@ -262,7 +262,7 @@ namespace AMDiS {
oppVertex_[0] = 2;
if (fill_opp_coords) {
- if (nb->getNewCoord(-1)) {
+ if (nb->isNewCoordSet()) {
oppCoord_[0] = *(nb->getNewCoord());
} else {
oppCoord_[0].setMidpoint(elInfoOld->neighbourCoord_[2][1],
@@ -312,7 +312,7 @@ namespace AMDiS {
oppVertex_[0] = 2;
if (fill_opp_coords) {
- if (elem->getFirstChild()->getNewCoord(-1)) {
+ if (elem->getFirstChild()->isNewCoordSet()) {
oppCoord_[0] = *(elem->getFirstChild()->getNewCoord());
} else {
oppCoord_[0].setMidpoint(elInfoOld->coord_[0],
@@ -348,7 +348,7 @@ namespace AMDiS {
oppVertex_[1] = 2;
if (fill_opp_coords) {
- if (nb->getNewCoord(-1)) {
+ if (nb->isNewCoordSet()) {
oppCoord_[1] = *(nb->getNewCoord());
} else {
oppCoord_[1].setMidpoint(elInfoOld->neighbourCoord_[2][0],
diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc
index a20ede18a9db741ec903f414f6fdfc6a9fb4bf66..069e0c7131d82953700ed23f805c51614db31e35 100644
--- a/AMDiS/src/Element.cc
+++ b/AMDiS/src/Element.cc
@@ -361,15 +361,6 @@ namespace AMDiS {
}
}
- double Element::getNewCoord(int j) const {
- if (j >= 0) {
- TEST_EXIT_DBG(newCoord)("newCoord = NULL\n");
- return (*newCoord)[j];
- } else {
- return (newCoord != NULL);
- }
- }
-
void Element::eraseNewCoord() {
if (newCoord != NULL) {
DELETE newCoord;
diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h
index bd7769a2d4d158dbf625ea5cab168c21df30b192..5a050181f7c52ba98d011931fa21b5fc43b393cb 100644
--- a/AMDiS/src/Element.h
+++ b/AMDiS/src/Element.h
@@ -246,21 +246,18 @@ namespace AMDiS {
return mark;
}
- /** \brief
- * Returns \ref newCoord[i]
- */
- double getNewCoord(int j) const;
+ /// Returns \ref newCoord[i]
+ inline double getNewCoord(int i) const {
+ TEST_EXIT_DBG(newCoord)("newCoord = NULL\n");
+ return (*newCoord)[i];
+ }
- /** \brief
- * Returns Element's \ref index
- */
+ /// Returns Element's \ref index
inline int getIndex() const {
return index;
}
- /** \brief
- * Returns \ref newCoord
- */
+ /// Returns \ref newCoord
inline WorldVector<double>* getNewCoord() const {
return newCoord;
}
diff --git a/AMDiS/src/ElementFunction.h b/AMDiS/src/ElementFunction.h
index 157696889a659125a9a1553b1b6279d36c16e069..a075e4175f4e89373c2b478fbaf23cbc3f8321c9 100644
--- a/AMDiS/src/ElementFunction.h
+++ b/AMDiS/src/ElementFunction.h
@@ -72,26 +72,21 @@ namespace AMDiS {
class ElementFunctionAnalytic : public ElementFunction<T>
{
public:
- /** \brief
- * constructor
- */
+ /// constructor
ElementFunctionAnalytic(const AbstractFunction<T, WorldVector<double> > *fct)
: ElementFunction<T>(),
fct_(fct)
- {};
+ {}
- /** \brief
- * evaluation at given coordinates.
- */
+ /// evaluation at given coordinates.
const T& operator()(const DimVec<double>& bary) const {
- const WorldVector<double> *worldCoords = this->elInfo_->coordToWorld(bary, NULL);
- return (*fct_)(*worldCoords);
- };
+ WorldVector<double> worldCoords;
+ this->elInfo_->coordToWorld(bary, &worldCoords);
+ return (*fct_)(worldCoords);
+ }
protected:
- /** \brief
- * function to be avaluated at world coordinates.
- */
+ /// function to be avaluated at world coordinates.
const AbstractFunction<T, WorldVector<double> > *fct_;
};
diff --git a/AMDiS/src/Error.hh b/AMDiS/src/Error.hh
index 4179336f01bde4c78b378657dfae2b71b05d0b39..3f0a2b38ad7b3f6c2276013b131ca51644240183 100644
--- a/AMDiS/src/Error.hh
+++ b/AMDiS/src/Error.hh
@@ -8,23 +8,23 @@ namespace AMDiS {
T Error<T>::errUFct(const DimVec<double>& lambda)
{
WorldVector<double> x;
- elinfo->coordToWorld(lambda, &x);
- return((*pU)(x));
+ elinfo->coordToWorld(lambda, x);
+ return ((*pU)(x));
}
template<typename T>
WorldVector<T> Error<T>::grdErrUFct(const DimVec<double>& lambda)
{
WorldVector<double> x;
- elinfo->coordToWorld(lambda, &x);
- return((*pGrdU)(x));
+ elinfo->coordToWorld(lambda, x);
+ return ((*pGrdU)(x));
}
template<typename T>
int Error<T>::maxErrAtQpFct(ElInfo* el_info)
{
- double err;
+ double err = 0.0;
const double *u_vec, *uh_vec;
elinfo = el_info;
diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc
index 1e189e574bda176807fec5fc8bfb97a22fcbd909..5c7e0796aadf5710277ebd4d2d6c02caa8d1410e 100644
--- a/AMDiS/src/FirstOrderAssembler.cc
+++ b/AMDiS/src/FirstOrderAssembler.cc
@@ -81,28 +81,28 @@ namespace AMDiS {
}
// create new assembler
- // if (!optimized) {
+ if (!optimized) {
newAssembler =
(type == GRD_PSI) ?
dynamic_cast<FirstOrderAssembler*>(NEW Stand10(op, assembler, quad)) :
dynamic_cast<FirstOrderAssembler*>(NEW Stand01(op, assembler, quad));
-// } else {
-// if (pwConst) {
-// newAssembler =
-// (type == GRD_PSI) ?
-// dynamic_cast<FirstOrderAssembler*>(NEW Pre10(op, assembler, quad)) :
-// dynamic_cast<FirstOrderAssembler*>(NEW Pre01(op, assembler, quad));
-// } else {
-// newAssembler =
-// (type == GRD_PSI) ?
-// dynamic_cast<FirstOrderAssembler*>( NEW Quad10(op, assembler, quad)) :
-// dynamic_cast<FirstOrderAssembler*>( NEW Quad01(op, assembler, quad));
-// }
-// }
+ } else {
+ if (pwConst) {
+ newAssembler =
+ (type == GRD_PSI) ?
+ dynamic_cast<FirstOrderAssembler*>(NEW Pre10(op, assembler, quad)) :
+ dynamic_cast<FirstOrderAssembler*>(NEW Pre01(op, assembler, quad));
+ } else {
+ newAssembler =
+ (type == GRD_PSI) ?
+ dynamic_cast<FirstOrderAssembler*>( NEW Quad10(op, assembler, quad)) :
+ dynamic_cast<FirstOrderAssembler*>( NEW Quad01(op, assembler, quad));
+ }
+ }
subAssemblers->push_back(newAssembler);
return newAssembler;
- };
+ }
Stand10::Stand10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index 4ebf8f8b40441dd4af1bfb851045393974459738..39006e79546821bf443075532f688d5938df60d3 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -849,7 +849,7 @@ namespace AMDiS {
if (b_no[i] < Global::getGeo(VERTEX, dim)) {
rvec[i] = (*f)(el_info->getCoord(b_no[i]));
} else {
- el_info->coordToWorld(*(*bary)[b_no[i]], &x);
+ el_info->coordToWorld(*(*bary)[b_no[i]], x);
rvec[i] = (*f)(x);
}
}
@@ -858,7 +858,7 @@ namespace AMDiS {
for (int i = 0; i < vertices; i++)
rvec[i] = (*f)(el_info->getCoord(i));
for (int i = vertices; i < nBasFcts; i++) {
- el_info->coordToWorld(*(*bary)[i], &x);
+ el_info->coordToWorld(*(*bary)[i], x);
rvec[i] = (*f)(x);
}
}
@@ -894,7 +894,7 @@ namespace AMDiS {
if (b_no[i] < Global::getGeo(VERTEX, dim)) {
rvec[i] = (*f)(el_info->getCoord(b_no[i]));
} else {
- el_info->coordToWorld(*(*bary)[b_no[i]], &x);
+ el_info->coordToWorld(*(*bary)[b_no[i]], x);
rvec[i] = (*f)(x);
}
}
@@ -902,7 +902,7 @@ namespace AMDiS {
for (int i = 0; i < vertices; i++)
rvec[i] = (*f)(el_info->getCoord(i));
for (int i = vertices; i < nBasFcts; i++) {
- el_info->coordToWorld(*(*bary)[i], &x);
+ el_info->coordToWorld(*(*bary)[i], x);
rvec[i] = (*f)(x);
}
}
@@ -1018,7 +1018,7 @@ namespace AMDiS {
det = el_info->getDet();
for (iq = 0; iq < numPoints; iq++) {
- el_info->coordToWorld(q->getLambda(iq), &x);
+ el_info->coordToWorld(q->getLambda(iq), x);
wdetf_qp[iq] = q->getWeight(iq)*det*((*f)(x));
}
diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc
index 996f6b88e953071b6b19f11c3f2b641c669dc4e2..b01b5e4399d41e7b8e9cbc10ffba06b6fa3bcaa6 100644
--- a/AMDiS/src/MacroReader.cc
+++ b/AMDiS/src/MacroReader.cc
@@ -4,7 +4,7 @@
#include "Boundary.h"
#include "FiniteElemSpace.h"
#include "Mesh.h"
-#include <string>
+#include <string.h>
#include "FixVec.h"
#include "FixVecConvert.h"
#include "PeriodicMap.h"
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index a153f20ee73e14ae9deffbabf10406c784e52d5f..5f7031808a39eaddf51f48176dce9768ad4a6fdc 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -284,8 +284,6 @@ namespace AMDiS {
return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
}
-
-
void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
{
FUNCNAME("Mesh::addDOFAdmin()");
diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h
index 449e1365a02b73960d3eaa6eb217973cbf18e5e6..ad4d3e47fba5acfdfd4131ec51c138b94e39edf5 100644
--- a/AMDiS/src/Mesh.h
+++ b/AMDiS/src/Mesh.h
@@ -493,7 +493,7 @@ namespace AMDiS {
* Adds a DOFAdmin to the mesh
*/
virtual void addDOFAdmin(DOFAdmin *admin);
-
+
/** \brief
* Traverses the mesh. The argument level specifies the element level if
* CALL_EL_LEVEL or CALL_LEAF_EL_LEVEL, or the multigrid level if
diff --git a/AMDiS/src/MultiGridSolver.cc b/AMDiS/src/MultiGridSolver.cc
index 689f67c2942526690b369867229ae47a13180039..a58e7d6904337c9352e63d1536bc87ab5a79d590 100644
--- a/AMDiS/src/MultiGridSolver.cc
+++ b/AMDiS/src/MultiGridSolver.cc
@@ -156,7 +156,10 @@ namespace AMDiS {
std::vector<double*>::iterator factorIt;
std::vector<double*>::iterator factorBegin = systemMatrix_->getOperatorFactorBegin();
+
ElementMatrix *elementMatrix = NULL;
+ ERROR_EXIT("Initialize elementMatrix with a non NULL value!\n");
+
BoundaryType *bound = GET_MEMORY(BoundaryType, numDOFs);
TraverseStack stack;
@@ -191,7 +194,7 @@ namespace AMDiS {
level = elInfo->getLevel();
basFcts->getLocalIndices(element, admin, dofs);
basFcts->getBound(elInfo, bound);
- elementMatrix = (*opBegin)->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
+ (*opBegin)->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
bool levelUsed = isMGLevel_[level];
@@ -1075,7 +1078,7 @@ namespace AMDiS {
level = elInfo->getLevel();
basFcts->getLocalIndices(element, admin, dofs);
basFcts->getBound(elInfo, bound);
- elementMatrix = (*opBegin)->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
+ (*opBegin)->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
bool levelUsed = isMGLevel_[level];
if (levelUsed) {
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 365781618c8abf44d5b1e5d1a5394bf91b3e2e14..a50211dd457266b3a60f8ce2ea0d4e98fbe4a582 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -372,7 +372,6 @@ namespace AMDiS {
Quadrature *quad1GrdPhi,
Quadrature *quad0)
{
- if (!assembler[rank])
#ifdef _OPENMP
#pragma omp critical (initAssembler)
#endif
diff --git a/AMDiS/src/ParMetisPartitioner.cc b/AMDiS/src/ParMetisPartitioner.cc
index 5381ca913cb3e496f95a705b2a433f7154931a69..327238dc3c1ca172453d541a5aa507b1b4fb0517 100644
--- a/AMDiS/src/ParMetisPartitioner.cc
+++ b/AMDiS/src/ParMetisPartitioner.cc
@@ -114,7 +114,7 @@ namespace AMDiS {
// write xyz element coordinates
if (ptr_xyz) {
- elInfo->coordToWorld(bary, &world);
+ elInfo->coordToWorld(bary, world);
for (int i = 0; i < dim_; i++) {
*ptr_xyz = static_cast<float>(world[i]);
ptr_xyz++;
diff --git a/AMDiS/src/ParMetisPartitioner.h b/AMDiS/src/ParMetisPartitioner.h
index b5633721e1b89816e3ddc2c44eb3ce8bd34f261a..62cc9067be620701a12b8627165b3324cc297f4e 100644
--- a/AMDiS/src/ParMetisPartitioner.h
+++ b/AMDiS/src/ParMetisPartitioner.h
@@ -52,49 +52,49 @@ namespace AMDiS {
inline void setParMetisIndex(int amdisIndex, int parMetisIndex) {
elem_a2p_[amdisIndex] = parMetisIndex + 1;
- };
+ }
inline int getParMetisIndex(int amdisIndex) {
int result = elem_a2p_[amdisIndex];
TEST_EXIT(result > 0)("invalid index\n");
return result - 1;
- };
+ }
inline void setAMDiSIndex(int parMetisIndex, int amdisIndex) {
elem_p2a_[parMetisIndex] = amdisIndex;
- };
+ }
inline int getAMDiSIndex(int parMetisIndex) {
return elem_p2a_[parMetisIndex];
- };
+ }
inline int *getAMDiSIndices() {
return elem_p2a_;
- };
+ }
inline int *getElementPtr() {
return eptr_;
- };
+ }
inline int *getElementInd() {
return eind_;
- };
+ }
inline int *getElementDist() {
return elmdist_;
- };
+ }
inline int getDim() {
return dim_;
- };
+ }
inline float *getXYZ() {
return xyz_;
- };
+ }
inline int getNumElements() {
return numElements_;
- };
+ }
protected:
int *eptr_;
@@ -113,9 +113,7 @@ namespace AMDiS {
int *elem_p2a_;
- /** \brief
- * The MPI communicator that should be used for mesh partition.
- */
+ /// The MPI communicator that should be used for mesh partition.
MPI::Comm *mpiComm;
};
diff --git a/AMDiS/src/ParallelProblem.cc b/AMDiS/src/ParallelProblem.cc
index ba5a51bbf92a45acdf346b68542bd15fdbce4a22..fe3075bd598eb1939ddc118f4d1e509ca3fa53db 100644
--- a/AMDiS/src/ParallelProblem.cc
+++ b/AMDiS/src/ParallelProblem.cc
@@ -957,7 +957,7 @@ namespace AMDiS {
if (!visited[fineDOFs[i]]) {
visited[fineDOFs[i]] = true;
- elInfo2->coordToWorld(*(basFcts->getCoords(i)), &worldCoord);
+ elInfo2->coordToWorld(*(basFcts->getCoords(i)), worldCoord);
elInfo1->worldToCoord(worldCoord, &baryCoord);
// for all coarse vertex DOFs
diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc
index 10bb8438eafc4d0a951784ad05fba0ad237afe11..b9fb370f3790c0b462026472e578c76f4abb2725 100644
--- a/AMDiS/src/ProblemScal.cc
+++ b/AMDiS/src/ProblemScal.cc
@@ -27,8 +27,6 @@
namespace AMDiS {
- ProblemScal *ProblemScal::traversePtr_ = NULL;
-
ProblemScal::~ProblemScal()
{
FUNCNAME("ProblemScal::~ProblemScal()");
@@ -586,15 +584,12 @@ namespace AMDiS {
systemMatrix->getAssembleFlag() |
rhs->getAssembleFlag();
- if (useGetBound_)
+ if (useGetBound)
assembleFlag |= Mesh::FILL_BOUND;
systemMatrix->clear();
rhs->set(0.0);
- traversePtr_ = this;
- mesh->traverse(-1, assembleFlag, &buildAfterCoarsenFct);
-
// fill boundary conditions
if (systemMatrix->getBoundaryManager())
systemMatrix->getBoundaryManager()->initMatrix(systemMatrix);
@@ -603,17 +598,25 @@ namespace AMDiS {
if (solution->getBoundaryManager())
solution->getBoundaryManager()->initVector(solution);
+ BoundaryType *bound;
TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL |
- Mesh::FILL_BOUND |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA |
- Mesh::FILL_NEIGH);
-
- // for all elements ...
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
+
while (elInfo) {
+ if (useGetBound) {
+ bound = GET_MEMORY(BoundaryType, feSpace_->getBasisFcts()->getNumber());
+ feSpace_->getBasisFcts()->getBound(elInfo, bound);
+ } else {
+ bound = NULL;
+ }
+
+ systemMatrix->assemble(1.0, elInfo, bound);
+ rhs->assemble(1.0, elInfo, bound);
+
+ if (useGetBound) {
+ FREE_MEMORY(bound, BoundaryType, feSpace_->getBasisFcts()->getNumber());
+ }
+
if (systemMatrix->getBoundaryManager())
systemMatrix->getBoundaryManager()->fillBoundaryConditions(elInfo, systemMatrix);
if (rhs->getBoundaryManager())
@@ -634,29 +637,6 @@ namespace AMDiS {
INFO(info_, 8)("buildAfterCoarsen needed %.5f seconds\n",
TIME_USED(first,clock()));
-
- return;
- }
-
- int ProblemScal::buildAfterCoarsenFct(ElInfo *elInfo)
- {
- BoundaryType *bound;
-
- if (traversePtr_->getBoundUsed()) {
- bound = GET_MEMORY(BoundaryType, traversePtr_->getFESpace()->getBasisFcts()->getNumber());
- traversePtr_->getFESpace()->getBasisFcts()->getBound(elInfo, bound);
- } else {
- bound = NULL;
- }
-
- traversePtr_->getSystemMatrix()->assemble(1.0, elInfo, bound);
- traversePtr_->getRHS()->assemble(1.0, elInfo, bound);
-
- if (traversePtr_->getBoundUsed()) {
- FREE_MEMORY(bound, BoundaryType, traversePtr_->getFESpace()->getBasisFcts()->getNumber());
- }
-
- return 0;
}
void ProblemScal::writeResidualMesh(AdaptInfo *adaptInfo, const std::string name)
diff --git a/AMDiS/src/ProblemScal.h b/AMDiS/src/ProblemScal.h
index 55ea24b9cd8c5e597e9aee5b95f84d74fe0dd16a..3571119f5747dd2766fa8be3c478b5d40151f262 100644
--- a/AMDiS/src/ProblemScal.h
+++ b/AMDiS/src/ProblemScal.h
@@ -70,7 +70,7 @@ namespace AMDiS {
matVec(NULL),
leftPrecon(NULL),
rightPrecon(NULL),
- useGetBound_(true),
+ useGetBound(true),
refinementManager_(NULL),
coarseningManager_(NULL),
info_(10)
@@ -310,48 +310,37 @@ namespace AMDiS {
return refinementManager_;
}
- /** \brief
- * Returns \ref solver.
- */
+ /// Returns \ref solver.
inline OEMSolver<DOFVector<double> >* getSolver() {
return solver;
}
- /** \brief
- * Returns \ref marker_.
- */
+ /// Returns \ref marker_.
inline Marker *getMarker() {
return marker;
}
- /** \brief
- * Returns \ref name.
- */
+ /// Returns \ref name.
virtual inline const std::string& getName() {
return name;
}
- /** \brief
- * Returns \ref useGetBound_.
- */
+ /// Returns \ref useGetBound.
inline bool getBoundUsed() {
- return useGetBound_;
+ return useGetBound;
}
- /** \brief
- * Returns \ref leftPrecon.
- */
+ /// Returns \ref leftPrecon.
inline Preconditioner<DOFVector<double> > *getLeftPrecon() {
return leftPrecon;
}
- /** \brief
- * Returns \ref rightPrecon.
- */
+ /// Returns \ref rightPrecon.
inline Preconditioner<DOFVector<double> > *getRightPrecon() {
return rightPrecon;
}
+ ///
inline CoarseningManager *getCoarseningManager() {
return coarseningManager_;
}
@@ -429,79 +418,47 @@ namespace AMDiS {
return fileWriters;
}
- protected:
- /** \brief
- * Used by \ref buildAfterCoarsen().
- */
- static int buildAfterCoarsenFct(ElInfo *elInfo);
-
protected:
/// Name of this problem.
std::string name;
- /** \brief
- * FiniteElemSpace of this problem.
- */
+ /// FiniteElemSpace of this problem.
FiniteElemSpace *feSpace_;
- /** \brief
- * Mesh of this problem.
- */
+ /// Mesh of this problem.
Mesh *mesh;
- /** \brief
- * Responsible for element marking.
- */
+ /// Responsible for element marking.
Marker *marker;
- /** \brief
- * Estimator of this problem. Used in \ref estimate().
- */
+ /// Estimator of this problem. Used in \ref estimate().
Estimator *estimator;
- /** \brief
- * Linear solver of this problem. Used in \ref solve().
- */
+ /// Linear solver of this problem. Used in \ref solve().
OEMSolver<DOFVector<double> > *solver;
- /** \brief
- * DOFVector storing the calculated solution of the problem.
- */
+ /// DOFVector storing the calculated solution of the problem.
DOFVector<double> *solution;
- /** \brief
- * DOFVector for the right hand side
- */
+ /// DOFVector for the right hand side
DOFVector<double> *rhs;
- /** \brief
- * System matrix
- */
+ /// System matrix
DOFMatrix *systemMatrix;
- /** \brief
- * Matrix-vector multiplication
- */
+ /// Matrix-vector multiplication
MatVecMultiplier<DOFVector<double> > *matVec;
- /** \brief
- * Left preconditioner. Used in \ref solver.
- */
+ /// Left preconditioner. Used in \ref solver.
Preconditioner<DOFVector<double> > *leftPrecon;
- /** \brief
- * Right preconditioner. Used in \ref solver.
- */
+ /// Right preconditioner. Used in \ref solver.
Preconditioner<DOFVector<double> > *rightPrecon;
- /** \brief
- * Determines whether domain boundaries should be considered at assembling.
- */
- bool useGetBound_;
+ /// Determines whether domain boundaries should be considered at assembling.
+ bool useGetBound;
- /** \brief
- * Writes the meshes and solution after the adaption loop.
- */
+ /// Writes the meshes and solution after the adaption loop.
std::vector<FileWriterInterface*> fileWriters;
/** \brief
@@ -518,15 +475,8 @@ namespace AMDiS {
*/
CoarseningManager *coarseningManager_;
- /** \brief
- * Info level.
- */
+ /// Info level.
int info_;
-
- /** \brief
- * Used for mesh traversal.
- */
- static ProblemScal *traversePtr_;
};
}
diff --git a/AMDiS/src/QPInfo.cc b/AMDiS/src/QPInfo.cc
index f7b0d1e9d6c910735255fc67073c85453fafc47f..f330a62d9738a1ac16d02271534ba965c17c35ae 100644
--- a/AMDiS/src/QPInfo.cc
+++ b/AMDiS/src/QPInfo.cc
@@ -101,7 +101,7 @@ namespace AMDiS {
for(i = 0; i < numPoints; i++) {
const DimVec<double>& lambda = quadrature_->getLambda(i);
TEST_EXIT_DBG(currentElInfo_)("currentElInfo_ not set\n");
- currentElInfo_->coordToWorld(lambda, &coordsAtQPs_[i]);
+ currentElInfo_->coordToWorld(lambda, coordsAtQPs_[i]);
}
coordsNumPointsValid_ = numPoints;
}
diff --git a/AMDiS/src/Recovery.cc b/AMDiS/src/Recovery.cc
index b3f691723d4e58fb19d122b5303b599342d23350..bfe4a3e7b6358273951249b7fe66c55d36c89e4b 100644
--- a/AMDiS/src/Recovery.cc
+++ b/AMDiS/src/Recovery.cc
@@ -290,7 +290,7 @@ void Recovery::compute_integrals(DOFVector<double> *uh, ElInfo *elInfo,
for (k=0; k<n_points; k++)
{
- elInfo->coordToWorld(quad->getLambda(k), &quad_pts);
+ elInfo->coordToWorld(quad->getLambda(k), quad_pts);
sum += quad->getWeight(k) *
(*(*matrix_fcts)[i][j])(quad_pts, *rec_struct->coords);
}
@@ -322,7 +322,7 @@ void Recovery::compute_integrals(DOFVector<double> *uh, ElInfo *elInfo,
vec_sum = 0.0;
for (k=0; k<n_points; k++)
{
- elInfo->coordToWorld(quad->getLambda(k), &quad_pts);
+ elInfo->coordToWorld(quad->getLambda(k), quad_pts);
if (f_vec)
fAtQP = (*f_vec)(quad_pts);
if (f_scal)
@@ -343,7 +343,7 @@ void Recovery::compute_integrals(DOFVector<double> *uh, ElInfo *elInfo,
sum = 0.0;
for (k=0; k<n_points; k++)
{
- elInfo->coordToWorld(quad->getLambda(k), &quad_pts);
+ elInfo->coordToWorld(quad->getLambda(k), quad_pts);
sum += uhAtQP[k] * quad->getWeight(k)
* (*(*matrix_fcts)[0][i])(quad_pts, *rec_struct->coords);
}
@@ -385,7 +385,7 @@ void Recovery::compute_interior_sums(DOFVector<double> *uh, ElInfo *elInfo,
sum = 0.0;
for (k=0; k<n_points; k++)
{
- elInfo->coordToWorld(quad->getLambda(k), &quad_pts);
+ elInfo->coordToWorld(quad->getLambda(k), quad_pts);
sum += (*(*matrix_fcts)[i][j])(quad_pts, *rec_struct->coords);
}
(*(rec_struct->A))[i][j] += sum;
@@ -406,7 +406,7 @@ void Recovery::compute_interior_sums(DOFVector<double> *uh, ElInfo *elInfo,
vec_sum = 0.0;
for (k=0; k<n_points; k++)
{
- elInfo->coordToWorld(quad->getLambda(k), &quad_pts);
+ elInfo->coordToWorld(quad->getLambda(k), quad_pts);
if (f_vec)
fAtQP = (*f_vec)(quad_pts);
if (f_scal)
@@ -680,7 +680,7 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
{
// Setting world coordinates of node.
el_info->coordToWorld(*basis_fcts->getCoords(n_count),
- &coordinates);
+ coordinates);
(*struct_vec)[k].coords = NEW WorldVector<double>;
*(*struct_vec)[k].coords = coordinates;
@@ -712,7 +712,7 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
{
// Setting world coordinates of node.
el_info->coordToWorld(*basis_fcts->getCoords(n_count),
- &coordinates);
+ coordinates);
(*struct_vec)[k].coords = NEW WorldVector<double>;
*(*struct_vec)[k].coords = coordinates;
@@ -741,7 +741,7 @@ void Recovery::fill_struct_vec(DOFVector<double> *uh,
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);
+ el_info->coordToWorld(*basis_fcts->getCoords(n_count), coordinates);
(*struct_vec)[k].coords = NEW WorldVector<double>;
*(*struct_vec)[k].coords = coordinates;
@@ -1023,7 +1023,7 @@ Recovery::recovery(DOFVector<double> *uh,
double fAtBary = 1.0;
if (f_vec) {
- elInfo->coordToWorld(bary, &barycenter);
+ elInfo->coordToWorld(bary, barycenter);
fAtBary = (*f_vec)(barycenter);
}
diff --git a/AMDiS/src/RecoveryEstimator.cc b/AMDiS/src/RecoveryEstimator.cc
index 81d5b535704041875a849fc18b299a8a04ff3dad..32c8d58bce59f22ea7ef80c324a59e261d076dae 100644
--- a/AMDiS/src/RecoveryEstimator.cc
+++ b/AMDiS/src/RecoveryEstimator.cc
@@ -146,7 +146,7 @@ double RecoveryEstimator::estimate(double timestep)
fAtQP = (*f_scal)(uhAtQP[i]);
if (f_vec)
{
- elInfo->coordToWorld(quad->getLambda(i), &quad_pt);
+ elInfo->coordToWorld(quad->getLambda(i), quad_pt);
fAtQP = (*f_vec)(quad_pt);
}
diff --git a/AMDiS/src/SMIAdapter.cc b/AMDiS/src/SMIAdapter.cc
index b1b15817843ffe47afa55bdd4921dc5e26267631..1c04881659af1071f084846984722a4911203cc7 100644
--- a/AMDiS/src/SMIAdapter.cc
+++ b/AMDiS/src/SMIAdapter.cc
@@ -194,7 +194,7 @@ namespace AMDiS {
int elementIndex;
DimVec<double> *bary = NULL;
- const WorldVector<double> *world = NULL;
+ WorldVector<double> world;
const DegreeOfFreedom *elementDofs = NULL;
DegreeOfFreedom dof;
@@ -259,9 +259,9 @@ namespace AMDiS {
if(!alreadyAdded[dof]) {
newNodeIndices[numNewNodes] = dof;
- world = elInfo->coordToWorld(*bary, NULL);
+ elInfo->coordToWorld(*bary, world);
for(j = 0; j < dim; j++) {
- nodeCoords[numNewNodes * dim + j] = (*world)[j];
+ nodeCoords[numNewNodes * dim + j] = world[j];
}
numNewNodes++;
alreadyAdded[dof] = 1;
@@ -312,9 +312,9 @@ namespace AMDiS {
if(!alreadyAdded[dof]) {
newNodeIndices[numNewNodes] = dof;
- world = elInfo->coordToWorld(*bary, NULL);
+ elInfo->coordToWorld(*bary, world);
for(j = 0; j < dim; j++) {
- nodeCoords[numNewNodes * dim + j] = (*world)[j];
+ nodeCoords[numNewNodes * dim + j] = world[j];
}
numNewNodes++;
alreadyAdded[dof] = 1;
diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc
index c1385c876a502505611ed1d9a46c0c0d55f251f5..8a023e9a6b5b221d53542e5b858fbf7104f14485 100644
--- a/AMDiS/src/SecondOrderAssembler.cc
+++ b/AMDiS/src/SecondOrderAssembler.cc
@@ -68,15 +68,15 @@ namespace AMDiS {
}
// create new assembler
- // if (!optimized) {
+ if (!optimized) {
newAssembler = NEW Stand2(op, assembler, quad);
-// } else {
-// if (pwConst) {
-// newAssembler = NEW Pre2(op, assembler, quad);
-// } else {
-// newAssembler = NEW Quad2(op, assembler, quad);
-// }
-// }
+ } else {
+ if (pwConst) {
+ newAssembler = NEW Pre2(op, assembler, quad);
+ } else {
+ newAssembler = NEW Quad2(op, assembler, quad);
+ }
+ }
subAssemblers->push_back(newAssembler);
return newAssembler;
diff --git a/AMDiS/src/SubAssembler.cc b/AMDiS/src/SubAssembler.cc
index d710c970c0d984be92f02867a02154df6a1d2480..0a40ba8ff811162dc5e8ece11213d5da98dd452c 100644
--- a/AMDiS/src/SubAssembler.cc
+++ b/AMDiS/src/SubAssembler.cc
@@ -135,7 +135,7 @@ namespace AMDiS {
// set new values
WorldVector<double>* k = &(coordsAtQPs[0]);
for (int l = 0; k < &(coordsAtQPs[nPoints]); ++k, ++l) {
- elInfo->coordToWorld(localQuad->getLambda(l), k);
+ elInfo->coordToWorld(localQuad->getLambda(l), *k);
}
// mark values as valid
diff --git a/AMDiS/src/SubElInfo.cc b/AMDiS/src/SubElInfo.cc
index a3fead77baf13a61ad662458dc6e8be1e0ee9ad9..9e7745f23eae5f6484c7eb5fae0e8e409dbc035f 100644
--- a/AMDiS/src/SubElInfo.cc
+++ b/AMDiS/src/SubElInfo.cc
@@ -29,7 +29,7 @@ namespace AMDiS {
*/
for (int i = 0; i < nPoints; i++) {
elInfo->coordToWorld((const DimVec<double>)((*lambda)[i]),
- &(worldCoords[i]));
+ worldCoords[i]);
}
det = ElInfo::calcDet(worldCoords);
diff --git a/AMDiS/src/SurfaceOperator.h b/AMDiS/src/SurfaceOperator.h
index f9322fd578dea93db515ba6ccf0f1f3dc292b16c..349d5efc217cd04b9f807235ed8837b30b972de6 100644
--- a/AMDiS/src/SurfaceOperator.h
+++ b/AMDiS/src/SurfaceOperator.h
@@ -145,7 +145,7 @@ namespace AMDiS {
// transform barycentric coords to world coords
for(i = 0; i < dim; i++) {
- elInfo->coordToWorld(coords_[i], &worldCoords[i]);
+ elInfo->coordToWorld(coords_[i], worldCoords[i]);
}
// set determinant for world coords of the side
@@ -175,7 +175,7 @@ namespace AMDiS {
// transform barycentric coords to world coords
for(i = 0; i < dim; i++) {
- elInfo->coordToWorld(coords_[i], &worldCoords[i]);
+ elInfo->coordToWorld(coords_[i], worldCoords[i]);
}
// set determinant for world coords of the side
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index 667b640e67f53261f05e69435e29f6811c55472e..d5b5f68b28b13964ef0f4b16df857bd28831c9d4 100644
--- a/AMDiS/src/Traverse.cc
+++ b/AMDiS/src/Traverse.cc
@@ -111,7 +111,6 @@ namespace AMDiS {
/* - all tree elements in pre-/in-/post-order */
/* depending on the traverse_info->level variable */
/****************************************************************************/
-
int Traverse::recursive(ElInfoStack *elInfoStack)
{
FUNCNAME("Traverse::recursive()");
@@ -363,10 +362,6 @@ namespace AMDiS {
((info_stack[stack_used] >= 2) || (el->getFirstChild() == NULL))) {
stack_used--;
el = elinfo_stack[stack_used]->getElement();
-
- if (calcDisplacement) {
- displacementStack.pop();
- }
}
/* goto next macro element */
@@ -403,15 +398,6 @@ namespace AMDiS {
elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
- if (calcDisplacement) {
- if (i == 0) {
- displacementStack.push(displacementStack.top() * 2);
- } else {
- displacementStack.push(displacementStack.top() * 2 + 1);
- }
- }
-
-
TEST_EXIT_DBG(stack_used < stack_size)
("stack_size=%d too small, level=(%d,%d)\n",
stack_size, elinfo_stack[stack_used]->getLevel());
@@ -419,10 +405,6 @@ namespace AMDiS {
info_stack[stack_used] = 0;
}
- if (calcDisplacement) {
- elinfo_stack[stack_used]->setDisplacement(displacementStack.top());
- }
-
return elinfo_stack[stack_used];
}
@@ -888,7 +870,7 @@ namespace AMDiS {
{
FUNCNAME("TraverseStack::traverseNeighbour2d()");
- Triangle *el, *el2, *sav_el;
+ Element *el, *el2, *sav_el;
ElInfo *old_elinfo, *elinfo, *elinfo2;
int i, nb, opp_vertex, stack2_used;
@@ -908,7 +890,7 @@ 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 = elinfo_stack[stack_used]->getElement();
sav_index = el->getIndex();
sav_el = el;
@@ -965,7 +947,7 @@ namespace AMDiS {
stack2_used = 1;
}
elinfo2 = save_elinfo_stack[stack2_used];
- el2 = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo2->getElement()));
+ el2 = elinfo2->getElement();
i = traverse_mel->getOppVertex(nb);
traverse_mel = traverse_mel->getNeighbour(nb);
@@ -984,17 +966,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 = 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;
+ elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
- nb = 1-i;
+ nb = 1 - i;
}
/****************************************************************************/
@@ -1004,15 +986,15 @@ namespace AMDiS {
/****************************************************************************/
elinfo = elinfo_stack[stack_used];
- el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
+ el = 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;
+ elinfo_stack[stack_used + 1]->fillElInfo(1 - nb, elinfo_stack[stack_used]);
+ info_stack[stack_used] = 2 - nb;
stack_used++;
nb = 2;
}
@@ -1023,25 +1005,24 @@ 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;
+ elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
nb = i;
elinfo = elinfo_stack[stack_used];
- el = dynamic_cast<Triangle*>(const_cast<Element*>( elinfo->getElement()));
+ el = elinfo->getElement();
stack2_used++;
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 = 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 = elinfo->getElement();
}
}
@@ -1077,7 +1058,7 @@ namespace AMDiS {
elinfo->fillDetGrdLambda();
}
- return(elinfo);
+ return elinfo;
}
void TraverseStack::update()
@@ -1103,17 +1084,4 @@ namespace AMDiS {
elinfo_stack[stack_used - levelDif + i]->getIChild());
}
}
-
- void TraverseStack::startDisplacementCalculation(int level)
- {
- calcDisplacement = true;
-
- displacementStack.empty();
- for (int i = 0; i <= elinfo_stack[stack_used]->getLevel() - level; i++) {
- displacementStack.push(0);
- }
-
- elinfo_stack[stack_used]->setDisplacement(0);
- }
-
}
diff --git a/AMDiS/src/Traverse.h b/AMDiS/src/Traverse.h
index 63dc4403e787e5376a0c86b2735d0caa09f10a74..7c15c2aa2d87f380f7640e00aa24197204f8440a 100644
--- a/AMDiS/src/Traverse.h
+++ b/AMDiS/src/Traverse.h
@@ -71,7 +71,6 @@ namespace AMDiS {
stack_size(0),
stack_used(0),
save_stack_used(0),
- calcDisplacement(false),
myThreadId_(0),
maxThreads_(1)
{
@@ -152,18 +151,6 @@ namespace AMDiS {
void getCoordsInElem(const ElInfo *upperElInfo,
DimMat<double> *coords);
- /** \brief
- * Starts the calculation of the displacement relative to given level.
- */
- void startDisplacementCalculation(int level);
-
- /** \brief
- * Stops the calculation of displacement information.
- */
- void stopDisplacementCalculation() {
- calcDisplacement = false;
- }
-
/** \brief
* Is used for parallel mesh traverse.
*/
@@ -225,14 +212,10 @@ namespace AMDiS {
}
private:
- /** \brief
- * Iterator to the current MacroElement
- */
+ /// Iterator to the current MacroElement
std::deque<MacroElement*>::const_iterator currentMacro;
- /** \brief
- * Mesh which is currently traversed
- */
+ /// Mesh which is currently traversed
Mesh* traverse_mesh;
/** \brief
@@ -247,24 +230,16 @@ namespace AMDiS {
*/
Flag traverse_fill_flag;
- /** \brief
- * current macro element
- */
+ /// current macro element
const MacroElement *traverse_mel;
- /** \brief
- * Current size of the stack
- */
+ /// Current size of the stack
int stack_size;
- /** \brief
- * Used size of the stack
- */
+ /// Used size of the stack
int stack_used;
- /** \brief
- *
- */
+ ///
std::vector<ElInfo*> elinfo_stack;
/** \brief
@@ -275,41 +250,21 @@ namespace AMDiS {
* visited.
*/
std::vector<int> info_stack;
-
- /** \brief
- *
- */
+
+ ///
const MacroElement *save_traverse_mel;
- /** \brief
- *
- */
+ ///
std::vector<ElInfo*> save_elinfo_stack;
- /** \brief
- *
- */
+ ///
std::vector<unsigned char> save_info_stack;
- /** \brief
- *
- */
+ ///
int save_stack_used;
-
- /** \brief
- *
- */
- int id;
- /** \brief
- * Stack for counting the deplacement information for dual traverse.
- */
- std::stack<int> displacementStack;
-
- /** \brief
- * True, if \ref displacementStack should be calculated.
- */
- bool calcDisplacement;
+ ///
+ int id;
/** \brief
* Is used for parallel mesh traverse. The thread with the id
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index 4b92dfe0e4f2be3a6a6ef0d53f6dbab1fe53c081..8fdb4605e68ea79872beab3abaa11dc147c5bdfd 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -66,15 +66,15 @@ namespace AMDiS {
}
// create new assembler
- // if (!optimized) {
- newAssembler = NEW StandardZOA(op, assembler, quad);
-// } else {
-// if (pwConst) {
-// newAssembler = NEW PrecalcZOA(op, assembler, quad);
-// } else {
-// newAssembler = NEW FastQuadZOA(op, assembler, quad);
-// }
-// }
+ if (!optimized) {
+ newAssembler = NEW StandardZOA(op, assembler, quad);
+ } else {
+ if (pwConst) {
+ newAssembler = NEW PrecalcZOA(op, assembler, quad);
+ } else {
+ newAssembler = NEW FastQuadZOA(op, assembler, quad);
+ }
+ }
subAssemblers->push_back(newAssembler);
return newAssembler;
@@ -360,11 +360,15 @@ namespace AMDiS {
void FastQuadZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
+ int myRank = omp_get_thread_num();
+ int nPoints = quadrature->getNumPoints();
+
if (firstCall) {
#ifdef _OPENMP
#pragma omp critical
#endif
{
+ cPtrs[myRank] = GET_MEMORY(double, nPoints);
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
@@ -373,14 +377,11 @@ namespace AMDiS {
}
}
- int nPoints = quadrature->getNumPoints();
- double *c = GET_MEMORY(double, nPoints);
-
+ double *c = cPtrs[myRank];
for (int iq = 0; iq < nPoints; iq++) {
c[iq] = 0.0;
}
- int myRank = omp_get_thread_num();
std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
@@ -394,7 +395,6 @@ namespace AMDiS {
(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
}
}
- FREE_MEMORY(c, double, nPoints);
}
PrecalcZOA::PrecalcZOA(Operator *op, Assembler *assembler, Quadrature *quad)