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

Generic reimplementation of getVectorAtQPs().

parent d5e25a2a
...@@ -107,8 +107,8 @@ namespace AMDiS { ...@@ -107,8 +107,8 @@ namespace AMDiS {
{ {
name = "standard first order assembler"; name = "standard first order assembler";
psi = owner->getRowFeSpace()->getBasisFcts(); psi = rowFeSpace->getBasisFcts();
phi = owner->getColFeSpace()->getBasisFcts(); phi = colFeSpace->getBasisFcts();
} }
...@@ -182,9 +182,9 @@ namespace AMDiS { ...@@ -182,9 +182,9 @@ namespace AMDiS {
#pragma omp critical #pragma omp critical
#endif #endif
{ {
const BasisFunction *basFcts = owner->getRowFeSpace()->getBasisFcts(); const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI); psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
basFcts = owner->getColFeSpace()->getBasisFcts(); basFcts = colFeSpace->getBasisFcts();
phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI); phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
firstCall = false; firstCall = false;
} }
...@@ -225,9 +225,9 @@ namespace AMDiS { ...@@ -225,9 +225,9 @@ namespace AMDiS {
#pragma omp critical #pragma omp critical
#endif #endif
{ {
const BasisFunction *basFcts = owner->getRowFeSpace()->getBasisFcts(); const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI); psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
basFcts = owner->getColFeSpace()->getBasisFcts(); basFcts = colFeSpace->getBasisFcts();
phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI); phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
firstCall = false; firstCall = false;
} }
...@@ -270,11 +270,10 @@ namespace AMDiS { ...@@ -270,11 +270,10 @@ namespace AMDiS {
#pragma omp critical #pragma omp critical
#endif #endif
{ {
q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFeSpace()->getBasisFcts(), q10 = Q10PsiPhi::provideQ10PsiPhi(rowFeSpace->getBasisFcts(),
owner->getColFeSpace()->getBasisFcts(), colFeSpace->getBasisFcts(),
quadrature); quadrature);
q1 = Q1Psi::provideQ1Psi(owner->getRowFeSpace()->getBasisFcts(), q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature);
quadrature);
firstCall = false; firstCall = false;
} }
} }
...@@ -309,8 +308,8 @@ namespace AMDiS { ...@@ -309,8 +308,8 @@ namespace AMDiS {
grdPhi(assembler->getRowFeSpace()->getMesh()->getDim(), nCol, NO_INIT); grdPhi(assembler->getRowFeSpace()->getMesh()->getDim(), nCol, NO_INIT);
tmpGrdPhi.resize(omp_get_overall_max_threads(), grdPhi); tmpGrdPhi.resize(omp_get_overall_max_threads(), grdPhi);
psi = owner->getRowFeSpace()->getBasisFcts(); psi = rowFeSpace->getBasisFcts();
phi = owner->getColFeSpace()->getBasisFcts(); phi = colFeSpace->getBasisFcts();
} }
void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat) void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
...@@ -356,9 +355,9 @@ namespace AMDiS { ...@@ -356,9 +355,9 @@ namespace AMDiS {
#pragma omp critical #pragma omp critical
#endif #endif
{ {
const BasisFunction *basFcts = owner->getRowFeSpace()->getBasisFcts(); const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI); psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
basFcts = owner->getColFeSpace()->getBasisFcts(); basFcts = colFeSpace->getBasisFcts();
phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI); phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
firstCall = false; firstCall = false;
} }
...@@ -406,11 +405,10 @@ namespace AMDiS { ...@@ -406,11 +405,10 @@ namespace AMDiS {
#pragma omp critical #pragma omp critical
#endif #endif
{ {
q01 = Q01PsiPhi::provideQ01PsiPhi(owner->getRowFeSpace()->getBasisFcts(), q01 = Q01PsiPhi::provideQ01PsiPhi(rowFeSpace->getBasisFcts(),
owner->getColFeSpace()->getBasisFcts(), colFeSpace->getBasisFcts(),
quadrature); quadrature);
q1 = Q1Psi::provideQ1Psi(owner->getRowFeSpace()->getBasisFcts(), q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature);
quadrature);
firstCall = false; firstCall = false;
} }
} }
...@@ -451,11 +449,10 @@ namespace AMDiS { ...@@ -451,11 +449,10 @@ namespace AMDiS {
#pragma omp critical #pragma omp critical
#endif #endif
{ {
q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFeSpace()->getBasisFcts(), q10 = Q10PsiPhi::provideQ10PsiPhi(rowFeSpace->getBasisFcts(),
owner->getColFeSpace()->getBasisFcts(), colFeSpace->getBasisFcts(),
quadrature); quadrature);
q1 = Q1Psi::provideQ1Psi(owner->getRowFeSpace()->getBasisFcts(), q1 = Q1Psi::provideQ1Psi(rowFeSpace->getBasisFcts(), quadrature);
quadrature);
firstCall = false; firstCall = false;
} }
} }
...@@ -466,9 +463,8 @@ namespace AMDiS { ...@@ -466,9 +463,8 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank]; VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb[0].set(0.0); Lb[0].set(0.0);
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) { for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
(static_cast<FirstOrderTerm*>(terms[myRank][i]))->getLb(elInfo, 1, Lb); (static_cast<FirstOrderTerm*>(terms[myRank][i]))->getLb(elInfo, 1, Lb);
}
Lb[0] *= elInfo->getDet(); Lb[0] *= elInfo->getDet();
......
...@@ -833,7 +833,7 @@ namespace AMDiS { ...@@ -833,7 +833,7 @@ namespace AMDiS {
if (no <= 0 || no > getNumber()) { if (no <= 0 || no > getNumber()) {
ERROR("something is wrong, doing nothing\n"); ERROR("something is wrong, doing nothing\n");
rvec[0] = 0.0; rvec[0] = 0.0;
return(const_cast<const double *>( rvec)); return(const_cast<const double *>(rvec));
} }
for (int i = 0; i < no; i++) { for (int i = 0; i < no; i++) {
...@@ -857,8 +857,10 @@ namespace AMDiS { ...@@ -857,8 +857,10 @@ namespace AMDiS {
return(const_cast<const double *>( rvec)); return(const_cast<const double *>( rvec));
} }
const WorldVector<double>* const WorldVector<double>*
Lagrange::interpol(const ElInfo *el_info, int no, Lagrange::interpol(const ElInfo *el_info,
int no,
const int *b_no, const int *b_no,
AbstractFunction<WorldVector<double>, WorldVector<double> > *f, AbstractFunction<WorldVector<double>, WorldVector<double> > *f,
WorldVector<double> *vec) WorldVector<double> *vec)
......
...@@ -23,55 +23,6 @@ namespace AMDiS { ...@@ -23,55 +23,6 @@ namespace AMDiS {
} }
double *OperatorTerm::getVectorAtQPs(DOFVectorBase<double>* vec,
const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
FUNCNAME("OperatorTerm::getVectorAtQPs()");
TEST_EXIT_DBG(elInfo->getMesh() == vec->getFeSpace()->getMesh())
("There is something wrong!\n");
return subAssembler->getVectorAtQPs(vec, elInfo, quad);
}
double *OperatorTerm::getVectorAtQPs(DOFVectorBase<double>* vec,
const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
FUNCNAME("OperatorTerm::getVectorAtQPs()");
TEST_EXIT(smallElInfo->getMesh() == vec->getFeSpace()->getMesh() ||
largeElInfo->getMesh() == vec->getFeSpace()->getMesh())
("There is something wrong!\n");
if (smallElInfo->getLevel() == largeElInfo->getLevel()) {
// Both elements have the same size, so we can use the simple procedure
// to determine the vecAtQPs.
if (vec->getFeSpace()->getMesh() == smallElInfo->getMesh())
return subAssembler->getVectorAtQPs(vec, smallElInfo, quad);
else
return subAssembler->getVectorAtQPs(vec, largeElInfo, quad);
} else {
// The two elements are different. If the vector is defined on the mesh of the
// small element, we can still use the simple procedure to determine the vecAtQPs.
if (vec->getFeSpace()->getMesh() == largeElInfo->getMesh())
return subAssembler->getVectorAtQPs(vec, smallElInfo, largeElInfo, quad);
else
return subAssembler->getVectorAtQPs(vec, smallElInfo, quad);
}
}
WorldVector<double>* OperatorTerm::getGradientsAtQPs(DOFVectorBase<double>* vec, WorldVector<double>* OperatorTerm::getGradientsAtQPs(DOFVectorBase<double>* vec,
const ElInfo* elInfo, const ElInfo* elInfo,
SubAssembler* subAssembler, SubAssembler* subAssembler,
......
...@@ -110,10 +110,11 @@ namespace AMDiS { ...@@ -110,10 +110,11 @@ namespace AMDiS {
* Determines the value of a dof vector at the quadrature points of a given * Determines the value of a dof vector at the quadrature points of a given
* element. It is used by all VecAtQP like operator terms. * element. It is used by all VecAtQP like operator terms.
*/ */
double *getVectorAtQPs(DOFVectorBase<double>* vec, template<typename T>
const ElInfo* elInfo, T* getVectorAtQPs(DOFVectorBase<T>* vec,
SubAssembler* subAssembler, const ElInfo* elInfo,
Quadrature *quad); SubAssembler* subAssembler,
Quadrature *quad);
/** \brief /** \brief
* Determines the value of a dof vector at the quadrature points of a given * Determines the value of a dof vector at the quadrature points of a given
...@@ -121,11 +122,12 @@ namespace AMDiS { ...@@ -121,11 +122,12 @@ namespace AMDiS {
* meshes using the dual traverse. The element, i.e. the small or the large one, * meshes using the dual traverse. The element, i.e. the small or the large one,
* is choosen, which corresponds to the mesh the dof vector is defined on. * is choosen, which corresponds to the mesh the dof vector is defined on.
*/ */
double *getVectorAtQPs(DOFVectorBase<double>* vec, template<typename T>
const ElInfo* smallElInfo, T* getVectorAtQPs(DOFVectorBase<T>* vec,
const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler* subAssembler, const ElInfo* largeElInfo,
Quadrature *quad); SubAssembler* subAssembler,
Quadrature *quad);
/// ///
WorldVector<double>* getGradientsAtQPs(DOFVectorBase<double>* vec, WorldVector<double>* getGradientsAtQPs(DOFVectorBase<double>* vec,
...@@ -180,4 +182,6 @@ namespace AMDiS { ...@@ -180,4 +182,6 @@ namespace AMDiS {
} }
#endif #include "OperatorTerm.hh"
#endif // AMDIS_OPERATORTERM_H
#include "ElInfo.h"
namespace AMDiS {
template<typename T>
T* OperatorTerm::getVectorAtQPs(DOFVectorBase<T>* vec,
const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
FUNCNAME("OperatorTerm::getVectorAtQPs()");
TEST_EXIT_DBG(elInfo->getMesh() == vec->getFeSpace()->getMesh())
("There is something wrong!\n");
return subAssembler->getVectorAtQPs(vec, elInfo, quad);
}
template<typename T>
T* OperatorTerm::getVectorAtQPs(DOFVectorBase<T>* vec,
const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
FUNCNAME("OperatorTerm::getVectorAtQPs()");
TEST_EXIT(smallElInfo->getMesh() == vec->getFeSpace()->getMesh() ||
largeElInfo->getMesh() == vec->getFeSpace()->getMesh())
("There is something wrong!\n");
if (smallElInfo->getLevel() == largeElInfo->getLevel()) {
// Both elements have the same size, so we can use the simple procedure
// to determine the vecAtQPs.
if (vec->getFeSpace()->getMesh() == smallElInfo->getMesh())
return subAssembler->getVectorAtQPs(vec, smallElInfo, quad);
else
return subAssembler->getVectorAtQPs(vec, largeElInfo, quad);
} else {
// The two elements are different. If the vector is defined on the mesh of the
// small element, we can still use the simple procedure to determine the vecAtQPs.
if (vec->getFeSpace()->getMesh() == largeElInfo->getMesh())
return subAssembler->getVectorAtQPs(vec, smallElInfo, largeElInfo, quad);
else
return subAssembler->getVectorAtQPs(vec, smallElInfo, quad);
}
}
}
...@@ -84,8 +84,8 @@ namespace AMDiS { ...@@ -84,8 +84,8 @@ namespace AMDiS {
{ {
name = "precalculated second order assembler"; name = "precalculated second order assembler";
q11 = Q11PsiPhi::provideQ11PsiPhi(owner->getRowFeSpace()->getBasisFcts(), q11 = Q11PsiPhi::provideQ11PsiPhi(rowFeSpace->getBasisFcts(),
owner->getColFeSpace()->getBasisFcts(), colFeSpace->getBasisFcts(),
quadrature); quadrature);
tmpLALt.resize(omp_get_overall_max_threads()); tmpLALt.resize(omp_get_overall_max_threads());
for (int i = 0; i < omp_get_overall_max_threads(); i++) { for (int i = 0; i < omp_get_overall_max_threads(); i++) {
...@@ -203,10 +203,10 @@ namespace AMDiS { ...@@ -203,10 +203,10 @@ namespace AMDiS {
#pragma omp critical #pragma omp critical
#endif #endif
{ {
psiFast = updateFastQuadrature(psiFast, owner->getRowFeSpace()->getBasisFcts(), psiFast =
INIT_GRD_PHI); updateFastQuadrature(psiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI);
phiFast = updateFastQuadrature(phiFast, owner->getRowFeSpace()->getBasisFcts(), phiFast =
INIT_GRD_PHI); updateFastQuadrature(phiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI);
} }
firstCall = false; firstCall = false;
} }
...@@ -270,8 +270,8 @@ namespace AMDiS { ...@@ -270,8 +270,8 @@ namespace AMDiS {
DimVec<double> grdPsi(dim, NO_INIT); DimVec<double> grdPsi(dim, NO_INIT);
VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT); VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);
const BasisFunction *psi = owner->getRowFeSpace()->getBasisFcts(); const BasisFunction *psi = rowFeSpace->getBasisFcts();
const BasisFunction *phi = owner->getColFeSpace()->getBasisFcts(); const BasisFunction *phi = colFeSpace->getBasisFcts();
int nPoints = quadrature->getNumPoints(); int nPoints = quadrature->getNumPoints();
......
...@@ -16,21 +16,27 @@ namespace AMDiS { ...@@ -16,21 +16,27 @@ namespace AMDiS {
int order, int order,
bool optimized, bool optimized,
FirstOrderType type) FirstOrderType type)
: nRow(0), : rowFeSpace(assembler->rowFeSpace),
colFeSpace(assembler->colFeSpace),
nRow(0),
nCol(0), nCol(0),
coordsAtQPs(NULL), coordsAtQPs(NULL),
coordsNumAllocated(0), coordsNumAllocated(0),
quadrature(quadrat), quadrature(quadrat),
psiFast(NULL), psiFast(NULL),
phiFast(NULL), phiFast(NULL),
owner(assembler),
symmetric(true), symmetric(true),
opt(optimized), opt(optimized),
firstCall(true), firstCall(true),
name("") name("")
{ {
const BasisFunction *psi = assembler->rowFeSpace->getBasisFcts(); FUNCNAME("SubAssembler::SubAssembler()");
const BasisFunction *phi = assembler->colFeSpace->getBasisFcts();
TEST_EXIT(rowFeSpace)("No row FE space defined!\n");
TEST_EXIT(colFeSpace)("No col FE space defined!\n");
const BasisFunction *psi = rowFeSpace->getBasisFcts();
const BasisFunction *phi = colFeSpace->getBasisFcts();
nRow = psi->getNumber(); nRow = psi->getNumber();
nCol = phi->getNumber(); nCol = phi->getNumber();
...@@ -69,7 +75,7 @@ namespace AMDiS { ...@@ -69,7 +75,7 @@ namespace AMDiS {
symmetric = false; symmetric = false;
} }
dim = assembler->rowFeSpace->getMesh()->getDim(); dim = rowFeSpace->getMesh()->getDim();
} }
...@@ -151,97 +157,13 @@ namespace AMDiS { ...@@ -151,97 +157,13 @@ namespace AMDiS {
} }
double* SubAssembler::getVectorAtQPs(DOFVectorBase<double>* dv, WorldVector<double>* SubAssembler::getGradientsAtQPs(DOFVectorBase<double>* vec,
const ElInfo* elInfo,
Quadrature *quad)
{
FUNCNAME("SubAssembler::getVectorAtQPs()");
const DOFVectorBase<double>* vec = dv ? dv : owner->operat->getUhOld();
TEST_EXIT_DBG(vec)("no dof vector!\n");
TEST_EXIT_DBG(elInfo->getMesh() == vec->getFeSpace()->getMesh())
("Vector and fe space do not fit together!\n");
Quadrature *localQuad = quad ? quad : quadrature;
if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid)
return &(valuesAtQPs[vec]->values[0]);
if (!valuesAtQPs[vec])
valuesAtQPs[vec] = new ValuesAtQPs;
valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());
double *values = &(valuesAtQPs[vec]->values[0]);
bool sameFeSpaces =
(vec->getFeSpace() == owner->rowFeSpace) ||
(vec->getFeSpace() == owner->colFeSpace);
if (opt && !quad && sameFeSpaces) {
const BasisFunction *psi = owner->rowFeSpace->getBasisFcts();
const BasisFunction *phi = owner->colFeSpace->getBasisFcts();
if (vec->getFeSpace()->getBasisFcts() == psi)
psiFast = updateFastQuadrature(psiFast, psi, INIT_PHI);
else if(vec->getFeSpace()->getBasisFcts() == phi)
phiFast = updateFastQuadrature(phiFast, phi, INIT_PHI);
}
// calculate new values
const BasisFunction *basFcts = vec->getFeSpace()->getBasisFcts();
if (opt && !quad && sameFeSpaces) {
if (psiFast->getBasisFunctions() == basFcts) {
vec->getVecAtQPs(elInfo, NULL, psiFast, values);
} else if (phiFast->getBasisFunctions() == basFcts) {
vec->getVecAtQPs(elInfo, NULL, phiFast, values);
} else {
vec->getVecAtQPs(elInfo, localQuad, NULL, values);
}
} else {
vec->getVecAtQPs(elInfo, localQuad, NULL, values);
}
valuesAtQPs[vec]->valid = true;
return values;
}
double* SubAssembler::getVectorAtQPs(DOFVectorBase<double>* dv,
const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
Quadrature *quad)
{
FUNCNAME("SubAssembler::getVectorAtQPs()");
const DOFVectorBase<double>* vec = dv ? dv : owner->operat->getUhOld();
TEST_EXIT_DBG(vec)("no dof vector!\n");
Quadrature *localQuad = quad ? quad : quadrature;
if (!valuesAtQPs[vec])
valuesAtQPs[vec] = new ValuesAtQPs;
valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());
double *values = &(valuesAtQPs[vec]->values[0]);
valuesAtQPs[vec]->valid = true;
vec->getVecAtQPs(smallElInfo, largeElInfo, localQuad, NULL, values);
return values;
}
WorldVector<double>* SubAssembler::getGradientsAtQPs(DOFVectorBase<double>* dv,
const ElInfo* elInfo, const ElInfo* elInfo,
Quadrature *quad) Quadrature *quad)
{ {
FUNCNAME("SubAssembler::getGradientsAtQPs()"); FUNCNAME("SubAssembler::getGradientsAtQPs()");
const DOFVectorBase<double>* vec = dv ? dv : owner->operat->getUhOld(); TEST_EXIT_DBG(vec)("No DOF vector!\n");
TEST_EXIT_DBG(vec)("no dof vector!\n");
if (gradientsAtQPs[vec] && gradientsAtQPs[vec]->valid) if (gradientsAtQPs[vec] && gradientsAtQPs[vec]->valid)
return &(gradientsAtQPs[vec]->values[0]); return &(gradientsAtQPs[vec]->values[0]);
...@@ -255,12 +177,11 @@ namespace AMDiS { ...@@ -255,12 +177,11 @@ namespace AMDiS {
WorldVector<double> *values = &(gradientsAtQPs[vec]->values[0]); WorldVector<double> *values = &(gradientsAtQPs[vec]->values[0]);
const BasisFunction *psi = owner->rowFeSpace->getBasisFcts(); const BasisFunction *psi = rowFeSpace->getBasisFcts();
const BasisFunction *phi = owner->colFeSpace->getBasisFcts(); const BasisFunction *phi = colFeSpace->getBasisFcts();
bool sameFeSpaces = bool sameFeSpaces =
(vec->getFeSpace() == owner->rowFeSpace) || vec->getFeSpace() == rowFeSpace || vec->getFeSpace() == colFeSpace;
(vec->getFeSpace() == owner->colFeSpace);
if (opt && !quad && sameFeSpaces) { if (opt && !quad && sameFeSpaces) {
if (vec->getFeSpace()->getBasisFcts() == psi) if (vec->getFeSpace()->getBasisFcts() == psi)
...@@ -287,16 +208,14 @@ namespace AMDiS { ...@@ -287,16 +208,14 @@ namespace AMDiS {
} }
WorldVector<double>* SubAssembler::getGradientsAtQPs(DOFVectorBase<double>* dv, WorldVector<double>* SubAssembler::getGradientsAtQPs(DOFVectorBase<double>* vec,
const ElInfo* smallElInfo, const ElInfo* smallElInfo,
const ElInfo* largeElInfo, const ElInfo* largeElInfo,
Quadrature *quad)