From f39c2d1ae56619c5181556ebf0a23702def04c2f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Wed, 22 Jan 2020 19:38:12 +0100
Subject: [PATCH] use quadrature rule for load vector b
---
dune/mmdg/dg.hh | 38 +++++++++++++++++---------------------
1 file changed, 17 insertions(+), 21 deletions(-)
diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 28f71d8..014257a 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -26,7 +26,7 @@ public:
using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
using VTKFunction = Dune::VTKFunction<GridView>;
using P1Function = Dune::NonconformingP1VTKFunction<GridView,
- Dune::DynamicVector<double>>;
+ Dune::DynamicVector<Scalar>>;
//constructor
DG (const GridView& gridView, const Mapper& mapper,
@@ -63,6 +63,11 @@ private:
//assemble stiffness matrix A and load vector b
void assembleSLE (const Scalar K, const Scalar mu)
{
+ //quadrature rule
+ const Dune::QuadratureRule<Scalar, dim>& rule =
+ Dune::QuadratureRules<Scalar, dim>::rule(
+ Dune::GeometryTypes::simplex(dim), problem_.quadratureOrder());
+
//we use the basis
// phi_elem,0 (x) = indicator(elem);
// phi_elem,i (x) = x[i]*indicator(elem);
@@ -72,43 +77,34 @@ private:
const int elemIdx = mapper_.index(elem);
const auto& geo = elem.geometry();
const double elemVol = geo.volume();
- const auto& center = geo.center();
//in the system of linear equations (SLE) Ad = b,
//the index elemIdxSLE refers to the basis function phi_elem,0
//and the indices elemIdxSLE + i + 1, i = 1,...,dim, refer to the
//basis function phi_elem,i
const int elemIdxSLE = (dim + 1)*elemIdx;
-/*
- //TODO: can be done outside of the loop?
- const Dune::QuadratureRule<double,dim>& rule =
- Dune::QuadratureRules<double,dim>::rule(geo.type(), problem_.quadratureOrder());
- // Dune::FieldVector<double, dim+1> update(0.0);
- //NOTE: how are quadrature rules in Dune applied correctly?
+ //fill load vector b using a quadrature rule
for (const auto& ip : rule)
{
- const auto& qp = ip.position();
- //NOTE: is the volume of the element taken into account
- //automatically?
- const double weight = ip.weight();
+ //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 qEvalution(problem_.q(qpGlobal));
//quadrature for int_elem q*phi_elem,0 dV
- b[elemIdxSLE] += weight * problem_.q(geo.global(qp)) * geo.integrationElement(qp);
+ b[elemIdxSLE] += weight * qEvalution * integrationElem;
//quadrature for int_elem q*phi_elem,i dV
for (int i = 0; i < dim; i++)
{
- b[elemIdxSLE + i + 1] += weight * qp[i] * problem_.q(geo.global(qp)) * geo.integrationElement(qp);
+ b[elemIdxSLE + i + 1] +=
+ weight * qpGlobal[i] * qEvalution * integrationElem;
}
}
-*/
- //NOTE: makeshift solution for source term q = -1
- b[elemIdxSLE] += -elemVol;
- for (int i = 0; i < dim; i++)
- {
- b[elemIdxSLE + i + 1] += -elemVol * center[i];
- }
for (int i = 0; i < dim; i++)
{
--
GitLab