]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add new tests 1635/head
authorAlexander Grayver <agrayver@erdw.ethz.ch>
Fri, 18 Sep 2015 13:46:08 +0000 (15:46 +0200)
committerAlexander Grayver <agrayver@erdw.ethz.ch>
Fri, 18 Sep 2015 14:08:49 +0000 (16:08 +0200)
tests/fe/fe_project_2d.cc [new file with mode: 0644]
tests/fe/fe_project_2d.output [new file with mode: 0644]

diff --git a/tests/fe/fe_project_2d.cc b/tests/fe/fe_project_2d.cc
new file mode 100644 (file)
index 0000000..861b613
--- /dev/null
@@ -0,0 +1,265 @@
+//\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 on a\r
+ * series of globally refined meshes and get output.\r
+ *\r
+ * Among FE spaces tested are: FE_ABF, FE_Nedelec, FE_RaviartThomas,\r
+ * FE_Q^dim (via FESystem)\r
+ *\r
+ * Alexander Grayver\r
+ */\r
+\r
+using namespace dealii;\r
+\r
+static const Point<2> vertices_nonaffine[] =\r
+{\r
+  Point<2> (-1., -1.),\r
+  Point<2> (0., -1.),\r
+  Point<2> (1., -1.),\r
+\r
+  Point<2> (-1., 0.),\r
+  Point<2> (0.3, 0.3),\r
+  Point<2> (1., 0.),\r
+\r
+  Point<2> (-1., 1.),\r
+  Point<2> (0., 1.),\r
+  Point<2> (1., 1.),\r
+};\r
+\r
+static const Point<2> vertices_affine[] =\r
+{\r
+  Point<2> (-1.4, -1.),\r
+  Point<2> (-0.4, -1.),\r
+  Point<2> (0.6, -1.),\r
+\r
+  Point<2> (-1.2, 0.),\r
+  Point<2> (-0.2, 0.),\r
+  Point<2> (0.8, 0.),\r
+\r
+  Point<2> (-1., 1.),\r
+  Point<2> (0., 1.),\r
+  Point<2> (1., 1.),\r
+};\r
+\r
+static const Point<2> vertices_rectangular[] =\r
+{\r
+  Point<2> (-1., -1.),\r
+  Point<2> (0., -1.),\r
+  Point<2> (1., -1.),\r
+\r
+  Point<2> (-1., 0.),\r
+  Point<2> (0., 0.),\r
+  Point<2> (1., 0.),\r
+\r
+  Point<2> (-1., 1.),\r
+  Point<2> (0., 1.),\r
+  Point<2> (1., 1.),\r
+};\r
+\r
+static const unsigned n_vertices = sizeof(vertices_rectangular) / sizeof(vertices_rectangular[0]);\r
+static const unsigned n_cells = 4;\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
+};\r
+\r
+template<>\r
+double VectorFunction<2>::value(const Point<2> &p, const unsigned int component) const\r
+{\r
+  Assert (component < 2,  ExcIndexRange (component, 0, 1));\r
+\r
+  const double PI = numbers::PI;\r
+  double val = 0.0;\r
+  switch (component)\r
+    {\r
+    case 0:\r
+      val = cos(PI*p(0))*sin(PI*p(1));\r
+      break;\r
+    case 1:\r
+      val = -sin(PI*p(0))*cos(PI*p(1));\r
+      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
+void create_tria(Triangulation<2> &triangulation, const Point<2> *vertices_parallelograms)\r
+{\r
+  const std::vector<Point<2> > vertices (&vertices_parallelograms[0],\r
+                                         &vertices_parallelograms[n_vertices]);\r
+\r
+  static const int cell_vertices[][GeometryInfo<2>::vertices_per_cell] =\r
+  {\r
+    {0, 1, 3, 4},\r
+    {1, 2, 4, 5},\r
+    {3, 4, 6, 7},\r
+    {4, 5, 7, 8}\r
+  };\r
+\r
+  std::vector<CellData<2> > cells (n_cells, CellData<2>());\r
+  for (unsigned i = 0; i<cells.size(); ++i)\r
+    {\r
+      for (unsigned int j=0; j<GeometryInfo<2>::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
+  triangulation.refine_global(1);\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||\tcurl(u_h)\tdiv(u_h)" << 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 unsigned int n_q_points = quadrature.size ();\r
+  MappingQ<dim> mapping(1);\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
+\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::L2_norm);\r
+\r
+      typename FEValuesViews::Vector<dim>::curl_type total_curl;\r
+      double total_div = 0;\r
+      total_curl *= 0.;\r
+\r
+      FEValues<dim> fe_values (mapping, fe, quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients);\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
+          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
+            }\r
+        }\r
+\r
+      deallog << dof_handler.n_dofs() << "\t\t"\r
+              << diff.l2_norm() << "\t"\r
+              << total_curl << "\t"\r
+              << total_div << 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 = 2;\r
+  unsigned order = 1;\r
+  unsigned n_cycles = 4;\r
+\r
+  deallog << "2d\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
+  test<dim>(FESystem<dim> (FE_Q<dim>(order), dim), n_cycles, true, vertices);\r
+  test<dim>(FE_ABF<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
+  test<dim>(FESystem<dim> (FE_Q<dim>(order), dim), n_cycles, true, vertices);\r
+  test<dim>(FE_ABF<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
+  test<dim>(FESystem<dim> (FE_Q<dim>(order), dim), n_cycles, true, vertices);\r
+  test<dim>(FE_ABF<dim> (order), n_cycles, true, vertices);\r
+\r
+  deallog << std::endl;\r
+}\r
diff --git a/tests/fe/fe_project_2d.output b/tests/fe/fe_project_2d.output
new file mode 100644 (file)
index 0000000..b30affb
--- /dev/null
@@ -0,0 +1,79 @@
+DEAL::2d
+Rectangular grid:
+dim: 2 FE_Nedelec<2>(1)
+DEAL::DoFs             ||u-u_h||       curl(u_h)       div(u_h)
+DEAL::144              0.127197        0.000000        0
+DEAL::544              0.032385        -0.000000       0
+DEAL::2112             0.008122        0.000000        0
+DEAL::8320             0.002032        -0.000000       0
+DEAL::dim: 2   FE_RaviartThomas<2>(1)
+DEAL::DoFs             ||u-u_h||       curl(u_h)       div(u_h)
+DEAL::144              0.127676        0.000000        0
+DEAL::544              0.032425        0.000000        0
+DEAL::2112             0.008124        0.000000        0
+DEAL::8320             0.002032        0.000000        0
+DEAL::dim: 2   FESystem<2>[FE_Q<2>(1)^2]
+DEAL::DoFs             ||u-u_h||       curl(u_h)       div(u_h)
+DEAL::50               0.229086        -0.000000       0
+DEAL::162              0.049159        -0.000000       0
+DEAL::578              0.011701        0.000000        0
+DEAL::2178             0.002887        0.000000        0
+DEAL::dim: 2   FE_ABF<2>(1)
+DEAL::DoFs             ||u-u_h||       curl(u_h)       div(u_h)
+DEAL::208              0.126166        0.000000        0
+DEAL::800              0.032829        -0.000000       0
+DEAL::3136             0.008676        -0.000000       0
+DEAL::12416            0.002544        -0.000000       0
+DEAL::
+Affine grid:
+dim: 2 FE_Nedelec<2>(1)
+DEAL::DoFs             ||u-u_h||       curl(u_h)       div(u_h)
+DEAL::144              0.126899        -0.010605       0
+DEAL::544              0.032013        -0.001829       0
+DEAL::2112             0.007983        -0.000255       0
+DEAL::8320             0.001993        -0.000033       0
+DEAL::dim: 2   FE_RaviartThomas<2>(1)
+DEAL::DoFs             ||u-u_h||       curl(u_h)       div(u_h)
+DEAL::144              0.142775        -0.000000       0.018305
+DEAL::544              0.036564        0.000000        0.003444
+DEAL::2112             0.009186        0.000000        0.000490
+DEAL::8320             0.002299        0.000000        0.000063
+DEAL::dim: 2   FESystem<2>[FE_Q<2>(1)^2]
+DEAL::DoFs             ||u-u_h||       curl(u_h)       div(u_h)
+DEAL::50               0.250113        -0.021467       -0.051342
+DEAL::162              0.052634        -0.002128       -0.005106
+DEAL::578              0.012425        -0.000250       -0.000600
+DEAL::2178             0.003058        -0.000031       -0.000074
+DEAL::dim: 2   FE_ABF<2>(1)
+DEAL::DoFs             ||u-u_h||       curl(u_h)       div(u_h)
+DEAL::208              0.150632        0.000000        0
+DEAL::800              0.038793        -0.000000       0
+DEAL::3136             0.010523        -0.000000       0
+DEAL::12416            0.003304        0.000000        0
+DEAL::
+Non-affine grid:
+dim: 2 FE_Nedelec<2>(1)
+DEAL::DoFs             ||u-u_h||       curl(u_h)       div(u_h)
+DEAL::144              0.153434        0.020872        0
+DEAL::544              0.037846        0.005462        0
+DEAL::2112             0.009393        0.000834        0
+DEAL::8320             0.002340        0.000110        0
+DEAL::dim: 2   FE_RaviartThomas<2>(1)
+DEAL::DoFs             ||u-u_h||       curl(u_h)       div(u_h)
+DEAL::144              0.163372        -0.033531       0
+DEAL::544              0.042792        -0.006665       0
+DEAL::2112             0.010774        -0.001430       0
+DEAL::8320             0.002695        -0.000330       0
+DEAL::dim: 2   FESystem<2>[FE_Q<2>(1)^2]
+DEAL::DoFs             ||u-u_h||       curl(u_h)       div(u_h)
+DEAL::50               0.303441        -0.149674       0
+DEAL::162              0.064869        -0.009370       0
+DEAL::578              0.014893        -0.000911       0
+DEAL::2178             0.003620        -0.000108       0
+DEAL::dim: 2   FE_ABF<2>(1)
+DEAL::DoFs             ||u-u_h||       curl(u_h)       div(u_h)
+DEAL::208              0.191291        -0.048789       0
+DEAL::800              0.047830        -0.008916       0
+DEAL::3136             0.012716        -0.002205       0
+DEAL::12416            0.003947        -0.000586       0
+DEAL::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.