From 5a72d71211fc85e9c5aa004d144eca314ddde516 Mon Sep 17 00:00:00 2001 From: Alexander Grayver Date: Thu, 24 Sep 2015 16:49:33 +0200 Subject: [PATCH] Add a test for checking correct behaviour of FE_PolyTensor elements and mapping on different types of meshes. --- tests/fe/fe_project_3d.cc | 421 ++++++++++++++++++++++++++++++++++ tests/fe/fe_project_3d.output | 37 +++ 2 files changed, 458 insertions(+) create mode 100644 tests/fe/fe_project_3d.cc create mode 100644 tests/fe/fe_project_3d.output diff --git a/tests/fe/fe_project_3d.cc b/tests/fe/fe_project_3d.cc new file mode 100644 index 0000000000..f058bfa3db --- /dev/null +++ b/tests/fe/fe_project_3d.cc @@ -0,0 +1,421 @@ +// +// 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 + +#include "../tests.h" + +/* + * This program projects a function into FE spaces defined on meshes + * consisting of: rectangular cells, affine cells, non-affine cells + * + * The error, curl and divergence are then numerically calculated + * on a series of globally refined meshes and output. + * + * Among FE spaces tested are FE_Nedelec and FE_RaviartThomas + * + * Alexander Grayver, Maien Hamed + */ + +using namespace dealii; + +static const Point<3> vertices_affine[] = { + Point<3> (-1., -1., -1.), + Point<3> (0., -1., -1.), + Point<3> (1., -1., -1.), + + Point<3> (-1.2, -1., 0.), + Point<3> (-0.2, -1., 0.), + Point<3> (0.8, -1., 0.), + + Point<3> (-1.4, -1., 1), + Point<3> (-0.4, -1., 1), + Point<3> (0.6, -1., 1), + + Point<3> (-1., 0., -1.), + Point<3> (0., 0., -1.), + Point<3> (1., 0., -1.), + + Point<3> (-1.2, 0., 0.), + Point<3> (-0.2, 0., 0.), + Point<3> (0.8, 0., 0.), + + Point<3> (-1.4, 0., 1), + Point<3> (-0.4, 0., 1), + Point<3> (0.6, 0., 1), + + Point<3> (-1., 1., -1.), + Point<3> (0., 1., -1.), + Point<3> (1., 1., -1.), + + Point<3> (-1.2, 1., 0.), + Point<3> (-0.2, 1., 0.), + Point<3> (0.8, 1., 0.), + + Point<3> (-1.4, 1., 1), + Point<3> (-0.4, 1., 1), + Point<3> (0.6, 1., 1) +}; + +static const Point<3> vertices_nonaffine[] = { + Point<3> (-1., -1., -1.), + Point<3> (0., -1., -1.), + Point<3> (1., -1., -1.), + + Point<3> (-1., -1., 0.), + Point<3> (0., -1., 0.), + Point<3> (1., -1., 0.), + + Point<3> (-1., -1., 1), + Point<3> (0., -1., 1), + Point<3> (1., -1., 1), + + Point<3> (-1., 0., -1.), + Point<3> (0., 0., -1.), + Point<3> (1., 0., -1.), + + Point<3> (-1., 0., 0.), + Point<3> (0.2, 0.3, 0.1), + Point<3> (1., 0., 0.), + + Point<3> (-1., 0., 1), + Point<3> (0., 0., 1), + Point<3> (1., 0., 1), + + Point<3> (-1., 1., -1.), + Point<3> (0., 1., -1.), + Point<3> (1., 1., -1.), + + Point<3> (-1., 1., 0.), + Point<3> (0., 1., 0.), + Point<3> (1., 1., 0.), + + Point<3> (-1., 1., 1.), + Point<3> (0., 1., 1.), + Point<3> (1., 1., 1.) +}; + +static const Point<3> vertices_rectangular[] = { + Point<3> (-1., -1., -1.), + Point<3> (0., -1., -1.), + Point<3> (1., -1., -1.), + + Point<3> (-1., -1., 0.), + Point<3> (0., -1., 0.), + Point<3> (1., -1., 0.), + + Point<3> (-1., -1., 1), + Point<3> (0., -1., 1), + Point<3> (1., -1., 1), + + Point<3> (-1., 0., -1.), + Point<3> (0., 0., -1.), + Point<3> (1., 0., -1.), + + Point<3> (-1., 0., 0.), + Point<3> (0., 0., 0.), + Point<3> (1., 0., 0.), + + Point<3> (-1., 0., 1), + Point<3> (0., 0., 1), + Point<3> (1., 0., 1), + + Point<3> (-1., 1., -1.), + Point<3> (0., 1., -1.), + Point<3> (1., 1., -1.), + + Point<3> (-1., 1., 0.), + Point<3> (0., 1., 0.), + Point<3> (1., 1., 0.), + + Point<3> (-1., 1., 1.), + Point<3> (0., 1., 1.), + Point<3> (1., 1., 1.) +}; + +static const unsigned n_vertices = sizeof(vertices_rectangular) / sizeof(vertices_rectangular[0]); +static const unsigned n_cells = 8; + +template +class VectorFunction : public Function +{ +public: + VectorFunction() : Function(dim) {} + virtual double value (const Point &p, const unsigned int component) const; + virtual void vector_value(const Point &p, Vector &values) const; + virtual Tensor< 1, dim > gradient (const Point< dim > &p, const unsigned int component=0) const; +}; + +template<> +double VectorFunction<3>::value(const Point<3> &p, const unsigned int component) const +{ + Assert (component < 3, ExcIndexRange (component, 0, 2)); + + const double PI = numbers::PI; + double val = 0.0; + switch(component) + { + case 0: val = -sin(PI*p(0))*cos(PI*p(1))*cos(PI*p(2)); break; + case 1: val = -cos(PI*p(0))*sin(PI*p(1))*cos(PI*p(2)); break; + case 2: val = 2*cos(PI*p(0))*cos(PI*p(1))*sin(PI*p(2)); break; + } + return val; +} + +template +void VectorFunction::vector_value(const Point &p, Vector &values) const +{ + for(int i = 0; i < dim; ++i) + values(i) = value(p, i); +} + +template<> +Tensor<1, 3> VectorFunction<3>::gradient(const Point<3> &p, const unsigned int component) const +{ + const double PI = numbers::PI; + Tensor<1, 3> val; + double x = p(0), y = p(1), z = p(2); + + switch(component) + { + case 0: + val[0] = -PI*cos(PI*x)*cos(PI*y)*cos(PI*z); + val[1] = PI*cos(PI*z)*sin(PI*x)*sin(PI*y); + val[2] = -2*PI*cos(PI*y)*sin(PI*x)*sin(PI*z); + break; + case 1: + val[0] = PI*cos(PI*z)*sin(PI*x)*sin(PI*y); + val[1] = -PI*cos(PI*x)*cos(PI*y)*cos(PI*z); + val[2] = -2*PI*cos(PI*x)*sin(PI*y)*sin(PI*z); + break; + case 2: + val[0] = PI*cos(PI*y)*sin(PI*x)*sin(PI*z); + val[1] = PI*cos(PI*x)*sin(PI*y)*sin(PI*z); + val[2] = 2*PI*cos(PI*x)*cos(PI*y)*cos(PI*z); + break; + } + return val; +} + +void create_tria(Triangulation<3>& triangulation, const Point<3> *vertices_parallelograms) +{ + const std::vector > vertices (&vertices_parallelograms[0], + &vertices_parallelograms[n_vertices]); + + // create grid with all possible combintations of face_flip, face_orientation and face_rotation flags + static const int cell_vertices[][GeometryInfo<3>::vertices_per_cell] = { + {0, 1, 9, 10, 3, 4, 12, 13}, // cell 1 standard + {1, 2, 10, 11, 4, 5, 13, 14}, // cell 2 standard + //{10, 11, 13, 14, 1, 2, 4, 5}, // cell 2 rotated by 270 deg + {9, 10, 18, 19, 12, 13, 21, 22}, // cell 3 standard + {10, 11, 19, 20, 13, 14, 22, 23}, // cell 4 standard + //{13, 14, 10, 11, 22, 23, 19, 20}, // cell 4 rotated by 90 deg + {3, 4, 12, 13, 6, 7, 15, 16}, // cell 5 standard + {4, 5, 13, 14, 7, 8, 16, 17}, // cell 6 standard + {12, 13, 21, 22, 15, 16, 24, 25}, // cell 7 standard + //{24, 25, 15, 16, 21, 22, 12, 13}, // cell 7 rotated by 180 deg + {13, 14, 22, 23, 16, 17, 25, 26} // cell 8 standard + }; + + std::vector > cells (n_cells, CellData<3>()); + for (unsigned i = 0; i::vertices_per_cell; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + } + + triangulation.create_triangulation (vertices, cells, SubCellData()); +} + +template +void test(const FiniteElement& fe, unsigned n_cycles, bool global, const Point *vertices_parallelograms) +{ + deallog << "dim: " << dim << "\t" << fe.get_name() << std::endl; + deallog << "DoFs\t\t||u-u_h||_1\tcurl(u_h)\ttangentials\tcurl(curl(u_h))\tcurl_curl_traces\tdiv(u_h)\tboundary_flux" << std::endl; + + Triangulation triangulation; + create_tria(triangulation, vertices_parallelograms); + + DoFHandler dof_handler(triangulation); + + VectorFunction fe_function; + const FEValuesExtractors::Vector vec (0); + const QGauss quadrature (fe.degree+1); + const QGauss face_quadrature (fe.degree+1); + const unsigned int n_q_points = quadrature.size (); + const unsigned int n_face_q_points = face_quadrature.size (); + //MappingQ mapping(2); + MappingQ1 mapping; + std::vector div_v(n_q_points); + std::vector::curl_type> curl_v(n_q_points); + std::vector > hessians(n_q_points); + + std::vector > face_values (n_face_q_points); + std::vector::curl_type> face_curls (n_face_q_points); + + for(unsigned cycle = 0; cycle < n_cycles; ++cycle) + { + dof_handler.distribute_dofs(fe); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + + Vector v(dof_handler.n_dofs()); + VectorTools::project(mapping, dof_handler, constraints, quadrature, fe_function, v); + + Vector diff(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, dof_handler, v, fe_function, diff, + QGauss(fe.degree + 2), VectorTools::L1_norm); + + typename FEValuesViews::Vector::curl_type total_curl, boundary_tangentials; + Tensor<1, dim> total_curl_curl, boundary_curl_curl_traces; + double total_div = 0; + double boundary_flux = 0; + total_curl *= 0.; + boundary_tangentials *= 0.; + + FEValues fe_values (mapping, fe, quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients | update_hessians); + FEFaceValues fe_face_values(mapping, fe, face_quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients | update_normal_vectors ); + unsigned int cell_index = 0; + + for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); + cell != dof_handler.end (); ++cell, ++cell_index) + { + fe_values.reinit (cell); + const std::vector& JxW_values = fe_values.get_JxW_values (); + fe_values[vec].get_function_divergences (v, div_v); + fe_values[vec].get_function_curls (v, curl_v); + fe_values[vec].get_function_hessians (v, hessians); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + total_div += JxW_values[q_point] * div_v[q_point]; + total_curl += JxW_values[q_point] * curl_v[q_point]; + if(dim == 3) + { + total_curl_curl[0] += JxW_values[q_point] * (hessians[q_point][1][0][1] + hessians[q_point][2][0][2] - hessians[q_point][0][1][1] - hessians[q_point][0][2][2]); + total_curl_curl[1] += JxW_values[q_point] * (hessians[q_point][2][1][2] + hessians[q_point][0][0][1] - hessians[q_point][1][2][2] - hessians[q_point][1][0][0]); + total_curl_curl[2] += JxW_values[q_point] * (hessians[q_point][0][0][2] + hessians[q_point][1][1][2] - hessians[q_point][2][0][0] - hessians[q_point][2][1][1]); + } + } + + for (unsigned int face=0; face::faces_per_cell; ++face) + { + fe_face_values.reinit(cell,face); + const std::vector& face_JxW_values = fe_face_values.get_JxW_values (); + fe_face_values[vec].get_function_values (v, face_values); + if(dim==3) fe_face_values[vec].get_function_curls (v, face_curls); + for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) + { + const Tensor<1,dim> &normal = fe_face_values.normal_vector(q_point); + + // boundary flux + if(cell->at_boundary(face)) + boundary_flux += face_JxW_values[q_point] * (face_values[q_point] * normal); + else + total_div -= face_JxW_values[q_point] * (face_values[q_point] * normal); + + // boundary tangentials (curl traces) + typename FEValuesViews::Vector::curl_type n_x_v; + if(dim==2) + n_x_v[0] = (-normal[1]*face_values[q_point][0] + normal[0]*face_values[q_point][1]); + else if(dim==3) + cross_product(*reinterpret_cast*>(&n_x_v), normal, face_values[q_point]); + + if(cell->at_boundary(face)) + boundary_tangentials += face_JxW_values[q_point] * n_x_v; + else + total_curl -= face_JxW_values[q_point] * n_x_v; + + // boundary curl curl traces + if(dim==3) + { + Tensor<1,dim> n_x_curl_u; + cross_product(n_x_curl_u, normal, *reinterpret_cast*>(&face_curls[q_point])); + if(cell->at_boundary(face)) + boundary_curl_curl_traces += face_JxW_values[q_point] * n_x_curl_u; + else + total_curl_curl -= face_JxW_values[q_point] * n_x_curl_u; + } + } + } + } + + deallog << dof_handler.n_dofs() << "\t\t" + << diff.l1_norm() << "\t" + << total_curl.norm() << "\t" + << boundary_tangentials.norm() << "\t" + << total_curl_curl.norm() << "\t" + << boundary_curl_curl_traces.norm() << "\t" + << total_div << "\t" + << boundary_flux << std::endl; + + if(global) + triangulation.refine_global(); + else + { + GridRefinement::refine_and_coarsen_fixed_number(triangulation, diff, 0.3, 0.0); + triangulation.prepare_coarsening_and_refinement(); + triangulation.execute_coarsening_and_refinement(); + } + } +} + +int main () +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(6); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double (1e-8); + + const static unsigned dim = 3; + unsigned order = 1; + unsigned n_cycles = 3; + + deallog << "3d\nRectangular grid:\n"; + + const Point *vertices = &vertices_rectangular[0]; + test(FE_Nedelec (order), n_cycles, true, vertices); + test(FE_RaviartThomas (order), n_cycles, true, vertices); + + deallog << "\nAffine grid:\n"; + + vertices = &vertices_affine[0]; + test(FE_Nedelec (order), n_cycles, true, vertices); + test(FE_RaviartThomas (order), n_cycles, true, vertices); + + deallog << "\nNon-affine grid:\n"; + + vertices = &vertices_nonaffine[0]; + test(FE_Nedelec (order), n_cycles, true, vertices); + test(FE_RaviartThomas (order), n_cycles, true, vertices); + + deallog << std::endl; +} diff --git a/tests/fe/fe_project_3d.output b/tests/fe/fe_project_3d.output new file mode 100644 index 0000000000..8b428b907c --- /dev/null +++ b/tests/fe/fe_project_3d.output @@ -0,0 +1,37 @@ +DEAL::3d +Rectangular grid: +dim: 3 FE_Nedelec<3>(1) +DEAL::DoFs ||u-u_h||_1 curl(u_h) tangentials curl(curl(u_h)) curl_curl_traces div(u_h) boundary_flux +DEAL::300 4.543286 0 0 0 0 0 0 +DEAL::1944 0.768808 0 0 0 0 0 0 +DEAL::13872 0.184452 0 0 0 0 0 0 +DEAL::dim: 3 FE_RaviartThomas<3>(1) +DEAL::DoFs ||u-u_h||_1 curl(u_h) tangentials curl(curl(u_h)) curl_curl_traces div(u_h) boundary_flux +DEAL::240 1.704443 0 0 0 0 0 0 +DEAL::1728 1.047194 0 0 0 0 0 0 +DEAL::13056 0.242682 0 0 0 0 0 0 +DEAL:: +Affine grid: +dim: 3 FE_Nedelec<3>(1) +DEAL::DoFs ||u-u_h||_1 curl(u_h) tangentials curl(curl(u_h)) curl_curl_traces div(u_h) boundary_flux +DEAL::300 4.043897 0 0 0.007018 0.007018 0 0 +DEAL::1944 0.785617 0 0 0.000527 0.000527 0 0 +DEAL::13872 0.194003 0 0 0.000014 0.000014 0 0 +DEAL::dim: 3 FE_RaviartThomas<3>(1) +DEAL::DoFs ||u-u_h||_1 curl(u_h) tangentials curl(curl(u_h)) curl_curl_traces div(u_h) boundary_flux +DEAL::240 2.853715 0 0 0 0 0 0 +DEAL::1728 1.108491 0 0 0 0 0 0 +DEAL::13056 0.270326 0 0 0 0 0 0 +DEAL:: +Non-affine grid: +dim: 3 FE_Nedelec<3>(1) +DEAL::DoFs ||u-u_h||_1 curl(u_h) tangentials curl(curl(u_h)) curl_curl_traces div(u_h) boundary_flux +DEAL::300 4.500380 0.052702 0.052702 0.893524 0.893783 -0.067389 -0.067366 +DEAL::1944 0.800197 0.006576 0.006576 0.160602 0.160603 -0.026230 -0.026230 +DEAL::13872 0.196793 0.002025 0.002025 0.065033 0.065033 -0.004489 -0.004489 +DEAL::dim: 3 FE_RaviartThomas<3>(1) +DEAL::DoFs ||u-u_h||_1 curl(u_h) tangentials curl(curl(u_h)) curl_curl_traces div(u_h) boundary_flux +DEAL::240 2.198354 0.020423 0.020435 1.920196 2.031656 0.049057 0.049057 +DEAL::1728 1.150110 0.022659 0.022658 0.638775 0.724501 0.016784 0.016784 +DEAL::13056 0.272206 0.004230 0.004230 0.112209 0.193524 0.003806 0.003806 +DEAL:: -- 2.39.5