]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add tests.
authorWolfgang Bangerth <bangerth@colostate.edu>
Fri, 28 Feb 2020 18:15:02 +0000 (11:15 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Fri, 28 Feb 2020 18:15:02 +0000 (11:15 -0700)
tests/numerics/nedelec_projection_mixed_01.cc [new file with mode: 0644]
tests/numerics/nedelec_projection_mixed_01.output [new file with mode: 0644]
tests/numerics/project_bv_curl_conf_04.cc [new file with mode: 0644]
tests/numerics/project_bv_curl_conf_04.output [new file with mode: 0644]

diff --git a/tests/numerics/nedelec_projection_mixed_01.cc b/tests/numerics/nedelec_projection_mixed_01.cc
new file mode 100644 (file)
index 0000000..cf76560
--- /dev/null
@@ -0,0 +1,67 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check VectorTools::project_boundary_values_curl_conforming_l2 for a
+// combination element that has edge DoFs other than those just from
+// the Nedelec base element
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+int
+main()
+{
+  initlog();
+
+  const unsigned int dim = 3;
+
+
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation, 0, 1, true);
+  MappingQ<dim>             mapping(1);
+  AffineConstraints<double> constraints;
+
+  FESystem<dim>   fe_system(FE_Q<dim>(2), 1, FE_Nedelec<dim>(0), 1);
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(fe_system);
+
+  constraints.clear();
+  VectorTools::project_boundary_values_curl_conforming_l2(
+    dof_handler,
+    1,
+    Functions::ZeroFunction<dim>(dim + 1),
+    0,
+    constraints,
+    mapping);
+  constraints.close();
+
+  constraints.print(deallog.get_file_stream());
+}
diff --git a/tests/numerics/nedelec_projection_mixed_01.output b/tests/numerics/nedelec_projection_mixed_01.output
new file mode 100644 (file)
index 0000000..8d00f4c
--- /dev/null
@@ -0,0 +1,5 @@
+
+    9 = 0
+    17 = 0
+    25 = 0
+    29 = 0
diff --git a/tests/numerics/project_bv_curl_conf_04.cc b/tests/numerics/project_bv_curl_conf_04.cc
new file mode 100644 (file)
index 0000000..8f93314
--- /dev/null
@@ -0,0 +1,419 @@
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/numerics/vector_tools.h>
+using namespace dealii;
+
+const unsigned int dim = 3;
+
+void
+apply_boundary_values(DoFHandler<dim> &          dof_handler,
+                      AffineConstraints<double> &constraints,
+                      unsigned int               n_comps,
+                      unsigned int               start_comp,
+                      Mapping<dim> &             mapping,
+                      Vector<double> &           dst)
+{
+  constraints.clear();
+  for (unsigned int i = 0; i < 6; i++)
+    VectorTools::project_boundary_values_curl_conforming_l2(
+      dof_handler,
+      start_comp,
+      Functions::ConstantFunction<dim>(1, n_comps),
+      i,
+      constraints,
+      mapping);
+  constraints.close();
+
+  constraints.distribute(dst);
+}
+
+bool
+test_boundary_values(DoFHandler<dim> &   dof_handler,
+                     Mapping<dim> &      mapping,
+                     FiniteElement<dim> &fe,
+                     unsigned int        n_comps,
+                     unsigned int        start_comp,
+                     Vector<double> &    vec)
+{
+  double ret = true;
+  // Initialize
+  QGaussLobatto<dim - 1>                         quadrature(3);
+  FEFaceValues<dim>                              fe_values(mapping,
+                              fe,
+                              quadrature,
+                              update_values | update_normal_vectors |
+                                update_quadrature_points | update_JxW_values);
+  typename DoFHandler<dim>::active_cell_iterator cell =
+                                                   dof_handler.begin_active(),
+                                                 endc = dof_handler.end();
+  unsigned                    num_cells               = 0;
+  std::vector<Vector<double>> local_averages(
+    6, Vector<double>(dof_handler.n_dofs()));
+  const FEValuesExtractors::Vector nedelec(start_comp);
+
+  // For testing we take only one cell!
+  assert(dof_handler.n_dofs() == cell->get_fe().dofs_per_cell);
+
+  for (; cell != endc; ++cell)
+    {
+      // For testing we take only one cell! Otherwise, local to global is
+      // missing
+      assert(num_cells == 0);
+      for (unsigned int face = 0;
+           face < dealii::GeometryInfo<dim>::faces_per_cell;
+           ++face)
+        {
+          if (cell->face(face)->at_boundary())
+            {
+              unsigned int   boundary_number = cell->face(face)->boundary_id();
+              Tensor<1, dim> expected_value;
+              expected_value[0] = 1;
+              expected_value[1] = 1;
+              expected_value[2] = 1;
+              assert(boundary_number < 6);
+              fe_values.reinit(cell, face);
+
+              std::vector<Vector<double>> local_values(
+                fe_values.n_quadrature_points, Vector<double>(n_comps));
+              fe_values.get_function_values(vec, local_values);
+
+              for (unsigned int q_point = 0;
+                   q_point < fe_values.n_quadrature_points;
+                   q_point++)
+                {
+                  // Compare n x v values
+                  Tensor<1, dim> nedelec_values;
+                  nedelec_values[0]     = local_values[q_point](start_comp);
+                  nedelec_values[1]     = local_values[q_point](start_comp + 1);
+                  nedelec_values[2]     = local_values[q_point](start_comp + 2);
+                  Tensor<1, dim> normal = fe_values.normal_vector(q_point);
+
+                  for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell;
+                       i++)
+                    {
+                      local_averages[boundary_number](i) +=
+                        ((cross_product_3d(normal, nedelec_values) -
+                          cross_product_3d(normal, expected_value)) *
+                         cross_product_3d(
+                           normal, fe_values[nedelec].value(i, q_point))) *
+                        fe_values.JxW(q_point);
+                    }
+                }
+            }
+        }
+      num_cells++;
+    }
+  // Test if ok
+  for (unsigned int bc = 0; bc < 6; bc++)
+    {
+      for (unsigned int basis = 0; basis < dof_handler.n_dofs(); basis++)
+        {
+          if (fabs(local_averages[bc](basis)) > 1.e-10)
+            {
+              std::cout << "Wrong values on boundary " << bc
+                        << " in basis function " << basis << " of size "
+                        << fabs(local_averages[bc](basis)) << std::endl;
+              ret = false;
+            }
+        }
+    }
+  return ret;
+}
+
+int
+main(int /*argc*/, char ** /*argv*/)
+{
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation, 0, 1, true);
+  MappingQ<dim>             mapping(1);
+  AffineConstraints<double> constraints;
+  Vector<double>            test_vec;
+  Vector<double>            test_vec2;
+
+  {
+    // Test 1: Only FE_Nedelec
+    std::cout << "Testing FE_Nedelec(0): ";
+    FE_Nedelec<dim> fe(0);
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe);
+    test_vec2.reinit(dof_handler.n_dofs());
+    apply_boundary_values(dof_handler, constraints, 3, 0, mapping, test_vec2);
+    if (test_boundary_values(dof_handler, mapping, fe, 3, 0, test_vec2))
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+  }
+  {
+    // Test 2: FESystem(FE_Nedelec(0))
+    std::cout << "Testing FESystem(FE_Nedelec(0)): ";
+    FESystem<dim>   fe_system(FE_Nedelec<dim>(0), 1);
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe_system);
+    test_vec.reinit(dof_handler.n_dofs());
+    apply_boundary_values(dof_handler, constraints, 3, 0, mapping, test_vec);
+    if (test_boundary_values(dof_handler, mapping, fe_system, 3, 0, test_vec))
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+    // Now, in the case of only one element both initializations should give
+    // identical vectors
+    std::cout << "Checking Consistency: ";
+    bool equal = (test_vec.size() == test_vec2.size());
+    if (equal)
+      {
+        for (unsigned int i = 0; i < test_vec.size(); i++)
+          {
+            equal = equal && (test_vec(i) == test_vec2(i));
+          }
+      }
+    if (equal)
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+  }
+  {
+    // Test 1b: Only FE_Nedelec
+    std::cout << "Testing FE_Nedelec(1): ";
+    FE_Nedelec<dim> fe(1);
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe);
+    test_vec2.reinit(dof_handler.n_dofs());
+    apply_boundary_values(dof_handler, constraints, 3, 0, mapping, test_vec2);
+    if (test_boundary_values(dof_handler, mapping, fe, 3, 0, test_vec2))
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+  }
+  {
+    // Test 2b: FESystem(FE_Nedelec(1))
+    std::cout << "Testing FESystem(FE_Nedelec(1)): ";
+    FESystem<dim>   fe_system(FE_Nedelec<dim>(1), 1);
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe_system);
+    test_vec.reinit(dof_handler.n_dofs());
+    apply_boundary_values(dof_handler, constraints, 3, 0, mapping, test_vec);
+    if (test_boundary_values(dof_handler, mapping, fe_system, 3, 0, test_vec))
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+    // Now, in the case of only one element both initializations should give
+    // identical vectors
+    std::cout << "Checking Consistency: ";
+    bool equal = (test_vec.size() == test_vec2.size());
+    if (equal)
+      {
+        for (unsigned int i = 0; i < test_vec.size(); i++)
+          {
+            equal = equal && (test_vec(i) == test_vec2(i));
+          }
+      }
+    if (equal)
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+  }
+  {
+    // Test 3: FESystem(FE_Nedelec(0), FE_Nedelec(0))
+    std::cout << "Testing FESystem(FE_Nedelec(0), FE_Nedelec(0)): ";
+    FESystem<dim>   fe_system(FE_Nedelec<dim>(0), 1, FE_Nedelec<dim>(0), 1);
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe_system);
+    test_vec.reinit(dof_handler.n_dofs());
+    apply_boundary_values(dof_handler, constraints, 6, 3, mapping, test_vec);
+    if (test_boundary_values(dof_handler, mapping, fe_system, 6, 3, test_vec))
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+  }
+  {
+    // Test 4: FESystem(FE_Nedelec(1), FE_Nedelec(0))
+    std::cout << "Testing FESystem(FE_Nedelec(1), FE_Nedelec(0)): ";
+    FESystem<dim>   fe_system(FE_Nedelec<dim>(1), 1, FE_Nedelec<dim>(0), 1);
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe_system);
+    test_vec.reinit(dof_handler.n_dofs());
+    apply_boundary_values(dof_handler, constraints, 6, 3, mapping, test_vec);
+    if (test_boundary_values(dof_handler, mapping, fe_system, 6, 3, test_vec))
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+  }
+  {
+    // Test 5: FESystem(FE_Nedelec(0), FE_Nedelec(1))
+    std::cout << "Testing FESystem(FE_Nedelec(0), FE_Nedelec(1)): ";
+    FESystem<dim>   fe_system(FE_Nedelec<dim>(0), 1, FE_Nedelec<dim>(1), 1);
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe_system);
+    test_vec.reinit(dof_handler.n_dofs());
+    apply_boundary_values(dof_handler, constraints, 6, 3, mapping, test_vec);
+    if (test_boundary_values(dof_handler, mapping, fe_system, 6, 3, test_vec))
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+  }
+  {
+    // Test 6: FESystem(FE_Nedelec(0), FE_Q(1))
+    std::cout << "Testing FESystem(FE_Nedelec(0), FE_Q(1)): ";
+    FESystem<dim>   fe_system(FE_Nedelec<dim>(0), 1, FE_Q<dim>(1), 1);
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe_system);
+    test_vec.reinit(dof_handler.n_dofs());
+    apply_boundary_values(dof_handler, constraints, 4, 0, mapping, test_vec);
+    if (test_boundary_values(dof_handler, mapping, fe_system, 4, 0, test_vec))
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+  }
+  {
+    // Test 7: FESystem(FE_Q(1), FE_Nedelec(0))
+    std::cout << "Testing FESystem(FE_Q(1), FE_Nedelec(0)): ";
+    FESystem<dim>   fe_system(FE_Q<dim>(1), 1, FE_Nedelec<dim>(0), 1);
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe_system);
+    test_vec.reinit(dof_handler.n_dofs());
+    apply_boundary_values(dof_handler, constraints, 4, 1, mapping, test_vec);
+    if (test_boundary_values(dof_handler, mapping, fe_system, 4, 1, test_vec))
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+  }
+  {
+    // Test 8: FESystem(FE_Nedelec(0), FE_Q(2))
+    std::cout << "Testing FESystem(FE_Nedelec(0), FE_Q(2)): ";
+    FESystem<dim>   fe_system(FE_Nedelec<dim>(0), 1, FE_Q<dim>(2), 1);
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe_system);
+    test_vec.reinit(dof_handler.n_dofs());
+    apply_boundary_values(dof_handler, constraints, 4, 0, mapping, test_vec);
+    if (test_boundary_values(dof_handler, mapping, fe_system, 4, 0, test_vec))
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+  }
+  {
+    // Test 9: FESystem(FE_Q(2), FE_Nedelec(0))
+    std::cout << "Testing FESystem(FE_Q(2), FE_Nedelec(0)): ";
+    FESystem<dim>   fe_system(FE_Q<dim>(2), 1, FE_Nedelec<dim>(0), 1);
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe_system);
+    test_vec.reinit(dof_handler.n_dofs());
+    apply_boundary_values(dof_handler, constraints, 4, 1, mapping, test_vec);
+    if (test_boundary_values(dof_handler, mapping, fe_system, 4, 1, test_vec))
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+  }
+  {
+    // Test 10: FESystem(FE_Nedelec(1), FE_Q(1))
+    std::cout << "Testing FESystem(FE_Nedelec(1), FE_Q(1)): ";
+    FESystem<dim>   fe_system(FE_Nedelec<dim>(1), 1, FE_Q<dim>(1), 1);
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe_system);
+    test_vec.reinit(dof_handler.n_dofs());
+    apply_boundary_values(dof_handler, constraints, 4, 0, mapping, test_vec);
+    if (test_boundary_values(dof_handler, mapping, fe_system, 4, 0, test_vec))
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+  }
+  {
+    // Test 11: FESystem(FE_Q(1), FE_Nedelec(1))
+    std::cout << "Testing FESystem(FE_Q(1), FE_Nedelec(1)): ";
+    FESystem<dim>   fe_system(FE_Q<dim>(1), 1, FE_Nedelec<dim>(1), 1);
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe_system);
+    test_vec.reinit(dof_handler.n_dofs());
+    apply_boundary_values(dof_handler, constraints, 4, 1, mapping, test_vec);
+    if (test_boundary_values(dof_handler, mapping, fe_system, 4, 1, test_vec))
+      {
+        std::cout << "OK" << std::endl;
+      }
+    else
+      {
+        std::cout << "Failed" << std::endl;
+        abort();
+      }
+  }
+
+  return 0;
+}
diff --git a/tests/numerics/project_bv_curl_conf_04.output b/tests/numerics/project_bv_curl_conf_04.output
new file mode 100644 (file)
index 0000000..b8d4974
--- /dev/null
@@ -0,0 +1,15 @@
+Testing FE_Nedelec(0): OK
+Testing FESystem(FE_Nedelec(0)): OK
+Checking Consistency: OK
+Testing FE_Nedelec(1): OK
+Testing FESystem(FE_Nedelec(1)): OK
+Checking Consistency: OK
+Testing FESystem(FE_Nedelec(0), FE_Nedelec(0)): OK
+Testing FESystem(FE_Nedelec(1), FE_Nedelec(0)): OK
+Testing FESystem(FE_Nedelec(0), FE_Nedelec(1)): OK
+Testing FESystem(FE_Nedelec(0), FE_Q(1)): OK
+Testing FESystem(FE_Q(1), FE_Nedelec(0)): OK
+Testing FESystem(FE_Nedelec(0), FE_Q(2)): OK
+Testing FESystem(FE_Q(2), FE_Nedelec(0)): OK
+Testing FESystem(FE_Nedelec(1), FE_Q(1)): OK
+Testing FESystem(FE_Q(1), FE_Nedelec(1)): OK

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.