From 45e4e7047e5c079ac8d03ae4ff1ad035a4498e7d Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Mon, 1 Dec 2008 08:22:45 +0000
Subject: [PATCH] * This and that
---
AMDiS/src/DOFMatrix.cc | 14 ++-
AMDiS/src/DOFMatrix.h | 18 ++--
AMDiS/src/Operator.h | 165 +++++++++++++++-----------------
AMDiS/src/ProblemVec.cc | 2 +-
AMDiS/src/ZeroOrderAssembler.cc | 68 ++++++++++++-
AMDiS/src/ZeroOrderAssembler.h | 5 +-
6 files changed, 160 insertions(+), 112 deletions(-)
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 28f2427b..59b329d0 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -632,17 +632,23 @@ namespace AMDiS {
void DOFMatrix::scal(double b)
{
- int i;
DOFMatrix::Iterator rowIterator(this, USED_DOFS);
- for(rowIterator.reset(); !rowIterator.end(); ++rowIterator) {
- for(i=0; i < static_cast<int>((*rowIterator).size()); i++) {
- if((*rowIterator)[i].col >= 0) {
+ for (rowIterator.reset(); !rowIterator.end(); ++rowIterator) {
+ for (int i = 0; i < static_cast<int>((*rowIterator).size()); i++) {
+ if ((*rowIterator)[i].col >= 0) {
(*rowIterator)[i].entry *= b;
}
}
}
}
+ void DOFMatrix::addOperator(Operator *op, double* factor, double* estFactor)
+ {
+ operators.push_back(op);
+ operatorFactor.push_back(factor);
+ operatorEstFactor.push_back(estFactor);
+ }
+
void DOFMatrix::copy(const DOFMatrix& rhs)
{
clear();
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index 0285814a..25cacc21 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -294,23 +294,19 @@ namespace AMDiS {
/** \brief
* a*x + y
*/
- void axpy(double a,
- const DOFMatrix& x,
- const DOFMatrix& y);
+ void axpy(double a, const DOFMatrix& x, const DOFMatrix& y);
/** \brief
* Multiplication with a scalar.
*/
void scal(double s);
- inline void addOperator(Operator *op,
- double* factor = NULL,
- double* estFactor = NULL)
- {
- operators.push_back(op);
- operatorFactor.push_back(factor);
- operatorEstFactor.push_back(estFactor);
- }
+ /** \brief
+ * Adds an operator to the DOFMatrix. A factor, that is multipled
+ * to the operator, and a multilier factor for the estimator may be
+ * also given.
+ */
+ void addOperator(Operator *op, double* factor = NULL, double* estFactor = NULL);
inline std::vector<double*>::iterator getOperatorFactorBegin() {
return operatorFactor.begin();
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 11210f56..fa863da0 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -674,7 +674,7 @@ namespace AMDiS {
: SecondOrderTerm(f_->getDegree()), vec(dv), f(f_)
{
setSymmetric(true);
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -743,7 +743,7 @@ namespace AMDiS {
f(f_)
{
setSymmetric(true);
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -811,7 +811,7 @@ namespace AMDiS {
: SecondOrderTerm(g_->getDegree()), g(g_)
{
setSymmetric(true);
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -877,7 +877,7 @@ namespace AMDiS {
: SecondOrderTerm(f_->getDegree()), vec(dv), f(f_), divFct(divFct_), symmetric(symm)
{
setSymmetric(symmetric);
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -949,7 +949,7 @@ namespace AMDiS {
AbstractFunction<double, WorldVector<double> > *f_)
: SecondOrderTerm(f_->getDegree()), vec(dv), f(f_)
{
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -1018,7 +1018,7 @@ namespace AMDiS {
WorldVector<double> > *f_)
: SecondOrderTerm(f_->getDegree()), vec(dv), f(f_)
{
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -1090,7 +1090,7 @@ namespace AMDiS {
divFct(divFct_), symmetric(symm)
{
setSymmetric(symmetric);
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -1163,7 +1163,7 @@ namespace AMDiS {
WorldVector<double> > *f_)
: SecondOrderTerm(f_->getDegree()), vec(dv), f(f_)
{
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -1233,7 +1233,7 @@ namespace AMDiS {
: SecondOrderTerm(f_->getDegree()), vec(dv), f(f_)
{
setSymmetric(true);
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -1303,7 +1303,7 @@ namespace AMDiS {
: SecondOrderTerm(f_->getDegree()), vec(dv), f(f_), divFct(divFct_), symmetric(symm)
{
setSymmetric(symmetric);
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -1383,7 +1383,7 @@ namespace AMDiS {
{
vecsAtQPs_.resize(vecs_.size());
gradsAtQPs_.resize(grads_.size());
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -1461,7 +1461,7 @@ namespace AMDiS {
{
vecsAtQPs_.resize(vecs_.size());
gradsAtQPs_.resize(grads_.size());
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -1535,12 +1535,12 @@ namespace AMDiS {
*/
FirstOrderTerm(int deg)
: OperatorTerm(deg)
- {};
+ {}
/** \brief
* Destructor.
*/
- virtual ~FirstOrderTerm() {};
+ virtual ~FirstOrderTerm() {}
/** \brief
* Evaluation of \f$ \Lambda b \f$.
@@ -1570,7 +1570,7 @@ namespace AMDiS {
result[iq] += resultQP * factor;
}
}
- };
+ }
};
// ============================================================================
@@ -1590,7 +1590,7 @@ namespace AMDiS {
*/
Simple_FOT()
: FirstOrderTerm(0)
- {};
+ {}
/** \brief
* Implements FirstOrderTerm::getLb().
@@ -1604,7 +1604,7 @@ namespace AMDiS {
for (int iq = 0; iq < numPoints; iq++) {
l1(Lambda, Lb[iq], 1.0);
}
- };
+ }
};
// ============================================================================
@@ -1627,12 +1627,12 @@ namespace AMDiS {
{
factor = NEW double;
*factor = f;
- };
+ }
FactorSimple_FOT(double *fptr)
: FirstOrderTerm(0),
factor(fptr)
- {};
+ {}
/** \brief
* Implements FirstOrderTerm::getLb().
@@ -1644,7 +1644,7 @@ namespace AMDiS {
for (int iq = 0; iq < numPoints; iq++) {
l1(Lambda, Lb[iq], (*factor));
}
- };
+ }
private:
/** \brief
@@ -1670,7 +1670,7 @@ namespace AMDiS {
Vector_FOT(WorldVector<double> b_)
: FirstOrderTerm(0), b(b_)
{
- };
+ }
/** \brief
* Implements FirstOrderTerm::getLb().
@@ -1684,7 +1684,7 @@ namespace AMDiS {
for (int iq = 0; iq < numPoints; iq++) {
lb(Lambda, b, Lb[iq], 1.0);
}
- };
+ }
/** \brief
* Implements FirstOrderTerm::eval().
@@ -1701,7 +1701,7 @@ namespace AMDiS {
result[iq] += b * grdUhAtQP[iq] * factor;
}
}
- };
+ }
protected:
/** \brief
@@ -1729,7 +1729,7 @@ namespace AMDiS {
WorldVector<double> *b_)
: FirstOrderTerm(f_->getDegree()), vec(dv), f(f_), b(b_)
{
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -1789,8 +1789,7 @@ namespace AMDiS {
*/
CoordsAtQP_FOT(AbstractFunction<double, WorldVector<double> > *g_)
: FirstOrderTerm(g_->getDegree()), g(g_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -1840,7 +1839,7 @@ namespace AMDiS {
VecCoordsAtQP_FOT(AbstractFunction<double, WorldVector<double> > *g_,
WorldVector<double> b_)
: FirstOrderTerm(g_->getDegree()), g(g_), b(b_)
- {};
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -1897,8 +1896,7 @@ namespace AMDiS {
VectorGradient_FOT(DOFVectorBase<double> *dv,
AbstractFunction<WorldVector<double>, WorldVector<double> > *f_)
: FirstOrderTerm(f_->getDegree()), vec(dv), f(f_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -1952,8 +1950,7 @@ namespace AMDiS {
VectorFct_FOT(DOFVectorBase<double> *dv,
AbstractFunction<WorldVector<double>, double> *vecFct_)
: FirstOrderTerm(vecFct_->getDegree()), vec(dv), vecFct(vecFct_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2010,8 +2007,7 @@ namespace AMDiS {
*/
VecFctAtQP_FOT(AbstractFunction<WorldVector<double>, WorldVector<double> > *g_)
: FirstOrderTerm(g_->getDegree()), g(g_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2063,7 +2059,7 @@ namespace AMDiS {
BinaryAbstractFunction<WorldVector<double>, double,
WorldVector<double> > *vecFct_)
: FirstOrderTerm(vecFct_->getDegree()), vec1(dv1), vec2(dv2), vecFct(vecFct_)
- {};
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2128,7 +2124,7 @@ namespace AMDiS {
{
vecsAtQPs_.resize(vecs_.size());
gradsAtQPs_.resize(grads_.size());
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2194,7 +2190,7 @@ namespace AMDiS {
{
vecsAtQPs_.resize(vecs_.size());
gradsAtQPs_.resize(grads_.size());
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2255,12 +2251,12 @@ namespace AMDiS {
/** \brief
* Constructor.
*/
- ZeroOrderTerm(int deg) : OperatorTerm(deg) {};
+ ZeroOrderTerm(int deg) : OperatorTerm(deg) {}
/** \brief
* Destructor.
*/
- virtual ~ZeroOrderTerm() {};
+ virtual ~ZeroOrderTerm() {}
/** \brief
* Evaluates \f$ c \f$
@@ -2285,7 +2281,7 @@ namespace AMDiS {
/** \brief
* Constructor.
*/
- Simple_ZOT(double f = 1.0) : ZeroOrderTerm(0), factor(f) {};
+ Simple_ZOT(double f = 1.0) : ZeroOrderTerm(0), factor(f) {}
/** \brief
* Implements ZeroOrderTerm::getC().
@@ -2295,7 +2291,7 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
C[iq] += factor;
}
- };
+ }
/** \brief
* Implements ZeroOrderTerm::eval().
@@ -2310,7 +2306,7 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
result[iq] += fac * factor * uhAtQP[iq];
}
- };
+ }
protected:
/** \brief
@@ -2336,8 +2332,7 @@ namespace AMDiS {
VecAtQP_ZOT(DOFVectorBase<double> *dv,
AbstractFunction<double, double> *f_)
: ZeroOrderTerm(f_ ? f_->getDegree() : 0), vec(dv), f(f_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2395,8 +2390,7 @@ namespace AMDiS {
AbstractFunction<double, double> *f2_)
: ZeroOrderTerm(f1->getDegree()+f2_->getDegree()),
vec1(dv1), vec2(dv2), f1(f1_), f2(f2_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2460,8 +2454,7 @@ namespace AMDiS {
BinaryAbstractFunction<double, double, double> *f_)
: ZeroOrderTerm(f_->getDegree()),
vec1(dv1), vec2(dv2), f(f_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2525,7 +2518,7 @@ namespace AMDiS {
vec2(dv2),
vec3(dv3),
f(f_)
- {};
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2587,8 +2580,7 @@ namespace AMDiS {
BinaryAbstractFunction<double,
WorldVector<double>, WorldVector<double> > *f_)
: ZeroOrderTerm(f_->getDegree()), vec(dv), f(f_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2645,7 +2637,7 @@ namespace AMDiS {
VecGradCoordsAtQP_ZOT(DOFVectorBase<double> *dv,
TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *f_)
: ZeroOrderTerm(f_->getDegree()), vec(dv), f(f_)
- {};
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2709,8 +2701,7 @@ namespace AMDiS {
VecAndCoordsAtQP_ZOT(DOFVectorBase<double> *dv,
BinaryAbstractFunction<double, double, WorldVector<double> > *f_)
: ZeroOrderTerm(f_->getDegree()), vec(dv), f(f_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2771,8 +2762,7 @@ namespace AMDiS {
TertiaryAbstractFunction<double, double, WorldVector<double>,
double > *f_)
: ZeroOrderTerm(f_->getDegree()), vec1(dv1), vec2(dv2), f(f_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2836,8 +2826,7 @@ namespace AMDiS {
FctGradient_ZOT(DOFVectorBase<double> *dv,
AbstractFunction<double, WorldVector<double> > *f_)
: ZeroOrderTerm(f_->getDegree()), vec(dv), f(f_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2892,8 +2881,7 @@ namespace AMDiS {
VecAndGradAtQP_ZOT(DOFVectorBase<double> *dv,
BinaryAbstractFunction<double, double, WorldVector<double> > *f_)
: ZeroOrderTerm(f_->getDegree()), vec(dv), f(f_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -2954,8 +2942,7 @@ namespace AMDiS {
VecAndGradAtQP_SOT(DOFVectorBase<double> *dv,
BinaryAbstractFunction<double, double, WorldVector<double> > *f_)
: SecondOrderTerm(f_->getDegree()), vec(dv), f(f_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -3026,7 +3013,7 @@ namespace AMDiS {
CoordsAtQP_ZOT(AbstractFunction<double, WorldVector<double> > *g_)
: ZeroOrderTerm(g_->getDegree()),
g(g_)
- {};
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -3083,7 +3070,7 @@ namespace AMDiS {
: SecondOrderTerm(g_->getDegree()), g(g_), xi(x_i), xj(x_j)
{
setSymmetric(xi == xj);
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -3153,7 +3140,7 @@ namespace AMDiS {
: SecondOrderTerm(f_->getDegree()), vec(dv), f(f_), xi(x_i), xj(x_j)
{
setSymmetric(xi == xj);
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -3219,8 +3206,7 @@ namespace AMDiS {
DOFVectorBase<double> *dGrd,
BinaryAbstractFunction<double, double, WorldVector<double> > *f_)
: ZeroOrderTerm(f_->getDegree()), vec(dv), vecGrd(dGrd), f(f_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -3283,8 +3269,7 @@ namespace AMDiS {
TertiaryAbstractFunction<double,
double, WorldVector<double>, WorldVector<double> > *f_)
: ZeroOrderTerm(f_->getDegree()), vec(dv), vecGrd1(dGrd1), vecGrd2(dGrd2), f(f_)
- {
- };
+ {}
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -3356,7 +3341,7 @@ namespace AMDiS {
: ZeroOrderTerm(f_->getDegree()), vecs(dv), f(f_)
{
vecsAtQPs.resize(vecs.size());
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -3407,7 +3392,7 @@ namespace AMDiS {
: ZeroOrderTerm(f_->getDegree()), vecs(dv), f(f_)
{
gradsAtQPs.resize(vecs.size());
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -3465,7 +3450,7 @@ namespace AMDiS {
vecs[0] = vec0;
vecs[1] = vec1;
vecs[2] = vec2;
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -3514,7 +3499,7 @@ namespace AMDiS {
: ZeroOrderTerm(f_->getDegree()), vec(vec_), vecs(dv), f(f_)
{
gradsAtQPs.resize(vecs.size());
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -3587,7 +3572,7 @@ namespace AMDiS {
{
vecsAtQPs_.resize(vecs_.size());
gradsAtQPs_.resize(grads_.size());
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -3650,7 +3635,7 @@ namespace AMDiS {
{
vecsAtQPs_.resize(vecs_.size());
gradsAtQPs_.resize(grads_.size());
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -3740,21 +3725,21 @@ namespace AMDiS {
/** \brief
* Destructor.
*/
- virtual ~Operator() {};
+ virtual ~Operator() {}
/** \brief
* Sets \ref optimized.
*/
inline void useOptimizedAssembler(bool opt) {
optimized = opt;
- };
+ }
/** \brief
* Returns \ref optimized.
*/
inline bool isOptimized() {
return optimized;
- };
+ }
/** \brief
* Adds a SecondOrderTerm to the Operator
@@ -3897,7 +3882,7 @@ namespace AMDiS {
++termIt) {
(*termIt)->eval(numPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
}
- };
+ }
/** \brief
@@ -3918,7 +3903,7 @@ namespace AMDiS {
++termIt) {
(*termIt)->eval(numPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
}
- };
+ }
/** \brief
* Evaluation of all terms in \ref firstOrderGrdPhi.
@@ -3938,7 +3923,7 @@ namespace AMDiS {
++termIt) {
(*termIt)->eval(numPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
}
- };
+ }
/** \brief
@@ -3959,7 +3944,7 @@ namespace AMDiS {
++termIt) {
(*termIt)->eval(numPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
}
- };
+ }
/** \brief
* Weak evaluation of all terms in \ref secondOrder.
@@ -3976,7 +3961,7 @@ namespace AMDiS {
++termIt) {
static_cast<SecondOrderTerm*>(*termIt)->weakEval(numPoints, grdUhAtQP, result);
}
- };
+ }
/** \brief
* Calls getLALt() for each term in \ref secondOrder
@@ -3992,7 +3977,7 @@ namespace AMDiS {
++termIt) {
static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, numPoints, LALt);
}
- };
+ }
/** \brief
* Calls getLb() for each term in \ref firstOrderGrdPsi
@@ -4008,7 +3993,7 @@ namespace AMDiS {
++termIt) {
static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, numPoints, Lb);
}
- };
+ }
/** \brief
* Calls getLb() for each term in \ref firstOrderGrdPhi
@@ -4024,7 +4009,7 @@ namespace AMDiS {
++termIt) {
static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, numPoints, Lb);
}
- };
+ }
/** \brief
* Calls getC() for each term in \ref zeroOrder
@@ -4040,14 +4025,14 @@ namespace AMDiS {
++termIt) {
static_cast<ZeroOrderTerm*>(*termIt)->getC(elInfo, numPoints, c);
}
- };
+ }
/** \brief
* Returns true, if there are second order terms. Returns false otherwise.
*/
inline bool secondOrderTerms() {
return secondOrder[omp_get_thread_num()].size() != 0;
- };
+ }
/** \brief
* Returns true, if there are first order terms (grdPsi).
@@ -4055,7 +4040,7 @@ namespace AMDiS {
*/
inline bool firstOrderTermsGrdPsi() {
return firstOrderGrdPsi[omp_get_thread_num()].size() != 0;
- };
+ }
/** \brief
* Returns true, if there are first order terms (grdPhi).
@@ -4063,7 +4048,7 @@ namespace AMDiS {
*/
inline bool firstOrderTermsGrdPhi() {
return firstOrderGrdPhi[omp_get_thread_num()].size() != 0;
- };
+ }
/** \brief
* Returns true, if there are zero order terms.
@@ -4071,7 +4056,7 @@ namespace AMDiS {
*/
inline bool zeroOrderTerms() {
return zeroOrder[omp_get_thread_num()].size() != 0;
- };
+ }
public:
/** \brief
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 1ba7847a..3d8817d3 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -833,7 +833,7 @@ namespace AMDiS {
double *estFactor)
{
FUNCNAME("ProblemVec::addMatrixOperator()");
-
+
if (!(*systemMatrix_)[i][j]) {
TEST_EXIT(i != j)("should have been created already\n");
(*systemMatrix_)[i][j] = NEW DOFMatrix(componentSpaces[i],
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index e3579a36..2d72f50d 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -175,14 +175,13 @@ namespace AMDiS {
// The constants are defined by the subElement Matrix.
+ int myRank = omp_get_thread_num();
int nPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, nPoints);
-
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(rowElInfo, nPoints, c);
@@ -205,6 +204,8 @@ namespace AMDiS {
}
}
}
+
+ FREE_MEMORY(c, double, nPoints);
}
void StandardZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
@@ -308,6 +309,69 @@ namespace AMDiS {
}
}
+ void FastQuadZOA::calculateElementMatrix(const ElInfo *rowElInfo,
+ const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
+ ElementMatrix *mat)
+ {
+ int nPoints = quadrature->getNumPoints();
+ int myRank = omp_get_thread_num();
+ DimMat<double> *m = smallElInfo->getSubElemCoordsMat();
+
+ 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();
+ phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
+ firstCall = false;
+ }
+ }
+
+ double *c = cPtrs[myRank];
+ for (int iq = 0; iq < nPoints; iq++) {
+ c[iq] = 0.0;
+ }
+
+ std::vector<OperatorTerm*>::iterator termIt;
+ for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
+ (static_cast<ZeroOrderTerm*>((*termIt)))->getC(rowElInfo, nPoints, c);
+ }
+
+ for (int iq = 0; iq < nPoints; iq++) {
+ c[iq] *= smallElInfo->getDet();
+
+ const double *psi = psiFast->getPhi(iq);
+ const double *phi = phiFast->getPhi(iq);
+
+ // for (int i = 0; i < nRow; i++) {
+ // for (int j = 0; j < nCol; j++) {
+ // (*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
+ // }
+ // }
+
+ for (int i = 0; i < nCol; i++) {
+ for (int j = 0; j < nRow; j++) {
+ double val = quadrature->getWeight(iq) * c[iq] * psi[i];
+
+ double tmpval = 0.0;
+ for (int k = 0; k < nCol; k++) {
+ tmpval += (*m)[j][k] * phi[k];
+ }
+ val *= tmpval;
+
+ (*mat)[j][i] += val;
+ }
+ }
+
+ }
+ }
+
void FastQuadZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
if (firstCall) {
diff --git a/AMDiS/src/ZeroOrderAssembler.h b/AMDiS/src/ZeroOrderAssembler.h
index aa33fcce..4c40f300 100644
--- a/AMDiS/src/ZeroOrderAssembler.h
+++ b/AMDiS/src/ZeroOrderAssembler.h
@@ -121,10 +121,7 @@ namespace AMDiS {
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
- ElementMatrix *mat)
- {
- ERROR_EXIT("CEM not yet\n");
- }
+ ElementMatrix *mat);
/** \brief
* Implements SubAssembler::calculateElementVector().
--
GitLab