Skip to content
GitLab
Commit bbd936aa authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Fixed a bug for FOT operators when using 2d meshes in 3d.

parent 005e656e
Show whitespace changes
Inline Side-by-side
AMDiS/src/FixVec.h
View file @ bbd936aa
......@@ -277,7 +277,8 @@ namespace AMDiS {
* VectorOfFixVecs constructors.
*/
MatrixOfFixVecs(int dim, int r, int c, InitType initType)
: rows(r), columns(c)
: rows(r),
columns(c)
{
TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or wrong initializer\n");
vec = new VectorOfFixVecs<FixVecType>*[rows];
......@@ -309,12 +310,14 @@ namespace AMDiS {
}
/// Returns \ref rows
inline int getNumberOfRows() const {
inline int getNumberOfRows() const
{
return rows;
}
/// Returns \ref columns
inline int getNumberOfColumns() const {
inline int getNumberOfColumns() const
{
return columns;
}
......@@ -337,7 +340,8 @@ namespace AMDiS {
* parts of an element.
*/
template<typename T>
class DimVec : public FixVec<T,PARTS> {
class DimVec : public FixVec<T,PARTS>
{
public:
DimVec() {}
......@@ -549,8 +553,7 @@ namespace AMDiS {
return result;
}
inline bool operator<(const WorldVector<double>& v1,
const WorldVector<double>& v2)
inline bool operator<(const WorldVector<double>& v1, const WorldVector<double>& v2)
{
int dow = Global::getGeo(WORLD);
for (int i = 0; i < dow; i++) {
......@@ -561,8 +564,7 @@ namespace AMDiS {
return false;
}
inline bool operator==(const WorldVector<double>& v1,
const WorldVector<double>& v2)
inline bool operator==(const WorldVector<double>& v1, const WorldVector<double>& v2)
{
int dow = Global::getGeo(WORLD);
for (int i = 0; i < dow; i++)
......
AMDiS/src/Operator.cc
View file @ bbd936aa
......@@ -53,12 +53,11 @@ namespace AMDiS {
void OperatorTerm::setSymmetric(bool symm)
{
if (symm) {
if (symm)
properties.setFlag(SYMMETRIC);
} else {
else
properties.unsetFlag(SYMMETRIC);
}
}
bool OperatorTerm::isSymmetric()
{
......@@ -145,24 +144,21 @@ namespace AMDiS {
// Both elements have the same size, so we can use the simple procedure
// to determine the gradients.
if (vec->getFESpace()->getMesh() == smallElInfo->getMesh()) {
if (vec->getFESpace()->getMesh() == smallElInfo->getMesh())
return subAssembler->getGradientsAtQPs(vec, smallElInfo, quad);
} else {
else
return subAssembler->getGradientsAtQPs(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 gradients.
if (vec->getFESpace()->getMesh() == largeElInfo->getMesh()) {
if (vec->getFESpace()->getMesh() == largeElInfo->getMesh())
return subAssembler->getGradientsAtQPs(vec, smallElInfo, largeElInfo, quad);
} else {
else
return subAssembler->getGradientsAtQPs(vec, smallElInfo, quad);
}
}
}
void OperatorTerm::lalt(const DimVec<WorldVector<double> >& Lambda,
......@@ -1423,82 +1419,102 @@ namespace AMDiS {
void CoordsAtQP_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad) {
Quadrature *quad)
{
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
void MatrixFct_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
void MatrixFct_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
lalt(Lambda, (*matrixFct)(vecAtQPs[iq]), *(LALt[iq]), symmetric, 1.0);
}
void VecAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
void VecAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq]));
}
void Vec2AtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
void Vec2AtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs1[iq], vecAtQPs2[iq]));
}
void CoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
void CoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, (*LALt[iq]), (*g)(coordsAtQPs[iq]));
}
void MatrixGradient_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
void MatrixGradient_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++)
lalt(Lambda, (*f)(gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
}
}
void FctGradient_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
void FctGradient_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*f)(gradAtQPs[iq]));
}
void VecAndGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const {
DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
}
void VecMatrixGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const {
DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
lalt(Lambda, (*f)(vecAtQPs[iq], gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
}
void VecAndCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt)
const { const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
void VecAndCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], coordsAtQPs[iq]));
}
void MatrixGradientAndCoords_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const {
DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
lalt(Lambda, (*f)(gradAtQPs[iq], coordsAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
}
void VecGradCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const {
DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]));
}
void VecAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
void VecAtQP_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)
......@@ -1508,60 +1524,61 @@ namespace AMDiS {
}
}
void CoordsAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
void CoordsAtQP_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++) {
for (int iq = 0; iq < nPoints; iq++)
l1(Lambda, Lb[iq], (*g)(coordsAtQPs[iq]));
}
}
void VecCoordsAtQP_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++) {
for (int iq = 0; iq < nPoints; iq++)
lb(Lambda, b, Lb[iq], (*g)(coordsAtQPs[iq]));
}
}
void VectorGradient_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
void VectorGradient_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
if (f) {
for (int iq = 0; iq < nPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++)
lb(Lambda, (*f)(gradAtQPs[iq]), Lb[iq], 1.0);
}
} else {
for (int iq = 0; iq < nPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++)
lb(Lambda, gradAtQPs[iq], Lb[iq], 1.0);
}
}
}
void VectorFct_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
void VectorFct_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++) {
for (int iq = 0; iq < nPoints; iq++)
lb(Lambda, (*vecFct)(vecAtQPs[iq]), Lb[iq], 1.0);
}
}
void VecFctAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
void VecFctAtQP_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++) {
for (int iq = 0; iq < nPoints; iq++)
lb(Lambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0);
}
}
void VecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const {
void VecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const
{
if (f) {
for (int iq = 0; iq < nPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vecAtQPs[iq]);
}
} else {
for (int iq = 0; iq < nPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++)
C[iq] += vecAtQPs[iq];
}
}
}
void VecAtQP_ZOT::eval(int nPoints,
const double *uhAtQP,
......@@ -1571,22 +1588,20 @@ namespace AMDiS {
double fac) const
{
if (f) {
for (int iq = 0; iq < nPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++)
result[iq] += fac * (*f)(vecAtQPs[iq]) * uhAtQP[iq];
}
} else {
for (int iq = 0; iq < nPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++)
result[iq] += fac * vecAtQPs[iq] * uhAtQP[iq];
}
}
}
void MultVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++) {
void MultVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]);
}
}
void MultVecAtQP_ZOT::eval(int nPoints,
const double *uhAtQP,
......@@ -1595,17 +1610,15 @@ namespace AMDiS {
double *result,
double fac) const
{
for (int iq = 0; iq < nPoints; iq++) {
result[iq] += fac * (*f1)(vecAtQPs1[iq]) *
(*f2)(vecAtQPs2[iq]) * uhAtQP[iq];
}
for (int iq = 0; iq < nPoints; iq++)
result[iq] += fac * (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]) * uhAtQP[iq];
}
void Vec2AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++) {
void Vec2AtQP_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]);
}
}
void Vec2AtQP_ZOT::eval(int nPoints,
const double *uhAtQP,
......@@ -1614,17 +1627,15 @@ namespace AMDiS {
double *result,
double fac) const
{
for (int iq = 0; iq < nPoints; iq++) {
result[iq] += fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq]) *
uhAtQP[iq];
}
for (int iq = 0; iq < nPoints; iq++)
result[iq] += fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq]) * uhAtQP[iq];
}
void Vec3AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++) {
void Vec3AtQP_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], vecAtQPs3[iq]);
}
}
void Vec3AtQP_ZOT::eval(int nPoints,
const double *uhAtQP,
......@@ -1633,21 +1644,18 @@ namespace AMDiS {
double *result,
double fac) const
{
for (int iq = 0; iq < nPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++)
result[iq] +=
fac *
(*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]) *
uhAtQP[iq];
}
fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]) * uhAtQP[iq];
}
void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints,
std::vector<double> &C) const {
for (int iq = 0; iq < nPoints; iq++) {
std::vector<double> &C) const
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
}
}
void VecAndCoordsAtQP_ZOT::eval(int nPoints,
const double *uhAtQP,
......@@ -2302,9 +2310,9 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*g)(coordsAtQPs[iq]);
double resultQP = 0.0;
for (int i = 0; i < dow; i++) {
for (int i = 0; i < dow; i++)
resultQP += grdUhAtQP[iq][i];
}
result[iq] += f * factor * resultQP;
}
}
......@@ -2324,12 +2332,11 @@ namespace AMDiS {
result[iq] += b * grdUhAtQP[iq] * factor;
}
} else {
for (int iq = 0; iq < nPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++)
result[iq] += gradAtQPs[iq] * grdUhAtQP[iq] * factor;
}
}
}
}
void VectorFct_FOT::eval(int nPoints,
const double *uhAtQP,
......@@ -2359,9 +2366,9 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
double resultQP = 0.0;
const WorldVector<double> &b = (*g)(coordsAtQPs[iq]);
for (int i = 0; i < dow; i++) {
for (int i = 0; i < dow; i++)
resultQP += b[i] * grdUhAtQP[iq][i];
}
result[iq] += f * resultQP;
}
}
......@@ -2369,9 +2376,8 @@ namespace AMDiS {
Assembler *Operator::getAssembler(int rank)
{
if (!assembler[rank]) {
if (!assembler[rank])
initAssembler(rank, NULL, NULL, NULL, NULL);
}
return assembler[rank];
}
......@@ -2393,10 +2399,9 @@ namespace AMDiS {
DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++)
lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*g)(coordsAtQPs[iq]));
}
}
void CoordsAtQP_IJ_SOT::eval(int nPoints,
const double *uhAtQP,
......@@ -2556,9 +2561,9 @@ namespace AMDiS {
std::vector<WorldVector<double>*> arg(size);
for (int iq = 0; iq < nPoints; iq++) {
for (int i = 0; i < size; i++) {
for (int i = 0; i < size; i++)
arg[i] = &(gradsAtQPs[i][iq]);
}
result[iq] += fac * (*f)(arg) * uhAtQP[iq];
}
}
......@@ -2624,7 +2629,9 @@ namespace AMDiS {
}
}
void VecAndGradAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
void VecAndGradAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
......@@ -2666,12 +2673,11 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
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];
}
lalt(Lambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg),
*(LALt[iq]), symmetric_, 1.0);
}
......@@ -2694,26 +2700,21 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
double resultQP = 0.0;
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];
}
WorldMatrix<double> A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
if (D2UhAtQP) {
for (int i = 0; i < dow; i++) {
for (int j = 0; j < dow; j++) {
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;
}
......@@ -2732,12 +2733,11 @@ namespace AMDiS {
if (grdUhAtQP) {
WorldMatrix<double> A;
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];
}
A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
result[iq] += A * grdUhAtQP[iq];
}
......@@ -2753,13 +2753,11 @@ 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_FOT::getLb(const ElInfo *elInfo,
int nPoints,
......@@ -2774,12 +2772,11 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
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];
}
lb(Lambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg),
result[iq], 1.0);
}
......@@ -2803,18 +2800,16 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
double resultQP = 0.0;
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];
}
const WorldVector<double> &b = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
for (int i = 0; i < dow; i++) {
for (int i = 0; i < dow; i++)
resultQP += grdUhAtQP[iq][i] * b[i];
}