Select Git revision
Hörl, Maximilian authored
dg.hh 18.25 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/io/file/vtk.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/mmdg/nonconformingp1vtkfunction.hh>
template<class GridView, class Mapper, class Problem>
class DG
{
public:
using Scalar = typename GridView::ctype;
static constexpr int dim = GridView::dimension;
using Matrix = Dune::DynamicMatrix<Scalar>; //NOTE: what is an appropriate sparse matrix type? -> BCRS
using Vector = Dune::DynamicVector<Scalar>;
//constructor
DG (const GridView& gridView, const Mapper& mapper,
const Problem& problem) :
gridView_(gridView), mapper_(mapper), problem_(problem),
dof((1 + dim) * gridView.size(0))
{
A = Matrix(dof, dof, 0.0); //initialize stiffness matrix
b = Vector(dof, 0.0); //initialize load vector
d = Vector(dof, 0.0); //initialize solution vector
}
const void operator() (const Scalar K, const Scalar mu)
{
//assemble stiffness matrix A and load vector b
assembleSLE(K, mu);
//NOTE: what would be an appropiate solver here?
A.solve(d, b);
std::cout << b << std::endl;
for (int i = 0; i < dof; i++)
for (int j = 0; j < i; j++)
/*assert*/if(std::abs(A[i][j] - A[j][i]) >
std::numeric_limits<Scalar>::epsilon())
std::cout << i << ", " << j << std::endl;
//NOTE: dim = 1
// auto exactSolution = [](const auto& pos) {return 0.5*pos[0]*(pos[0] - 1);};
//NOTE: dim = 2
auto exactSolution = [](const auto& pos)
{
double solution = 0.0;
for (int m = 1; m <= 21; m = m + 2)
{
for (int n = 1; n <= 21; n = n + 2)
{
solution += sin(M_PI * m * pos[0]) * sin(M_PI * n * pos[1])
/ (m * n * (m*m + n*n));
}
}
solution *= -16 / pow(M_PI, 4);
return solution;
};
//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((dim + 1)*gridView_.size(0), 0.0);
Dune::DynamicVector<double>
exactPressure((dim + 1)*gridView_.size(0), 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++)
{
exactPressure[elemIdxSLE + k] = 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);
using VTKFunction = Dune::VTKFunction<GridView>;
using P1Function = Dune::NonconformingP1VTKFunction<GridView,
Dune::DynamicVector<double>>;
vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
new P1Function(gridView_, pressure, "pressure")));
vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
new P1Function(gridView_, exactPressure, "exactPressure")));
vtkWriter.pwrite("pressureData", "", "data");
}
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
private:
//assemble stiffness matrix A and load vector b
void assembleSLE (const Scalar K, const Scalar mu)
{
//NOTE:
//const int order = 3; //order of the quadrature rule
for (const auto& elem : elements(gridView_))
{
const auto& geo = elem.geometry();
std::cout << mapper_.index(elem) << "\t";
for (int k = 0; k < geo.corners(); k++)
{
std::cout << geo.corner(k) << "\t";
} std::cout << "\n\n";
}
//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();
const double elemVol = geo.volume();
const auto& center = geo.center();
std::cout << "===============================\nElement " << elemIdx;
for (int k = 0; k < geo.corners(); k++)
{
std::cout << "\t" << geo.corner(k);
}
std::cout << "\n\n";
//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(),order);
//NOTE: how are quadrature rules in Dune applied correctly?
for (const auto& ip : rule)
{
const auto qp = ip.position();
//NOTE: is the volume of the element taken into account
//automatically?
const auto weight = ip.weight();
//quadrature for int_elem q*phi_elem,0 dV
b[elemIdxSLE] += weight * q(qp);
//quadrature for int_elem q*phi_elem,i dV
for (int i = 0; i < dim; i++)
{
b[elemIdxSLE + i + 1] += weight * qp[i] * q(qp);
}
}
*/
//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++)
{
//exact evaluation of
// int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV
A[elemIdxSLE + i + 1][elemIdxSLE + i + 1] += K * elemVol;
}
//iterate over all intersection with the boundary of elem
for (const auto& intersct : intersections(gridView_, elem))
{
const auto& normal = intersct.centerUnitOuterNormal();
const auto& intersctGeo = intersct.geometry();
const double intersctVol = intersctGeo.volume(); const auto& intersctCenter = intersctGeo.center();
std::cout << "+++++++++++++++++++++++++++\nIntersection ";
for (int k = 0; k < intersctGeo.corners(); k++)
{
std::cout << "\t" << intersctGeo.corner(k);
}
std::cout << "\tnormal: " << normal << "\n\n";
//TODO: quadrature rule cannot be used for dim = 1!
// const Dune::QuadratureRule<double,dim-1>& secondOrderRule =
// Dune::QuadratureRules<double,dim-1>::rule(
// intersct.type(), 2);
//storage for multiply used integral values,
//note that only the lower diagonal and diagonal entries of
//quadraticIntregrals are used
Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);
Dune::FieldMatrix<Scalar, dim, dim> quadraticIntregrals(0.0);
//NOTE: is there a better type for symmetric matrix?
for (int i = 0; i < dim; i++)
{
//use midpoint rule for exact evaluation of
// int_intersct x_i ds
linearIntegrals[i] = intersctVol * intersctCenter[i];
for (int j = 0; j <= i; j++)
{
const auto& leftCorner = intersctGeo.corner(0);
const auto& rightCorner = intersctGeo.corner(1);
std::cout << "left corner: " << leftCorner << "\tright corner: " << rightCorner << "\n";
quadraticIntregrals[i][j] = intersctVol / 3 *
( leftCorner[i] * leftCorner[j] +
rightCorner[i] * rightCorner[j]
+ 0.5 * (leftCorner[i] * rightCorner[j] +
leftCorner[j] * rightCorner[i]) );
std::cout << i << ", " << j << "\t" << leftCorner[i] * leftCorner[j] << "\t" <<
rightCorner[i] * rightCorner[j] << "\t" <<
0.5 * (leftCorner[i] * rightCorner[j] +
leftCorner[j] * rightCorner[i]) << "\t" << intersctVol / 3 *
( leftCorner[i] * leftCorner[j] +
rightCorner[i] * rightCorner[j]
+ 0.5 * (leftCorner[i] * rightCorner[j] +
leftCorner[j] * rightCorner[i]) ) << "\t" << quadraticIntregrals[i][j] << "\n\n";
/*
//use second order quadrature rule for exact evaluation of
// int_intersct x_i*x_j ds
//NOTE: use Simpson's rule instead manually?
for (const auto& ip : secondOrderRule)
{
const auto& qp = ip.position();
quadraticIntregrals[i][j] += //NOTE: volume necessary?
ip.weight() * qp[i] * qp[j] * intersctVol;
}
*/
//NOTE: probably unnecessary
quadraticIntregrals[j][i] = quadraticIntregrals[i][j];
}
}
std::cout << "linearIntegrals:\n" << linearIntegrals << "\n\n";
std::cout << "quadraticIntregrals:\n" << quadraticIntregrals << "\n\n";
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
A[elemIdxSLE][elemIdxSLE] += mu * intersctVol;
std::cout << elemIdxSLE << ", " << elemIdxSLE <<"\t" << intersctCenter<< ":\n"
<< mu * intersctVol << "\t" << A[elemIdx][elemIdx] << "\n\n";
if (intersct.neighbor()) //intersct has neighboring element
{
std::cout << "------------------------------\nNeighbor\t";
for (int k = 0; k < intersct.outside().geometry().corners(); k++)
{
std::cout << "\t" << intersct.outside().geometry().corner(k);
}
std::cout << "\n\n";
//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 * K * normal[i] * vol(intersct)
A[elemIdxSLE + i + 1][elemIdxSLE] +=
mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE << ":\n"
<< mu * linearIntegrals[i] << "\t" <<
- 0.5 * K * normal[i] * intersctVol << "\t"
<< A[elemIdxSLE + i + 1][elemIdxSLE] << "\n\n";
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 * K * normal[i] * int_intersct x_j ds
A[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
+= mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]);
std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE + j + 1 << ":\n"
<< mu * quadraticIntregrals[i][j] << "\t" <<
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]) << "\t" <<
A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] << "\n\n";
}
}
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[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;
//stiffness matrix A is symmetric
A[neighborIdxSLE][elemIdxSLE] += A[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 * K * normal[i] * vol(intersct)
A[elemIdxSLE + i + 1][neighborIdxSLE] +=
-mu * linearIntegrals[i] + 0.5 * K * normal[i] * intersctVol;
std::cout << elemIdxSLE + i + 1 << ", " << neighborIdxSLE << ":\n"
<< -mu * linearIntegrals[i] << "\t" <<
0.5 * K * normal[i] * intersctVol << "\t" <<
A[elemIdxSLE + i + 1][neighborIdxSLE] << "\n\n";
//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 * K * normal[i] * vol(intersct)
A[elemIdxSLE][neighborIdxSLE + i + 1] +=
-mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
std::cout << elemIdxSLE << ", " << neighborIdxSLE + i + 1<< ":\n"
<< -mu * linearIntegrals[i] << "\t" <<
-0.5 * K * normal[i] * intersctVol << "\t" <<
A[elemIdxSLE][neighborIdxSLE + i + 1] << "\n\n";
//stiffness matrix A is symmetric
A[neighborIdxSLE][elemIdxSLE + i + 1] +=
A[elemIdxSLE + i + 1][neighborIdxSLE];
A[neighborIdxSLE + i + 1][elemIdxSLE] +=
A[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 * K * normal[j] * int_intersct x_i ds
//as well as
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_neighbor,j) ds
// = -0.5 * K * normal[i] * int_intersct x_j ds
A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
-mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[j] * linearIntegrals[i]
- normal[i] * linearIntegrals[j]);
std::cout << elemIdxSLE + i + 1 << ", " << neighborIdxSLE + j + 1<< ":\n"
<< -mu * quadraticIntregrals[i][j] << "\t" <<
- 0.5 * K * (normal[j] * linearIntegrals[i]
- normal[i] * linearIntegrals[j]) << "\t" <<
A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] << "\n\n";
//stiffness matrix A is symmetric
A[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1];
if (i != j)
{
// 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 * K * normal[i] * int_intersct x_j ds
//as well as
// int_intersct avg(K*grad(phi_elem,j))
// *jump(phi_neighbor,i) ds
// = -0.5 * K * normal[j] * int_intersct x_i ds A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] +=
-mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
- normal[j] * linearIntegrals[i]);
std::cout << elemIdxSLE + j + 1 << ", " << neighborIdxSLE + i + 1<< ":\n"
<< -mu * quadraticIntregrals[i][j] << "\t" <<
- 0.5 * K * (normal[i] * linearIntegrals[j]
- normal[j] * linearIntegrals[i]) << "\t" <<
A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] << "\n\n";
//stiffness matrix A is symmetric
A[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] +=
A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1];
}
}
}
}
else //boundary facet
{
std::cout << "------------------------------\nBoundary\n\n";
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 for boundary facets
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_elem,0) ds
// = K * normal[i] * vol(intersct)
A[elemIdxSLE + i + 1][elemIdxSLE] +=
mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE << "\t"
<< intersctCenter << ":\n" << mu * linearIntegrals[i] << "\t" <<
- 0.5 * K * normal[i] * intersctVol << "\t"
<< A[elemIdxSLE + i + 1][elemIdxSLE] << "\n\n";
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 for boundary facets
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_elem,j) ds
// = 0.5 * K * normal[i] * int_intersct x_j ds
A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]);
std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE + j+ 1 << "\t"
<< intersctCenter << ":\n" << mu * quadraticIntregrals[i][j] << "\t" <<
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]) << "\t"
<< A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] << "\n\n";
}
}
}
}
//stiffness matrix A is symmetric
for (int i = 0; i < dim; i++)
{
A[elemIdxSLE][elemIdxSLE + i + 1] = A[elemIdxSLE + i + 1][elemIdxSLE];
for(int j = 0; j < i; j++)
{
A[elemIdxSLE + j + 1][elemIdxSLE + i + 1] =
A[elemIdxSLE + i + 1][elemIdxSLE + j + 1];
}
} }
std::cout << A << std::endl;
}
Matrix A; //stiffness matrix
Vector b; //load vector
Vector d; //solution vector
};
#endif