Skip to content

GitLab

Commit 23e3a357 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files
Options

* Faster assembling for 2D-vecphase

parent 6eda4c4f
Hide whitespace changes
Inline Side-by-side
Showing with 311 additions and 300 deletions
AMDiS/src/Assembler.cc
View file @ 23e3a357
......@@ -498,7 +498,7 @@ namespace AMDiS {
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
double psival;
double psival;
double *phival = GET_MEMORY(double, nCol);
int nPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, nPoints);
......@@ -544,7 +544,7 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
for (int j = 0; j < nCol; j++) {
(*mat)[i][j] += quadrature->getWeight(iq)*c[iq]*psival*phival[j];
(*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psival * phival[j];
}
}
}
......@@ -1191,12 +1191,14 @@ namespace AMDiS {
Quad2::~Quad2()
{
int nPoints = quadrature->getNumPoints();
for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
for (int j = 0; j < nPoints; j++) {
DELETE tmpLALt[i][j];
if (!firstCall) {
int nPoints = quadrature->getNumPoints();
for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
for (int j = 0; j < nPoints; j++) {
DELETE tmpLALt[i][j];
}
DELETE [] tmpLALt[i];
}
DELETE [] tmpLALt[i];
}
}
......@@ -1219,7 +1221,7 @@ namespace AMDiS {
}
DimMat<double> **LALt = tmpLALt[myRank];
for (int i = 0; i < nPoints; i++) {
LALt[i]->set(0.0);
}
......@@ -1227,8 +1229,8 @@ namespace AMDiS {
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
}
VectorOfFixVecs<DimVec<double> > *grdPsi, *grdPhi;
VectorOfFixVecs< DimVec<double> > *grdPsi, *grdPhi;
if (symmetric) {
for (int iq = 0; iq < nPoints; iq++) {
......@@ -1239,10 +1241,11 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
(*mat)[i][i] += quadrature->getWeight(iq) *
((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[i]));
xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[i]);
for (int j = i + 1; j < nCol; j++) {
double val = quadrature->getWeight(iq) * ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[j]));
double val = quadrature->getWeight(iq) *
xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
(*mat)[i][j] += val;
(*mat)[j][i] += val;
}
......@@ -1257,8 +1260,8 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) {
(*mat)[i][j] += quadrature->getWeight(iq) *
((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[j]));
(*mat)[i][j] += quadrature->getWeight(iq) *
xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
}
}
}
......
AMDiS/src/DOFMatrix.cc
View file @ 23e3a357
......@@ -440,9 +440,6 @@ namespace AMDiS {
}
Operator *operat = op ? op : operators[0];
// !!! Do no combine the next two lines into one !!!
// ElementMatrix *elementMatrix = NULL;
operat->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
if (op) {
......@@ -452,17 +449,14 @@ namespace AMDiS {
::std::vector<double*>::iterator factorIt;
for (it = operators.begin(), factorIt = operatorFactor.begin();
it != operators.end();
++it, ++factorIt)
{
(*it)->getElementMatrix(elInfo,
elementMatrix,
*factorIt ? **factorIt : 1.0);
}
++it, ++factorIt) {
(*it)->getElementMatrix(elInfo,
elementMatrix,
*factorIt ? **factorIt : 1.0);
}
}
addElementMatrix(factor, *elementMatrix, bound);
// DELETE elementMatrix;
}
Flag DOFMatrix::getAssembleFlag()
......
AMDiS/src/DOFVector.cc
View file @ 23e3a357
......@@ -71,9 +71,9 @@ namespace AMDiS {
DOFAdmin *admin = feSpace->getAdmin();
// count number of nodes and dofs per node
::std::vector<int> numNodeDOFs;
::std::vector<int> numNodePreDOFs;
::std::vector<DimVec<double>*> bary;
std::vector<int> numNodeDOFs;
std::vector<int> numNodePreDOFs;
std::vector<DimVec<double>*> bary;
int numNodes = 0;
int numDOFs = 0;
......@@ -102,7 +102,7 @@ namespace AMDiS {
}
// traverse mesh
::std::vector<bool> visited(getUsedSize(), false);
std::vector<bool> visited(getUsedSize(), false);
TraverseStack stack;
......@@ -213,10 +213,10 @@ namespace AMDiS {
}
template<>
const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *elInfo,
const Quadrature *quad,
const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *elInfo,
const Quadrature *quad,
const FastQuadrature *quadFast,
WorldVector<double> *grdAtQPs) const
WorldVector<double> *grdAtQPs) const
{
FUNCNAME("DOFVector<double>::getGrdAtQPs()");
......@@ -231,17 +231,12 @@ namespace AMDiS {
("invalid basis functions");
Element *el = elInfo->getElement();
int dim = elInfo->getMesh()->getDim();
int dow = Global::getGeo(WORLD);
const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
const BasisFunction *basFcts = feSpace->getBasisFcts();
int numPoints = quadrature->getNumPoints();
static WorldVector<double> *grd = NULL;
WorldVector<double> *result;
if (grdAtQPs) {
......@@ -589,9 +584,9 @@ namespace AMDiS {
}
// count number of nodes and dofs per node
::std::vector<int> numNodeDOFs;
::std::vector<int> numNodePreDOFs;
::std::vector<DimVec<double>*> bary;
std::vector<int> numNodeDOFs;
std::vector<int> numNodePreDOFs;
std::vector<DimVec<double>*> bary;
int numNodes = 0;
int numDOFs = 0;
......@@ -619,7 +614,7 @@ namespace AMDiS {
}
// traverse mesh
::std::vector<bool> visited(getUsedSize(), false);
std::vector<bool> visited(getUsedSize(), false);
TraverseStack stack;
Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
......@@ -661,8 +656,8 @@ namespace AMDiS {
void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) {
::std::vector<DegreeOfFreedom>::iterator it;
::std::vector<DegreeOfFreedom>::iterator end = vec.end();
std::vector<DegreeOfFreedom>::iterator it;
std::vector<DegreeOfFreedom>::iterator end = vec.end();
DegreeOfFreedom pos = 0;
for (it = vec.begin(); it != end; ++it, ++pos) {
if (*it == dof) *it = pos;
......
AMDiS/src/DOFVector.h
View file @ 23e3a357
......@@ -78,7 +78,7 @@ namespace AMDiS {
nBasFcts(0)
{};
DOFVectorBase(const FiniteElemSpace *f, ::std::string n);
DOFVectorBase(const FiniteElemSpace *f, std::string n);
virtual ~DOFVectorBase();
......@@ -121,28 +121,28 @@ namespace AMDiS {
operatorEstFactor.push_back(estFactor);
};
inline ::std::vector<double*>::iterator getOperatorFactorBegin() {
inline std::vector<double*>::iterator getOperatorFactorBegin() {
return operatorFactor.begin();
};
inline ::std::vector<double*>::iterator getOperatorFactorEnd() {
inline std::vector<double*>::iterator getOperatorFactorEnd() {
return operatorFactor.end();
};
inline ::std::vector<double*>::iterator getOperatorEstFactorBegin() {
inline std::vector<double*>::iterator getOperatorEstFactorBegin() {
return operatorEstFactor.begin();
};
inline ::std::vector<double*>::iterator getOperatorEstFactorEnd() {
inline std::vector<double*>::iterator getOperatorEstFactorEnd() {
return operatorEstFactor.end();
};
inline ::std::vector<Operator*>::iterator getOperatorsBegin() {
inline std::vector<Operator*>::iterator getOperatorsBegin() {
return operators.begin();
};
inline ::std::vector<Operator*>::iterator getOperatorsEnd() {
inline std::vector<Operator*>::iterator getOperatorsEnd() {
return operators.end();
};
......@@ -156,22 +156,22 @@ namespace AMDiS {
*/
T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);
inline ::std::vector<Operator*>& getOperators() {
inline std::vector<Operator*>& getOperators() {
return operators;
};
inline ::std::vector<double*>& getOperatorFactor() {
inline std::vector<double*>& getOperatorFactor() {
return operatorFactor;
};
inline ::std::vector<double*>& getOperatorEstFactor() {
inline std::vector<double*>& getOperatorEstFactor() {
return operatorEstFactor;
};
/** \brief
* Returns \ref name
*/
inline const ::std::string& getName() const {
inline const std::string& getName() const {
return name;
};
......@@ -192,7 +192,7 @@ namespace AMDiS {
/** \brief
*
*/
::std::string name;
std::string name;
/** \brief
*
......@@ -202,17 +202,17 @@ namespace AMDiS {
/** \brief
*
*/
::std::vector<Operator*> operators;
std::vector<Operator*> operators;
/** \brief
*
*/
::std::vector<double*> operatorFactor;
std::vector<double*> operatorFactor;
/** \brief
*
*/
::std::vector<double*> operatorEstFactor;
std::vector<double*> operatorEstFactor;
/** \brief
*
......@@ -228,7 +228,7 @@ namespace AMDiS {
* Are used to store temporary local dofs of an element.
* Thread safe.
*/
::std::vector<DegreeOfFreedom* > localIndices;
std::vector<DegreeOfFreedom* > localIndices;
};
......@@ -314,12 +314,12 @@ namespace AMDiS {
/** \brief
* Constructs a DOFVector with name n belonging to FiniteElemSpace f
*/
DOFVector(const FiniteElemSpace* f, ::std::string n);
DOFVector(const FiniteElemSpace* f, std::string n);
/** \brief
* Initialization.
*/
void init(const FiniteElemSpace* f, ::std::string n);
void init(const FiniteElemSpace* f, std::string n);
/** \brief
* Copy Constructor
......@@ -340,14 +340,14 @@ namespace AMDiS {
/** \brief
* Returns iterator to the begin of \ref vec
*/
typename ::std::vector<T>::iterator begin() {
typename std::vector<T>::iterator begin() {
return vec.begin();
};
/** \brief
* Returns iterator to the end of \ref vec
*/
typename ::std::vector<T>::iterator end() {
typename std::vector<T>::iterator end() {
return vec.end();
};
......@@ -356,7 +356,7 @@ namespace AMDiS {
* DOFIndexedBase::compress()
*/
virtual void compressDOFIndexed(int first, int last,
::std::vector<DegreeOfFreedom> &newDof);
std::vector<DegreeOfFreedom> &newDof);
/** \brief
* Sets \ref refineInter to b
......@@ -399,7 +399,7 @@ namespace AMDiS {
/** \brief
* Returns \ref vec
*/
::std::vector<T>& getVector() {
std::vector<T>& getVector() {
return vec;
};
......@@ -580,13 +580,13 @@ namespace AMDiS {
// ===== Serializable implementation =====
void serialize(::std::ostream &out) {
void serialize(std::ostream &out) {
unsigned int size = vec.size();
out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
};
void deserialize(::std::istream &in) {
void deserialize(std::istream &in) {
unsigned int size;
in.read(reinterpret_cast<char*>(&size), sizeof(unsigned int));
vec.resize(size);
......@@ -626,7 +626,7 @@ namespace AMDiS {
/** \brief
* Name of this DOFVector
*/
::std::string name;
std::string name;
/** \brief
* FiniteElemSpace of the vector
......@@ -636,7 +636,7 @@ namespace AMDiS {
/** \brief
* Data container
*/
::std::vector<T> vec;
std::vector<T> vec;
/** \brief
* Specifies whether interpolation should be performed after refinement
......@@ -709,7 +709,7 @@ namespace AMDiS {
* Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
* as DOFContainer at DOFAdmin
*/
DOFVectorDOF(const FiniteElemSpace* feSpace_, ::std::string name_)
DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
: DOFVector<DegreeOfFreedom>(feSpace_, name_)
{
feSpace->getAdmin()->addDOFContainer(this);
......
AMDiS/src/DOFVector.hh
View file @ 23e3a357
......@@ -18,7 +18,7 @@
namespace AMDiS {
template<typename T>
DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, ::std::string n)
DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
: feSpace(f),
name(n),
elementVector(NULL),
......@@ -42,7 +42,7 @@ namespace AMDiS {
}
template<typename T>
DOFVector<T>::DOFVector(const FiniteElemSpace* f, ::std::string n)
DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
: DOFVectorBase<T>(f, n),
refineInter(false),
coarsenOperation(NO_OPERATION)
......@@ -51,7 +51,7 @@ namespace AMDiS {
}
template<typename T>
void DOFVector<T>::init(const FiniteElemSpace* f, ::std::string n)
void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
{
name = n;
feSpace = f;
......@@ -244,7 +244,7 @@ namespace AMDiS {
T m;
Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator) {
m = ::std::min(m, *vecIterator);
m = std::min(m, *vecIterator);
}
return m;
......@@ -264,7 +264,7 @@ namespace AMDiS {
T m;
Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator) {
m = ::std::max(m, *vecIterator);
m = std::max(m, *vecIterator);
}
return m;
......@@ -338,7 +338,7 @@ namespace AMDiS {
++rowIterator, ++vecIterator) {
sum = 0;
row = &(a[rowIterator.getDOFIndex()]);
for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
for(std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
colIterator != rowIterator->end();
colIterator++) {
jcol = colIterator->col;
......@@ -367,7 +367,7 @@ namespace AMDiS {
++rowIterator, ++xIterator) {
ax = alpha * (*xIterator);
row = &(a[rowIterator.getDOFIndex()]);
for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
for(std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
colIterator != rowIterator->end();
colIterator++) {
jcol = colIterator->col;
......@@ -726,7 +726,7 @@ namespace AMDiS {
template<typename T>
void DOFVector<T>::compressDOFIndexed(int first, int last,
::std::vector<DegreeOfFreedom> &newDOF)
std::vector<DegreeOfFreedom> &newDOF)
{
int i, j;
for(i = first; i <= last; i++) {
......@@ -754,8 +754,8 @@ namespace AMDiS {
if (op) {
op->getElementVector(elInfo, this->elementVector);
} else {
::std::vector<Operator*>::iterator it;
::std::vector<double*>::iterator factorIt;
std::vector<Operator*>::iterator it;
std::vector<double*>::iterator factorIt;
for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();
it != this->operators.end();
++it, ++factorIt) {
......@@ -776,7 +776,7 @@ namespace AMDiS {
Flag DOFVectorBase<T>::getAssembleFlag()
{
Flag fillFlag(0);
::std::vector<Operator*>::iterator op;
std::vector<Operator*>::iterator op;
for(op = this->operators.begin(); op != this->operators.end(); ++op) {
fillFlag |= (*op)->getFillFlag();
}
......@@ -962,7 +962,7 @@ namespace AMDiS {
double sum = 0;
if (!add)
*vecIterator = 0.0;
for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
for(std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
colIterator != rowIterator->end();
colIterator++) {
......@@ -1002,7 +1002,7 @@ namespace AMDiS {
!rowIterator.end();
++rowIterator, ++xIterator) {
T ax = (*xIterator);
for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
for (std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
colIterator != rowIterator->end();
colIterator++) {
int jcol = colIterator->col;
......
AMDiS/src/ElInfo2d.cc
View file @ 23e3a357
......@@ -14,6 +14,21 @@
namespace AMDiS {
ElInfo2d::ElInfo2d(Mesh *aMesh)
: ElInfo(aMesh)
{
e1 = NEW WorldVector<double>;
e2 = NEW WorldVector<double>;
normal = NEW WorldVector<double>;
}
ElInfo2d::~ElInfo2d()
{
DELETE e1;
DELETE e2;
DELETE normal;
}
void ElInfo2d::fillMacroInfo(const MacroElement * mel)
{
FUNCNAME("ElInfo::fillMacroInfo()");
......@@ -394,21 +409,19 @@ namespace AMDiS {
double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
{
FUNCNAME("ElInfo2d::calcGrdLambda()");
WorldVector<double> e1, e2;
double adet = 0.0;
const WorldVector<double> v0 = coord_[0];
testFlag(Mesh::FILL_COORDS);
double adet = 0.0;
int dim = mesh_->getDim();
for (int i = 0; i < dimOfWorld; i++) {
e1[i] = coord_[1][i] - v0[i];
e2[i] = coord_[2][i] - v0[i];
(*e1)[i] = coord_[1][i] - coord_[0][i];
(*e2)[i] = coord_[2][i] - coord_[0][i];
}
if (dimOfWorld == 2) {
double sdet = e1[0] * e2[1] - e1[1] * e2[0];
double sdet = (*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0];
adet = abs(sdet);
if (adet < 1.0E-25) {
......@@ -418,24 +431,18 @@ namespace AMDiS {
grd_lam[i][j] = 0.0;
} else {
double det1 = 1.0 / sdet;
double a11 = e2[1] * det1; /* (a_ij) = A^{-T} */
double a21 = -e2[0] * det1;
double a12 = -e1[1] * det1;
double a22 = e1[0] * det1;
grd_lam[1][0] = a11;
grd_lam[1][1] = a21;
grd_lam[2][0] = a12;
grd_lam[2][1] = a22;
grd_lam[1][0] = (*e2)[1] * det1; // a11: (a_ij) = A^{-T}
grd_lam[1][1] = -(*e2)[0] * det1; // a21
grd_lam[2][0] = -(*e1)[1] * det1; // a12
grd_lam[2][1] = (*e1)[0] * det1; // a22
grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
}
} else {
WorldVector<double> normal;
vectorProduct(e1, e2, normal);
} else {
vectorProduct(*e1, *e2, *normal);
adet = norm(&normal);
adet = norm(normal);
if (adet < 1.0E-15) {
MSG("abs(det) = %lf\n", adet);
......@@ -443,8 +450,8 @@ namespace AMDiS {
for (int j = 0; j < dimOfWorld; j++)
grd_lam[i][j] = 0.0;
} else {
vectorProduct(e2, normal, grd_lam[1]);
vectorProduct(normal, e1, grd_lam[2]);
vectorProduct(*e2, *normal, grd_lam[1]);
vectorProduct(*normal, *e1, grd_lam[2]);
double adet2 = 1.0 / (adet * adet);
......
AMDiS/src/ElInfo2d.h
View file @ 23e3a357
......@@ -43,7 +43,9 @@ namespace AMDiS {
/** \brief
* Constructor. Calls ElInfo's protected Constructor.
*/
ElInfo2d(Mesh* aMesh) : ElInfo(aMesh) {};
ElInfo2d(Mesh* aMesh);
~ElInfo2d();
/** \brief
* 2-dimensional realisation of ElInfo's fillElInfo method.
......@@ -74,6 +76,12 @@ namespace AMDiS {
* 2-dimensional realisation of ElInfo's getElementNormal method.
*/
double getElementNormal( WorldVector<double> &normal) const;
protected:
/** \brief
* Temp vectors for function \ref calcGrdLambda.
*/
WorldVector<double> *e1, *e2, *normal;
};