--- /dev/null
+//\r
+// This file is part of the deal.II library.\r
+//\r
+// The deal.II library is free software; you can use it, redistribute\r
+// it, and/or modify it under the terms of the GNU Lesser General\r
+// Public License as published by the Free Software Foundation; either\r
+// version 2.1 of the License, or (at your option) any later version.\r
+// The full text of the license can be found in the file LICENSE at\r
+// the top level of the deal.II distribution.\r
+//\r
+// ---------------------------------------------------------------------\r
+#include <deal.II/base/function.h>\r
+#include <deal.II/dofs/dof_handler.h>\r
+#include <deal.II/dofs/dof_accessor.h>\r
+#include <deal.II/dofs/dof_tools.h>\r
+#include <deal.II/grid/grid_generator.h>\r
+#include <deal.II/grid/tria.h>\r
+#include <deal.II/grid/grid_tools.h>\r
+#include <deal.II/grid/grid_refinement.h>\r
+#include <deal.II/fe/fe_nedelec.h>\r
+#include <deal.II/fe/fe_q.h>\r
+#include <deal.II/fe/fe_system.h>\r
+#include <deal.II/fe/fe_raviart_thomas.h>\r
+#include <deal.II/fe/fe_abf.h>\r
+#include <deal.II/fe/fe_values.h>\r
+#include <deal.II/fe/mapping_q.h>\r
+#include <deal.II/lac/constraint_matrix.h>\r
+#include <deal.II/numerics/vector_tools.h>\r
+\r
+#include "../tests.h"\r
+\r
+/*\r
+ * This program projects a function into FE spaces defined on meshes\r
+ * consisting of: rectangular cells, affine cells, non-affine cells\r
+ * \r
+ * The error, curl and divergence are then numerically calculated\r
+ * on a series of globally refined meshes and output.\r
+ * \r
+ * Among FE spaces tested are FE_Nedelec and FE_RaviartThomas\r
+ * \r
+ * Alexander Grayver, Maien Hamed\r
+ */\r
+\r
+using namespace dealii;\r
+\r
+static const Point<3> vertices_affine[] = {\r
+ Point<3> (-1., -1., -1.),\r
+ Point<3> (0., -1., -1.),\r
+ Point<3> (1., -1., -1.),\r
+\r
+ Point<3> (-1.2, -1., 0.),\r
+ Point<3> (-0.2, -1., 0.),\r
+ Point<3> (0.8, -1., 0.),\r
+\r
+ Point<3> (-1.4, -1., 1),\r
+ Point<3> (-0.4, -1., 1),\r
+ Point<3> (0.6, -1., 1),\r
+\r
+ Point<3> (-1., 0., -1.),\r
+ Point<3> (0., 0., -1.),\r
+ Point<3> (1., 0., -1.),\r
+\r
+ Point<3> (-1.2, 0., 0.),\r
+ Point<3> (-0.2, 0., 0.),\r
+ Point<3> (0.8, 0., 0.),\r
+\r
+ Point<3> (-1.4, 0., 1),\r
+ Point<3> (-0.4, 0., 1),\r
+ Point<3> (0.6, 0., 1),\r
+\r
+ Point<3> (-1., 1., -1.),\r
+ Point<3> (0., 1., -1.),\r
+ Point<3> (1., 1., -1.),\r
+\r
+ Point<3> (-1.2, 1., 0.),\r
+ Point<3> (-0.2, 1., 0.),\r
+ Point<3> (0.8, 1., 0.),\r
+\r
+ Point<3> (-1.4, 1., 1),\r
+ Point<3> (-0.4, 1., 1),\r
+ Point<3> (0.6, 1., 1)\r
+};\r
+\r
+static const Point<3> vertices_nonaffine[] = {\r
+ Point<3> (-1., -1., -1.),\r
+ Point<3> (0., -1., -1.),\r
+ Point<3> (1., -1., -1.),\r
+\r
+ Point<3> (-1., -1., 0.),\r
+ Point<3> (0., -1., 0.),\r
+ Point<3> (1., -1., 0.),\r
+\r
+ Point<3> (-1., -1., 1),\r
+ Point<3> (0., -1., 1),\r
+ Point<3> (1., -1., 1),\r
+\r
+ Point<3> (-1., 0., -1.),\r
+ Point<3> (0., 0., -1.),\r
+ Point<3> (1., 0., -1.),\r
+\r
+ Point<3> (-1., 0., 0.),\r
+ Point<3> (0.2, 0.3, 0.1),\r
+ Point<3> (1., 0., 0.),\r
+\r
+ Point<3> (-1., 0., 1),\r
+ Point<3> (0., 0., 1),\r
+ Point<3> (1., 0., 1),\r
+\r
+ Point<3> (-1., 1., -1.),\r
+ Point<3> (0., 1., -1.),\r
+ Point<3> (1., 1., -1.),\r
+\r
+ Point<3> (-1., 1., 0.),\r
+ Point<3> (0., 1., 0.),\r
+ Point<3> (1., 1., 0.),\r
+\r
+ Point<3> (-1., 1., 1.),\r
+ Point<3> (0., 1., 1.),\r
+ Point<3> (1., 1., 1.)\r
+};\r
+\r
+static const Point<3> vertices_rectangular[] = {\r
+ Point<3> (-1., -1., -1.),\r
+ Point<3> (0., -1., -1.),\r
+ Point<3> (1., -1., -1.),\r
+\r
+ Point<3> (-1., -1., 0.),\r
+ Point<3> (0., -1., 0.),\r
+ Point<3> (1., -1., 0.),\r
+\r
+ Point<3> (-1., -1., 1),\r
+ Point<3> (0., -1., 1),\r
+ Point<3> (1., -1., 1),\r
+\r
+ Point<3> (-1., 0., -1.),\r
+ Point<3> (0., 0., -1.),\r
+ Point<3> (1., 0., -1.),\r
+\r
+ Point<3> (-1., 0., 0.),\r
+ Point<3> (0., 0., 0.),\r
+ Point<3> (1., 0., 0.),\r
+\r
+ Point<3> (-1., 0., 1),\r
+ Point<3> (0., 0., 1),\r
+ Point<3> (1., 0., 1),\r
+\r
+ Point<3> (-1., 1., -1.),\r
+ Point<3> (0., 1., -1.),\r
+ Point<3> (1., 1., -1.),\r
+\r
+ Point<3> (-1., 1., 0.),\r
+ Point<3> (0., 1., 0.),\r
+ Point<3> (1., 1., 0.),\r
+\r
+ Point<3> (-1., 1., 1.),\r
+ Point<3> (0., 1., 1.),\r
+ Point<3> (1., 1., 1.)\r
+};\r
+\r
+static const unsigned n_vertices = sizeof(vertices_rectangular) / sizeof(vertices_rectangular[0]);\r
+static const unsigned n_cells = 8;\r
+\r
+template<int dim>\r
+class VectorFunction : public Function<dim>\r
+{\r
+public:\r
+ VectorFunction() : Function<dim>(dim) {}\r
+ virtual double value (const Point<dim> &p, const unsigned int component) const;\r
+ virtual void vector_value(const Point<dim> &p, Vector<double> &values) const;\r
+ virtual Tensor< 1, dim > gradient (const Point< dim > &p, const unsigned int component=0) const;\r
+};\r
+\r
+template<>\r
+double VectorFunction<3>::value(const Point<3> &p, const unsigned int component) const\r
+{\r
+ Assert (component < 3, ExcIndexRange (component, 0, 2));\r
+\r
+ const double PI = numbers::PI;\r
+ double val = 0.0;\r
+ switch(component)\r
+ {\r
+ case 0: val = -sin(PI*p(0))*cos(PI*p(1))*cos(PI*p(2)); break;\r
+ case 1: val = -cos(PI*p(0))*sin(PI*p(1))*cos(PI*p(2)); break;\r
+ case 2: val = 2*cos(PI*p(0))*cos(PI*p(1))*sin(PI*p(2)); break;\r
+ }\r
+ return val;\r
+}\r
+\r
+template<int dim>\r
+void VectorFunction<dim>::vector_value(const Point<dim> &p, Vector<double> &values) const\r
+{\r
+ for(int i = 0; i < dim; ++i)\r
+ values(i) = value(p, i);\r
+}\r
+\r
+template<>\r
+Tensor<1, 3> VectorFunction<3>::gradient(const Point<3> &p, const unsigned int component) const\r
+{\r
+ const double PI = numbers::PI;\r
+ Tensor<1, 3> val;\r
+ double x = p(0), y = p(1), z = p(2);\r
+\r
+ switch(component)\r
+ {\r
+ case 0:\r
+ val[0] = -PI*cos(PI*x)*cos(PI*y)*cos(PI*z);\r
+ val[1] = PI*cos(PI*z)*sin(PI*x)*sin(PI*y);\r
+ val[2] = -2*PI*cos(PI*y)*sin(PI*x)*sin(PI*z);\r
+ break;\r
+ case 1:\r
+ val[0] = PI*cos(PI*z)*sin(PI*x)*sin(PI*y);\r
+ val[1] = -PI*cos(PI*x)*cos(PI*y)*cos(PI*z);\r
+ val[2] = -2*PI*cos(PI*x)*sin(PI*y)*sin(PI*z);\r
+ break;\r
+ case 2:\r
+ val[0] = PI*cos(PI*y)*sin(PI*x)*sin(PI*z);\r
+ val[1] = PI*cos(PI*x)*sin(PI*y)*sin(PI*z);\r
+ val[2] = 2*PI*cos(PI*x)*cos(PI*y)*cos(PI*z);\r
+ break;\r
+ }\r
+ return val;\r
+}\r
+\r
+void create_tria(Triangulation<3>& triangulation, const Point<3> *vertices_parallelograms)\r
+{\r
+ const std::vector<Point<3> > vertices (&vertices_parallelograms[0],\r
+ &vertices_parallelograms[n_vertices]);\r
+\r
+ // create grid with all possible combintations of face_flip, face_orientation and face_rotation flags\r
+ static const int cell_vertices[][GeometryInfo<3>::vertices_per_cell] = {\r
+ {0, 1, 9, 10, 3, 4, 12, 13}, // cell 1 standard\r
+ {1, 2, 10, 11, 4, 5, 13, 14}, // cell 2 standard\r
+ //{10, 11, 13, 14, 1, 2, 4, 5}, // cell 2 rotated by 270 deg\r
+ {9, 10, 18, 19, 12, 13, 21, 22}, // cell 3 standard\r
+ {10, 11, 19, 20, 13, 14, 22, 23}, // cell 4 standard\r
+ //{13, 14, 10, 11, 22, 23, 19, 20}, // cell 4 rotated by 90 deg\r
+ {3, 4, 12, 13, 6, 7, 15, 16}, // cell 5 standard\r
+ {4, 5, 13, 14, 7, 8, 16, 17}, // cell 6 standard\r
+ {12, 13, 21, 22, 15, 16, 24, 25}, // cell 7 standard\r
+ //{24, 25, 15, 16, 21, 22, 12, 13}, // cell 7 rotated by 180 deg\r
+ {13, 14, 22, 23, 16, 17, 25, 26} // cell 8 standard\r
+ };\r
+\r
+ std::vector<CellData<3> > cells (n_cells, CellData<3>());\r
+ for (unsigned i = 0; i<cells.size(); ++i)\r
+ {\r
+ for (unsigned int j=0; j<GeometryInfo<3>::vertices_per_cell; ++j)\r
+ cells[i].vertices[j] = cell_vertices[i][j];\r
+ cells[i].material_id = 0;\r
+ }\r
+\r
+ triangulation.create_triangulation (vertices, cells, SubCellData());\r
+}\r
+\r
+template <int dim>\r
+void test(const FiniteElement<dim>& fe, unsigned n_cycles, bool global, const Point<dim> *vertices_parallelograms)\r
+{\r
+ deallog << "dim: " << dim << "\t" << fe.get_name() << std::endl;\r
+ 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;\r
+\r
+ Triangulation<dim> triangulation;\r
+ create_tria(triangulation, vertices_parallelograms);\r
+\r
+ DoFHandler<dim> dof_handler(triangulation);\r
+\r
+ VectorFunction<dim> fe_function;\r
+ const FEValuesExtractors::Vector vec (0);\r
+ const QGauss<dim> quadrature (fe.degree+1);\r
+ const QGauss<dim-1> face_quadrature (fe.degree+1);\r
+ const unsigned int n_q_points = quadrature.size ();\r
+ const unsigned int n_face_q_points = face_quadrature.size ();\r
+ //MappingQ<dim> mapping(2);\r
+ MappingQ1<dim> mapping;\r
+ std::vector<double> div_v(n_q_points);\r
+ std::vector<typename FEValuesViews::Vector<dim>::curl_type> curl_v(n_q_points);\r
+ std::vector<Tensor<3,dim> > hessians(n_q_points);\r
+\r
+ std::vector<Tensor<1,dim> > face_values (n_face_q_points);\r
+ std::vector<typename FEValuesViews::Vector<dim>::curl_type> face_curls (n_face_q_points);\r
+\r
+ for(unsigned cycle = 0; cycle < n_cycles; ++cycle)\r
+ {\r
+ dof_handler.distribute_dofs(fe);\r
+\r
+ ConstraintMatrix constraints;\r
+ DoFTools::make_hanging_node_constraints(dof_handler, constraints);\r
+ constraints.close();\r
+\r
+ Vector<double> v(dof_handler.n_dofs());\r
+ VectorTools::project(mapping, dof_handler, constraints, quadrature, fe_function, v);\r
+\r
+ Vector<float> diff(triangulation.n_active_cells());\r
+ VectorTools::integrate_difference(mapping, dof_handler, v, fe_function, diff,\r
+ QGauss<dim>(fe.degree + 2), VectorTools::L1_norm);\r
+\r
+ typename FEValuesViews::Vector<dim>::curl_type total_curl, boundary_tangentials;\r
+ Tensor<1, dim> total_curl_curl, boundary_curl_curl_traces;\r
+ double total_div = 0;\r
+ double boundary_flux = 0;\r
+ total_curl *= 0.;\r
+ boundary_tangentials *= 0.;\r
+\r
+ FEValues<dim> fe_values (mapping, fe, quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients | update_hessians);\r
+ FEFaceValues<dim> fe_face_values(mapping, fe, face_quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients | update_normal_vectors );\r
+ unsigned int cell_index = 0;\r
+\r
+ for (typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active ();\r
+ cell != dof_handler.end (); ++cell, ++cell_index)\r
+ {\r
+ fe_values.reinit (cell);\r
+ const std::vector<double>& JxW_values = fe_values.get_JxW_values ();\r
+ fe_values[vec].get_function_divergences (v, div_v);\r
+ fe_values[vec].get_function_curls (v, curl_v);\r
+ fe_values[vec].get_function_hessians (v, hessians);\r
+ for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)\r
+ {\r
+ total_div += JxW_values[q_point] * div_v[q_point];\r
+ total_curl += JxW_values[q_point] * curl_v[q_point];\r
+ if(dim == 3)\r
+ {\r
+ 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]);\r
+ 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]);\r
+ 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]);\r
+ }\r
+ }\r
+\r
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)\r
+ {\r
+ fe_face_values.reinit(cell,face);\r
+ const std::vector<double>& face_JxW_values = fe_face_values.get_JxW_values ();\r
+ fe_face_values[vec].get_function_values (v, face_values);\r
+ if(dim==3) fe_face_values[vec].get_function_curls (v, face_curls);\r
+ for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)\r
+ {\r
+ const Tensor<1,dim> &normal = fe_face_values.normal_vector(q_point);\r
+\r
+ // boundary flux\r
+ if(cell->at_boundary(face))\r
+ boundary_flux += face_JxW_values[q_point] * (face_values[q_point] * normal);\r
+ else\r
+ total_div -= face_JxW_values[q_point] * (face_values[q_point] * normal);\r
+\r
+ // boundary tangentials (curl traces)\r
+ typename FEValuesViews::Vector<dim>::curl_type n_x_v;\r
+ if(dim==2)\r
+ n_x_v[0] = (-normal[1]*face_values[q_point][0] + normal[0]*face_values[q_point][1]);\r
+ else if(dim==3)\r
+ cross_product(*reinterpret_cast<Tensor<1,dim>*>(&n_x_v), normal, face_values[q_point]);\r
+\r
+ if(cell->at_boundary(face))\r
+ boundary_tangentials += face_JxW_values[q_point] * n_x_v;\r
+ else\r
+ total_curl -= face_JxW_values[q_point] * n_x_v;\r
+\r
+ // boundary curl curl traces\r
+ if(dim==3)\r
+ {\r
+ Tensor<1,dim> n_x_curl_u;\r
+ cross_product(n_x_curl_u, normal, *reinterpret_cast<Tensor<1,dim>*>(&face_curls[q_point]));\r
+ if(cell->at_boundary(face))\r
+ boundary_curl_curl_traces += face_JxW_values[q_point] * n_x_curl_u;\r
+ else\r
+ total_curl_curl -= face_JxW_values[q_point] * n_x_curl_u;\r
+ }\r
+ }\r
+ }\r
+ }\r
+\r
+ deallog << dof_handler.n_dofs() << "\t\t"\r
+ << diff.l1_norm() << "\t"\r
+ << total_curl.norm() << "\t"\r
+ << boundary_tangentials.norm() << "\t"\r
+ << total_curl_curl.norm() << "\t"\r
+ << boundary_curl_curl_traces.norm() << "\t"\r
+ << total_div << "\t"\r
+ << boundary_flux << std::endl;\r
+\r
+ if(global)\r
+ triangulation.refine_global();\r
+ else\r
+ {\r
+ GridRefinement::refine_and_coarsen_fixed_number(triangulation, diff, 0.3, 0.0);\r
+ triangulation.prepare_coarsening_and_refinement();\r
+ triangulation.execute_coarsening_and_refinement();\r
+ }\r
+ }\r
+}\r
+\r
+int main ()\r
+{\r
+ std::ofstream logfile ("output");\r
+ deallog << std::setprecision(6);\r
+ deallog << std::fixed;\r
+ deallog.attach(logfile);\r
+ deallog.depth_console(0);\r
+ deallog.threshold_double (1e-8);\r
+\r
+ const static unsigned dim = 3;\r
+ unsigned order = 1;\r
+ unsigned n_cycles = 3;\r
+\r
+ deallog << "3d\nRectangular grid:\n";\r
+\r
+ const Point<dim> *vertices = &vertices_rectangular[0];\r
+ test<dim>(FE_Nedelec<dim> (order), n_cycles, true, vertices);\r
+ test<dim>(FE_RaviartThomas<dim> (order), n_cycles, true, vertices);\r
+\r
+ deallog << "\nAffine grid:\n";\r
+\r
+ vertices = &vertices_affine[0];\r
+ test<dim>(FE_Nedelec<dim> (order), n_cycles, true, vertices);\r
+ test<dim>(FE_RaviartThomas<dim> (order), n_cycles, true, vertices);\r
+\r
+ deallog << "\nNon-affine grid:\n";\r
+\r
+ vertices = &vertices_nonaffine[0];\r
+ test<dim>(FE_Nedelec<dim> (order), n_cycles, true, vertices);\r
+ test<dim>(FE_RaviartThomas<dim> (order), n_cycles, true, vertices);\r
+\r
+ deallog << std::endl;\r
+}\r