]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix VectorTools::project_boundary_values_div_conforming for more than one FE_RT element
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Thu, 13 Sep 2018 20:53:23 +0000 (22:53 +0200)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Thu, 13 Sep 2018 21:00:58 +0000 (23:00 +0200)
include/deal.II/numerics/vector_tools.templates.h
tests/numerics/project_bv_div_conf_02.cc [new file with mode: 0644]
tests/numerics/project_bv_div_conf_02.output [new file with mode: 0644]

index 755f61f47f3bbefa8cb256985759d5df6070fe17..2a06dd02d0b123e0b0b29e45a2b8ea185ab26289 100644 (file)
@@ -6411,9 +6411,14 @@ namespace VectorTools
 
           for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
             dof_values(i) +=
-              tmp *
-              (normals[q_point] *
-               fe_values[vec].value(fe.face_to_cell_index(i, face), q_point));
+              tmp * (normals[q_point] *
+                     fe_values[vec].value(
+                       fe.face_to_cell_index(i,
+                                             face,
+                                             cell->face_orientation(face),
+                                             cell->face_flip(face),
+                                             cell->face_rotation(face)),
+                       q_point));
         }
 
       std::vector<types::global_dof_index> face_dof_indices(fe.dofs_per_face);
@@ -6424,7 +6429,13 @@ namespace VectorTools
       // Copy the computed values in the AffineConstraints only, if the degree
       // of freedom is not already constrained.
       for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
-        if (!(constraints.is_constrained(face_dof_indices[i])))
+        if (!(constraints.is_constrained(face_dof_indices[i])) &&
+            fe.get_nonzero_components(fe.face_to_cell_index(
+              i,
+              face,
+              cell->face_orientation(face),
+              cell->face_flip(face),
+              cell->face_rotation(face)))[first_vector_component])
           {
             constraints.add_line(face_dof_indices[i]);
 
@@ -6509,9 +6520,14 @@ namespace VectorTools
 
           for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
             dof_values_local(i) +=
-              tmp *
-              (normals[q_point] *
-               fe_values[vec].value(fe.face_to_cell_index(i, face), q_point));
+              tmp * (normals[q_point] *
+                     fe_values[vec].value(
+                       fe.face_to_cell_index(i,
+                                             face,
+                                             cell->face_orientation(face),
+                                             cell->face_flip(face),
+                                             cell->face_rotation(face)),
+                       q_point));
         }
 
       std::vector<types::global_dof_index> face_dof_indices(fe.dofs_per_face);
@@ -6520,7 +6536,13 @@ namespace VectorTools
                                         cell->active_fe_index());
 
       for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
-        if (projected_dofs[face_dof_indices[i]] < fe.degree)
+        if (projected_dofs[face_dof_indices[i]] < fe.degree &&
+            fe.get_nonzero_components(fe.face_to_cell_index(
+              i,
+              face,
+              cell->face_orientation(face),
+              cell->face_flip(face),
+              cell->face_rotation(face)))[first_vector_component])
           {
             dof_values[face_dof_indices[i]]     = dof_values_local(i);
             projected_dofs[face_dof_indices[i]] = fe.degree;
diff --git a/tests/numerics/project_bv_div_conf_02.cc b/tests/numerics/project_bv_div_conf_02.cc
new file mode 100644 (file)
index 0000000..83e61bd
--- /dev/null
@@ -0,0 +1,132 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Test VectorTools::project_boundary_values_div_conforming for the
+// case that the DoFHandler constains is more than one FE_RaviartThomas element.
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_raviart_thomas.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/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+template <int dim>
+class BoundaryFunctionDisp : public dealii::Function<dim>
+{
+public:
+  BoundaryFunctionDisp()
+    : Function<dim>(dim)
+  {}
+
+  virtual double
+  value(const Point<dim> & /*point*/,
+        const unsigned int /*component*/) const override
+  {
+    return 1.;
+  }
+};
+
+template <int dim>
+class BoundaryFunctionVelo : public dealii::Function<dim>
+{
+public:
+  BoundaryFunctionVelo()
+    : Function<dim>(dim)
+  {}
+
+  virtual double
+  value(const Point<dim> & /*point*/,
+        const unsigned int /* component */) const override
+  {
+    return -1.;
+  }
+};
+
+template <int dim>
+void
+test_boundary_values(const FiniteElement<dim> &fe)
+{
+  Triangulation<dim> triangulation;
+
+  GridGenerator::hyper_cube(triangulation, -1., 1., true);
+  triangulation.refine_global(1);
+
+  DoFHandler<dim> dof_handler(triangulation);
+
+  dof_handler.distribute_dofs(fe);
+
+  BoundaryFunctionDisp<dim> boundary_function_disp;
+  BoundaryFunctionVelo<dim> boundary_function_velo;
+
+  AffineConstraints<double> constraints;
+
+  {
+    constraints.clear();
+    VectorTools::project_boundary_values_div_conforming(
+      dof_handler,
+      0, /*first_vector_component*/
+      boundary_function_disp,
+      0, /*bdry_id*/
+      constraints);
+    VectorTools::project_boundary_values_div_conforming(
+      dof_handler,
+      dim, /*first_vector_component*/
+      boundary_function_velo,
+      0, /*bdry_id*/
+      constraints);
+    constraints.close();
+  }
+
+  constraints.print(deallog.get_file_stream());
+}
+
+int
+main()
+{
+  initlog();
+  {
+    constexpr unsigned int dim = 2;
+    FE_RaviartThomas<dim>  u(1);
+    FE_DGQ<dim>            p(1);
+    FESystem<dim>          fesys(u, 2, p, 1);
+    test_boundary_values(fesys);
+  }
+
+  {
+    constexpr unsigned int dim = 3;
+    FE_RaviartThomas<dim>  u(1);
+    FE_DGQ<dim>            p(1);
+    FESystem<dim>          fesys(u, 2, p, 1);
+    test_boundary_values(fesys);
+  }
+}
diff --git a/tests/numerics/project_bv_div_conf_02.output b/tests/numerics/project_bv_div_conf_02.output
new file mode 100644 (file)
index 0000000..96f2a10
--- /dev/null
@@ -0,0 +1,41 @@
+
+    0 = 1.00000
+    1 = 0
+    2 = -1.00000
+    3 = 0
+    52 = 1.00000
+    53 = 0
+    54 = -1.00000
+    55 = 0
+    0 = 1.00000
+    1 = 0
+    2 = 0
+    3 = 0
+    4 = -1.00000
+    5 = 0
+    6 = 0
+    7 = 0
+    152 = 1.00000
+    153 = 0
+    154 = 0
+    155 = 0
+    156 = -1.00000
+    157 = 0
+    158 = 0
+    159 = 0
+    288 = 1.00000
+    289 = 0
+    290 = 0
+    291 = 0
+    292 = -1.00000
+    293 = 0
+    294 = 0
+    295 = 0
+    424 = 1.00000
+    425 = 0
+    426 = 0
+    427 = 0
+    428 = -1.00000
+    429 = 0
+    430 = 0
+    431 = 0

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.