From c110fb104f33888120c742f95b156379749d399f Mon Sep 17 00:00:00 2001 From: David Wells Date: Sun, 8 Jan 2017 12:49:21 -0500 Subject: [PATCH] Add a test checking ChartManifold/MappingQ convergence rates. --- tests/mappings/mapping_q_convergence.cc | 257 ++++++++++++++++++++ tests/mappings/mapping_q_convergence.output | 76 ++++++ 2 files changed, 333 insertions(+) create mode 100644 tests/mappings/mapping_q_convergence.cc create mode 100644 tests/mappings/mapping_q_convergence.output diff --git a/tests/mappings/mapping_q_convergence.cc b/tests/mappings/mapping_q_convergence.cc new file mode 100644 index 0000000000..27230c557e --- /dev/null +++ b/tests/mappings/mapping_q_convergence.cc @@ -0,0 +1,257 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include "../tests.h" + +// Test that we achieve the correct rates of convergence with MappingQ when +// all cells are curved. Credit for this test case goes to Alexander +// Grayver. This test verifies that several issues with the combination of +// MappingQ and a ChartManifold have been fixed: more precisely, the computed +// errors for this program were, for a time, much lower for the 8.2 release +// than for the development branch. At the time of writing this (on 8.5pre) +// the results of this test are identical when compiled against either 8.2 or +// 8.5pre. +// +// In addition to checking the L2 errors of a projection, this test also +// verifies that there is no loss in convergence order. + +using namespace dealii; +// Like FESeries::linear_regression. Included here so that this test can also +// be run with older versions of deal.II. +double regression_slope(const std::vector &x, + const std::vector &y) +{ + FullMatrix K(2,2), invK(2,2); + Vector X(2), B(2); + + Assert (x.size() == y.size(), + ExcMessage("x and y are expected to have the same size")); + + Assert (x.size() >= 2, + dealii::ExcMessage("at least two points are required for linear regression fit")); + + double sum_1 = 0.0, + sum_x = 0.0, + sum_x2= 0.0, + sum_y = 0.0, + sum_xy= 0.0; + + for (unsigned int i = 0; i < x.size(); i++) + { + sum_1 += 1.0; + sum_x += x[i]; + sum_x2 += x[i]*x[i]; + sum_y += y[i]; + sum_xy += x[i]*y[i]; + } + + K(0,0) = sum_1; + K(0,1) = sum_x; + K(1,0) = sum_x; + K(1,1) = sum_x2; + + B(0) = sum_y; + B(1) = sum_xy; + + invK.invert(K); + invK.vmult(X,B,false); + + return X(1); +} + + +double zvalue (const double x, const double y) +{ + double xh = x * 5., yh = y * 5.; + return (xh * exp(-xh*xh - yh*yh)) / 10.; +} + +template +class Geometry: public ChartManifold +{ +public: + virtual Point pull_back(const Point &space_point) const; + virtual Point push_forward(const Point &chart_point) const; +}; + +template +Point Geometry::pull_back(const Point &space_point) const +{ + const double d = space_point[dim - 1]; + const double z = zvalue(space_point[0], dim == 3 ? space_point[1] : 0); + + double d_hat = 0.; + if ((d - z) <= 0) + d_hat = (d - z) / (1. + z); + else + d_hat = (d - z) / (1. - z); + + Point p; + for (unsigned i = 0; i < dim - 1; ++i) + p[i] = space_point[i]; + p[dim - 1] = d_hat; + + return p; +} + +template +Point Geometry::push_forward(const Point &chart_point) const +{ + const double d_hat = chart_point[dim - 1]; + const double z = zvalue(chart_point[0], dim == 3 ? chart_point[1] : 0); + + double d = 0.; + if (d_hat <= 0) + d = d_hat + (d_hat + 1.) * z; + else + d = d_hat - (d_hat - 1.) * z; + + Point p; + for (unsigned i = 0; i < dim - 1; ++i) + p[i] = chart_point[i]; + p[dim - 1] = d; + + return p; +} + + + +template +class TranscendentalManufacturedSolution : public Function +{ +public: + virtual double value (const Point &p, const unsigned int /*component*/) const + { + return std::cos(p[0]) + 2.0*std::sin(2*p[1]); + } +}; + + + + +template +void create_tria(Triangulation &triangulation, const Geometry &geometry) +{ + GridGenerator::hyper_cube (triangulation, -1.0, 1.0); + GridTools::transform (std_cxx11::bind(&Geometry::push_forward, + std_cxx11::cref(geometry), + std_cxx11::_1), + triangulation); + + triangulation.set_manifold(0, geometry); + for (Triangulation<3>::active_cell_iterator cell=triangulation.begin_active(); + cell!=triangulation.end(); ++cell) + cell->set_all_manifold_ids(0); + + triangulation.refine_global(1); +} + + + +template +void test(const FiniteElement &fe) +{ + Geometry geometry; + TranscendentalManufacturedSolution fe_function; + ConstraintMatrix constraints; + constraints.close(); + + deallog << "FE degree: " << fe.degree << std::endl; + for (unsigned mapping_p = 1; mapping_p < fe.degree + 3; ++mapping_p) + { + deallog << "mapping order: " << mapping_p << std::endl; + Triangulation triangulation; + create_tria(triangulation, geometry); + DoFHandler dof_handler; + + std::vector log_refinements; + std::vector log_l2_errors; + + MappingQ mapping(mapping_p, true); + for (unsigned int refinement_n = 1; refinement_n < 4; ++refinement_n) + { + triangulation.refine_global(1); + dof_handler.clear(); + dof_handler.initialize(triangulation, fe); + + Vector v(dof_handler.n_dofs()); + VectorTools::project(mapping, dof_handler, constraints, + QGauss(fe.degree + (mapping_p + 2)/2), fe_function, v); + + Vector diff(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, dof_handler, v, fe_function, diff, + // superconvergence with QGauss(k + 1) + // QGauss(fe.degree + 2), + // normal convergence with QIterated + QIterated(QGauss<1>(fe.degree + 1), 2), + VectorTools::L2_norm); + log_refinements.push_back(std::log10(std::pow(2.0, -double(refinement_n)))); + log_l2_errors.push_back(std::log10(diff.l2_norm())); + if (log_refinements.size() > 1) + { + deallog << "current slope: " + << (*(log_l2_errors.end() - 1) - *(log_l2_errors.end() - 2)) + /(*(log_refinements.end() - 1) - *(log_refinements.end() - 2)) + << std::endl; + } + } + deallog << "Last number of DoFs: " + << dof_handler.n_dofs() + << std::endl; + deallog << "Last L2 error: " + << std::pow(10.0, log_l2_errors.back()) + << std::endl; + deallog << "regression slope: " + << regression_slope(log_refinements, log_l2_errors) + << std::endl; + } +} + +int main () +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog << std::setprecision (5); + deallog.depth_console (0); + + const static unsigned dim = 3; + + for (unsigned p = 1; p < 4; ++p) + { + test(FE_Q(QGaussLobatto<1>(p+1))); + } +} diff --git a/tests/mappings/mapping_q_convergence.output b/tests/mappings/mapping_q_convergence.output new file mode 100644 index 0000000000..c37bdbc4a8 --- /dev/null +++ b/tests/mappings/mapping_q_convergence.output @@ -0,0 +1,76 @@ + +DEAL::FE degree: 1 +DEAL::mapping order: 1 +DEAL::current slope: 2.0935 +DEAL::current slope: 2.0244 +DEAL::Last number of DoFs: 4913 +DEAL::Last L2 error: 0.0099885 +DEAL::regression slope: 2.0589 +DEAL::mapping order: 2 +DEAL::current slope: 2.0933 +DEAL::current slope: 2.0244 +DEAL::Last number of DoFs: 4913 +DEAL::Last L2 error: 0.0099885 +DEAL::regression slope: 2.0588 +DEAL::mapping order: 3 +DEAL::current slope: 2.0933 +DEAL::current slope: 2.0244 +DEAL::Last number of DoFs: 4913 +DEAL::Last L2 error: 0.0099885 +DEAL::regression slope: 2.0588 +DEAL::FE degree: 2 +DEAL::mapping order: 1 +DEAL::current slope: 2.8710 +DEAL::current slope: 2.9558 +DEAL::Last number of DoFs: 35937 +DEAL::Last L2 error: 0.00031701 +DEAL::regression slope: 2.9134 +DEAL::mapping order: 2 +DEAL::current slope: 2.8700 +DEAL::current slope: 2.9556 +DEAL::Last number of DoFs: 35937 +DEAL::Last L2 error: 0.00031699 +DEAL::regression slope: 2.9128 +DEAL::mapping order: 3 +DEAL::current slope: 2.8700 +DEAL::current slope: 2.9556 +DEAL::Last number of DoFs: 35937 +DEAL::Last L2 error: 0.00031699 +DEAL::regression slope: 2.9128 +DEAL::mapping order: 4 +DEAL::current slope: 2.8700 +DEAL::current slope: 2.9556 +DEAL::Last number of DoFs: 35937 +DEAL::Last L2 error: 0.00031699 +DEAL::regression slope: 2.9128 +DEAL::FE degree: 3 +DEAL::mapping order: 1 +DEAL::current slope: 4.0630 +DEAL::current slope: 4.0165 +DEAL::Last number of DoFs: 117649 +DEAL::Last L2 error: 3.3894e-06 +DEAL::regression slope: 4.0398 +DEAL::mapping order: 2 +DEAL::current slope: 4.0627 +DEAL::current slope: 4.0164 +DEAL::Last number of DoFs: 117649 +DEAL::Last L2 error: 3.3894e-06 +DEAL::regression slope: 4.0396 +DEAL::mapping order: 3 +DEAL::current slope: 4.0627 +DEAL::current slope: 4.0164 +DEAL::Last number of DoFs: 117649 +DEAL::Last L2 error: 3.3894e-06 +DEAL::regression slope: 4.0396 +DEAL::mapping order: 4 +DEAL::current slope: 4.0627 +DEAL::current slope: 4.0164 +DEAL::Last number of DoFs: 117649 +DEAL::Last L2 error: 3.3894e-06 +DEAL::regression slope: 4.0396 +DEAL::mapping order: 5 +DEAL::current slope: 4.0627 +DEAL::current slope: 4.0164 +DEAL::Last number of DoFs: 117649 +DEAL::Last L2 error: 3.3894e-06 +DEAL::regression slope: 4.0396 -- 2.39.5