...
 
Commits (33)
This diff is collapsed.
......@@ -75,6 +75,32 @@ struct FlatVectorBackend<typename Dune::FieldMatrix<K, n, m> >
};
template<class K, int n, int m>
struct FlatVectorBackend<typename Dune::FieldVector< typename Dune::FieldVector<K, n>, m> >
{
// Note: The ordering of first and second index is reverse compared to the FieldMatrix
// specialization. This is due to the 'wrong' definition of the JacobianRange in
// discreteglobalbasisfunction.hh. When fixing that issue, the order here has to be changed
// as well.
template<class VV, class Index>
static auto getEntry(VV&& v, const Index& i) -> decltype(v[i%m][i/m])
{
return v[i%m][i/m];
}
template<class VV>
static std::size_t size(VV&& v)
{
auto size = 0;
for (std::size_t i=0; i<v.size(); ++i)
size += v[i].size();
return size;
}
};
} // namespace Dune::Functions
} // namespace Dune
......
......@@ -110,6 +110,7 @@ public:
{
using FiniteElement = typename Node::FiniteElement;
using FiniteElementRange = typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
using FiniteElementRangeField = typename FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
using FunctionBaseClass = typename Dune::LocalFiniteElementFunctionBase<FiniteElement>::type;
// Note that we capture j by reference. Hence we can switch
......@@ -126,7 +127,7 @@ public:
using FunctionFromCallable = typename Dune::Functions::FunctionFromCallable<FiniteElementRange(LocalDomain), decltype(localFj), FunctionBaseClass>;
auto interpolationValues = std::vector<FiniteElementRange>();
auto interpolationCoefficients = std::vector<FiniteElementRangeField>();
auto&& fe = node.finiteElement();
......@@ -139,7 +140,7 @@ public:
{
// We cannot use localFj directly because interpolate requires a Dune::VirtualFunction like interface
fe.localInterpolation().interpolate(FunctionFromCallable(localFj), interpolationValues);
fe.localInterpolation().interpolate(FunctionFromCallable(localFj), interpolationCoefficients);
for (size_t i=0; i<fe.localBasis().size(); ++i)
{
auto multiIndex = localIndexSet_.index(node.localIndex(i));
......@@ -149,7 +150,7 @@ public:
if (interpolateHere)
{
auto&& vectorBlock = vector_[multiIndex];
FlatVectorBackend<CoefficientBlock>::getEntry(vectorBlock, j) = interpolationValues[i];
FlatVectorBackend<CoefficientBlock>::getEntry(vectorBlock, j) = interpolationCoefficients[i];
}
}
}
......
This diff is collapsed.
......@@ -6,3 +6,12 @@ target_link_dune_default_libraries("poisson-pq2")
add_executable("stokes-taylorhood" stokes-taylorhood.cc)
target_link_dune_default_libraries("stokes-taylorhood")
add_executable("poisson-mfem" poisson-mfem.cc)
target_link_dune_default_libraries("poisson-mfem")
add_executable("rt0-basistest" rt0-basistest.cc)
target_link_dune_default_libraries("rt0-basistest")
add_executable("derivative-test" derivative-test.cc)
target_link_dune_default_libraries("derivative-test")
#include <config.h>
#include <dune/common/parallel/mpihelper.hh>
#include<dune/geometry/type.hh>
#include<dune/geometry/quadraturerules.hh>
#include <dune/istl/bvector.hh>
#include <dune/grid/common/gridview.hh>
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/grid/yaspgrid.hh>
#include "dune/functions/functionspacebases/interpolate.hh"
#include <dune/functions/functionspacebases/pqknodalbasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <math.h>
using namespace Dune;
template <class GridView, class FunctionA, class FunctionB>
double l2difference (const GridView& gridView, FunctionA& functionA, FunctionB& functionB)
{
double difference = 0.0;
for (const auto& element: elements(gridView))
{
auto geometry = element.geometry();
functionA.bind(element);
functionB.bind(element);
const auto& quad = QuadratureRules<double, GridView::dimension>::rule(element.type(), 2);
for (const auto& quadPoint : quad)
{
auto quadPos = quadPoint.position();
auto diff = functionA(quadPos);
diff -= functionB(quadPos);
auto norm = diff.two_norm2();
auto integrationElement = geometry.integrationElement(quadPos);
difference += norm * quadPoint.weight() * integrationElement;
}
}
difference = sqrt(difference);
return difference;
}
template <class GridView, class FunctionA, class FunctionB>
double h1difference (const GridView& gridView, FunctionA& functionA, FunctionB& functionB)
{
double difference = 0.0;
for (const auto& element: elements(gridView))
{
auto geometry = element.geometry();
functionA.bind(element);
functionB.bind(element);
const auto& quad = QuadratureRules<double, GridView::dimension>::rule(element.type(), 2);
for (const auto& quadPoint : quad)
{
auto quadPos = quadPoint.position();
auto diff = functionB(quadPos);
auto grad = functionA.derivative(quadPos);
for (size_t i=0; i<diff.N(); ++i)
for (size_t j=0; j<diff.M(); ++j)
diff[i][j] -= grad[i][j];
auto norm = diff.frobenius_norm2();
auto integrationElement = geometry.integrationElement(quadPos);
difference += norm * quadPoint.weight() * integrationElement;
}
}
difference = sqrt(difference);
return difference;
}
int main (int argc, char *argv[]) try
{
// Set up MPI, if available
MPIHelper::instance(argc, argv);
//////////////////////////////////////////////////////////////////////////////////////////////
// Generate the grid
//////////////////////////////////////////////////////////////////////////////////////////////
const int dim = 3;
typedef YaspGrid<dim> GridType;
FieldVector<double,dim> ur(1.0);
std::array<int,dim> numElements = {10,10,10};
GridType grid(ur,numElements);
typedef GridType::LeafGridView GridView;
GridView gridView = grid.leafGridView();
for (size_t refineLevel = 0; refineLevel < 3; ++refineLevel)
{
//////////////////////////////////////////////////////////////////////////////////////////////
// Choose a finite element space
//////////////////////////////////////////////////////////////////////////////////////////////
using namespace Functions::BasisBuilder;
auto basis = makeBasis(gridView, power<2>(pq<1>(),flatInterleaved()));
//////////////////////////////////////////////////////////////////////////////////////////////
// Define ranges
//////////////////////////////////////////////////////////////////////////////////////////////
using FunctionRange = FieldVector<double,2>;
using DerivativeRange = FieldMatrix<double,2,dim>;
//////////////////////////////////////////////////////////////////////////////////////////////
// Array type formats
//////////////////////////////////////////////////////////////////////////////////////////////
typedef BlockVector<FieldVector<double,1> > VectorType;
typedef Dune::Functions::HierarchicVectorWrapper<VectorType, double> HierarchicVectorView;
//////////////////////////////////////////////////////////////////////////////////////////////
// Define previous time step and current iterate (pressure, flux), exact solution
//////////////////////////////////////////////////////////////////////////////////////////////
VectorType x;
HierarchicVectorView(x).resize(basis);
x = 0;
//////////////////////////////////////////////////////////////////////////////////////////////
// Make a discrete function from the FE basis and the coefficient vector
//////////////////////////////////////////////////////////////////////////////////////////////
auto discreteFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FunctionRange>(basis, TypeTree::hybridTreePath(), HierarchicVectorView(x));
//////////////////////////////////////////////////////////////////////////////////////////////
// Transform discrete functions to local functions and store in a tuple
//////////////////////////////////////////////////////////////////////////////////////////////
auto localDiscreteFunction = localFunction(discreteFunction);
//////////////////////////////////////////////////////////////////////////////////////////////
// Define continuous function and its derivative
//////////////////////////////////////////////////////////////////////////////////////////////
auto pi = std::acos(-1.0);
auto function = [&pi] (const auto& x)
{
FieldVector<double, 2> result(0.0);
result[0] = sin(2.*pi*x[0]) * sin(2.*pi*x[1]) * sin(2.*pi*x[2]);
result[1] = pow(x[0],2) * exp(x[1]) - x[2];
return result;
};
auto derivative = [&pi] (const auto& x)
{
FieldMatrix<double, 2, 3> result(0.0);
result[0][0] = 2.*pi * cos(2.*pi*x[0]) * sin(2.*pi*x[1]) * sin(2.*pi*x[2]);
result[0][1] = 2.*pi * sin(2.*pi*x[0]) * cos(2.*pi*x[1]) * sin(2.*pi*x[2]);
result[0][2] = 2.*pi * sin(2.*pi*x[0]) * sin(2.*pi*x[1]) * cos(2.*pi*x[2]);
result[1][0] = 2. * x[0] * exp(x[1]);
result[1][1] = pow(x[0],2) * exp(x[1]);
result[1][2] = - 1.;
return result;
};
auto gridFunction = localFunction(Functions::makeGridViewFunction(function, gridView));
auto gridDerivative = localFunction(Functions::makeGridViewFunction(derivative, gridView));
//////////////////////////////////////////////////////////////////////////////////////////////
// Interpolate function
//////////////////////////////////////////////////////////////////////////////////////////////
interpolate(basis, Dune::TypeTree::hybridTreePath(), HierarchicVectorView(x), function);
//////////////////////////////////////////////////////////////////////////////////////////////
// Define VTK writer and output
//////////////////////////////////////////////////////////////////////////////////////////////
// SubsamplingVTKWriter<GridView> vtkWriter(gridView,1);
// vtkWriter.addVertexData(discreteFunction, VTK::FieldInfo("function", VTK::FieldInfo::Type::vector, 2));
// vtkWriter.write("derivative-test");
//////////////////////////////////////////////////////////////////////////////////////////////
// Check l2 difference for function and its derivative
//////////////////////////////////////////////////////////////////////////////////////////////
auto l2diff = l2difference(gridView, localDiscreteFunction, gridFunction);
auto h1diff = h1difference(gridView, localDiscreteFunction, gridDerivative);
std::cout << "Refine level: " << refineLevel << std::endl;
std::cout << "L2 interpolation error: " << l2diff << std::endl;
std::cout << "H1-seminorm interpolation error: " << h1diff << std::endl << std::endl;
grid.globalRefine(1);
}
}
// Error handling
catch (Exception e) {
std::cout << e << std::endl;
}
This diff is collapsed.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include <config.h>
#include <vector>
#include <cmath>
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/raviartthomasbasis.hh>
#include <dune/functions/functionspacebases/pqknodalbasis.hh>
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
using namespace Dune;
//////////////////////////////////////////////////////////
//
// Test single RT0 basis functions by plotting.
// Use BasisBuilder.
//
//////////////////////////////////////////////////////////
int main (int argc, char *argv[]) try
{
// Set up MPI, if available
MPIHelper::instance(argc, argv);
///////////////////////////////////
// Generate the grid
///////////////////////////////////
const int dim = 2;
typedef YaspGrid<dim> GridType;
FieldVector<double,dim> l(1);
std::array<int,dim> elements = {3, 3};
GridType grid(l,elements);
typedef GridType::LeafGridView GridView;
GridView gridView = grid.leafGridView();
/////////////////////////////////////////////////////////
// Choose a finite element space
/////////////////////////////////////////////////////////
using namespace Functions::BasisBuilder;
auto basis = makeBasis(
gridView,
composite(
rt<0, GeometryType::BasicType::cube>(),
pq<0>()
));
/////////////////////////////////////////////////////////
// Stiffness matrix and right hand side vector
/////////////////////////////////////////////////////////
typedef BlockVector<BlockVector<FieldVector<double,1> > > VectorType;
typedef Dune::Functions::HierarchicVectorWrapper<VectorType, double> HierarchicVectorView;
VectorType u;
HierarchicVectorView(u).resize(basis);
// using namespace TypeTree::Indices;
for (size_t i=0; i<u[0].size(); i++) {
u = 0;
u[0][i][0] = 1.0;
////////////////////////////////////////////////////////////////////////////
// Make a discrete function from the FE basis and the coefficient vector
////////////////////////////////////////////////////////////////////////////
using FluxRange = FieldVector<double,dim>;
auto fluxFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FluxRange>(basis, TypeTree::hybridTreePath(TypeTree::Indices::_0), HierarchicVectorView(u));
//////////////////////////////////////////////////////////////////////////////////////////////
// Write result to VTK file
//////////////////////////////////////////////////////////////////////////////////////////////
SubsamplingVTKWriter<GridView> vtkWriter(gridView,2);
vtkWriter.addVertexData(fluxFunction, VTK::FieldInfo("flux", VTK::FieldInfo::Type::vector, dim));
vtkWriter.write("rt0-basistest-result_"+std::to_string(i+1));
}
}
// Error handling
catch (Exception e) {
std::cout << e << std::endl;
}