Skip to content

GitLab

Commit 3918b92b authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Bugfix in IJ_SOT operators.

parent f1e48be6
Hide whitespace changes
Inline Side-by-side
AMDiS/src/Operator.cc
View file @ 3918b92b
......@@ -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++) {
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment