Commit 7742f089 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

test for curvedsurfacegrid added

parent 91bba18c
......@@ -7,5 +7,5 @@ Module: dune-curvedsurfacegrid
Version: 0.1.0
Maintainer: florian.stenger@tu-dresden.de
#depending on
Depends: dune-common dune-geometry dune-grid dune-curvilineargeometry
Suggests: dune-alugrid
Depends: dune-common dune-geometry dune-grid dune-localfunctions dune-curvilineargeometry
Suggests: dune-alugrid dune-foamgrid
add_subdirectory(surfacedistance)
add_subdirectory(test)
set(HEADERS
backuprestore.hh
......
......@@ -201,6 +201,8 @@ namespace Dune
return l;
}
GlobalCoordinate normal ( const LocalCoordinate &local ) const { return mapping_->normal( local ); }
ctype integrationElement ( const LocalCoordinate &local ) const { return mapping_->integrationElement( local ); }
ctype volume () const { return mapping_->volume(std::numeric_limits<ctype>::epsilon()); }
......
if (dune-foamgrid_FOUND)
dune_add_test(SOURCES test_curvedsurfacegrid.cc
COMPILE_DEFINITIONS DUNE_GRID_PATH=${PROJECT_SOURCE_DIR}/example/)
add_dune_alberta_flags(test_curvedsurfacegrid WORLDDIM 3)
endif ()
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#include <iostream>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/curvedsurfacegrid/curvedsurfacegrid.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/localfunctions/lagrange/pk.hh>
//#include <dune/grid/albertagrid.hh>
#include <dune/foamgrid/foamgrid.hh>
#define STR(s) STR_HELPER(s)
#define STR_HELPER(s) #s
const int order = 3;
const int quad_order = order+5;
const int num_levels = 5;
// Test Hausdorf-distance
template <class Grid, class Projection>
typename Grid::ctype inf_error(Grid const& grid, Projection const& projection)
{
using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
typename Grid::ctype dist = 0;
for (auto const& element : elements(grid.leafGridView()))
{
auto geometry = element.geometry();
auto const& quadRule = QuadProvider::rule(element.type(), quad_order);
for (auto const& qp : quadRule) {
auto x = geometry.global(qp.position());
auto y = projection(x);
dist = std::max(dist, (x - y).two_norm());
}
}
return dist;
}
// Test integrated Hausdorf-distance
template <class Grid, class Projection>
typename Grid::ctype L2_error(Grid const& grid, Projection const& projection)
{
using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
typename Grid::ctype dist = 0;
for (auto const& element : elements(grid.leafGridView()))
{
auto geometry = element.geometry();
auto const& quadRule = QuadProvider::rule(element.type(), quad_order);
for (auto const& qp : quadRule) {
auto x = geometry.global(qp.position());
auto y = projection(x);
dist += (x - y).two_norm2() * qp.weight() * geometry.integrationElement(qp.position());
}
}
return std::sqrt(dist);
}
// Test integrated error in normal vectors
template <class Grid, class Projection>
typename Grid::ctype normal_error(Grid const& grid, Projection const& projection)
{
using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
typename Grid::ctype dist = 0;
for (auto const& element : elements(grid.leafGridView()))
{
auto geometry = element.geometry();
auto const& quadRule = QuadProvider::rule(element.type(), quad_order);
for (auto const& qp : quadRule) {
auto x = geometry.impl().normal(qp.position());
auto y = projection(geometry.global(qp.position()));
if ((y + x).two_norm() < 0.9)
x *= -1;
dist += (x - y).two_norm2() * qp.weight() * geometry.integrationElement(qp.position());
}
}
return std::sqrt(dist);
}
// Test integrated error in mean curvature
template <class Grid, class Projection>
typename Grid::ctype curvature_error(Grid const& grid, Projection const& projection)
{
using LocalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::LocalCoordinate;
using GlobalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::GlobalCoordinate;
using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
using FiniteElementType = Dune::PkLocalFiniteElement<typename Grid::ctype, typename Grid::ctype, 2, order+2>;
FiniteElementType fe;
using LBTraits = typename FiniteElementType::Traits::LocalBasisType::Traits;
std::vector<typename LBTraits::JacobianType> shapeGradients;
std::vector<GlobalCoordinate> gradients;
std::vector<GlobalCoordinate> normals;
typename Grid::ctype dist = 0;
for (auto const& element : elements(grid.leafGridView()))
{
auto geometry = element.geometry();
auto const& localBasis = fe.localBasis();
auto const& localInterpolation = fe.localInterpolation();
auto f = [&](auto const& local) -> GlobalCoordinate
{
auto x = geometry.impl().normal(local);
auto y = projection(geometry.global(local));
if ((y + x).two_norm() < 0.9)
x *= -1;
return x;
};
localInterpolation.interpolate(f, normals);
shapeGradients.resize(localBasis.size());
gradients.resize(localBasis.size());
auto const& quadRule = QuadProvider::rule(element.type(), quad_order);
for (auto const& qp : quadRule) {
auto jInvT = geometry.jacobianInverseTransposed(qp.position());
localBasis.evaluateJacobian(qp.position(), shapeGradients);
for (std::size_t i = 0; i < gradients.size(); ++i)
jInvT.mv(shapeGradients[i][0], gradients[i]);
typename Grid::ctype H = 0;
for (std::size_t i = 0; i < gradients.size(); ++i)
H += normals[i].dot(gradients[i]);
dist += std::abs(H - 2) * qp.weight() * geometry.integrationElement(qp.position());
}
}
return dist;
}
// longest edge in the grid
template <class Grid>
typename Grid::ctype edge_length(Grid const& grid)
{
typename Grid::ctype h = 0;
for (auto const& e : edges(grid.leafGridView()))
h = std::max(h, e.geometry().volume());
return h;
}
int main(int argc, char** argv)
{
using namespace Dune;
MPIHelper& helper = MPIHelper::instance(argc, argv);
// using HostGrid = AlbertaGrid<2,3>;
using HostGrid = FoamGrid<2,3>;
std::unique_ptr<HostGrid> hostGrid = GmshReader<HostGrid>::read( STR(DUNE_GRID_PATH) "sphere_coarse.msh");
SphereCoordFunction sphere({0.0, 0.0, 0.0}, 1.0);
using Grid = CurvedSurfaceGrid<HostGrid, SphereCoordFunction, order>;
Grid grid(*hostGrid, sphere);
std::vector<typename Grid::ctype> inf_errors, L2_errors, normal_errors, curvature_errors, edge_lengths;
for (int i = 0; i < num_levels; ++i) {
if (i > 0)
grid.globalRefine(1);
inf_errors.push_back( inf_error(grid, sphere) );
L2_errors.push_back( L2_error(grid, sphere) );
normal_errors.push_back( normal_error(grid, sphere) );
curvature_errors.push_back( curvature_error(grid, sphere) );
edge_lengths.push_back( edge_length(grid) );
}
std::vector<typename Grid::ctype> eocInf(num_levels, 0), eocL2(num_levels, 0), eocNormal(num_levels, 0), eocCurvature(num_levels, 0);
for (int i = 1; i < num_levels; ++i) {
eocInf[i] = std::log(inf_errors[i]/inf_errors[i-1]) / std::log(edge_lengths[i]/edge_lengths[i-1]);
eocL2[i] = std::log(L2_errors[i]/L2_errors[i-1]) / std::log(edge_lengths[i]/edge_lengths[i-1]);
eocNormal[i] = std::log(normal_errors[i]/normal_errors[i-1]) / std::log(edge_lengths[i]/edge_lengths[i-1]);
eocCurvature[i] = std::log(curvature_errors[i]/curvature_errors[i-1]) / std::log(edge_lengths[i]/edge_lengths[i-1]);
}
auto print_line = [](auto i, auto const& data) {
std::cout << i;
for (auto const& d : data)
std::cout << '\t' << "| " << data;
std::cout << std::endl;
};
auto print_break = [] {
std::cout.width(8 + 8*16);
std::cout.fill('-');
std::cout << '-' << std::endl;
std::cout.fill(' ');
};
std::cout.setf(std::ios::scientific);
print_line("level", std::vector<std::string>{"err_inf", "eoc_inf", "err_L2", "eoc_L2", "err_norm", "eoc_norm", "err_curv", "eoc_curv"});
print_break();
for (int i = 0; i < num_levels; ++i)
print_line(i, std::vector<double>{inf_errors[i], eocInf[i], L2_errors[i], eocL2[i], normal_errors[i], eocNormal[i], curvature_errors[i], eocCurvature[i]});
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment