Select Git revision
dg.hh 19.28 KiB
#ifndef DUNE_MMDG_DG_HH
#define DUNE_MMDG_DG_HH
#include <cmath>
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/solver.hh>
#include <dune/istl/umfpack.hh>
#include <dune/mmdg/nonconformingp1vtkfunction.hh>
template<class GridView, class Mapper, class Problem>
class DG
{
public:
static constexpr int dim = GridView::dimension;
static constexpr auto simplex = Dune::GeometryTypes::simplex(dim);
static constexpr auto edge = Dune::GeometryTypes::simplex(dim-1);
using Scalar = typename GridView::ctype;
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
using VTKFunction = Dune::VTKFunction<GridView>;
using P1Function = Dune::NonconformingP1VTKFunction<GridView,
Dune::DynamicVector<Scalar>>;
using QRuleElem = Dune::QuadratureRule<Scalar, dim>;
using QRulesElem = Dune::QuadratureRules<Scalar, dim>;
using QRuleEdge = Dune::QuadratureRule<Scalar, dim-1>;
using QRulesEdge = Dune::QuadratureRules<Scalar, dim-1>;
//constructor
DG (const GridView& gridView, const Mapper& mapper, const Problem& problem) :
gridView_(gridView), mapper_(mapper), problem_(problem),
dof_((1 + dim) * gridView.size(0))
{
//initialize stiffness matrix A, load vector b and solution vector d
A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
b = Vector(dof_);
d = Vector(dof_);
}
void operator() (const Scalar mu)
{
//assemble stiffness matrix A and load vector b
assembleSLE(mu);
A->compress();
Dune::InverseOperatorResult result;
Dune::UMFPack<Matrix> solver(*A);
solver.apply(d, b, result);
//write solution to a vtk file
writeVTKoutput();
}
//compute L2 error using a quadrature rule
Scalar computeL2error(const int quadratureOrder = 5) const
{
if (!problem_.hasExactSolution())
{
DUNE_THROW(Dune::Exception, "ERROR: Problem has no exact solution or "
"exact solution is not implemented"); }
Scalar error(0.0);
for (const auto& elem : elements(gridView_))
{
const int elemIdxSLE = (dim + 1) * mapper_.index(elem);
const auto& geo = elem.geometry();
const QRuleElem& qrule = QRulesElem::rule(simplex, quadratureOrder);
//determine L2-error on element elem using a quadrature rule
for (const auto& ip : qrule)
{
const auto& qpLocal = ip.position();
const auto& qpGlobal = geo.global(qpLocal);
const Scalar quadratureFactor =
ip.weight() * geo.integrationElement(qpLocal);
Scalar numericalSolutionAtQP = d[elemIdxSLE];
for (int i = 0; i < dim; i++)
{
numericalSolutionAtQP += d[elemIdxSLE + i + 1] * qpGlobal[i];
}
const Scalar errorEvaluation =
numericalSolutionAtQP - problem_.exactSolution(qpGlobal);
error += quadratureFactor * errorEvaluation * errorEvaluation;
}
}
return std::sqrt(error);
}
const GridView& gridView_;
const Mapper& mapper_; //element mapper for gridView_
const Problem& problem_; //the DG problem to be solved
const int dof_; //degrees of freedom
protected:
DG (const GridView& gridView, const Mapper& mapper, const Problem& problem,
const int dof) :
gridView_(gridView), mapper_(mapper), problem_(problem), dof_(dof)
{
//initialize stiffness matrix A, load vector b and solution vector d
A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
b = Vector(dof_);
d = Vector(dof_);
}
//assemble stiffness matrix A and load vector b
void assembleSLE (const Scalar mu)
{
//quadrature rules
const QRuleElem& rule_q =
QRulesElem::rule(simplex, problem_.quadratureOrder_q());
const QRuleElem& rule_K_Elem =
QRulesElem::rule(simplex, problem_.quadratureOrder_K());
const QRuleEdge& rule_K_Edge =
QRulesEdge::rule(edge, problem_.quadratureOrder_K());
const QRuleEdge& rule_Kplus1_Edge =
QRulesEdge::rule(edge, problem_.quadratureOrder_K() + 1);
const QRuleEdge& rule_2ndOrder = QRulesEdge::rule(edge, 2);
const QRuleEdge& rule_boundary =
QRulesEdge::rule(edge, problem_.quadratureOrderBoundary());
//we use the basis // phi_elem,0 (x) = indicator(elem);
// phi_elem,i (x) = x[i]*indicator(elem);
//for elem in elements(gridView) and i = 1,...,dim
for (const auto& elem : elements(gridView_))
{
const int elemIdx = mapper_.index(elem);
const auto& geo = elem.geometry();
//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;
//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 quadratureFactor =
ip.weight() * geo.integrationElement(qpLocal);
const Scalar qEvaluation(problem_.q(qpGlobal));
//quadrature for int_elem q*phi_elem,0 dV
b[elemIdxSLE] += quadratureFactor * qEvaluation;
//quadrature for int_elem q*phi_elem,i dV
for (int i = 0; i < dim; i++)
{
b[elemIdxSLE + i + 1] += quadratureFactor * qpGlobal[i] * qEvaluation;
}
}
//evaluate
// int_elem K*grad(phi_elem,i)*grad(phi_elem,j) dV
// = int_elem K[j][i] dV
//using a quadrature rule
for (const auto& ip : rule_K_Elem)
{
//quadrature point in local and global coordinates
const auto& qpLocal = ip.position();
const auto& qpGlobal = geo.global(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) +=
quadratureFactor * problem_.K(qpGlobal)[j][i];
}
}
}
//iterate over all intersection with the boundary of elem
for (const auto& intersct : intersections(gridView_, elem))
{
if constexpr (dim != 1)
{
if (gridView_.grid().isInterface(intersct))
{
continue;
}
}
const auto& normal = intersct.centerUnitOuterNormal(); const auto& intersctGeo = intersct.geometry();
const double intersctVol = intersctGeo.volume();
const auto& intersctCenter = intersctGeo.center();
//create storage for multiply used integral values
//linearIntegrals[i] = int_intersct x[i] ds
Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);
//quadraticIntregrals[i][j] = int_intersct x[i]*x[j] ds
//NOTE: is there a better type for a symmetric matrix?
//note that only the lower diagonal and diagonal entries of
//quadraticIntregrals are used
Dune::FieldMatrix<Scalar, dim, dim> quadraticIntregrals(0.0);
//KIntegrals[i] = int_intersct K[:,i]*normal ds
Dune::FieldVector<Scalar, dim> KIntegrals(0.0);
//KIntegralsX[i][j] = int_intersct x[j]*(K[:,i]*normal) ds
Dune::FieldMatrix<Scalar, dim, dim> KIntegralsX(0.0);
//fill vector linearIntegrals and matrix quadraticIntregrals
for (const auto& ip : rule_2ndOrder)
{
//quadrature point in local and global coordinates
const auto& qpLocal = ip.position();
const auto& qpGlobal = intersctGeo.global(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 //TODO
linearIntegrals[i] = intersctVol * intersctCenter[i];
for (int j = 0; j <= i; j++)
{
quadraticIntregrals[i][j] +=
quadratureFactor * qpGlobal[i] * qpGlobal[j];
}
}
}
//quadraticIntregrals is symmetric
for (int i = 0; i < dim; i++)
{
for (int j = 0; j < i; j++)
{
quadraticIntregrals[j][i] = quadraticIntregrals[i][j];
}
}
//fill vector KIntegrals using a quadrature rule
for (const auto& ip : rule_K_Edge)
{
//quadrature point in local and global coordinates
const auto& qpLocal = ip.position();
const auto& qpGlobal = intersctGeo.global(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] += quadratureFactor * (KEvaluation[i] * normal);
}
}
//fill vector KIntegralsX using a quadrature rule
for (const auto& ip : rule_Kplus1_Edge)
{
//quadrature point in local and global coordinates
const auto& qpLocal = ip.position();
const auto& qpGlobal = intersctGeo.global(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++)
{
for (int j = 0; j < dim; j++)
{
KIntegralsX[i][j] += quadratureFactor *
qpGlobal[j] * (KEvaluation[i] * normal);
}
}
}
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
A->entry(elemIdxSLE, elemIdxSLE) += mu * intersctVol;
if (intersct.neighbor()) //intersct has neighboring element
{
//index of the neighboring element
const int neighborIdx = mapper_.index(intersct.outside());
const int neighborIdxSLE = (dim + 1) * neighborIdx;
for (int i = 0; i < dim; i++)
{ //we use the relations
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
// = mu * int_intersct x[i] ds
//and
// int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds
// = 0.5 * int_intersct K[:,i] * normal ds
A->entry(elemIdxSLE + i + 1, elemIdxSLE) +=
mu * linearIntegrals[i] - 0.5 * KIntegrals[i];
for (int j = 0; j <= i; j++)
{
//we use the relations
// int_intersct mu*jump(phi_elem,i)*jump(phi_elem,j) ds
// = mu * int_intersct x[i] * x[j] ds
//and
// int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,j) ds
// = 0.5 * int_intersct x[j] * ( K[:,i] * normal ) ds
A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1)
+= mu * quadraticIntregrals[i][j]
- 0.5 * ( KIntegralsX[i][j] + KIntegralsX[j][i] );
}
}
if (neighborIdx > elemIdx)
{ //make sure that each facet is considered only once
continue;
}
//exact evaluation of
// int_intersct mu * jump(phi_elem,0) * jump(phi_neighbor,0) ds
A->entry(elemIdxSLE, neighborIdxSLE) += -mu * intersctVol;
//stiffness matrix A is symmetric
A->entry(neighborIdxSLE, elemIdxSLE) +=
A->entry(elemIdxSLE, neighborIdxSLE);
for (int i = 0; i < dim; i++)
{
//we use the relations
// int_intersct mu * jump(phi_elem,i) * jump(phi_neighbor,0) ds
// = -mu * int_intersct x[i] ds
//and
// int_intersct avg(K*grad(phi_neighbor,i)) * jump(phi_elem,0) ds
// = -0.5 * int_intersct K[:,i] * normal ds
A->entry(elemIdxSLE + i + 1, neighborIdxSLE) +=
-mu * linearIntegrals[i] + 0.5 * KIntegrals[i];
//we use the relations
// int_intersct mu * jump(phi_elem,0) * jump(phi_neighbor,i) ds
// = -mu * int_intersct x[i] ds
//and
// int_intersct avg(K*grad(phi_neighbor,i)) * jump(phi_elem,0) ds
// = 0.5 * int_intersct K[:,i] * normal ds
A->entry(elemIdxSLE, neighborIdxSLE + i + 1) +=
-mu * linearIntegrals[i] - 0.5 * KIntegrals[i];
//stiffness matrix A is symmetric
A->entry(neighborIdxSLE, elemIdxSLE + i + 1) +=
A->entry(elemIdxSLE + i + 1, neighborIdxSLE);
A->entry(neighborIdxSLE + i + 1, elemIdxSLE) +=
A->entry(elemIdxSLE, neighborIdxSLE + i + 1);
for (int j = 0; j <= i; j++)
{
//we use the relations
// int_intersct mu * jump(phi_elem,i) * jump(phi_neighbor,j) ds
// = -mu * int_intersct x[i] * x[j] ds
//and
// int_intersct avg(K*grad(phi_neighbor,j)) * jump(phi_elem,i) ds
// = 0.5 * int_intersct x[i] * ( K[:,j] * normal ) ds
//as well as
// int_intersct avg(K*grad(phi_elem,i)) * jump(phi_neighbor,j) ds
// = -0.5 * int_intersct x[j] * ( K[:,i] * normal ) ds
A->entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1) +=
-mu * quadraticIntregrals[i][j]
- 0.5 * ( KIntegralsX[j][i] - KIntegralsX[i][j] );
//stiffness matrix A is symmetric
A->entry(neighborIdxSLE + j + 1, elemIdxSLE + i + 1) +=
A->entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1);
if (i != j)
{
//we use the relations
// int_intersct mu * jump(phi_elem,j) * jump(phi_neighbor,i) ds
// = -mu * int_intersct x[i] * x[j] ds
//and
// int_intersct avg(K*grad(phi_neighbor,i))*jump(phi_elem,j) ds
// = 0.5 * int_intersct x[j] * ( K[:,i] * normal ) ds
//as well as
// int_intersct avg(K*grad(phi_elem,j))*jump(phi_neighbor,i) ds
// = -0.5 * int_intersct x[i] * ( K[:,j] * normal ) ds
A->entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1) +=
-mu * quadraticIntregrals[i][j]
- 0.5 * ( KIntegralsX[i][j] - KIntegralsX[j][i] );
//stiffness matrix A is symmetric
A->entry(neighborIdxSLE + i + 1, elemIdxSLE + j + 1) +=
A->entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1);
}
}
}
}
else //boundary facet
{
//fill load vector b using a quadrature rule (boundary terms) for (const auto& ip : rule_boundary)
{
//quadrature point in local and global coordinates
const auto& qpLocal = ip.position();
const auto& qpGlobal = intersctGeo.global(qpLocal);
const Scalar quadratureFactor =
ip.weight() * intersctGeo.integrationElement(qpLocal);
const Scalar boundaryEvaluation(problem_.boundary(qpGlobal));
const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
problem_.K(qpGlobal);
//quadrature for int_intersct mu * g * phi_elem,0 ds
b[elemIdxSLE] += quadratureFactor * mu * boundaryEvaluation;
//quadrature for
// int_intersct (mu * phi_elem,i - K[:,i] * normal) * g ds
//with K*grad(phi_elem,i) = K[:,i]
for (int i = 0; i < dim; i++)
{
b[elemIdxSLE + i + 1] +=
quadratureFactor * (mu * qpGlobal[i] -
(KEvaluation[i] * normal)) * boundaryEvaluation;
}
}
for (int i = 0; i < dim; i++)
{
//we evaluate the integrals
// int_intersct mu * phi_elem,0 * phi_elem,i ds
//and
// int_intersct phi_elem,0 * K*grad(phi_elem,i) * normal ds
//with K*grad(phi_elem,i) = K[:,i]
A->entry(elemIdxSLE + i + 1, elemIdxSLE) +=
mu * linearIntegrals[i] - KIntegrals[i];
for (int j = 0; j <= i; j++)
{
//we evaluate the integrals
// int_intersct mu * phi_elem,i * phi_elem,j ds
//and
// int_intersct (phi_elem,i * K[:,j]
// + phi_elem,j * K[:,i]) * normal ds
//with K*grad(phi_elem,i) = K[:,i]
A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
mu * quadraticIntregrals[i][j]
- ( KIntegralsX[i][j] + KIntegralsX[j][i] );
}
}
}
}
//stiffness matrix A is symmetric
for (int i = 0; i < dim; i++)
{
A->entry(elemIdxSLE, elemIdxSLE + i + 1) =
A->entry(elemIdxSLE + i + 1, elemIdxSLE);
for(int j = 0; j < i; j++)
{
A->entry(elemIdxSLE + j + 1, elemIdxSLE + i + 1) =
A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1);
}
}
}
}
//writes the solution to a vtk file
void writeVTKoutput (const int bulkDOF) const {
//storage for pressure data, for each element we store the
//pressure at the corners of the element, the VTKFunction will
//create the output using a linear interpolation for each element
Dune::DynamicVector<Scalar> pressure(bulkDOF, 0.0);
Dune::DynamicVector<Scalar> exactPressure(bulkDOF, 0.0);
for (const auto& elem : elements(gridView_))
{
const int elemIdxSLE = (dim + 1)*mapper_.index(elem);
const auto& geo = elem.geometry();
for (int k = 0; k < geo.corners(); k++)
{
if (problem_.hasExactSolution())
{
exactPressure[elemIdxSLE + k] = problem_.exactSolution(geo.corner(k));
}
//contribution of the basis function
// phi_elem,0 (x) = indicator(elem);
//at the kth corner of elem
pressure[elemIdxSLE + k] = d[elemIdxSLE];
for (int i = 0; i < dim; i++)
{
//contribution of the basis function
// phi_elem,i (x) = x[i]*indicator(elem);
//at the kth corner of elem
pressure[elemIdxSLE + k] += d[elemIdxSLE + i + 1] * geo.corner(k)[i];
}
}
}
//create output directory if necessary
mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
Dune::VTKWriter<GridView> vtkWriter(gridView_, Dune::VTK::nonconforming);
vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
new P1Function(gridView_, pressure, "pressure")));
if (problem_.hasExactSolution())
{
vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
new P1Function(gridView_, exactPressure, "exactPressure")));
}
vtkWriter.pwrite("pressureData", "", "data");
}
std::shared_ptr<Matrix> A; //stiffness matrix
Vector b; //load vector
Vector d; //solution vector
private:
//writes the solution to a vtk file
void writeVTKoutput () const
{
writeVTKoutput(dof_);
}
};
#endif