diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh new file mode 100644 index 0000000000000000000000000000000000000000..e5b4adbb0a8c0524605cb99e9e83ff3f4f3ea16c --- /dev/null +++ b/dune/mmdg/dg.hh @@ -0,0 +1,374 @@ +#ifndef DUNE_MMDG_DG_HH +#define DUNE_MMDG_DG_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; + + //constructor + DG (const GridView& gridView, const Mapper& mapper, + const Problem& problem) : + gridView_(gridView), mapper_(mapper), problem_(problem) {} + + const void operator() (const Scalar K, const Scalar mu) const + { + //NOTE: what is an appropriate sparse matrix type? + const int dof = (1 + dim)*gridView_.size(0); //degrees of freedom + using Matrix = Dune::DynamicMatrix<Scalar>; + using Vector = Dune::DynamicVector<Scalar>; + Matrix A(dof, dof, 0.0); //stiffness matrix + Vector b(dof, 0.0); //load vector + Vector d; + + //NOTE: +// const int order = 3; //order of the quadrature rule + + //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(); + + /* //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[(1 + dim)*elemIdx] += weight * q(qp); + + //quadrature for int_elem q*phi_elem,i dV + for (int i = 0; i < dim; i++) + { + b[(1 + dim)*elemIdx + i + 1] += weight * qp[i] * q(qp); + } + } + */ + + //NOTE: makeshift solution for source term q = -1 + b[(1 + dim)*elemIdx] = -elemVol; + for (int i = 0; i < dim; i++) + { + b[(1 + dim)*elemIdx + i + 1] = -elemVol * geo.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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + 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(); + + //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); + quadraticIntregrals[i][j] = intersctVol / 3 * + ( leftCorner[i] * leftCorner[j] + + rightCorner[i] * rightCorner[j] + + 0.5 * (leftCorner[i] * rightCorner[j] + + leftCorner[j] * rightCorner[i]) ); + /* + //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; + } +*/ + } + } + + //exact evaluation of + // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds + A[(1 + dim)*elemIdx][(1 + dim)*elemIdx] += + mu * intersctGeo.volume(); + + if (intersct.neighbor()) //intersct has neighboring element + { + //index of the neighboring element + const int neighborIdx = mapper_.index(intersct.outside()); + + 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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] += + mu * linearIntegrals[i] + - 0.5 * K * normal[i] * intersctVol; + + //stiffness matrix A is symmetric + A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] += + A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx]; + + 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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1] + += mu * quadraticIntregrals[i][j] + - 0.5 * K * (normal[i] * linearIntegrals[j] + + normal[j] * linearIntegrals[i]); + + //stiffness matrix A is symmetric + A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1] + += A[(1 + dim)*elemIdx + i + 1] + [(1 + dim)*elemIdx + j + 1]; + } + } + + 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[(1 + dim)*elemIdx][(1 + dim)*neighborIdx] += + -mu * intersctVol; + + //stiffness matrix A is symmetric + A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx] += + A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx]; + + 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[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx] += + -mu * linearIntegrals[i] + - 0.5 * K * normal[i] * intersctVol; + + //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[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + i + 1] += + -mu * linearIntegrals[i] + - 0.5 * K * normal[i] * intersctVol; + + //stiffness matrix A is symmetric + A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx + i + 1] += + A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx]; + A[(1 + dim)*neighborIdx + i + 1][(1 + dim)*elemIdx] += + A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + 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[(1 + dim)*elemIdx + i + 1] + [(1 + dim)*neighborIdx + j + 1] += + -mu * quadraticIntregrals[i][j] + - 0.5 * K * (normal[j] * linearIntegrals[i] + - normal[i] * linearIntegrals[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[(1 + dim)*elemIdx + j + 1] + [(1 + dim)*neighborIdx + i + 1] += + -mu * quadraticIntregrals[i][j] + - 0.5 * K * (normal[i] * linearIntegrals[j] + - normal[j] * linearIntegrals[i]); + + //stiffness matrix A is symmetric + A[(1 + dim)*neighborIdx + j + 1] + [(1 + dim)*elemIdx + i + 1] += + A[(1 + dim)*elemIdx + i + 1] + [(1 + dim)*neighborIdx + j + 1]; + A[(1 + dim)*neighborIdx + i + 1] + [(1 + dim)*elemIdx + j + 1] += + A[(1 + dim)*elemIdx + j + 1] + [(1 + dim)*neighborIdx + i + 1]; + } + } + } + else //boundary facet + { + 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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] += + mu * linearIntegrals[i] + - 0.5 * K * normal[i] * intersctVol; + + //stiffness matrix A is symmetric + A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] += + A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx]; + + 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[(1 + dim)*elemIdx + i + 1] + [(1 + dim)*elemIdx + j + 1] += + mu * quadraticIntregrals[i][j] + - 0.5 * K * (normal[i] * linearIntegrals[j] + + normal[j] * linearIntegrals[i]); + + //stiffness matrix A is symmetric + A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1] + += A[(1 + dim)*elemIdx + i + 1] + [(1 + dim)*elemIdx + j + 1]; + } + } + } + } + } + + std::cout << A << std::endl; + + //NOTE: what would be an appropiate solver here? + A.solve(d, b); + + //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); + + for (const auto& elem : elements(gridView_)) + { + const int elemIdx = mapper_.index(elem); + const auto& geo = elem.geometry(); + + for (int k = 0; k < geo.corners(); k++) + { + //contribution of the basis function + // phi_elem,0 (x) = indicator(elem); + //at the kth corner of elem + pressure[(dim + 1)*elemIdx + k] = d[(1 + dim)*elemIdx]; + + 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[(dim + 1)*elemIdx + k] += + d[(1 + dim)*elemIdx + i + 1] * geo.corner(k)[i]; + } + } + } + + //create output directory if necessary + mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); + + auto vtkWriter = + std::make_shared<Dune::VTKWriter<GridView>>(gridView_); + Dune::VTKSequenceWriter<GridView> + vtkSqWriter(vtkWriter, "pressure", "data", "data"); + + using VTKFunction = Dune::VTKFunction<GridView>; + using P1Function = Dune::NonconformingP1VTKFunction<GridView, + Dune::DynamicVector<double>>; + VTKFunction* vtkOperator = + new P1Function(gridView_, pressure, "pressure"); + + //NOTE: why does this only work with a sequence writer? + vtkSqWriter.addVertexData( + std::shared_ptr<const VTKFunction>(vtkOperator) ); + vtkSqWriter.write(0.0); + } + + const GridView& gridView_; + + const Mapper& mapper_; //element mapper for gridView_ + + const Problem& problem_; //the DG problem to be solved +}; + +#endif diff --git a/dune/mmdg/problems/dgproblem.hh b/dune/mmdg/problems/dgproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..8f62e78740576005225a8ad1997e44b657afb999 --- /dev/null +++ b/dune/mmdg/problems/dgproblem.hh @@ -0,0 +1,38 @@ +#ifndef __DUNE_MMDG_DGPROBLEM_HH__ +#define __DUNE_MMDG_DGPROBLEM_HH__ + +//abstract class providing an interface for a problem that is to be solved with +//a discontinuous Galerkin scheme +template<class Coordinate, class ctype> +class DGProblem +{ + public: + using Scalar = ctype; + using Vector = Coordinate; + + //the exact solution at position pos + virtual Scalar exactSolution (const Vector& pos) const + { + return Scalar(0.0); + } + + //indicates whether an exact solution is implemented for the problem + virtual bool hasExactSolution() const + { + return false; + }; + + //source term + virtual Scalar q (const Vector& pos) const //= 0; + { + return Scalar(0.0); + } + + //boundary condition at position pos + virtual Scalar boundary (const Vector& pos) const //= 0; + { + return Scalar(0.0); + } +}; + +#endif diff --git a/src/dg.cc b/src/dg.cc index b51f281fc12aa94c172438f3c9e2abc21ec208a4..58bfe30efc8fd91c216d9bd53118880f795912a9 100644 --- a/src/dg.cc +++ b/src/dg.cc @@ -4,6 +4,7 @@ #include <sys/stat.h> #include <sys/types.h> +#include <type_traits> #include <dune/common/parallel/mpihelper.hh> #include <dune/common/exceptions.hh> @@ -16,13 +17,15 @@ #include <dune/grid/io/file/dgfparser/dgfyasp.hh> #include <dune/grid/io/file/dgfparser/dgfug.hh> #include <dune/grid/io/file/vtk.hh> +#include <dune/grid/yaspgrid.hh> #include <dune/geometry/quadraturerules.hh> #include <dune/mmesh/mmesh.hh> +#include <dune/mmdg/dg.hh> #include <dune/mmdg/clockhelper.hh> -#include <dune/mmdg/nonconformingp1vtkfunction.hh> +#include <dune/mmdg/problems/dgproblem.hh> int main(int argc, char** argv) { @@ -35,7 +38,8 @@ int main(int argc, char** argv) static constexpr int dim = GRIDDIM; - using Grid = Dune::MovingMesh<dim>; + using Grid = std::conditional<dim == 1, Dune::YaspGrid<dim>, + Dune::MovingMesh<dim>>::type; using GridFactory = Dune::DGFGridFactory<Grid>; using GridView = typename Grid::LeafGridView; using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>; @@ -63,303 +67,16 @@ int main(int argc, char** argv) //element mapper for the leaf grid view const Mapper mapper(gridView, Dune::mcmgElementLayout()); - const int dof = (1 + dim)*gridView.size(0); //degrees of freedom - using Matrix = Dune::DynamicMatrix<double>; //Dune::FieldMatrix<double, dof, dof>; - using Vector = Dune::DynamicVector<double>; //Dune::FieldVector<double, dof>; - Matrix A(dof, dof, 0.0); //stiffness matrix - Vector b(dof, 0.0); - Vector d; + const DGProblem<Dune::FieldVector<double, dim>, double> dummyProblem; - const auto q = [](const Coordinate& x) {return x.two_norm();}; //source term + const DG dg(gridView, mapper, dummyProblem); - const int order = 3; //order of the quadrature rule - - //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(); - - //TODO: can be done outside of the loop? - const Dune::QuadratureRule<double,dim>& rule = - Dune::QuadratureRules<double,dim>::rule(geo.type(),order); - - 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[(1 + dim)*elemIdx] += weight * q(qp); - - //quadrature for int_elem q*phi_elem,i dV - for (int i = 0; i < dim; i++) - { - b[(1 + dim)*elemIdx + i + 1] += weight * qp[i] * q(qp); - } - } - - for (int i = 0; i < dim; i++) - { - //exact evaluation of - // int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV - A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + 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(); - - //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<double, dim> linearIntegrals(0.0); - Dune::FieldMatrix<double, dim, dim> quadraticIntregrals(0.0); - - for (int i = 0; i < dim; i++) - { - //use midpoint rule for exact evaluation of - // int_intersct x_i ds - linearIntegrals[i] = mu * intersctVol * 0.5 * intersctCenter[i]; - - for (int j = 0; j <= i; j++) - { - //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; - } - } - } - - //exact evaluation of - // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds - A[(1 + dim)*elemIdx][(1 + dim)*elemIdx] += mu * intersctGeo.volume(); - - if (intersct.neighbor()) //intersct has neighboring element - { - //index of the neighboring element - const int neighborIdx = mapper.index(intersct.outside()); - - 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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] = - mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol; - - //stiffness matrix A is symmetric - A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] = - A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx]; - - 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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1] = - mu * quadraticIntregrals[i][j] - 0.5 * K * (normal[i] * - linearIntegrals[j] + normal[j] * linearIntegrals[i]); - - //stiffness matrix A is symmetric - A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1] = - A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1]; - } - } - - 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[(1 + dim)*elemIdx][(1 + dim)*neighborIdx] = -mu * intersctVol; - - //stiffness matrix A is symmetric - A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx] = - A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx]; - - 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[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx] = - -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol; - - //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[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + i + 1] = - -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol; - - //stiffness matrix A is symmetric - A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx + i + 1] = - A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx]; - A[(1 + dim)*neighborIdx + i + 1][(1 + dim)*elemIdx] = - A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + 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[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx + j + 1] = - -mu * quadraticIntregrals[i][j] - - 0.5 * K * (normal[j] * linearIntegrals[i] - - normal[i] * linearIntegrals[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[(1 + dim)*elemIdx + j + 1][(1 + dim)*neighborIdx + i + 1] = - -mu * quadraticIntregrals[i][j] - - 0.5 * K * (normal[i] * linearIntegrals[j] - - normal[j] * linearIntegrals[i]); - - //stiffness matrix A is symmetric - A[(1 + dim)*neighborIdx + j + 1][(1 + dim)*elemIdx + i + 1] = - A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx + j + 1]; - A[(1 + dim)*neighborIdx + i + 1][(1 + dim)*elemIdx + j + 1] = - A[(1 + dim)*elemIdx + j + 1][(1 + dim)*neighborIdx + i + 1]; - } - } - } - else //boundary facet - { - 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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] = - mu * linearIntegrals[i] - K * normal[i] * intersctVol; - - //stiffness matrix A is symmetric - A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] = - A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx]; - - 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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1] = - mu * quadraticIntregrals[i][j] - K * (normal[i] * - linearIntegrals[j] + normal[j] * linearIntegrals[i]); - - //stiffness matrix A is symmetric - A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1] = - A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1]; - } - } - } - } - } - - A.solve(d, b); - - //storage for pressure data, for each element we store the pressure at the - //corners of the element, the VTKOperatorP1 will create the output using a - //linear interpolation for each element - Dune::DynamicVector<double> pressure(3*gridView.size(0), 0.0); - - for (const auto& elem : elements(gridView)) - { - const int elemIdx = mapper.index(elem); - const auto& geo = elem.geometry(); - - for (int k = 0; k < geo.corners(); k++) - { - //contribution of the basis function - // phi_elem,0 (x) = indicator(elem); - //at the kth corner of elem - pressure[3*elemIdx + k] = d[(1 + dim)*elemIdx]; - - 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[3*elemIdx + k] += - d[(1 + dim)*elemIdx + i + 1] * geo.corner(k)[i]; - } - } - } - -/* - typedef Dune::VTKFunction< GridView > VTKFunction; - typedef Dune::NonconformingP1VTKFunction<GridView, - Dune::DynamicVector<double>> P1Function; - VTKFunction* p1 = new P1Function(gridView, lp, "linearpressure"); - vtkwriter.addVertexData( std::shared_ptr< const VTKFunction >(p1) ); - */ - - //create output directory if necessary - mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); - - auto vtkWriter = std::make_shared<Dune::VTKWriter<GridView>>(gridView); - Dune::VTKSequenceWriter<GridView> - vtkSqWriter(vtkWriter, "pressure", "data", "data"); - - using VTKOperatorP1 = Dune::VTKFunction<GridView>; - using P1Function = Dune::NonconformingP1VTKFunction<GridView, - Dune::DynamicVector<double>>; - VTKOperatorP1* vtkOperator = new P1Function(gridView, pressure, "pressure"); - - - vtkSqWriter.addVertexData( std::shared_ptr<const VTKOperatorP1>(vtkOperator) ); - vtkSqWriter.write(0.0); + dg(K, mu); //print total run time const double timeTotal = ClockHelper::elapsedRuntime(timeInitial); - std::cout << "Finished with a total run time of " << timeTotal << ".\n"; + std::cout << "Finished with a total run time of " << timeTotal << + " Seconds.\n"; return 0; } diff --git a/src/grids/cube1d.dgf b/src/grids/cube1d.dgf index 10e02bb3085de7d407728863beadfb65b9f38200..61ec2a60ca44174b2e00266cd1840897b617d454 100644 --- a/src/grids/cube1d.dgf +++ b/src/grids/cube1d.dgf @@ -3,7 +3,7 @@ DGF Interval 0 1 -5 +20 # Simplex diff --git a/src/grids/cube2d.dgf b/src/grids/cube2d.dgf index f1c490a1ab48f79e39c668b00f8f8845fde53420..4ef9b6a3b250bd9dc9feb9ca1a4f0bce49e533a7 100644 --- a/src/grids/cube2d.dgf +++ b/src/grids/cube2d.dgf @@ -3,7 +3,7 @@ DGF Interval 0 0 1 1 -2 2 +20 20 # Simplex