Skip to content

GitLab

Commit 06dbd114 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Made some operators a little bit faster.

parent cf933621
Hide whitespace changes
Inline Side-by-side
Showing with 436 additions and 138 deletions
AMDiS/src/DOFVector.hh
View file @ 06dbd114
......@@ -981,13 +981,13 @@ namespace AMDiS {
{
FUNCNAME("DOFVector<T>::getVecAtQPs()");
TEST_EXIT(quad || quadFast)("neither quad nor quadFast defined\n");
TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
if (quad && quadFast)
TEST_EXIT(quad == quadFast->getQuadrature())
TEST_EXIT_DBG(quad == quadFast->getQuadrature())
("quad != quadFast->quadrature\n");
TEST_EXIT(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
Element *el = elInfo->getElement();
......@@ -1012,8 +1012,8 @@ namespace AMDiS {
result = localvec;
}
T *localVec = new T[nBasFcts];
T *localVec = localVectors[omp_get_thread_num()];
getLocalVector(el, localVec);
for (int i = 0; i < nPoints; i++) {
......@@ -1025,8 +1025,6 @@ namespace AMDiS {
((*(basFcts->getPhi(j)))(quad->getLambda(i))));
}
delete [] localVec;
return const_cast<const T*>(result);
}
......
AMDiS/src/FirstOrderAssembler.cc
View file @ 06dbd114
......@@ -338,8 +338,6 @@ namespace AMDiS {
void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
VectorOfFixVecs<DimVec<double> > *grdPhi;
if (firstCall) {
#ifdef _OPENMP
#pragma omp critical
......@@ -363,25 +361,23 @@ namespace AMDiS {
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
const double *psi = psiFast->getPhi(iq);
grdPhi = phiFast->getGradient(iq);
for (int i = 0; i < nRow; i++) {
double factor = quadrature->getWeight(iq) * psi[i];
for (int j = 0; j < nCol; j++)
mat[i][j] += factor * (Lb[iq] * (*grdPhi)[j]);
for (int i = 0; i < nCol; i++) {
double factor =
quadrature->getWeight(iq) * (Lb[iq] * phiFast->getGradient(iq, i));
for (int j = 0; j < nRow; j++)
mat[j][i] += factor * psi[j];
}
}
}
Pre01::Pre01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
{
}
{}
void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
......
AMDiS/src/Operator.cc
View file @ 06dbd114
......@@ -697,6 +697,17 @@ namespace AMDiS {
auxFESpaces.push_back(dv->getFESpace());
}
VecAtQP_FOT::VecAtQP_FOT(DOFVectorBase<double> *dv,
AbstractFunction<double, double> *af,
int bIdx)
: FirstOrderTerm(af->getDegree()), vec(dv), f(af)
{
TEST_EXIT(dv)("No vector!\n");
bOne = bIdx;
auxFESpaces.push_back(dv->getFESpace());
}
VectorGradient_FOT::VectorGradient_FOT(DOFVectorBase<double> *dv,
AbstractFunction<WorldVector<double>, WorldVector<double> > *af)
: FirstOrderTerm(af->getDegree()), vec(dv), f(af)
......@@ -824,6 +835,15 @@ namespace AMDiS {
auxFESpaces.push_back(dv2->getFESpace());
}
Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
int bIdx)
: FirstOrderTerm(8), vec1(dv1), vec2(dv2)
{
bOne = bIdx;
auxFESpaces.push_back(dv1->getFESpace());
auxFESpaces.push_back(dv2->getFESpace());
}
void Vec2AtQP_FOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
......@@ -836,16 +856,21 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> >& Lb) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
if (b)
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
lb_one(Lambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
} else if (b) {
for (int iq = 0; iq < nPoints; iq++)
lb(Lambda, *b, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
else
} else {
for (int iq = 0; iq < nPoints; iq++)
l1(Lambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
}
}
void Vec2AtQP_FOT::eval(int nPoints,
const double *uhAtQP,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
......@@ -858,14 +883,34 @@ namespace AMDiS {
Vec3FctAtQP_FOT::Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,DOFVectorBase<double> *dv3,
BinaryAbstractFunction<double, double, double> *f_,
WorldVector<double> *b_)
: FirstOrderTerm(10), vec1(dv1), vec2(dv2), vec3(dv3), f(f_), b(b_)
WorldVector<double> *bvec)
: FirstOrderTerm(10),
vec1(dv1),
vec2(dv2),
vec3(dv3),
f(f_),
b(bvec)
{
auxFESpaces.push_back(dv1->getFESpace());
auxFESpaces.push_back(dv2->getFESpace());
auxFESpaces.push_back(dv3->getFESpace());
}
Vec3FctAtQP_FOT::Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,DOFVectorBase<double> *dv3,
BinaryAbstractFunction<double, double, double> *f_,
int b)
: FirstOrderTerm(10),
vec1(dv1),
vec2(dv2),
vec3(dv3),
f(f_)
{
bOne = b;
auxFESpaces.push_back(dv1->getFESpace());
auxFESpaces.push_back(dv2->getFESpace());
auxFESpaces.push_back(dv3->getFESpace());
}
void Vec3FctAtQP_FOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
......@@ -876,19 +921,25 @@ namespace AMDiS {
}
void Vec3FctAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const
void Vec3FctAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
if (b)
lb(Lambda, *b, Lb[iq], (vec1AtQPs[iq])*(*f)(vec2AtQPs[iq],vec3AtQPs[iq]));
else
l1(Lambda, Lb[iq], (vec1AtQPs[iq])*(*f)(vec2AtQPs[iq],vec3AtQPs[iq]));
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
lb_one(Lambda, Lb[iq], vec1AtQPs[iq] * (*f)(vec2AtQPs[iq], vec3AtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++) {
if (b)
lb(Lambda, *b, Lb[iq], vec1AtQPs[iq] * (*f)(vec2AtQPs[iq], vec3AtQPs[iq]));
else
l1(Lambda, Lb[iq], vec1AtQPs[iq] * (*f)(vec2AtQPs[iq], vec3AtQPs[iq]));
}
}
}
void Vec3FctAtQP_FOT::eval(int nPoints,
const double *uhAtQP,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
......@@ -1582,10 +1633,15 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> >& Lb) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
if (b)
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
lb_one(Lambda, Lb[iq], (*f)(vecAtQPs[iq]));
} else if (b) {
for (int iq = 0; iq < nPoints; iq++)
lb(Lambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
else
} else {
for (int iq = 0; iq < nPoints; iq++)
l1(Lambda, Lb[iq], (*f)(vecAtQPs[iq]));
}
}
......
AMDiS/src/Operator.h
View file @ 06dbd114
......@@ -54,7 +54,9 @@ namespace AMDiS {
OperatorTerm(int deg)
: properties(0),
degree(deg),
auxFESpaces(0)
dimOfWorld(Global::getGeo(WORLD)),
auxFESpaces(0),
bOne(-1)
{}
/// Destructor.
......@@ -65,14 +67,11 @@ namespace AMDiS {
* each OperatorTerm belonging to this SubAssembler. E.g., vectors
* and coordinates at quadrature points can be calculated here.
*/
virtual void initElement(const ElInfo*,
SubAssembler*,
Quadrature *quad = NULL)
virtual void initElement(const ElInfo*, SubAssembler*, Quadrature *quad = NULL)
{}
virtual void initElement(const ElInfo* largeElInfo,
const ElInfo* smallElInfo,
SubAssembler*,
virtual void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler*,
Quadrature *quad = NULL)
{}
......@@ -100,6 +99,12 @@ namespace AMDiS {
return degree;
}
/// Sets one component of the b vector to be one. See \ref bOne.
void setB(int b)
{
bOne = b;
}
/// Evaluation of the OperatorTerm at all quadrature points.
virtual void eval(int nPoints,
const double *uhAtQP,
......@@ -166,12 +171,11 @@ namespace AMDiS {
double factor) const;
/// Evaluation of \f$ \Lambda \cdot b\f$.
static inline void lb(const DimVec<WorldVector<double> >& Lambda,
const WorldVector<double>& b,
DimVec<double>& Lb,
double factor)
inline void lb(const DimVec<WorldVector<double> >& Lambda,
const WorldVector<double>& b,
DimVec<double>& Lb,
double factor) const
{
const int dimOfWorld = Global::getGeo(WORLD);
const int dim = Lambda.size() - 1;
for (int i = 0; i <= dim; i++) {
......@@ -184,15 +188,25 @@ namespace AMDiS {
}
}
/// Evaluation of \f$ \Lambda \cdot b\f$.
inline void lb_one(const DimVec<WorldVector<double> >& Lambda,
DimVec<double>& Lb,
double factor) const
{
const int dim = Lambda.size();
for (int i = 0; i < dim; i++)
Lb[i] += Lambda[i][bOne] * factor;
}
/** \brief
* Evaluation of \f$ \Lambda \cdot b\f$ if b contains the value 1.0 in
* each component.
*/
static inline void l1(const DimVec<WorldVector<double> >& Lambda,
DimVec<double>& Lb,
double factor)
inline void l1(const DimVec<WorldVector<double> >& Lambda,
DimVec<double>& Lb,
double factor) const
{
const int dimOfWorld = Global::getGeo(WORLD);
const int dim = Lambda.size() - 1;
for (int i = 0; i <= dim; i++) {
......@@ -212,12 +226,24 @@ namespace AMDiS {
/// Polynomial degree of the term. Used to detemine the degree of the quadrature.
int degree;
/// Stores the dimension of the world.
int dimOfWorld;
/// List off all fe spaces, the operator term makes use off.
std::vector<const FiniteElemSpace*> auxFESpaces;
/// Pointer to the Operator this OperatorTerm belongs to.
Operator* operat;
/** \brief
* In many cases, the vector b in the evaluation \f$ \Lambda \cdot b\f$ has zeros
* in all components expect one that is set to one. Using the function \ref lb is
* then unnecessary time consuming. Instead, this variable defines the component
* of the vector b to be one. The function \ref lb_one is used if this variable is
* not -1.
*/
int bOne;
/// Flag for piecewise constant terms
static const Flag PW_CONST;
......@@ -1466,17 +1492,27 @@ namespace AMDiS {
: FirstOrderTerm(0), b(wv)
{}
/** \brief
* Implements FirstOrderTerm::getLb().
*/
/// Constructor.
Vector_FOT(int bIdx)
: FirstOrderTerm(0)
{
bOne = bIdx;
}
/// Implements FirstOrderTerm::getLb().
inline void getLb(const ElInfo *elInfo,
int nPoints,
VectorOfFixVecs<DimVec<double> >&Lb) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
lb(Lambda, b, Lb[iq], 1.0);
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
lb_one(Lambda, Lb[iq], 1.0);
} else {
for (int iq = 0; iq < nPoints; iq++)
lb(Lambda, b, Lb[iq], 1.0);
}
}
/// Implements FirstOrderTerm::eval().
......@@ -1507,10 +1543,15 @@ namespace AMDiS {
{
public:
/// Constructor.
VecAtQP_FOT(DOFVectorBase<double> *dv,
VecAtQP_FOT(DOFVectorBase<double> *dv,
AbstractFunction<double, double> *af,
WorldVector<double> *wv);
/// Constructor.
VecAtQP_FOT(DOFVectorBase<double> *dv,
AbstractFunction<double, double> *af,
int bIdx);
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
......@@ -1931,22 +1972,22 @@ namespace AMDiS {
class FctVecAtQP_FOT : public FirstOrderTerm
{
public:
FctVecAtQP_FOT(DOFVectorBase<double> *dv,
AbstractFunction<double, WorldVector<double> > *f_,
WorldVector<double> *b_);
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const;
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
double factor) const;
void getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const;
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
double factor) const;
protected:
DOFVectorBase<double>* vec;
......@@ -1960,22 +2001,25 @@ namespace AMDiS {
class Vec2AtQP_FOT : public FirstOrderTerm
{
public:
Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
WorldVector<double> *b_);
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
int bIdx);
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const;
void eval(int nPoints,
const double *uhAtQP,
void getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const;
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
double factor) const;
protected:
DOFVectorBase<double>* vec1;
DOFVectorBase<double>* vec2;
......@@ -1988,18 +2032,22 @@ namespace AMDiS {
class Vec3FctAtQP_FOT : public FirstOrderTerm
{
public:
Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3,
BinaryAbstractFunction<double, double, double> *f,
WorldVector<double> *b);
Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,DOFVectorBase<double> *dv3,
BinaryAbstractFunction<double, double, double> *f_,
WorldVector<double> *b_);
Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3,
BinaryAbstractFunction<double, double, double> *f,
int b);
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const;
void getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const;
void eval(int nPoints,
const double *uhAtQP,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
......@@ -3513,42 +3561,30 @@ namespace AMDiS {
const FiniteElemSpace *rowFESpace,
const FiniteElemSpace *colFESpace = NULL);
/** \brief
* Destructor.
*/
/// Destructor.
virtual ~Operator() {}
/** \brief
* Sets \ref optimized.
*/
/// Sets \ref optimized.
inline void useOptimizedAssembler(bool opt)
{
optimized = opt;
}
/** \brief
* Returns \ref optimized.
*/
/// Returns \ref optimized.
inline bool isOptimized()
{
return optimized;
}
/** \brief
* Adds a SecondOrderTerm to the Operator
*/
/// Adds a SecondOrderTerm to the Operator
template <typename T>
void addSecondOrderTerm(T *term);
/** \brief
* Adds a FirstOrderTerm to the Operator
*/
/// Adds a FirstOrderTerm to the Operator
template <typename T>
void addFirstOrderTerm(T *term,
FirstOrderType type = GRD_PHI);
/** \brief
* Adds a ZeroOrderTerm to the Operator
*/
void addFirstOrderTerm(T *term, FirstOrderType type = GRD_PHI);
/// Adds a ZeroOrderTerm to the Operator
template <typename T>
void addZeroOrderTerm(T *term);
......@@ -3604,14 +3640,10 @@ namespace AMDiS {
return auxFESpaces;
}
/** \brief
* Sets \ref uhOld.
*/
/// Sets \ref uhOld.
void setUhOld(const DOFVectorBase<double> *uhOld);
/** \brief
* Returns \ref uhOld.
*/
/// Returns \ref uhOld.
inline const DOFVectorBase<double> *getUhOld()
{
return uhOld;
......
AMDiS/src/ProblemVec.cc
View file @ 06dbd114
......@@ -619,8 +619,9 @@ namespace AMDiS {
{
FUNCNAME("ProblemVec::buildAfterCoarsen()");
clock_t first = clock();
// buildAfterCoarsen_sebastianMode(adaptInfo, flag);
clock_t first = clock();
#ifdef _OPENMP
double wtime = omp_get_wtime();
#endif
......@@ -648,13 +649,6 @@ namespace AMDiS {
for (int i = 0; i < nComponents; i++) {
#if 0
std::vector<int> nnzInRow;
int overallNnz, averageNnz;
componentSpaces[i]->getMesh()->computeMatrixFillin(componentSpaces[i],
nnzInRow, overallNnz, averageNnz);
#endif
MSG("%d DOFs for %s\n",
componentSpaces[i]->getAdmin()->getUsedSize(),
componentSpaces[i]->getName().c_str());
......@@ -795,6 +789,201 @@ namespace AMDiS {
#endif
}
void ProblemVec::buildAfterCoarsen_sebastianMode(AdaptInfo *adaptInfo, Flag flag)
{
FUNCNAME("ProblemVec::buildAfterCoarsen()");
clock_t first = clock();
#ifdef _OPENMP
double wtime = omp_get_wtime();
#endif
for (int i = 0; i < static_cast<int>(meshes.size()); i++)
meshes[i]->dofCompress();
Flag assembleFlag =
flag |
(*systemMatrix)[0][0]->getAssembleFlag() |
rhs->getDOFVector(0)->getAssembleFlag() |
Mesh::CALL_LEAF_EL |
Mesh::FILL_COORDS |
Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA |
Mesh::FILL_NEIGH;
if (useGetBound)
assembleFlag |= Mesh::FILL_BOUND;
traverseInfo.updateStatus();
// Used to calculate the overall number of non zero entries.
int nnz = 0;
/// === INITIALIZE ===