Skip to content
Snippets Groups Projects
Commit 523c9fb7 authored by Hörl, Maximilian's avatar Hörl, Maximilian
Browse files

use quadratureFactor instead of weight*intElem in dg

parent 07a740ae
No related branches found
No related tags found
No related merge requests found
......@@ -111,7 +111,7 @@ public:
const Mapper& mapper_; //element mapper for gridView_
const Problem& problem_; //the DG problem to be solved
const int dof_;
const int dof_; //degrees of freedom
protected:
DG (const Grid& grid, const GridView& gridView, const Mapper& mapper,
......@@ -156,25 +156,24 @@ protected:
//basis function phi_elem,i
const int elemIdxSLE = (dim + 1)*elemIdx;
//fill load vector b using a quadrature rule (bulk term)
//fill load vector b using a quadrature rule
for (const auto& ip : rule_q)
{
//quadrature point in local and global coordinates
const auto& qpLocal = ip.position();
const auto& qpGlobal = geo.global(qpLocal);
const Scalar weight(ip.weight());
const Scalar integrationElem(geo.integrationElement(qpLocal));
const Scalar quadratureFactor =
ip.weight() * geo.integrationElement(qpLocal);
const Scalar qEvaluation(problem_.q(qpGlobal));
//quadrature for int_elem q*phi_elem,0 dV
b[elemIdxSLE] += weight * integrationElem * qEvaluation;
b[elemIdxSLE] += quadratureFactor * qEvaluation;
//quadrature for int_elem q*phi_elem,i dV
for (int i = 0; i < dim; i++)
{
b[elemIdxSLE + i + 1] +=
weight * integrationElem * qpGlobal[i] * qEvaluation;
b[elemIdxSLE + i + 1] += quadratureFactor * qpGlobal[i] * qEvaluation;
}
}
......@@ -188,15 +187,15 @@ protected:
const auto& qpLocal = ip.position();
const auto& qpGlobal = geo.global(qpLocal);
const Scalar weight(ip.weight());
const Scalar integrationElem(geo.integrationElement(qpLocal));
const Scalar quadratureFactor =
ip.weight() * geo.integrationElement(qpLocal);
for (int i = 0; i < dim; i++)
{
for (int j = 0; j <= i; j++)
{
A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
weight * integrationElem * problem_.K(qpGlobal)[j][i];
quadratureFactor * problem_.K(qpGlobal)[j][i];
}
}
}
......@@ -240,19 +239,19 @@ protected:
const auto& qpLocal = ip.position();
const auto& qpGlobal = intersctGeo.global(qpLocal);
const Scalar weight(ip.weight());
const Scalar integrationElem(intersctGeo.integrationElement(qpLocal));
const Scalar quadratureFactor =
ip.weight() * intersctGeo.integrationElement(qpLocal);
for (int i = 0; i < dim; i++)
{
//use midpoint rule for exact evaluation of
// int_intersct x[i] ds
// int_intersct x[i] ds //TODO
linearIntegrals[i] = intersctVol * intersctCenter[i];
for (int j = 0; j <= i; j++)
{
quadraticIntregrals[i][j] +=
weight * integrationElem * qpGlobal[i] * qpGlobal[j];
quadratureFactor * qpGlobal[i] * qpGlobal[j];
}
}
}
......@@ -273,15 +272,14 @@ protected:
const auto& qpLocal = ip.position();
const auto& qpGlobal = intersctGeo.global(qpLocal);
const Scalar weight(ip.weight());
const Scalar integrationElem(intersctGeo.integrationElement(qpLocal));
const Scalar quadratureFactor =
ip.weight() * intersctGeo.integrationElement(qpLocal);
const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
problem_.K(qpGlobal);
for (int i = 0; i < dim; i++)
{
KIntegrals[i] +=
weight * integrationElem * (KEvaluation[i] * normal);
KIntegrals[i] += quadratureFactor * (KEvaluation[i] * normal);
}
}
......@@ -292,8 +290,8 @@ protected:
const auto& qpLocal = ip.position();
const auto& qpGlobal = intersctGeo.global(qpLocal);
const Scalar weight(ip.weight());
const Scalar integrationElem(intersctGeo.integrationElement(qpLocal));
const Scalar quadratureFactor =
ip.weight() * intersctGeo.integrationElement(qpLocal);
const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
problem_.K(qpGlobal);
......@@ -301,7 +299,7 @@ protected:
{
for (int j = 0; j < dim; j++)
{
KIntegralsX[i][j] += weight * integrationElem *
KIntegralsX[i][j] += quadratureFactor *
qpGlobal[j] * (KEvaluation[i] * normal);
}
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment