Commit 44677772 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* small update to increase performance of assembling procedure

parent dd8ffdf7
...@@ -279,12 +279,12 @@ namespace AMDiS { ...@@ -279,12 +279,12 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> > &surfVert) VectorOfFixVecs<DimVec<double> > &surfVert)
{ {
double surfDet; double surfDet;
int dim = surfVert[0].getSize()-1; int dim = surfVert[0].getSize() - 1;
FixVec<WorldVector<double>, VERTEX> worldCoords(dim-1, NO_INIT); FixVec<WorldVector<double>, VERTEX> worldCoords(dim - 1, NO_INIT);
// transform barycentric coords to world coords // transform barycentric coords to world coords
for(int i = 0; i < dim; i++) { for (int i = 0; i < dim; i++) {
loc_elInfo->coordToWorld(surfVert[i], &worldCoords[i]); loc_elInfo->coordToWorld(surfVert[i], worldCoords[i]);
} }
// calculate determinant for surface // calculate determinant for surface
......
...@@ -17,11 +17,11 @@ namespace AMDiS { ...@@ -17,11 +17,11 @@ namespace AMDiS {
const double &fac) const double &fac)
{ {
double val = 0.0; double val = 0.0;
const WorldVector<double> *worldCoordsAtQP; WorldVector<double> worldCoordsAtQP;
for (int iq = 0; iq < nQPts; ++iq) { for (int iq = 0; iq < nQPts; ++iq) {
worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL); elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
val += q->getWeight(iq) * fabs((*f)(*worldCoordsAtQP)); val += q->getWeight(iq) * fabs((*f)(worldCoordsAtQP));
} }
double nrm = det * val; double nrm = det * val;
...@@ -34,11 +34,11 @@ namespace AMDiS { ...@@ -34,11 +34,11 @@ namespace AMDiS {
const double &fac) const double &fac)
{ {
double val = 0.0; double val = 0.0;
const WorldVector<double> *worldCoordsAtQP; WorldVector<double> worldCoordsAtQP;
for (int iq = 0; iq < nQPts; ++iq) { for (int iq = 0; iq < nQPts; ++iq) {
worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL); elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
val += q->getWeight(iq) * sqr((*f)(*worldCoordsAtQP)); val += q->getWeight(iq) * sqr((*f)(worldCoordsAtQP));
} }
double nrm = det * val; double nrm = det * val;
...@@ -51,15 +51,15 @@ namespace AMDiS { ...@@ -51,15 +51,15 @@ namespace AMDiS {
const double &fac) const double &fac)
{ {
double val = 0.0; double val = 0.0;
double norm_grd2; double norm_grd2 = 0.0;
const WorldVector<double> *worldCoordsAtQP; WorldVector<double> worldCoordsAtQP;
for (int iq = 0; iq < nQPts; ++iq) { for (int iq = 0; iq < nQPts; ++iq) {
worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL); elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
norm_grd2 = 0.0; norm_grd2 = 0.0;
for (int j = 0; j < dim; j++) 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; val += q->getWeight(iq) * norm_grd2;
} }
...@@ -142,14 +142,14 @@ namespace AMDiS { ...@@ -142,14 +142,14 @@ namespace AMDiS {
q, q,
NULL, NULL,
NULL); NULL);
const WorldVector<double> *worldCoordsAtQP; WorldVector<double> worldCoordsAtQP;
for (int iq = 0; iq < nQPts; ++iq) { for (int iq = 0; iq < nQPts; ++iq) {
worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL); elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
val += q->getWeight(iq) * sqr((*u)(*worldCoordsAtQP) - uhAtQPs[iq]); val += q->getWeight(iq) * sqr((*u)(worldCoordsAtQP) - uhAtQPs[iq]);
if (relErr) if (relErr)
val_nrm += q->getWeight(iq) * sqr((*u)(*worldCoordsAtQP)); val_nrm += q->getWeight(iq) * sqr((*u)(worldCoordsAtQP));
} }
double err = det * val; double err = det * val;
...@@ -173,23 +173,22 @@ namespace AMDiS { ...@@ -173,23 +173,22 @@ namespace AMDiS {
q, q,
NULL, NULL,
NULL); NULL);
const WorldVector<double> *worldCoordsAtQP; WorldVector<double> worldCoordsAtQP;
for (int iq = 0; iq < nQPts; ++iq) { for (int iq = 0; iq < nQPts; ++iq) {
elInfo->coordToWorld(q->getLambda(iq), worldCoordsAtQP);
worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
norm_err_grd2 = 0.0; norm_err_grd2 = 0.0;
for (int j = 0; j < dim; ++j) for (int j = 0; j < dim; ++j)
norm_err_grd2 += norm_err_grd2 +=
sqr(((*grdu)(*worldCoordsAtQP))[j] - grdUhAtQPs[iq][j]); sqr(((*grdu)(worldCoordsAtQP))[j] - grdUhAtQPs[iq][j]);
val += q->getWeight(iq) * norm_err_grd2; val += q->getWeight(iq) * norm_err_grd2;
if (relErr) { if (relErr) {
norm_grd2 = 0.0; norm_grd2 = 0.0;
for (int j = 0; j < dim; ++j) 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; val_nrm += q->getWeight(iq) * norm_grd2;
} }
......
...@@ -249,8 +249,8 @@ ElementLevelSet::calcIntersecNormal_2d(WorldVector<double> &normalVec) ...@@ -249,8 +249,8 @@ ElementLevelSet::calcIntersecNormal_2d(WorldVector<double> &normalVec)
// ===== Get world coordinates of intersection points. ===== // ===== Get world coordinates of intersection points. =====
WorldVector<double> sP1; WorldVector<double> sP1;
WorldVector<double> sP2; WorldVector<double> sP2;
elInfo->coordToWorld((*elIntersecPoints)[0], &sP1); elInfo->coordToWorld((*elIntersecPoints)[0], sP1);
elInfo->coordToWorld((*elIntersecPoints)[1], &sP2); elInfo->coordToWorld((*elIntersecPoints)[1], sP2);
// ===== Calculate normal vector. ===== // ===== Calculate normal vector. =====
double norm2 = 0.0; double norm2 = 0.0;
...@@ -326,9 +326,9 @@ ElementLevelSet::calcIntersecNormal_3d(WorldVector<double> &normalVec) ...@@ -326,9 +326,9 @@ ElementLevelSet::calcIntersecNormal_3d(WorldVector<double> &normalVec)
WorldVector<double> sP1; WorldVector<double> sP1;
WorldVector<double> sP2; WorldVector<double> sP2;
WorldVector<double> sP3; WorldVector<double> sP3;
elInfo->coordToWorld((*elIntersecPoints)[0], &sP1); elInfo->coordToWorld((*elIntersecPoints)[0], sP1);
elInfo->coordToWorld((*elIntersecPoints)[1], &sP2); elInfo->coordToWorld((*elIntersecPoints)[1], sP2);
elInfo->coordToWorld((*elIntersecPoints)[2], &sP3); elInfo->coordToWorld((*elIntersecPoints)[2], sP3);
// ===== Calculate normal vector. ===== // ===== Calculate normal vector. =====
WorldVector<double> A = sP2-sP1; WorldVector<double> A = sP2-sP1;
......
...@@ -406,13 +406,11 @@ namespace AMDiS { ...@@ -406,13 +406,11 @@ namespace AMDiS {
ZeroOrderAssembler::getSubAssembler(op, this, quad0, false); ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);
} }
ElementMatrix *Assembler::initElementMatrix(ElementMatrix *elMat, void Assembler::initElementMatrix(ElementMatrix *elMat,
const ElInfo *rowElInfo, const ElInfo *rowElInfo,
const ElInfo *colElInfo) const ElInfo *colElInfo)
{ {
if (!elMat) { TEST_EXIT_DBG(elMat)("No ElementMatrix!\n");
elMat = NEW ElementMatrix(nRow, nCol);
}
elMat->set(0.0); elMat->set(0.0);
...@@ -434,25 +432,17 @@ namespace AMDiS { ...@@ -434,25 +432,17 @@ namespace AMDiS {
&(elMat->colIndices)); &(elMat->colIndices));
} }
} }
return elMat;
} }
ElementVector *Assembler::initElementVector(ElementVector *elVec, void Assembler::initElementVector(ElementVector *elVec, const ElInfo *elInfo)
const ElInfo *elInfo)
{ {
if (!elVec) { TEST_EXIT_DBG(elVec)("No ElementVector!\n");
elVec = NEW ElementVector(nRow);
}
elVec->set(0.0); elVec->set(0.0);
DOFAdmin *rowAdmin = rowFESpace->getAdmin(); rowFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
Element *element = elInfo->getElement(); rowFESpace->getAdmin(),
&(elVec->dofIndices));
rowFESpace->getBasisFcts()->getLocalIndicesVec(element, rowAdmin, &(elVec->dofIndices));
return elVec;
} }
void Assembler::checkQuadratures() void Assembler::checkQuadratures()
......
...@@ -69,13 +69,13 @@ namespace AMDiS { ...@@ -69,13 +69,13 @@ namespace AMDiS {
virtual ~Assembler() {} virtual ~Assembler() {}
ElementMatrix *initElementMatrix(ElementMatrix *elMat, void initElementMatrix(ElementMatrix *elMat,
const ElInfo *rowElInfo, const ElInfo *rowElInfo,
const ElInfo *colElInfo = NULL); const ElInfo *colElInfo = NULL);
ElementVector *initElementVector(ElementVector *elVec, void initElementVector(ElementVector *elVec,
const ElInfo *elInfo); const ElInfo *elInfo);
/** \brief /** \brief
* Assembles the element matrix for the given ElInfo * Assembles the element matrix for the given ElInfo
...@@ -230,9 +230,7 @@ namespace AMDiS { ...@@ -230,9 +230,7 @@ namespace AMDiS {
const ElInfo *smallElInfo, const ElInfo *largeElInfo, const ElInfo *smallElInfo, const ElInfo *largeElInfo,
ElementVector *vec); ElementVector *vec);
/** \brief /// Checks whether this is a new travese.
* Checks whether this is a new travese.
*/
inline void checkForNewTraverse() { inline void checkForNewTraverse() {
if (lastTraverseId != ElInfo::traverseId[omp_get_thread_num()]) { if (lastTraverseId != ElInfo::traverseId[omp_get_thread_num()]) {
lastVecEl = lastMatEl = NULL; lastVecEl = lastMatEl = NULL;
......
...@@ -12,6 +12,7 @@ namespace AMDiS { ...@@ -12,6 +12,7 @@ namespace AMDiS {
BoundaryManager::BoundaryManager(const FiniteElemSpace *feSpace) BoundaryManager::BoundaryManager(const FiniteElemSpace *feSpace)
{ {
localBounds.resize(omp_get_overall_max_threads()); localBounds.resize(omp_get_overall_max_threads());
dofIndices.resize(omp_get_overall_max_threads());
allocatedMemoryLocalBounds = feSpace->getBasisFcts()->getNumber(); allocatedMemoryLocalBounds = feSpace->getBasisFcts()->getNumber();
for (int i = 0; i < static_cast<int>(localBounds.size()); i++) { for (int i = 0; i < static_cast<int>(localBounds.size()); i++) {
localBounds[i] = GET_MEMORY(BoundaryType, allocatedMemoryLocalBounds); localBounds[i] = GET_MEMORY(BoundaryType, allocatedMemoryLocalBounds);
...@@ -23,6 +24,7 @@ namespace AMDiS { ...@@ -23,6 +24,7 @@ namespace AMDiS {
localBCs = bm.localBCs; localBCs = bm.localBCs;
allocatedMemoryLocalBounds = bm.allocatedMemoryLocalBounds; allocatedMemoryLocalBounds = bm.allocatedMemoryLocalBounds;
localBounds.resize(bm.localBounds.size()); localBounds.resize(bm.localBounds.size());
dofIndices.resize(bm.localBounds.size());
for (int i = 0; i < static_cast<int>(localBounds.size()); i++) { for (int i = 0; i < static_cast<int>(localBounds.size()); i++) {
localBounds[i] = GET_MEMORY(BoundaryType, allocatedMemoryLocalBounds); localBounds[i] = GET_MEMORY(BoundaryType, allocatedMemoryLocalBounds);
} }
...@@ -51,29 +53,26 @@ namespace AMDiS { ...@@ -51,29 +53,26 @@ namespace AMDiS {
void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo, void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo,
DOFVectorBase<double> *vec) 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) { 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 // get boundaries of all DOFs
BoundaryType *localBound = localBounds[omp_get_thread_num()]; BoundaryType *localBound = localBounds[omp_get_thread_num()];
basisFcts->getBound(elInfo, localBound); basisFcts->getBound(elInfo, localBound);
// get dof indices // get dof indices
basisFcts->getLocalIndicesVec(elInfo->getElement(), admin, &dofIndices); basisFcts->getLocalIndicesVec(elInfo->getElement(),
feSpace->getAdmin(), &dofVec);
// apply non dirichlet boundary conditions // apply non dirichlet boundary conditions
for (it = localBCs.begin(); it != localBCs.end(); ++it) { for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) { if ((*it).second) {
if (!(*it).second->isDirichlet()) { 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 { ...@@ -82,7 +81,7 @@ namespace AMDiS {
for (it = localBCs.begin(); it != localBCs.end(); ++it) { for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) { if ((*it).second) {
if ((*it).second->isDirichlet()) { 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 { ...@@ -92,28 +91,26 @@ namespace AMDiS {
void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo, void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo,
DOFMatrix *mat) 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) { 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 // get boundaries of all DOFs
BoundaryType *localBound = localBounds[omp_get_thread_num()]; BoundaryType *localBound = localBounds[omp_get_thread_num()];
basisFcts->getBound(elInfo, localBound); basisFcts->getBound(elInfo, localBound);
// get dof indices // get dof indices
basisFcts->getLocalIndicesVec(elInfo->getElement(), admin, &dofIndices); basisFcts->getLocalIndicesVec(elInfo->getElement(),
feSpace->getAdmin(), &dofVec);
// apply non dirichlet boundary conditions // apply non dirichlet boundary conditions
for (it = localBCs.begin(); it != localBCs.end(); ++it) { for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) { if ((*it).second) {
if (!(*it).second->isDirichlet()) { 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 { ...@@ -122,7 +119,7 @@ namespace AMDiS {
for (it = localBCs.begin(); it != localBCs.end(); ++it) { for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) { if ((*it).second) {
if ((*it).second->isDirichlet()) { if ((*it).second->isDirichlet()) {
(*it).second->fillBoundaryCondition(mat, elInfo, &dofIndices[0], localBound, nBasFcts); (*it).second->fillBoundaryCondition(mat, elInfo, &dofVec[0], localBound, nBasFcts);
} }
} }
} }
......
...@@ -31,6 +31,7 @@ namespace AMDiS { ...@@ -31,6 +31,7 @@ namespace AMDiS {
class DOFMatrix; class DOFMatrix;
class FiniteElemSpace; class FiniteElemSpace;
template<typename T> class Vector;
template<typename T> class DOFVectorBase; template<typename T> class DOFVectorBase;
// ============================================================================ // ============================================================================
...@@ -107,16 +108,15 @@ namespace AMDiS { ...@@ -107,16 +108,15 @@ namespace AMDiS {
}; };
protected: protected:
/** \brief /// Map of managed local boundary conditions.
* Map of managed local boundary conditions.
*/
std::map<BoundaryType, BoundaryCondition*> localBCs; std::map<BoundaryType, BoundaryCondition*> localBCs;
/** \brief /// Temporary thread-safe variable for functions fillBoundaryconditions.
* Temporary thread-safe variable for functions fillBoundaryconditions.
*/
std::vector<BoundaryType*> localBounds; std::vector<BoundaryType*> localBounds;
/// Temporary thread-safe variable for functions fillBoundaryconditions.
std::vector<Vector<DegreeOfFreedom> > dofIndices;
/** \brief /** \brief
* Stores the number of byte that were allocated in the constructor for * Stores the number of byte that were allocated in the constructor for
* each localBounds value. Is used to free the memory in the destructor. * each localBounds value. Is used to free the memory in the destructor.
......
...@@ -249,10 +249,6 @@ namespace AMDiS { ...@@ -249,10 +249,6 @@ namespace AMDiS {
int nRow = elMat.rowIndices.getSize(); int nRow = elMat.rowIndices.getSize();
int nCol = elMat.colIndices.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 for (int i = 0; i < nRow; i++) { // for all rows of element matrix
row = elMat.rowIndices[i]; row = elMat.rowIndices[i];
BoundaryCondition *condition = BoundaryCondition *condition =
...@@ -433,37 +429,40 @@ namespace AMDiS { ...@@ -433,37 +429,40 @@ namespace AMDiS {
} }
void DOFMatrix::assemble(double factor, ElInfo *elInfo, void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
const BoundaryType *bound, Operator *op)
{ {
FUNCNAME("DOFMatrix::assemble()"); FUNCNAME("DOFMatrix::assemble()");
if (!op && operators.size() == 0) { if (operators.size() == 0)
return; return;
}
Operator *operat = op ? op : operators[0]; operators[0]->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
operat->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
if (op) { std::vector<Operator*>::iterator it = operators.begin();
op->getElementMatrix(elInfo, elementMatrix); std::vector<double*>::iterator factorIt = operatorFactor.begin();
} else { for (; it != operators.end(); ++it, ++factorIt) {
std::vector<Operator*>::iterator it;
std::vector<double*>::iterator factorIt;
for (it = operators.begin(), factorIt = operatorFactor.begin();
it != operators.end();
++it, ++factorIt) {
if ((*it)->getNeedDualTraverse() == false) { if ((*it)->getNeedDualTraverse() == false) {
(*it)->getElementMatrix(elInfo, (*it)->getElementMatrix(elInfo,
elementMatrix, elementMatrix,
*factorIt ? **factorIt : 1.0); *factorIt ? **factorIt : 1.0);
} }
} }
}
addElementMatrix(factor, *elementMatrix, bound); 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, void DOFMatrix::assemble(double factor,
ElInfo *rowElInfo, ElInfo *colElInfo, ElInfo *rowElInfo, ElInfo *colElInfo,
ElInfo *smallElInfo, ElInfo *largeElInfo, ElInfo *smallElInfo, ElInfo *largeElInfo,
......
...@@ -357,9 +357,10 @@ namespace AMDiS { ...@@ -357,9 +357,10 @@ namespace AMDiS {
* matrix should be cleared by calling clear dof matrix(). * matrix should be cleared by calling clear dof matrix().
*/ */