DUNE-FEM (unstable)

recipe-integration.cc

See explanation at Integrate a function on a grid

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file COPYING in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// C++ includes
#include<math.h>
#include<iostream>
// dune-common includes
// dune-geometry includes
#include<dune/geometry/referenceelements.hh>
// dune-grid includes
int main(int argc, char** argv)
{
// Maybe initialize Mpi
[[maybe_unused]] Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
// [set up grid]
const int dim = 4;
using Grid = Dune::YaspGrid<dim>;
Dune::FieldVector<double,dim> len; for (auto& l:len) l=1.0;
std::array<int,dim> cells; for (auto& c : cells) c=5;
Grid grid(len,cells);
// [small vectors and matrices]
Dune::FieldVector<double,4> x({1,2,3,4}); // make a vector
auto y(x); // copy constructor
y *= 1.0/3.0; // scaling
[[maybe_unused]] auto s = x*y; // scalar product
[[maybe_unused]] auto norm = x.two_norm(); // Euclidean norm
Dune::FieldMatrix<double,4,4> A({{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}}); // make a matrix
A.mv(x,y); // matvec: y = Ax
A.usmv(0.5,x,y); // axpy: y += 0.5*Ax
// [a function to integrate]
auto u = [](const auto& x){return std::exp(x.two_norm());};
// [integration with midpoint rule]
double integral=0.0;
auto gv = grid.leafGridView(); // extract the grid view
for (const auto& e : elements(gv))
integral += u(e.geometry().center())*e.geometry().volume();
std::cout << "integral = " << integral << std::endl;
// [integration with quadrature rule]
double integral2 = 0.0;
for (const auto& e : elements(gv))
{
auto geo = e.geometry();
auto quadrature = QR::rule(geo.type(),5);
for (const auto& qp : quadrature)
integral2 += u(geo.global(qp.position()))
*geo.integrationElement(qp.position())*qp.weight();
}
std::cout << "integral2 = " << integral2 << std::endl;
// [integrating a flux]
auto f = [](const auto& x){return x;};
double divergence=0.0;
for (const auto& i : elements(gv)) {
for (const auto& I : intersections(gv,i))
if (!I.neighbor())
{
auto geoI = I.geometry();
divergence += f(geoI.center())*I.centerUnitOuterNormal()*geoI.volume();
}
}
std::cout << "divergence = " << divergence << std::endl;
}
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:641
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
A real mpi helper.
Definition: mpihelper.hh:181
static DUNE_EXPORT MPIHelper & instance(int &argc, char **&argv)
Get the singleton instance of the helper.
Definition: mpihelper.hh:252
A container for all quadrature rules of dimension dim
Definition: quadraturerules.hh:260
[ provides Dune::Grid ]
Definition: yaspgrid.hh:166
IteratorRange<... > intersections(const GV &gv, const Entity &e)
Iterates over all Intersections of an Entity with respect to the given GridView.
IteratorRange<... > elements(const GV &gv)
Iterates over all elements / cells (entities with codimension 0) of a GridView.
Helpers for dealing with MPI.
Various parser methods to get data into a ParameterTree object.
A simple timing class.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 20, 23:30, 2024)