diff --git a/AMDiS/src/AdaptInstationary.cc b/AMDiS/src/AdaptInstationary.cc
index e57ba42d52cdae8c025167e37e724e48b9dd1691..3ed06ba55e2b3e2ef967c14d72c7cce3b6951f09 100644
--- a/AMDiS/src/AdaptInstationary.cc
+++ b/AMDiS/src/AdaptInstationary.cc
@@ -80,7 +80,8 @@ namespace AMDiS {
// estimate before first adaption
if (adaptInfo->getTime() <= adaptInfo->getStartTime())
problemIteration_->oneIteration(adaptInfo, ESTIMATE);
-
+
+
// increment time
adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep());
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 1b2b1fa54bba85dfa4bef45d948cc91ffaad8a95..2ba2e35df10619c7a37a8acfb12114cdc5662603 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -72,11 +72,13 @@ namespace AMDiS {
DOFMatrix::~DOFMatrix()
{
FUNCNAME("DOFMatrix::~DOFMatrix()");
+
if (rowFESpace && rowFESpace->getAdmin())
(const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this);
-
- if (boundaryManager) delete boundaryManager;
- if (inserter) delete inserter;
+ if (boundaryManager)
+ delete boundaryManager;
+ if (inserter)
+ delete inserter;
}
void DOFMatrix::print() const
@@ -135,11 +137,10 @@ namespace AMDiS {
int non_symmetric = !symmetric();
- if (non_symmetric) {
+ if (non_symmetric)
MSG("matrix `%s' not symmetric.\n", name.data());
- } else {
+ else
MSG("matrix `%s' is symmetric.\n", name.data());
- }
}
DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
@@ -154,12 +155,11 @@ namespace AMDiS {
if (rhs.inserter == 0 && inserter == 0)
matrix = rhs.matrix;
- if (rhs.boundaryManager) {
+ if (rhs.boundaryManager)
boundaryManager = new BoundaryManager(*rhs.boundaryManager);
- } else {
+ else
boundaryManager = NULL;
- }
-
+
nRow = rhs.nRow;
nCol = rhs.nCol;
elementMatrix.change_dim(nRow, nCol);
@@ -302,9 +302,8 @@ namespace AMDiS {
{
FUNCNAME("DOFMatrix::assemble2()");
- if (!op && operators.size() == 0) {
+ if (!op && operators.size() == 0)
return;
- }
set_to_zero(elementMatrix);
@@ -335,8 +334,7 @@ namespace AMDiS {
{
// call the operatos cleanup procedures
for (std::vector<Operator*>::iterator it = operators.begin();
- it != operators.end();
- ++it)
+ it != operators.end(); ++it)
(*it)->finishAssembling();
}
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index 0515b6b2823004161203100e3c849af451e83d23..731536f13682026afe3dd3515edc7a0555674c73 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -382,6 +382,23 @@ namespace AMDiS {
boundaryManager = bm;
}
+ /// Calculates the average of non zero entries per row in matrix.
+ void calculateNnz()
+ {
+ nnzPerRow = 0;
+
+ if (num_rows(matrix) != 0)
+ nnzPerRow = int(double(matrix.nnz()) / num_rows(matrix) * 1.2);
+ if (nnzPerRow < 5)
+ nnzPerRow= 5;
+ }
+
+ /// Returns \ref nnzPerRow.
+ int getNnz()
+ {
+ return nnzPerRow;
+ }
+
private:
template <typename T>
void s_write(::std::ostream &out, const T& value)
@@ -536,6 +553,13 @@ namespace AMDiS {
*/
std::set<int> applyDBCs;
+ /* \brief
+ * Number of non zero entries per row (average). For instationary problems this
+ * information may be used in the next timestep to accelerate insertion of
+ * elemnts during assembling.
+ */
+ int nnzPerRow;
+
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
std::map<DegreeOfFreedom, bool> isRankDOF;
#endif
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 69e2b75e797b2fbc4668f51d7d432b755a945103..fce9869225b8550e916d847fc0abb412a030d74b 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -970,11 +970,9 @@ namespace AMDiS {
#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 25000) default(shared) private(i)
#endif
- for (i = 0; i < maxI; i++) {
- if (!admin->isDOFFree(i)) {
+ for (i = 0; i < maxI; i++)
+ if (!admin->isDOFFree(i))
y[i] = alpha * y[i] + x[i];
- }
- }
}
template<typename T>
@@ -984,12 +982,9 @@ namespace AMDiS {
{
typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
- for(vIterator.reset(), rIterator.reset();
- !vIterator.end();
- ++vIterator, ++rIterator)
- {
- *rIterator = scal * (*vIterator);
- };
+ for (vIterator.reset(), rIterator.reset();
+ !vIterator.end(); ++vIterator, ++rIterator)
+ *rIterator = scal * (*vIterator);
return result;
}
@@ -1001,12 +996,10 @@ namespace AMDiS {
{
typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
- for(vIterator.reset(), rIterator.reset();
- !vIterator.end();
- ++vIterator, ++rIterator)
- {
- *rIterator = (*vIterator) + scal;
- };
+ for (vIterator.reset(), rIterator.reset();
+ !vIterator.end(); ++vIterator, ++rIterator)
+ *rIterator = (*vIterator) + scal;
+
return result;
}
@@ -1018,20 +1011,18 @@ namespace AMDiS {
typename DOFVector<T>::Iterator v1Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v1)), USED_DOFS);
typename DOFVector<T>::Iterator v2Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v2)), USED_DOFS);
typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
- for(v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
- !v1Iterator.end();
- ++v1Iterator, ++v2Iterator, ++rIterator)
- {
- *rIterator = (*v1Iterator) + (*v2Iterator);
- };
- return result;
+ for (v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
+ !v1Iterator.end(); ++v1Iterator, ++v2Iterator, ++rIterator)
+ *rIterator = (*v1Iterator) + (*v2Iterator);
+ return result;
}
template<typename T>
const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
{
static T* localVec = NULL;
+ static int localVecSize = 0;
const DOFAdmin* admin = feSpace->getAdmin();
T *result;
@@ -1039,8 +1030,14 @@ namespace AMDiS {
if (d) {
result = d;
} else {
- if(localVec) delete [] localVec;
- localVec = new T[nBasFcts];
+ if (localVec && nBasFcts > localVecSize) {
+ delete [] localVec;
+ localVec = new T[nBasFcts];
+ }
+ if (!localVec)
+ localVec = new T[nBasFcts];
+
+ localVecSize = nBasFcts;
result = localVec;
}
diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc
index 29356cc5e00a27455bfe8e05a667b84c5a197d08..ce14cf4d8df1d677bb40a0e209e5174fedb68e95 100644
--- a/AMDiS/src/FirstOrderAssembler.cc
+++ b/AMDiS/src/FirstOrderAssembler.cc
@@ -47,8 +47,8 @@ namespace AMDiS {
&standardSubAssemblersGrdPhi);
int myRank = omp_get_thread_num();
- std::vector<OperatorTerm*> opTerms
- = (type == GRD_PSI) ? op->firstOrderGrdPsi[myRank] : op->firstOrderGrdPhi[myRank];
+ std::vector<OperatorTerm*> opTerms =
+ (type == GRD_PSI) ? op->firstOrderGrdPsi[myRank] : op->firstOrderGrdPhi[myRank];
// check if a assembler is needed at all
if (opTerms.size() == 0)
@@ -103,15 +103,15 @@ namespace AMDiS {
Stand10::Stand10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
- {}
+ {
+ psi = owner->getRowFESpace()->getBasisFcts();
+ phi = owner->getColFESpace()->getBasisFcts();
+ }
void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
- const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
- const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
-
int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
@@ -132,7 +132,7 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
for (int j = 0; j < nCol; j++)
- mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j];
+ mat[i][j] += quadrature->getWeight(iq) * phival[j] * (Lb[iq] * grdPsi);
}
}
}
@@ -140,7 +140,6 @@ namespace AMDiS {
void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
- const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
@@ -189,13 +188,11 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
Lb[iq].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, nPoints, Lb);
- }
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
@@ -203,10 +200,9 @@ namespace AMDiS {
const double *phi = phiFast->getPhi(iq);
grdPsi = psiFast->getGradient(iq);
- for (int i = 0; i < nRow; i++) {
+ for (int i = 0; i < nRow; i++)
for (int j = 0; j < nCol; j++)
mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j];
- }
}
}
@@ -232,21 +228,17 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
- for (int iq = 0; iq < nPoints; iq++) {
- Lb[iq].set(0.0);
- }
- for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
+ for (int iq = 0; iq < nPoints; iq++)
+ Lb[iq].set(0.0);
+ 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();
grdPsi = psiFast->getGradient(iq);
- for (int i = 0; i < nRow; i++) {
+ for (int i = 0; i < nRow; i++)
vec[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
- }
}
}
@@ -305,13 +297,13 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> >
grdPhi(assembler->getRowFESpace()->getMesh()->getDim(), nCol, NO_INIT);
tmpGrdPhi.resize(omp_get_overall_max_threads(), grdPhi);
+
+ psi = owner->getRowFESpace()->getBasisFcts();
+ phi = owner->getColFESpace()->getBasisFcts();
}
void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
- const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
- const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
-
int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> >& Lb = tmpLb[myRank];
@@ -320,6 +312,7 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++)
Lb[iq].set(0.0);
+
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
@@ -332,7 +325,7 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
for (int j = 0; j < nCol; j++)
- mat[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * psival) * grdPhi[j]);
+ mat[i][j] += quadrature->getWeight(iq) * psival * (Lb[iq] * grdPhi[j]);
}
}
}
diff --git a/AMDiS/src/FirstOrderAssembler.h b/AMDiS/src/FirstOrderAssembler.h
index d6e06b6d9336ffcecef1978a01e7dea9c4b9db7d..b84b4146531c5a22439e40fd6bbbaf50622fb74a 100644
--- a/AMDiS/src/FirstOrderAssembler.h
+++ b/AMDiS/src/FirstOrderAssembler.h
@@ -78,6 +78,9 @@ namespace AMDiS {
/// Implements SubAssembler::calculateElementVector().
void calculateElementVector(const ElInfo *, ElementVector&);
+
+ protected:
+ const BasisFunction *psi, *phi;
};
@@ -104,6 +107,8 @@ namespace AMDiS {
protected:
std::vector<VectorOfFixVecs<DimVec<double> > > tmpGrdPhi;
+
+ const BasisFunction *psi, *phi;
};
diff --git a/AMDiS/src/FixVec.h b/AMDiS/src/FixVec.h
index 14058f778eb9b8c8cb8fc317aecc44b9a8a782db..a785ac338472c4a7a286376664e8fb8b86efa959 100644
--- a/AMDiS/src/FixVec.h
+++ b/AMDiS/src/FixVec.h
@@ -233,17 +233,20 @@ namespace AMDiS {
}
/// Returns the \ref size of this VectorOfFixVecs
- inline int getSize() const {
+ inline int getSize() const
+ {
return size;
}
/// Returns \ref dim
- inline int getDim() const {
+ inline int getDim() const
+ {
return dim;
}
/// Returns the size of the contained FixVecs
- inline int getSizeOfFixVec() const {
+ inline int getSizeOfFixVec() const
+ {
return vec[0]->getSize();
}
@@ -515,7 +518,8 @@ namespace AMDiS {
/// creates and inits and double array
double *createAndInitArray(int size, ...);
- inline WorldVector<double> operator*(const WorldVector<double>& v, double d) {
+ inline WorldVector<double> operator*(const WorldVector<double>& v, double d)
+ {
WorldVector<double> result = v;
result *= d;
return result;
diff --git a/AMDiS/src/FixVec.hh b/AMDiS/src/FixVec.hh
index c46166b2cfb85df4ad9ea19408300be4b0ed749a..7ccf9b0f60b2dca615b246921cf473b5caf1cfc2 100644
--- a/AMDiS/src/FixVec.hh
+++ b/AMDiS/src/FixVec.hh
@@ -8,9 +8,8 @@ namespace AMDiS {
thisIt != this->end();
thisIt++) {
*thisIt = 0;
- for (T* vIt = v.begin(); vIt != v.end(); vIt++, mIt++) {
+ for (T* vIt = v.begin(); vIt != v.end(); vIt++, mIt++)
*thisIt += *vIt * *mIt;
- }
}
}
@@ -19,9 +18,8 @@ namespace AMDiS {
{
double erg = 0.0;
- for (int i = 0; i < a.getSize() ; i++) {
+ for (int i = 0; i < a.getSize() ; i++)
erg = erg + ((a[i] - b[i]) * (a[i] - b[i]));
- }
return sqrt(erg);
}
@@ -39,9 +37,8 @@ namespace AMDiS {
thisIt++) {
*thisIt = 0;
- for (T* vIt = v.begin(); vIt != v.end(); vIt++, mIt++) {
- *thisIt += *vIt * *mIt;
- }
+ for (T* vIt = v.begin(); vIt != v.end(); vIt++, mIt++)
+ *thisIt += *vIt * *mIt;
}
}
@@ -70,9 +67,8 @@ namespace AMDiS {
template<typename T>
void WorldMatrix<T>::setDiag(T value)
{
- for (int i = 0; i < this->rows; i++) {
+ for (int i = 0; i < this->rows; i++)
this->valArray[i * this->cols + i] = value;
- }
}
template<typename T>
@@ -85,19 +81,16 @@ namespace AMDiS {
T *thisIt = this->begin();
- for (T* v1It = v1.begin(); v1It != v1.end(); v1It++) {
- for (T* v2It = v2.begin(); v2It != v2.end(); v2It++, thisIt++) {
+ for (T* v1It = v1.begin(); v1It != v1.end(); v1It++)
+ for (T* v2It = v2.begin(); v2It != v2.end(); v2It++, thisIt++)
*thisIt = *v1It * *v2It;
- }
- }
}
template<typename T>
const WorldMatrix<T>& operator*=(WorldMatrix<T>& m, T scal)
{
- for (T* mIt = m.begin(); mIt != m.end(); mIt++) {
+ for (T* mIt = m.begin(); mIt != m.end(); mIt++)
*mIt *= scal;
- }
return m;
}
@@ -108,9 +101,8 @@ namespace AMDiS {
T *m1It, *m2It;
for (m1It = m1.begin(), m2It = m2.begin();
m1It != m1.end();
- m1It++, m2It++) {
+ m1It++, m2It++)
*m1It -= *m2It;
- }
return m1;
}
@@ -121,9 +113,8 @@ namespace AMDiS {
T* m1It, *m2It;
for (m1It = m1.begin(), m2It = m2.begin();
m1It != m1.end();
- m1It++, m2It++) {
+ m1It++, m2It++)
*m1It += *m2It;
- }
return m1;
}
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index 034a17783f4a537bdd521d83159dddf37c0653a3..99a7f352f766255e60fde7412444b853c875bada 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -921,9 +921,13 @@ namespace AMDiS {
if (indices) {
result = indices;
} else {
- if (localVec)
+ if (localVec && nBasFcts > localVecSize) {
delete [] localVec;
- localVec = new DegreeOfFreedom[nBasFcts];
+ localVec = new DegreeOfFreedom[nBasFcts];
+ }
+ if (!localVec)
+ localVec = new DegreeOfFreedom[nBasFcts];
+
localVecSize = nBasFcts;
result = localVec;
}
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 20f676dce03dd08e2ac664a4ff0c2a6b3ebb1518..537664e0ce16b35051538e6a07bce74de1d0613b 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -578,6 +578,94 @@ namespace AMDiS {
}
}
+
+ Vec2AndGrad2AtQP_ZOT::Vec2AndGrad2AtQP_ZOT(DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ QuartAbstractFunction<double, double, double, WorldVector<double>,WorldVector<double> > *af)
+ : ZeroOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), f(af)
+ {
+ TEST_EXIT(dv1)("No first vector!\n");
+ TEST_EXIT(dv2)("No second vector!\n");
+
+ auxFESpaces.push_back(dv1->getFESpace());
+ auxFESpaces.push_back(dv2->getFESpace());
+ }
+
+ void Vec2AndGrad2AtQP_ZOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ vecAtQPs1 = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
+ vecAtQPs2 = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
+ gradAtQPs1 = getGradientsAtQPs(vec1, elInfo, subAssembler, quad);
+ gradAtQPs2 = getGradientsAtQPs(vec2, elInfo, subAssembler, quad);
+ }
+
+
+ void Vec2AndGrad2AtQP_ZOT::eval(int nPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double fac) const
+ {
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] +=
+ fac *
+ (*f)(vecAtQPs1[iq], vecAtQPs2[iq], gradAtQPs1[iq], gradAtQPs2[iq]) *
+ uhAtQP[iq];
+ }
+
+ void Vec2AndGrad2AtQP_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const
+ {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], gradAtQPs1[iq], gradAtQPs2[iq]);
+ }
+
+ Vec2AndGradVecAtQP_ZOT::Vec2AndGradVecAtQP_ZOT(DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ DOFVectorBase<double> *dGrd,
+ TertiaryAbstractFunction<double, double,double, WorldVector<double> > *af)
+ : ZeroOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), vecGrd(dGrd), f(af)
+ {
+ TEST_EXIT(dv1)("No vector!\n");
+ TEST_EXIT(dv2)("No vector!\n");
+ TEST_EXIT(dGrd)("No gradient vector!\n");
+
+ auxFESpaces.push_back(dv1->getFESpace());
+ auxFESpaces.push_back(dv2->getFESpace());
+ auxFESpaces.push_back(dGrd->getFESpace());
+ }
+
+ void Vec2AndGradVecAtQP_ZOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ vec1AtQPs = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
+ vec2AtQPs = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
+ gradAtQPs = getGradientsAtQPs(vecGrd, elInfo, subAssembler, quad);
+ }
+
+ void Vec2AndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints,
+ std::vector<double> &C) const
+ {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs[iq]);
+ }
+
+ void Vec2AndGradVecAtQP_ZOT::eval(int nPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double fac) const
+ {
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] +=
+ fac * (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs[iq]) * uhAtQP[iq];
+ }
+
General_ZOT::General_ZOT(std::vector<DOFVectorBase<double>*> vecs,
std::vector<DOFVectorBase<double>*> grads,
TertiaryAbstractFunction<double, WorldVector<double>, std::vector<double>, std::vector<WorldVector<double> > > *af)
@@ -663,6 +751,181 @@ namespace AMDiS {
auxFESpaces.push_back(dv2->getFESpace());
}
+ Vec2AndGradAtQP_FOT::Vec2AndGradAtQP_FOT(DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f_, WorldVector<double> *b_)
+ : FirstOrderTerm(8),
+ vec1(dv1),
+ vec2(dv2),
+ f(f_),
+ b(b_)
+ {
+ auxFESpaces.push_back(dv1->getFESpace());
+ auxFESpaces.push_back(dv2->getFESpace());
+ }
+
+ void Vec2AndGradAtQP_FOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ vec1AtQPs = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
+ vec2AtQPs = subAssembler->getVectorAtQPs(vec2, elInfo, quad);
+ gradAtQPs1 = getGradientsAtQPs(vec1, elInfo, subAssembler, quad);
+ }
+
+
+ void Vec2AndGradAtQP_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], (*f)(vec1AtQPs[iq], vec2AtQPs[iq],gradAtQPs1[iq]));
+ else
+ l1(Lambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs1[iq]));
+ }
+ }
+
+ void Vec2AndGradAtQP_FOT::eval(int nPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double fac) const
+ {
+ if (grdUhAtQP)
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += fac *
+ (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs1[iq]) *
+ ((*b) * grdUhAtQP[iq]);
+ }
+
+ FctVecAtQP_FOT::FctVecAtQP_FOT(DOFVectorBase<double> *dv,
+ AbstractFunction<double, WorldVector<double> > *f_,
+ WorldVector<double> *b_)
+ : FirstOrderTerm(f_->getDegree()), vec(dv), f(f_), b(b_)
+ {
+ auxFESpaces.push_back(dv->getFESpace());
+ }
+
+
+ void FctVecAtQP_FOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ vecAtQPs = subAssembler->getVectorAtQPs(vec, elInfo, quad);
+ coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
+ }
+
+
+ void FctVecAtQP_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], (*f)(coordsAtQPs[iq]) * (vecAtQPs[iq]));
+ else
+ l1(Lambda, Lb[iq], (*f)(coordsAtQPs[iq]) * (vecAtQPs[iq]));
+ }
+ }
+
+ void FctVecAtQP_FOT::eval(int nPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double fac) const
+ {
+ if (grdUhAtQP)
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] +=
+ fac * (*f)(coordsAtQPs[iq]) * (vecAtQPs[iq]) * ((*b) * grdUhAtQP[iq]);
+ }
+
+ Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ WorldVector<double> *b_)
+ : FirstOrderTerm(8), vec1(dv1), vec2(dv2), b(b_)
+ {
+ auxFESpaces.push_back(dv1->getFESpace());
+ auxFESpaces.push_back(dv2->getFESpace());
+ }
+
+ void Vec2AtQP_FOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ vec1AtQPs = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
+ vec2AtQPs = subAssembler->getVectorAtQPs(vec2, elInfo, quad);
+ }
+
+ void Vec2AtQP_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])*(vec2AtQPs[iq]));
+ else
+ l1(Lambda, Lb[iq], (vec1AtQPs[iq])*(vec2AtQPs[iq]));
+ }
+ }
+
+ void Vec2AtQP_FOT::eval(int nPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double fac) const
+ {
+ if (grdUhAtQP)
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += fac * (vec1AtQPs[iq]) * (vec2AtQPs[iq]) * ((*b) * grdUhAtQP[iq]);
+ }
+
+ 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_)
+ {
+ 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)
+ {
+ vec1AtQPs = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
+ vec2AtQPs = subAssembler->getVectorAtQPs(vec2, elInfo, quad);
+ vec3AtQPs = subAssembler->getVectorAtQPs(vec3, elInfo, quad);
+ }
+
+
+ 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]));
+ }
+ }
+
+ void Vec3FctAtQP_FOT::eval(int nPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double fac) const
+ {
+ if (grdUhAtQP)
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += fac *
+ (vec1AtQPs[iq]) * (*f)(vec2AtQPs[iq] ,vec3AtQPs[iq]) *
+ ((*b) * grdUhAtQP[iq]);
+ }
+
+
General_FOT::General_FOT(std::vector<DOFVectorBase<double>*> vecs,
std::vector<DOFVectorBase<double>*> grads,
TertiaryAbstractFunction<WorldVector<double>,
@@ -1580,17 +1843,13 @@ namespace AMDiS {
WorldMatrix<double> A = (*matrixFct)(vecAtQPs[iq]);
- if (D2UhAtQP) {
- for (int i = 0; i < dow; i++) {
- for (int j = 0; j < dow; j++) {
+ if (D2UhAtQP)
+ for (int i = 0; i < dow; i++)
+ for (int j = 0; j < dow; j++)
resultQP += A[i][j] * D2UhAtQP[iq][j][i];
- }
- }
- }
- if (grdUhAtQP) {
+ if (grdUhAtQP)
resultQP += (*divFct)(A) * grdUhAtQP[iq];
- }
result[iq] += resultQP * factor;
}
@@ -2594,12 +2853,10 @@ namespace AMDiS {
coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);
- for (int i = 0; i < nVecs; i++) {
+ for (int i = 0; i < nVecs; i++)
vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
- }
- for (int i = 0; i < nGrads; i++) {
+ for (int i = 0; i < nGrads; i++)
gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
- }
}
void General_ZOT::getC(const ElInfo *, int nPoints,
@@ -2612,12 +2869,12 @@ namespace AMDiS {
std::vector<WorldVector<double> > gradsArg(nGrads);
for (int iq = 0; iq < nPoints; iq++) {
- for (int i = 0; i < nVecs; i++) {
+ for (int i = 0; i < nVecs; i++)
vecsArg[i] = vecsAtQPs_[i][iq];
- }
- for (int i = 0; i < nGrads; i++) {
+ for (int i = 0; i < nGrads; i++)
+
gradsArg[i] = gradsAtQPs_[i][iq];
- }
+
C[iq] += (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
}
}
@@ -2976,7 +3233,7 @@ namespace AMDiS {
}
}
- void VecGrad_SOT::eval(int numPoints,
+ void VecGrad_SOT::eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
@@ -2986,7 +3243,7 @@ namespace AMDiS {
int dow = Global::getGeo(WORLD);
if (D2UhAtQP) {
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
double resultQP = 0.0;
for (int i = 0; i < dow; i++) {
@@ -2997,28 +3254,28 @@ namespace AMDiS {
}
}
- void VecGrad_SOT::weakEval(int numPoints,
+ void VecGrad_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
{
if (grdUhAtQP) {
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
}
- void VecGrad_SOT::getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const {
+ void VecGrad_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
}
}
void VecGrad_SOT::initElement(const ElInfo* elInfo,
- SubAssembler* subAssembler,
- Quadrature *quad)
+ SubAssembler* subAssembler,
+ Quadrature *quad)
{
vecAtQPs = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
gradAtQPs = subAssembler->getGradientsAtQPs(vec2, elInfo, quad);
@@ -3032,14 +3289,13 @@ namespace AMDiS {
grad2AtQPs = subAssembler->getGradientsAtQPs(vec2, elInfo, quad);
}
- void FctGrad2_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
+ void FctGrad2_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
lb(Lambda, (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]), Lb[iq], 1.0);
- }
}
- void FctGrad2_FOT::eval(int numPoints,
+ void FctGrad2_FOT::eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
@@ -3047,7 +3303,7 @@ namespace AMDiS {
double factor) const
{
if (grdUhAtQP) {
- for (int iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
WorldVector<double> b = (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]);
result[iq] += b * grdUhAtQP[iq] * factor;
}
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 5d4cd4555c14dce9012d6eea748b44c8987524c0..6dfc72c57820578e78c372d0a5169a1391de094d 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -1872,7 +1872,7 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> >& Lb) const;
/// Implements FirstOrderTerm::eval().
- void eval(int numPoints,
+ void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
@@ -1897,6 +1897,128 @@ namespace AMDiS {
};
+ class Vec2AndGradAtQP_FOT : public FirstOrderTerm
+ {
+ public:
+
+ Vec2AndGradAtQP_FOT(DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f_, WorldVector<double> *b_);
+
+ 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;
+
+ protected:
+ DOFVectorBase<double>* vec1;
+ DOFVectorBase<double>* vec2;
+
+ double *vec1AtQPs;
+ double *vec2AtQPs;
+ WorldVector<double> *gradAtQPs1;
+
+ TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f;
+ WorldVector<double> *b;
+ };
+
+
+ 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,
+ 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;
+
+ protected:
+ DOFVectorBase<double>* vec;
+ double *vecAtQPs;
+ WorldVector<double> *coordsAtQPs;
+ AbstractFunction<double, WorldVector<double> > *f;
+ WorldVector<double> *b;
+ };
+
+
+ class Vec2AtQP_FOT : public FirstOrderTerm
+ {
+ public:
+
+ Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ WorldVector<double> *b_);
+
+ 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;
+
+ protected:
+ DOFVectorBase<double>* vec1;
+ DOFVectorBase<double>* vec2;
+ double *vec1AtQPs;
+ double *vec2AtQPs;
+ WorldVector<double> *b;
+ };
+
+
+ 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_);
+
+ 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;
+
+ protected:
+ DOFVectorBase<double>* vec1;
+ DOFVectorBase<double>* vec2;
+ DOFVectorBase<double>* vec3;
+ BinaryAbstractFunction<double, double, double>* f;
+ double *vec1AtQPs;
+ double *vec2AtQPs;
+ double *vec3AtQPs;
+ WorldVector<double> *b;
+ };
+
+
class General_FOT : public FirstOrderTerm
{
public:
@@ -1908,22 +2030,16 @@ namespace AMDiS {
std::vector<double>,
std::vector<WorldVector<double> > > *f);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo*,
SubAssembler* ,
Quadrature *quad= NULL);
- /** \brief
- * Implements FirstOrderTerm::getLb().
- */
+ /// Implements FirstOrderTerm::getLb().
void getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& result) const;
- /** \brief
- * Implenetation of FirstOrderTerm::eval().
- */
+ /// Implenetation of FirstOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -2051,9 +2167,7 @@ namespace AMDiS {
C[iq] += factor;
}
- /** \brief
- * Implements ZeroOrderTerm::eval().
- */
+ /// Implements ZeroOrderTerm::eval().
inline void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *,
@@ -2066,9 +2180,7 @@ namespace AMDiS {
}
protected:
- /** \brief
- * Constant factor of zero order term.
- */
+ /// Constant factor of zero order term.
double factor;
};
@@ -2704,18 +2816,18 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- inline void getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const;
+ inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
/// Implements SecondOrderTerm::eval().
- void eval(int numPoints,
- const double *uhAtQP,
+ void eval(int nPoints,
+ const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
double factor) const;
/// Implements SecondOrderTerm::weakEval().
- void weakEval(int numPoints,
+ void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const;
@@ -3204,32 +3316,90 @@ namespace AMDiS {
double fac) const;
protected:
- /** \brief
- * DOFVector to be evaluated at quadrature points.
- */
+ /// DOFVector to be evaluated at quadrature points.
DOFVector<double>* vec;
- /** \brief
- * Vector v at quadrature points.
- */
+ /// Vector v at quadrature points.
double *vecAtQPs;
- /** \brief
- * Vector of DOFVectors to be evaluated at quadrature points.
- */
+ /// Vector of DOFVectors to be evaluated at quadrature points.
std::vector<DOFVector<double>*> vecs;
- /** \brief
- * Vectors at quadrature points.
- */
+ /// Vectors at quadrature points.
std::vector<WorldVector<double>*> gradsAtQPs;
- /** \brief
- * Function for c.
- */
+ /// Function for c.
BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *f;
};
+
+ class Vec2AndGrad2AtQP_ZOT : public ZeroOrderTerm
+ {
+ public:
+ Vec2AndGrad2AtQP_ZOT(DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ QuartAbstractFunction<double, double, double, WorldVector<double>, WorldVector<double> > *af);
+
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+
+ void eval(int nPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double fac) const;
+
+ protected:
+ DOFVectorBase<double>* vec1;
+ DOFVectorBase<double>* vec2;
+
+ double *vecAtQPs1;
+ double *vecAtQPs2;
+
+ WorldVector<double> *gradAtQPs1;
+ WorldVector<double> *gradAtQPs2;
+
+ QuartAbstractFunction<double, double, double, WorldVector<double>,WorldVector<double> > *f;
+ };
+
+
+ class Vec2AndGradVecAtQP_ZOT : public ZeroOrderTerm
+ {
+ public:
+ Vec2AndGradVecAtQP_ZOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ DOFVectorBase<double> *dGrd,
+ TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f);
+
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+
+ void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
+
+ void eval(int nPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double fac) const;
+
+ protected:
+ DOFVectorBase<double>* vec1;
+ DOFVectorBase<double>* vec2;
+
+ double *vec1AtQPs;
+ double *vec2AtQPs;
+
+ DOFVectorBase<double>* vecGrd;
+ WorldVector<double> *gradAtQPs;
+
+ TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f;
+ };
+
+
+
class General_ZOT : public ZeroOrderTerm
{
public:
@@ -3238,20 +3408,14 @@ namespace AMDiS {
std::vector<DOFVectorBase<double>*> grads,
TertiaryAbstractFunction<double, WorldVector<double>, std::vector<double>, std::vector<WorldVector<double> > > *f);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements ZeroOrderTerm::getC().
- */
+ /// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
- /** \brief
- * Implements ZeroOrderTerm::eval().
- */
+ /// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 1127a4cbe19fe51ceffaee25a8a8b0368885a1cf..b0ba62000086c0c94f487f6a64aaa4221b52bc1f 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -574,8 +574,6 @@ namespace AMDiS {
{
FUNCNAME("ProblemVec::refineMesh()");
- // if (adaptInfo->getTimestepNumber() == 9) return 0;
-
int nMeshes = static_cast<int>(meshes.size());
Flag refineFlag = 0;
for (int i = 0; i < nMeshes; i++)
@@ -603,18 +601,16 @@ namespace AMDiS {
FUNCNAME("ProblemVec::oneIteration()");
if (allowFirstRef) {
- for (int i = 0; i < nComponents; i++) {
+ for (int i = 0; i < nComponents; i++)
adaptInfo->allowRefinement(true, i);
- }
+
allowFirstRef = false;
} else {
- for (int i = 0; i < nComponents; i++) {
- if (adaptInfo->spaceToleranceReached(i)) {
+ for (int i = 0; i < nComponents; i++)
+ if (adaptInfo->spaceToleranceReached(i))
adaptInfo->allowRefinement(false, i);
- } else {
+ else
adaptInfo->allowRefinement(true, i);
- }
- }
}
return StandardProblemIteration::oneIteration(adaptInfo, toDo);
@@ -658,20 +654,13 @@ namespace AMDiS {
DOFMatrix* dof_matrix= (*systemMatrix)[i][j];
DOFMatrix::base_matrix_type& base_matrix= dof_matrix->getBaseMatrix();
- int nnz_per_row= 0;
- if (num_rows(base_matrix) != 0)
- nnz_per_row= int(double(base_matrix.nnz()) / num_rows(base_matrix) * 1.2);
- if (nnz_per_row < 5)
- nnz_per_row= 5;
+ dof_matrix->calculateNnz();
// Correct dimensionality of matrix
base_matrix.change_dim(componentSpaces[i]->getAdmin()->getUsedSize(),
componentSpaces[j]->getAdmin()->getUsedSize());
set_to_zero(base_matrix);
-
- // Reuse old sparsity information (if available) or default
- // dof_matrix->startInsertion(nnz_per_row);
}
}
}
@@ -993,6 +982,9 @@ namespace AMDiS {
TraverseStack stack;
#endif
+ // == Initialize matrix and vector. If we have to assemble in parallel, ===
+ // == temporary thread owned matrix and vector must be created. ===
+
#ifdef _OPENMP
#pragma omp parallel
#endif
@@ -1005,6 +997,7 @@ namespace AMDiS {
DOFMatrix *tmpMatrix = NULL;
DOFVector<double> *tmpVector = NULL;
+#ifdef _OPENMP
if (matrix) {
tmpMatrix = new DOFMatrix(matrix->getRowFESpace(), matrix->getColFESpace(), "tmp");
@@ -1016,7 +1009,7 @@ namespace AMDiS {
tmpMatrix->getBaseMatrix().change_dim(matrix->getRowFESpace()->getAdmin()->getUsedSize(),
matrix->getColFESpace()->getAdmin()->getUsedSize());
- tmpMatrix->startInsertion();
+ tmpMatrix->startInsertion(10);
}
if (vector) {
@@ -1028,6 +1021,20 @@ namespace AMDiS {
*tmpVector = *vector;
tmpVector->set(0.0);
}
+#else
+ if (matrix) {
+ tmpMatrix = matrix;
+ tmpMatrix->startInsertion(matrix->getNnz());
+ }
+ if (vector) {
+ tmpVector = vector;
+ tmpVector->set(0.0);
+ }
+#endif
+
+
+ // == Traverse mesh (either sequentially or in parallel) and assemble. ==
+
// Because we are using the parallel traverse stack, each thread will
// traverse only a part of the mesh.
@@ -1059,6 +1066,11 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
+
+ // == Finally, if we have assembled in parallel, we have to add the thread ==
+ // == private matrix and vector to the global one. ==
+
+#ifdef _OPENMP
tmpMatrix->finishInsertion();
// After mesh traverse, all thread have to added their private matrices and
@@ -1067,19 +1079,13 @@ namespace AMDiS {
// the same time.
if (matrix) {
-#ifdef _OPENMP
#pragma omp critical
-#endif
matrix->getBaseMatrix() += tmpMatrix->getBaseMatrix();
}
-#ifdef _OPENMP
#pragma omp barrier
-#endif
-#ifdef _OPENMP
#pragma omp master
-#endif
{
if (matrix)
matrix->startInsertion();
@@ -1087,25 +1093,27 @@ namespace AMDiS {
if (matrix) {
// Remove rows corresponding to DOFs on a Dirichlet boundary.
-#ifdef _OPENMP
#pragma omp critical
-#endif
matrix->removeRowsWithDBC(tmpMatrix->getApplyDBCs());
delete tmpMatrix;
}
if (vector) {
-#ifdef _OPENMP
#pragma omp critical
-#endif
*vector += *tmpVector;
delete tmpVector;
}
+#else
+ if (matrix)
+ matrix->removeRowsWithDBC(matrix->getApplyDBCs());
+#endif
+
if (useGetBound)
delete [] bound;
+
} // pragma omp parallel
}
diff --git a/AMDiS/src/Quadrature.h b/AMDiS/src/Quadrature.h
index de6b738137be3028968cd86029cf3978e04b15e8..74090b2685c4ce6645207d3b412bd3ba4939f00a 100644
--- a/AMDiS/src/Quadrature.h
+++ b/AMDiS/src/Quadrature.h
@@ -46,9 +46,7 @@ namespace AMDiS {
class Quadrature
{
protected:
- /** \brief
- * Avoids call of default constructor
- */
+ /// Avoids call of default constructor
Quadrature();
/** \brief
@@ -93,43 +91,39 @@ namespace AMDiS {
*/
double integrateStdSimplex(AbstractFunction<double, DimVec<double> > *f);
- /** \brief
- * Returns \ref name
- */
- inline const std::string& getName() {
+ /// Returns \ref name
+ inline const std::string& getName()
+ {
return name;
}
/// Returns \ref n_points
- inline int getNumPoints() const {
+ inline int getNumPoints() const
+ {
return n_points;
}
- /** \brief
- * Returns \ref w[p]
- */
- inline double getWeight(int p) const {
+ /// Returns \ref w[p]
+ inline double getWeight(int p) const
+ {
return w[p];
}
- /** \brief
- * Returns \ref w.
- */
- inline double* getWeight() const {
+ /// Returns \ref w.
+ inline double* getWeight() const
+ {
return w;
}
- /** \brief
- * Returns \ref dim
- */
- inline int getDim() const {
+ /// Returns \ref dim
+ inline int getDim() const
+ {
return dim;
}
- /** \brief
- * Returns \ref degree
- */
- inline int getDegree() const {
+ /// Returns \ref degree
+ inline int getDegree() const
+ {
return degree;
}
@@ -162,7 +156,8 @@ namespace AMDiS {
* Returns \ref lambda[a][b] which is the b-th coordinate entry of the a-th
* quadrature point
*/
- inline double getLambda(int a, int b) const {
+ inline double getLambda(int a, int b) const
+ {
return (lambda ? (*lambda)[a][b] : 0.0);
}
@@ -170,51 +165,39 @@ namespace AMDiS {
* Returns \ref lambda[a] which is a DimVec<double> containing the
* coordiantes of the a-th quadrature point
*/
- inline const DimVec<double>& getLambda(int a) const {
+ inline const DimVec<double>& getLambda(int a) const
+ {
return (*lambda)[a];
}
- /** \brief
- * Returns \ref lambda which is a VectorOfFixvecs<DimVec<double> >.
- */
- VectorOfFixVecs<DimVec<double> > *getLambda() const {
+ /// Returns \ref lambda which is a VectorOfFixvecs<DimVec<double> >.
+ VectorOfFixVecs<DimVec<double> > *getLambda() const
+ {
return lambda;
}
public:
- /** \brief
- * Maximal number of quadrature points for the different dimensions
- */
+ /// Maximal number of quadrature points for the different dimensions
static const int maxNQuadPoints[4];
protected:
- /** \brief
- * Name of this Quadrature
- */
+ /// Name of this Quadrature
std::string name;
- /** \brief
- * Quadrature is exact of this degree
- */
+ /// Quadrature is exact of this degree
int degree;
- /** \brief
- * Quadrature for dimension dim
- */
+ /// Quadrature for dimension dim
int dim;
/// Number of quadrature points
int n_points;
- /** \brief
- * Vector of quadrature points given in barycentric coordinates
- */
+ /// Vector of quadrature points given in barycentric coordinates
VectorOfFixVecs<DimVec<double> > *lambda;
- /** \brief
- * Vector of quadrature weights
- */
+ /// Vector of quadrature weights
double *w;
protected:
@@ -298,9 +281,7 @@ namespace AMDiS {
/** \} */
};
- // ============================================================================
- // ===== class FastQuadrature =================================================
- // ============================================================================
+
/** \brief
* Pre-compute the values of all basis functions at all quadrature nodes;
@@ -319,7 +300,6 @@ namespace AMDiS {
*/
const Flag INIT_D2_PHI=4;
- // ============================================================================
/**
* \ingroup Integration
@@ -355,142 +335,112 @@ namespace AMDiS {
basisFunctions(basFcts)
{}
- /** \brief
- * Copy constructor
- */
+ /// Copy constructor
FastQuadrature(const FastQuadrature&);
- /** \brief
- * Extended copy constructor
- */
+ /// Extended copy constructor
FastQuadrature(const FastQuadrature&, const Flag);
- /** \brief
- * Destructor
- */
+ /// Destructor
~FastQuadrature();
public:
- /** \brief
- * Returns a FastQuadrature for the given BasisFunction, Quadrature, and
- * flags
- */
+ /// Returns a FastQuadrature for the given BasisFunction, Quadrature, and flags.
static FastQuadrature* provideFastQuadrature(const BasisFunction*,
const Quadrature&,
Flag);
- /** \brief
- * inits FastQuadrature like speciefied in flag
- */
+ /// inits FastQuadrature like speciefied in flag
void init(Flag init_flag);
- inline bool initialized(Flag flag) {
- if (flag == INIT_PHI) {
+ inline bool initialized(Flag flag)
+ {
+ if (flag == INIT_PHI)
return (phi != NULL);
- }
- if (flag == INIT_GRD_PHI) {
+ if (flag == INIT_GRD_PHI)
return (grdPhi != NULL);
- }
- if (flag == INIT_D2_PHI) {
+ if (flag == INIT_D2_PHI)
return (D2Phi != NULL);
- }
ERROR_EXIT("invalid flag\n");
return false;
}
- /** \brief
- * Returns \ref quadrature
- */
- inline const Quadrature* getQuadrature() const {
+ /// Returns \ref quadrature
+ inline const Quadrature* getQuadrature() const
+ {
return quadrature;
}
- /** \brief
- * Returns \ref max_points
- */
- inline int getMaxQuadPoints() {
+ /// Returns \ref max_points
+ inline int getMaxQuadPoints()
+ {
return max_points;
}
- /** \brief
- * Returns (*\ref D2Phi)[q][i][j][m]
- */
+ /// Returns (*\ref D2Phi)[q][i][j][m]
const double getSecDer(int q, int i, int j, int m) const;
- /** \brief
- * Returns (*\ref D2Phi)[q]
- */
+ /// Returns (*\ref D2Phi)[q]
const VectorOfFixVecs<DimMat<double> > *getSecDer(int q) const;
- /** \brief
- * Returns (*\ref grdPhi)[q][i][j]
- */
- inline const double getGradient(int q, int i ,int j) const {
+ /// Returns (*\ref grdPhi)[q][i][j]
+ inline const double getGradient(int q, int i ,int j) const
+ {
return (grdPhi) ? (*grdPhi)[q][i][j] : 0.0;
}
- /** \brief
- * Returns (*\ref grdPhi)[q]
- */
- inline VectorOfFixVecs<DimVec<double> >* getGradient(int q) const {
+ /// Returns (*\ref grdPhi)[q]
+ inline VectorOfFixVecs<DimVec<double> >* getGradient(int q) const
+ {
return (grdPhi) ? &((*grdPhi)[q]) : NULL;
}
- /** \brief
- * Returns \ref phi[q][i]
- */
- inline const double getPhi(int q, int i) const {
+ /// Returns \ref phi[q][i]
+ inline const double getPhi(int q, int i) const
+ {
return (phi) ? phi[q][i] : 0.0;
}
- /** \brief
- * Returns \ref phi[q]
- */
- inline const double *getPhi(int q) const {
+ /// Returns \ref phi[q]
+ inline const double *getPhi(int q) const
+ {
return (phi) ? phi[q] : NULL;
}
- /** \brief
- * Returns \ref quadrature ->integrateStdSimplex(f)
- */
- inline double integrateStdSimplex(AbstractFunction<double, DimVec<double> > *f) {
+ /// Returns \ref quadrature ->integrateStdSimplex(f)
+ inline double integrateStdSimplex(AbstractFunction<double, DimVec<double> > *f)
+ {
return quadrature->integrateStdSimplex(f);
}
- /** \brief
- * Returns \ref quadrature ->getNumPoints()
- */
- inline int getNumPoints() const {
+ /// Returns \ref quadrature ->getNumPoints()
+ inline int getNumPoints() const
+ {
return quadrature->getNumPoints();
}
- /** \brief
- * Returns \ref quadrature ->getWeight(p)
- */
- inline double getWeight(int p) const {
+ /// Returns \ref quadrature ->getWeight(p)
+ inline double getWeight(int p) const
+ {
return quadrature->getWeight(p);
}
- /** \brief
- * Returns \ref quadrature ->getDim()
- */
- inline int getDim() const {
+ /// Returns \ref quadrature ->getDim()
+ inline int getDim() const
+ {
return quadrature->getDim();
}
- /** \brief
- * Returns \ref quadrature ->getDegree()
- */
- inline int getDegree() const {
+ /// Returns \ref quadrature ->getDegree()
+ inline int getDegree() const
+ {
return quadrature->getDegree();
}
- /** \brief
- * Returns \ref quadrature ->grdFAtQp(f, vec)
- */
+ /// Returns \ref quadrature ->grdFAtQp(f, vec)
inline const WorldVector<double>
*grdFAtQp(const AbstractFunction<WorldVector<double>,
DimVec<double> >& f,
@@ -499,33 +449,28 @@ namespace AMDiS {
return quadrature->grdFAtQp(f, vec);
}
- /** \brief
- * Returns \ref quadrature ->fAtQp(f, vec)
- */
+ /// Returns \ref quadrature ->fAtQp(f, vec)
inline const double *fAtQp(const AbstractFunction<double,
DimVec<double> >& f, double *vec) const
{
return quadrature->fAtQp(f, vec);
}
- /** \brief
- * Returns \ref quadrature ->getLambda(a,b)
- */
- inline double getLambda(int a,int b) const {
+ /// Returns \ref quadrature ->getLambda(a,b)
+ inline double getLambda(int a,int b) const
+ {
return quadrature->getLambda(a,b);
}
- /** \brief
- * Returns \ref quadrature ->getLambda(a)
- */
- inline const DimVec<double>& getLambda(int a) const {
+ /// Returns \ref quadrature ->getLambda(a)
+ inline const DimVec<double>& getLambda(int a) const
+ {
return quadrature->getLambda(a);
}
- /** \brief
- * Returns \ref basisFunctions
- */
- inline BasisFunction* getBasisFunctions() const {
+ /// Returns \ref basisFunctions
+ inline BasisFunction* getBasisFunctions() const
+ {
return basisFunctions;
}
@@ -574,14 +519,10 @@ namespace AMDiS {
*/
static int max_points;
- /** \brief
- * This FastQuadrature stores values for Quadrature quadrature
- */
+ /// This FastQuadrature stores values for Quadrature quadrature
Quadrature* quadrature;
- /** \brief
- * Values stored for basis functions basisFunctions
- */
+ /// Values stored for basis functions basisFunctions
BasisFunction* basisFunctions;
};
diff --git a/AMDiS/src/RCNeighbourList.h b/AMDiS/src/RCNeighbourList.h
index 797c96065d70c7ec4cca18bc55a5b6aa64bf4ff5..c5a9cac0ac9f7752bfd6d57a194cf7c83de39a7d 100644
--- a/AMDiS/src/RCNeighbourList.h
+++ b/AMDiS/src/RCNeighbourList.h
@@ -57,26 +57,29 @@ namespace AMDiS {
virtual ~RCNeighbourList();
/// Sets flag of \ref rclist[i] = true
- inline void setCoarsePatch(int i) {
+ inline void setCoarsePatch(int i)
+ {
rclist[i]->flag = true;
}
/// Sets flag of \ref rclist[i] = f
- inline void setCoarsePatch(int i, bool f) {
+ inline void setCoarsePatch(int i, bool f)
+ {
rclist[i]->flag = f;
}
/// Returns \ref rclist[i].flag
- inline const bool isPatchCoarse(int i) const {
+ inline const bool isPatchCoarse(int i) const
+ {
return rclist[i]->flag;
}
/** \brief
* If \ref rclist[i].neigh[j] is not a NULL pointer
- * \ref rclist[i].neigh[j]->no will be returned. Otherwise the return value
- * is -1
+ * \ref rclist[i].neigh[j]->no will be returned. Otherwise the return value is -1
*/
- inline int getNeighbourNr(int i, int j) const {
+ inline int getNeighbourNr(int i, int j) const
+ {
return (rclist[i]->neigh[j])?rclist[i]->neigh[j]->no:-1;
}
@@ -85,30 +88,35 @@ namespace AMDiS {
* \ref rclist[i].neigh[j]->el will be returned. Otherwise the return value
* is NULL
*/
- inline Element* getNeighbourElement(int i, int j) const {
+ inline Element* getNeighbourElement(int i, int j) const
+ {
return rclist[i]->neigh[j]? rclist[i]->neigh[j]->el : NULL;
}
/// Returns \ref rclist[i].el
- inline Element* getElement(int i) const {
+ inline Element* getElement(int i) const
+ {
if (static_cast<int>(rclist.size()) <= i)
return NULL;
return rclist[i]->el;
}
/// Sets \ref rclist[i].el to el and \ref rclist[i].flag to cp.
- inline const Element* setElement(int i, const Element* el, bool cp = false) {
+ inline const Element* setElement(int i, const Element* el, bool cp = false)
+ {
rclist[i]->el = const_cast<Element*>(el);
rclist[i]->flag = cp;
return el;
}
/// Returns \ref rclist[i].elType
- inline int getType(int i) const {
+ inline int getType(int i) const
+ {
return rclist[i]->elType;
}
- inline void setType(int i, int type) const {
+ inline void setType(int i, int type) const
+ {
rclist[i]->elType = type;
}
@@ -116,22 +124,26 @@ namespace AMDiS {
virtual bool doCoarsePatch(int n_neigh);
/// Sets \ref rclist[i].oppVertex[j] = k
- inline void setOppVertex(int i,int j,int k) {
+ inline void setOppVertex(int i,int j,int k)
+ {
rclist[i]->oppVertex[j] = k;
}
/// Returns \ref rclist[i].oppVertex[j]
- inline int getOppVertex(int i, int j) {
+ inline int getOppVertex(int i, int j)
+ {
return rclist[i]->oppVertex[j];
}
/// Sets \ref rclist[i].elType = t
- inline void setElType(int i,unsigned char t) {
+ inline void setElType(int i,unsigned char t)
+ {
rclist[i]->elType = t;
}
/// Sets \ref coarseningManager = cm
- inline void setCoarseningManager(CoarseningManager *cm) {
+ inline void setCoarseningManager(CoarseningManager *cm)
+ {
coarseningManager = cm;
}
diff --git a/AMDiS/src/SubAssembler.cc b/AMDiS/src/SubAssembler.cc
index 16273256269b32d5767c0937acc708d793516936..ba8d370653f11c99dc95ebc85f5a0fba8af4474c 100644
--- a/AMDiS/src/SubAssembler.cc
+++ b/AMDiS/src/SubAssembler.cc
@@ -114,7 +114,7 @@ namespace AMDiS {
}
}
- WorldVector<double>* SubAssembler::getCoordsAtQPs(const ElInfo* elInfo,
+ WorldVector<double>* SubAssembler::getCoordsAtQPs(const ElInfo* elInfo,
Quadrature *quad)
{
Quadrature *localQuad = quad ? quad : quadrature;
@@ -122,9 +122,8 @@ namespace AMDiS {
const int nPoints = localQuad->getNumPoints();
// already calculated for this element ?
- if (coordsValid) {
+ if (coordsValid)
return coordsAtQPs;
- }
if (coordsAtQPs) {
if (coordsNumAllocated != nPoints) {
@@ -139,9 +138,8 @@ namespace AMDiS {
// set new values
WorldVector<double>* k = &(coordsAtQPs[0]);
- for (int l = 0; k < &(coordsAtQPs[nPoints]); ++k, ++l) {
+ for (int l = 0; k < &(coordsAtQPs[nPoints]); ++k, ++l)
elInfo->coordToWorld(localQuad->getLambda(l), *k);
- }
// mark values as valid
coordsValid = true;
@@ -166,9 +164,9 @@ namespace AMDiS {
if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid)
return valuesAtQPs[vec]->values.getValArray();
- if (!valuesAtQPs[vec]) {
+ if (!valuesAtQPs[vec])
valuesAtQPs[vec] = new ValuesAtQPs;
- }
+
valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());
double *values = valuesAtQPs[vec]->values.getValArray();
@@ -219,14 +217,11 @@ namespace AMDiS {
TEST_EXIT_DBG(vec)("no dof vector!\n");
- // if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid)
- // return valuesAtQPs[vec]->values.getValArray();
-
Quadrature *localQuad = quad ? quad : quadrature;
- if (!valuesAtQPs[vec]) {
+ if (!valuesAtQPs[vec])
valuesAtQPs[vec] = new ValuesAtQPs;
- }
+
valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());
double *values = valuesAtQPs[vec]->values.getValArray();
diff --git a/AMDiS/src/SubAssembler.h b/AMDiS/src/SubAssembler.h
index 97a3e0f3c19eb6103025baaf70afce507f5eb883..4f5ae9bedf832350b821858d2e403d4a876f7f63 100644
--- a/AMDiS/src/SubAssembler.h
+++ b/AMDiS/src/SubAssembler.h
@@ -185,7 +185,7 @@ namespace AMDiS {
/// Used for \ref getCoordsAtQP(). Stores the number of allocated WorldVectors.
int coordsNumAllocated;
- /// Needed Quadrature. Constructed in the constructor of SubAssembler
+ /// Quadrature object to be used for assembling.
Quadrature *quadrature;
/// FastQuadrature for row basis functions