]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a test. 16884/head
authorDavid Wells <drwells@email.unc.edu>
Sat, 13 Apr 2024 11:43:08 +0000 (07:43 -0400)
committerDavid Wells <drwells@email.unc.edu>
Sat, 13 Apr 2024 11:52:50 +0000 (07:52 -0400)
This fails on master due to twisted elements.

tests/simplex/orientation_06.cc [new file with mode: 0644]
tests/simplex/orientation_06.output [new file with mode: 0644]

diff --git a/tests/simplex/orientation_06.cc b/tests/simplex/orientation_06.cc
new file mode 100644 (file)
index 0000000..b71172c
--- /dev/null
@@ -0,0 +1,223 @@
+/* ------------------------------------------------------------------------
+ *
+ * SPDX-License-Identifier: LGPL-2.1-or-later
+ * Copyright (C) 2023 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * Part of the source code is dual licensed under Apache-2.0 WITH
+ * LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+ * governing the source code and code contributions can be found in
+ * LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+ *
+ * ------------------------------------------------------------------------
+ */
+
+// Make sure that we get the correct convergence order for a combination of P3
+// and Q3 elements and and two reversed quadrilateral line orientations.
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools_topology.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools_boundary.h>
+#include <deal.II/numerics/vector_tools_integrate_difference.h>
+#include <deal.II/numerics/vector_tools_interpolate.h>
+
+#include "../tests.h"
+
+// u
+class Exact : public Function<2>
+{
+public:
+  double
+  value(const Point<2> &p, unsigned int /*component*/ = 0) const
+  {
+    return std::sin(p[0]) * std::sin(p[1]);
+  }
+};
+
+// -Laplacian(u) + u
+class Forcing : public Function<2>
+{
+public:
+  double
+  value(const Point<2> &p, unsigned int /*component*/ = 0) const
+  {
+    return 3.0 * std::sin(p[0]) * std::sin(p[1]);
+  }
+};
+
+int
+main()
+{
+  initlog();
+
+  const std::vector<Point<2>> vertices{{0.0, 0.0},
+                                       {1.0, 0.0},
+                                       {0.0, 2.0},
+                                       {3.0, 1.0},
+                                       {3.0, 3.0},
+                                       {2.0, 3.0},
+                                       {0.0, 3.0},
+                                       {3.0, 0.0}};
+
+  // faces 1 and 3 of the quadrilateral will have reversed orientations
+  std::vector<CellData<2>> cells(5);
+  cells[0].vertices = {0u, 1u, 2u};
+  cells[1].vertices = {3u, 4u, 5u};
+  cells[2].vertices = {1u, 3u, 7u};
+  cells[3].vertices = {5u, 2u, 6u};
+  cells[4].vertices = {1u, 3u, 2u, 5u};
+
+  const hp::MappingCollection<2> mapping(MappingFE<2>(FE_SimplexP<2>(1)),
+                                         MappingFE<2>(FE_Q<2>(1)));
+
+  const hp::FECollection<2> fe(FE_SimplexP<2>(3),
+                               FE_Q<2>(QIterated<1>(QTrapezoid<1>(), 3)));
+  const hp::QCollection<2>  quadrature_formula(QGaussSimplex<2>(4),
+                                              QGauss<2>(4));
+
+  Triangulation<2> tria;
+  GridTools::invert_cells_with_negative_measure(vertices, cells);
+  tria.create_triangulation(vertices, cells, SubCellData());
+
+  for (unsigned int r = 0; r < 5; ++r)
+    {
+      tria.refine_global(1);
+      DoFHandler<2> dof_handler(tria);
+      for (const auto &cell : dof_handler.active_cell_iterators())
+        {
+          if (cell->reference_cell() == ReferenceCells::Triangle)
+            cell->set_active_fe_index(0);
+          else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
+            cell->set_active_fe_index(1);
+        }
+      dof_handler.distribute_dofs(fe);
+
+      // set up matrices and vectors:
+      AffineConstraints<double> constraints(dof_handler.locally_owned_dofs(),
+                                            dof_handler.locally_owned_dofs());
+      VectorTools::interpolate_boundary_values(
+        mapping, dof_handler, 0, Exact(), constraints);
+      constraints.close();
+      SparsityPattern        sparsity_pattern;
+      DynamicSparsityPattern dsp(dof_handler.n_dofs());
+      DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
+      sparsity_pattern.copy_from(dsp);
+
+      SparseMatrix<double> system_matrix(sparsity_pattern);
+      Vector<double>       system_rhs(dof_handler.n_dofs());
+      Vector<double>       solution(dof_handler.n_dofs());
+
+      hp::FEValues<2> hp_fe_values(mapping,
+                                   fe,
+                                   quadrature_formula,
+                                   update_values | update_gradients |
+                                     update_quadrature_points |
+                                     update_JxW_values);
+
+      FullMatrix<double>                   cell_matrix;
+      Vector<double>                       cell_rhs;
+      std::vector<types::global_dof_index> local_dof_indices;
+
+      Forcing forcing;
+      for (const auto &cell : dof_handler.active_cell_iterators())
+        {
+          hp_fe_values.reinit(cell);
+          const auto &fe_values = hp_fe_values.get_present_fe_values();
+
+          const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
+          cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+          cell_rhs.reinit(dofs_per_cell);
+          local_dof_indices.resize(dofs_per_cell);
+          cell->get_dof_indices(local_dof_indices);
+
+          cell_matrix = 0;
+          cell_rhs    = 0;
+
+          for (const unsigned int q_index :
+               fe_values.quadrature_point_indices())
+            {
+              for (const unsigned int i : fe_values.dof_indices())
+                for (const unsigned int j : fe_values.dof_indices())
+                  cell_matrix(i, j) += (fe_values.shape_grad(i, q_index) *
+                                          fe_values.shape_grad(j, q_index) +
+                                        fe_values.shape_value(i, q_index) *
+                                          fe_values.shape_value(j, q_index)) *
+                                       fe_values.JxW(q_index);
+
+              for (const unsigned int i : fe_values.dof_indices())
+                cell_rhs(i) +=
+                  fe_values.shape_value(i, q_index) *
+                  forcing.value(fe_values.quadrature_point(q_index)) *
+                  fe_values.JxW(q_index);
+            }
+
+          constraints.distribute_local_to_global(cell_matrix,
+                                                 cell_rhs,
+                                                 local_dof_indices,
+                                                 system_matrix,
+                                                 system_rhs);
+        }
+
+      SolverControl            solver_control(1000, 1e-12);
+      SolverCG<Vector<double>> solver(solver_control);
+
+      // Cut down on a lot of solver iterations by interpolating the solution as
+      // the initial guess:
+      VectorTools::interpolate(mapping, dof_handler, Exact(), solution);
+
+      PreconditionSSOR<> precondition;
+      precondition.initialize(system_matrix,
+                              PreconditionSSOR<>::AdditionalData());
+      deallog << "rhs l2 = " << system_rhs.l2_norm() << std::endl;
+      solver.solve(system_matrix, solution, system_rhs, precondition);
+      constraints.distribute(solution);
+
+      Vector<double> cell_errors(tria.n_active_cells());
+      VectorTools::integrate_difference(mapping,
+                                        dof_handler,
+                                        solution,
+                                        Exact(),
+                                        cell_errors,
+                                        quadrature_formula,
+                                        VectorTools::L2_norm);
+
+      deallog << "global error = "
+              << VectorTools::compute_global_error(tria,
+                                                   cell_errors,
+                                                   VectorTools::L2_norm)
+              << std::endl;
+
+      std::ofstream out("out-" + std::to_string(r) + ".vtu");
+      DataOut<2>    data_out;
+      data_out.attach_dof_handler(dof_handler);
+      data_out.add_data_vector(solution, "solution");
+      data_out.add_data_vector(cell_errors, "cell_errors");
+      data_out.build_patches(mapping);
+      data_out.write_vtu(out);
+    }
+}
diff --git a/tests/simplex/orientation_06.output b/tests/simplex/orientation_06.output
new file mode 100644 (file)
index 0000000..c921f99
--- /dev/null
@@ -0,0 +1,21 @@
+
+DEAL::rhs l2 = 2.27397
+DEAL:cg::Starting value 1.73734
+DEAL:cg::Convergence step 23 value 8.54649e-13
+DEAL::global error = 0.00166288
+DEAL::rhs l2 = 2.08466
+DEAL:cg::Starting value 2.45428
+DEAL:cg::Convergence step 39 value 8.26837e-13
+DEAL::global error = 0.000105423
+DEAL::rhs l2 = 2.69354
+DEAL:cg::Starting value 3.47449
+DEAL:cg::Convergence step 67 value 7.08130e-13
+DEAL::global error = 6.61142e-06
+DEAL::rhs l2 = 3.77050
+DEAL:cg::Starting value 4.91796
+DEAL:cg::Convergence step 105 value 9.99712e-13
+DEAL::global error = 4.13368e-07
+DEAL::rhs l2 = 5.33304
+DEAL:cg::Starting value 6.95866
+DEAL:cg::Convergence step 169 value 9.77791e-13
+DEAL::global error = 2.58314e-08

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.