diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc index 190881742e717054a100ebed5400dc1b442d5621..239ce2ee5d2218dd719a415669ca384f7a788afd 100644 --- a/AMDiS/src/Operator.cc +++ b/AMDiS/src/Operator.cc @@ -2101,9 +2101,8 @@ 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++) { - resultQP += D2UhAtQP[iq][i][i]; - } + for (int i = 0; i < dow; i++) + resultQP += D2UhAtQP[iq][i][i]; result[iq] += factor * f * resultQP; } } @@ -2134,9 +2133,8 @@ namespace AMDiS { for (int iq = 0; iq < nPoints; iq++) { double factor = (*f)(gradAtQPs[iq]); double resultQP = 0.0; - for (int i = 0; i < dow; i++) { - resultQP += D2UhAtQP[iq][i][i]; - } + for (int i = 0; i < dow; i++) + resultQP += D2UhAtQP[iq][i][i]; result[iq] += fac * factor * resultQP; } } @@ -2167,9 +2165,8 @@ namespace AMDiS { 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++) { - resultQP += D2UhAtQP[iq][i][i]; - } + for (int i = 0; i < dow; i++) + resultQP += D2UhAtQP[iq][i][i]; result[iq] += fac * factor * resultQP; } } @@ -2201,17 +2198,13 @@ namespace AMDiS { WorldMatrix<double> A = (*f)(vecAtQPs[iq], gradAtQPs[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) { - resultQP += (*divFct)(A) * grdUhAtQP[iq]; - } + if (grdUhAtQP) + resultQP += (*divFct)(A) * grdUhAtQP[iq]; result[iq] += resultQP * factor; } @@ -2243,9 +2236,8 @@ namespace AMDiS { for (int iq = 0; iq < nPoints; iq++) { double factor = (*f)(vecAtQPs[iq], coordsAtQPs[iq]); double resultQP = 0.0; - for (int i = 0; i < dow; i++) { - resultQP += D2UhAtQP[iq][i][i]; - } + for (int i = 0; i < dow; i++) + resultQP += D2UhAtQP[iq][i][i]; result[iq] += fac * factor * resultQP; } } @@ -2277,18 +2269,14 @@ namespace AMDiS { WorldMatrix<double> A = (*f)(gradAtQPs[iq], coordsAtQPs[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; } } @@ -2319,9 +2307,8 @@ namespace AMDiS { for (int iq = 0; iq < nPoints; iq++) { double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]); double resultQP = 0.0; - for (int i = 0; i < dow; i++) { - resultQP += D2UhAtQP[iq][i][i]; - } + for (int i = 0; i < dow; i++) + resultQP += D2UhAtQP[iq][i][i]; result[iq] += fac * factor * resultQP; } } @@ -2352,9 +2339,8 @@ namespace AMDiS { for (int iq = 0; iq < nPoints; iq++) { double factor = (*f)(vecAtQPs[iq]); double resultQP = 0.0; - for (int i = 0; i < dow; i++) { - resultQP += grdUhAtQP[iq][i]; - } + for (int i = 0; i < dow; i++) + resultQP += grdUhAtQP[iq][i]; result[iq] += fac * factor * resultQP; } } @@ -2373,9 +2359,8 @@ 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++) { - resultQP += grdUhAtQP[iq][i]; - } + for (int i = 0; i < dow; i++) + resultQP += grdUhAtQP[iq][i]; result[iq] += f * factor * resultQP; } } @@ -2451,8 +2436,7 @@ namespace AMDiS { double resultQP = 0.0; const WorldVector<double> &b = (*g)(coordsAtQPs[iq]); for (int i = 0; i < dow; i++) - resultQP += b[i] * grdUhAtQP[iq][i]; - + resultQP += b[i] * grdUhAtQP[iq][i]; result[iq] += f * resultQP; } } @@ -2497,7 +2481,7 @@ namespace AMDiS { if (D2UhAtQP) { for (int iq = 0; iq < nPoints; iq++) { double factor = (*g)(coordsAtQPs[iq]); - result[iq] = D2UhAtQP[iq][xi][xj] * factor * f; + result[iq] += D2UhAtQP[iq][xi][xj] * factor * f; } } } @@ -2541,7 +2525,7 @@ namespace AMDiS { if (D2UhAtQP) { for (int iq = 0; iq < nPoints; iq++) { double factor = (*f)(vecAtQPs[iq]); - result[iq] = D2UhAtQP[iq][xi][xj] * factor * fac; + result[iq] += D2UhAtQP[iq][xi][xj] * factor * fac; } } } @@ -3171,13 +3155,12 @@ namespace AMDiS { std::vector<WorldVector<double> > gradsArg(nGrads); for (int iq = 0; iq < nPoints; iq++) { - for (int i = 0; i < nVecs; i++) { - vecsArg[i] = vecsAtQPs_[i][iq]; - } - for (int i = 0; i < nGrads; i++) { - gradsArg[i] = gradsAtQPs_[i][iq]; - } - result[iq] += fac * (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg) * uhAtQP[iq]; + for (int i = 0; i < nVecs; i++) + vecsArg[i] = vecsAtQPs_[i][iq]; + for (int i = 0; i < nGrads; i++) + gradsArg[i] = gradsAtQPs_[i][iq]; + result[iq] += + fac * (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg) * uhAtQP[iq]; } } @@ -3264,11 +3247,11 @@ namespace AMDiS { } void VecGrad_SOT::eval(int nPoints, - const double *uhAtQP, - const WorldVector<double> *grdUhAtQP, - const WorldMatrix<double> *D2UhAtQP, - double *result, - double fac) const + const double *uhAtQP, + const WorldVector<double> *grdUhAtQP, + const WorldMatrix<double> *D2UhAtQP, + double *result, + double fac) const { int dow = Global::getGeo(WORLD); @@ -3276,17 +3259,16 @@ namespace AMDiS { 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++) { - resultQP += D2UhAtQP[iq][i][i]; - } + for (int i = 0; i < dow; i++) + resultQP += D2UhAtQP[iq][i][i]; result[iq] += fac * resultQP * factor; } } } void VecGrad_SOT::weakEval(int nPoints, - const WorldVector<double> *grdUhAtQP, - WorldVector<double> *result) const + const WorldVector<double> *grdUhAtQP, + WorldVector<double> *result) const { if (grdUhAtQP) { for (int iq = 0; iq < nPoints; iq++) { @@ -3296,11 +3278,12 @@ namespace AMDiS { } } - void VecGrad_SOT::getLALt(const ElInfo *elInfo, int nPoints, 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 < nPoints; iq++) { - l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq])); - } + for (int iq = 0; iq < nPoints; iq++) + l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq])); } void VecGrad_SOT::initElement(const ElInfo* elInfo, @@ -3319,18 +3302,20 @@ namespace AMDiS { grad2AtQPs = subAssembler->getGradientsAtQPs(vec2, elInfo, quad); } - void FctGrad2_FOT::getLb(const ElInfo *elInfo, int nPoints, 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 < nPoints; iq++) lb(Lambda, (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]), Lb[iq], 1.0); } void FctGrad2_FOT::eval(int nPoints, - const double *uhAtQP, - const WorldVector<double> *grdUhAtQP, - const WorldMatrix<double> *D2UhAtQP, - double *result, - double factor) const + const double *uhAtQP, + const WorldVector<double> *grdUhAtQP, + const WorldMatrix<double> *D2UhAtQP, + double *result, + double factor) const { if (grdUhAtQP) { for (int iq = 0; iq < nPoints; iq++) {