From: Luca Heltai Date: Tue, 1 Mar 2016 10:48:50 +0000 (+0100) Subject: Test for bug #1646 X-Git-Tag: v8.5.0-rc1~1131^2~35 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=08d2d30fd40c32740eab02e749a38da501a5ee45;p=dealii.git Test for bug #1646 --- diff --git a/tests/fe/mapping_q_manifold_01.cc b/tests/fe/mapping_q_manifold_01.cc new file mode 100644 index 0000000000..898de99926 --- /dev/null +++ b/tests/fe/mapping_q_manifold_01.cc @@ -0,0 +1,193 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + +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 VectorFunction : public Function +{ +public: + VectorFunction(unsigned n_components): + Function(n_components) + {} + virtual double value (const Point &p, const unsigned int component) const; + virtual void value (const Point &p, Vector &values) const; + virtual void vector_value(const Point &p, Vector &values) const; +}; + +template<> +double VectorFunction<3>::value(const Point<3> &p, const unsigned int component) const +{ + return (1 - p(0)*p(0)) * (1 - p(1)*p(1)) * (1 - p(2)*p(2)); +} + +template<> +double VectorFunction<2>::value(const Point<2> &p, const unsigned int component) const +{ + return (1 - p(0)*p(0)) * (1 - p(1)*p(1)); +} + +template +void VectorFunction::value(const Point &p, Vector &values) const +{ + for (unsigned i = 0; i < values.size(); ++i) + values[i] = value(p, i); +} + +template +void VectorFunction::vector_value(const Point &p, Vector &values) const +{ + for (int i = 0; i < dim; ++i) + values(i) = value(p, i); +} + +template +void create_tria(Triangulation &triangulation, const Geometry &geometry) +{ + GridGenerator::hyper_cube (triangulation, -1., 1.); + triangulation.refine_global(3); + + 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); +} + +template +void test(const FiniteElement &fe) +{ + deallog << "dim: " << dim << "\t" << fe.get_name() << std::endl; + deallog << "Mapping degree\t||u-u_h||" << std::endl; + + Geometry geometry; + Triangulation triangulation; + create_tria(triangulation, geometry); + + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + + VectorFunction fe_function(fe.n_components()); + + for (unsigned mapping_p = 1; mapping_p <= 5; ++mapping_p) + { + MappingQ mapping(mapping_p, true); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + + Vector v(dof_handler.n_dofs()); + VectorTools::project(mapping, dof_handler, constraints, + QGauss(fe.degree + 1), fe_function, v); + + Vector diff(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, dof_handler, v, fe_function, diff, + QGauss(fe.degree + 2), VectorTools::L2_norm); + + deallog << mapping.get_degree() << "\t" + << diff.l2_norm() << std::endl; + + DataOut data_out; + std::ofstream output("output_" + Utilities::int_to_string(mapping_p) + ".vtk"); + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector(v, "v"); + data_out.add_data_vector(diff, "e"); + data_out.build_patches(mapping, fe.degree); + data_out.write_vtk(output); + output.close(); + } +} + +int main () +{ + deallog << std::setprecision (5); + deallog.attach (std::cout); + deallog.depth_console (0); + deallog.threshold_double (1e-12); + + const static unsigned dim = 3; + + for (unsigned p = 1; p < 5; ++p) + { + test(FE_Q(QGaussLobatto<1>(p+1))); + + //test(FE_RaviartThomas (p)); + //test(FE_Nedelec (p)); + //test(FESystem (FE_Q(p), dim)); + } +} diff --git a/tests/fe/mapping_q_manifold_01.output b/tests/fe/mapping_q_manifold_01.output new file mode 100644 index 0000000000..b39f7eeef1 --- /dev/null +++ b/tests/fe/mapping_q_manifold_01.output @@ -0,0 +1,8 @@ + +DEAL::dim: 3 FE_Q<3>(QGaussLobatto(5)) +DEAL::Mapping degree ||u-u_h|| +DEAL::1 6.8494e-13 +DEAL::2 6.8495e-13 +DEAL::3 6.8497e-13 +DEAL::4 6.8491e-13 +DEAL::5 6.8492e-13