]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix a bug in compute_no_normal_flux_constraints.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 10 Aug 2011 22:53:44 +0000 (22:53 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 10 Aug 2011 22:53:44 +0000 (22:53 +0000)
git-svn-id: https://svn.dealii.org/trunk@24044 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/numerics/vectors.templates.h
tests/deal.II/no_flux_07.cc [new file with mode: 0644]
tests/deal.II/no_flux_07/cmp/generic [new file with mode: 0644]

index f27eaee2dc60e4c5d5ea815e47f70c537294ef39..6ef4cab8f353b255c2982e4ee48b46ce7fb0e926 100644 (file)
@@ -252,6 +252,12 @@ and DoF handlers embedded in higher dimensional space have been added.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Fixed: The function VectorTools::compute_no_normal_flux_constraints had
+a bug that led to an exception whenever we were computing constraints for
+vector fields located on edges shared between two faces of a 3d cell if those
+faces were not perpendicular. This is now fixed.
+<br>
+(Wolfgang Bangerth, Thomas Geenen, Timo Heister, 2011/08/10)
 
 <li> New: The function FullMatrix::triple_product() adds triple products
 like Schur complements to existing matrices.
index 8c06e646370f3f6119654f1232c346752be9ccca..6c12bd2f8916de76defbfb700a180efb53079141 100644 (file)
@@ -2286,7 +2286,7 @@ namespace internal
                        constraints.add_entry (dof_indices.dof_indices[1],
                                               dof_indices.dof_indices[0],
                                               -constraining_vector[0]/constraining_vector[1]);
-                     
+
                      if (std::fabs (constraining_vector[2]/constraining_vector[1])
                          > std::numeric_limits<double>::epsilon())
                        constraints.add_entry (dof_indices.dof_indices[1],
@@ -3177,7 +3177,7 @@ project_boundary_values_curl_conforming (const DoFHandler<dim>& dof_handler,
                                                       // do
                   if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
                     return;
-                  
+
                                                    // this is only
                                                    // implemented, if the
                                                    // FE is a Nedelec
@@ -3208,7 +3208,7 @@ project_boundary_values_curl_conforming (const DoFHandler<dim>& dof_handler,
                                                    // freedom is not
                                                    // already constrained.
                   const double tol = 1e-13;
-                  
+
                   for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
                     if (!(constraints.is_constrained (face_dof_indices[dof])))
                       {
@@ -3261,7 +3261,7 @@ project_boundary_values_curl_conforming (const DoFHandler<dim>& dof_handler,
                                                       // do
                   if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
                     return;
-                  
+
                                                    // this is only
                                                    // implemented, if the
                                                    // FE is a Nedelec
@@ -3358,7 +3358,7 @@ project_boundary_values_curl_conforming (const DoFHandler<dim>& dof_handler,
                                                    // values in the global
                                                    // vector.
                   const double tol = 0.5 * superdegree * 1e-13 / cell->face (face)->diameter ();
-                  
+
                   for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
                     if (std::abs (computed_constraints[face_dof_indices[dof]] - dof_values[dof]) > tol)
                       computed_constraints[face_dof_indices[dof]] = dof_values[dof];
@@ -3433,7 +3433,7 @@ project_boundary_values_curl_conforming (const hp::DoFHandler<dim>& dof_handler,
                                                       // do
                   if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
                     return;
-                  
+
                                                    // This is only
                                                    // implemented, if the
                                                    // FE is a Nedelec
@@ -3465,7 +3465,7 @@ project_boundary_values_curl_conforming (const hp::DoFHandler<dim>& dof_handler,
                                                       cell->active_fe_index ());
 
                   const double tol = 0.5 * cell->get_fe ().degree * 1e-13  / cell->face (face)->diameter ();
-                  
+
                   for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
                     if (!(constraints.is_constrained (face_dof_indices[dof])))
                       {
@@ -3522,7 +3522,7 @@ project_boundary_values_curl_conforming (const hp::DoFHandler<dim>& dof_handler,
                                                       // do
                   if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
                     return;
-                  
+
                                                    // This is only
                                                    // implemented, if the
                                                    // FE is a Nedelec
@@ -3589,7 +3589,7 @@ project_boundary_values_curl_conforming (const hp::DoFHandler<dim>& dof_handler,
                     }
 
                   const double tol = 0.5 * superdegree * 1e-13  / cell->face (face)->diameter ();
-                  
+
                   for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
                     if (std::abs (computed_constraints[face_dof_indices[dof]] - dof_values[dof]) > tol)
                       computed_constraints[face_dof_indices[dof]] = dof_values[dof];
@@ -3859,7 +3859,7 @@ VectorTools::project_boundary_values_div_conforming (const DoFHandler<dim>& dof_
                                                       // do
                   if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
                     return;
-                  
+
                                                    // This is only
                                                    // implemented, if the
                                                    // FE is a Raviart-Thomas
@@ -4417,7 +4417,7 @@ compute_no_normal_flux_constraints (const DH<dim,spacedim>         &dof_handler,
                                             // then construct constraints
                                             // from this:
            const internal::VectorTools::VectorDoFTuple<dim> &
-             dof_indices = same_dof_range[0]->first;         
+             dof_indices = same_dof_range[0]->first;
            internal::VectorTools::add_constraint (dof_indices, normal,
                                                   constraints);
 
@@ -4587,7 +4587,28 @@ compute_no_normal_flux_constraints (const DH<dim,spacedim>         &dof_handler,
                                                 // calculate the
                                                 // tangent as the
                                                 // outer product of
-                                                // the normal vectors
+                                                // the normal
+                                                // vectors. since
+                                                // these vectors do
+                                                // not need to be
+                                                // orthogonal (think,
+                                                // for example, the
+                                                // case of the
+                                                // deal.II/no_flux_07
+                                                // test: a sheared
+                                                // cube in 3d, with
+                                                // Q2 elements, where
+                                                // we have
+                                                // constraints from
+                                                // the two normal
+                                                // vectors of two
+                                                // faces of the
+                                                // sheared cube that
+                                                // are not
+                                                // perpendicular to
+                                                // each other), we
+                                                // have to normalize
+                                                // the outer product
                Tensor<1,dim> tangent;
                switch (dim)
                  {
@@ -4602,8 +4623,10 @@ compute_no_normal_flux_constraints (const DH<dim,spacedim>         &dof_handler,
                                                           // it in
                                                           // the
                                                           // current
-                                                          // form to
-                                                          // make
+                                                          // form
+                                                          // (with
+                                                          // [dim-2])
+                                                          // to make
                                                           // sure
                                                           // that
                                                           // compilers
@@ -4642,8 +4665,10 @@ compute_no_normal_flux_constraints (const DH<dim,spacedim>         &dof_handler,
                          Assert (false, ExcNotImplemented());
                  }
 
-               Assert (std::fabs (tangent.norm()-1) < 1e-12,
-                       ExcInternalError());
+               Assert (std::fabs (tangent.norm()) > 1e-12,
+                       ExcMessage("Two normal vectors from adjacent faces are almost "
+                                  "parallel."));
+               tangent /= tangent.norm();
 
                tangential_vectors.push_back (tangent);
              }
diff --git a/tests/deal.II/no_flux_07.cc b/tests/deal.II/no_flux_07.cc
new file mode 100644 (file)
index 0000000..9c9263f
--- /dev/null
@@ -0,0 +1,89 @@
+//----------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2007, 2008, 2011 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------------------------------------------------
+
+
+// The function VectorTools::compute_no_normal_flux_constraints had a
+// bug that led to an exception whenever we were computing constraints
+// for vector fields located on edges shared between two faces of a 3d
+// cell if those faces were not perpendicular: in a case like that, we
+// computed the unconstrained tangent field as the cross product of
+// the two face normals, but the resulting tangent did not have unit
+// length
+//
+// we can test this using a 3d cube that we shear and by computing
+// constraints for a Q2^3 field for which there are DoFs on edge
+// midpoints.
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/numerics/vectors.h>
+
+#include <fstream>
+
+
+template<int dim>
+void test (const Triangulation<dim>& tr,
+                     const FiniteElement<dim>& fe)
+{
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+
+  std::set<unsigned char> boundary_ids;
+  boundary_ids.insert (0);
+
+  ConstraintMatrix cm;
+  VectorTools::compute_no_normal_flux_constraints (dof, 0, boundary_ids, cm);
+
+  cm.print (deallog.get_file_stream ());
+}
+
+
+template<int dim>
+void test_hyper_cube()
+{
+  Triangulation<dim> tr;
+
+                                  // create a hypercube, then shear
+                                  // its top surface to make the
+                                  // result non-square
+  GridGenerator::hyper_cube(tr);
+  Point<dim> shift;
+  for (unsigned int i=GeometryInfo<dim>::vertices_per_cell/2;
+       i < GeometryInfo<dim>::vertices_per_cell; ++i)
+    tr.begin_active()->vertex(i)[0] += 0.5;
+
+  FESystem<dim> fe (FE_Q<dim>(2), dim);
+  test(tr, fe);
+}
+
+
+int main()
+{
+  std::ofstream logfile ("no_flux_07/output");
+  deallog << std::setprecision (2);
+  deallog << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+  deallog.threshold_double(1.e-12);
+
+  test_hyper_cube<3>();
+}
diff --git a/tests/deal.II/no_flux_07/cmp/generic b/tests/deal.II/no_flux_07/cmp/generic
new file mode 100644 (file)
index 0000000..1109ca2
--- /dev/null
@@ -0,0 +1,55 @@
+
+    0 = 0
+    1 = 0
+    2 = 0
+    3 = 0
+    4 = 0
+    5 = 0
+    6 = 0
+    7 = 0
+    8 = 0
+    9 = 0
+    10 = 0
+    11 = 0
+    12 = 0
+    13 = 0
+    14 = 0
+    15 = 0
+    16 = 0
+    17 = 0
+    18 = 0
+    19 = 0
+    20 = 0
+    21 = 0
+    22 = 0
+    23 = 0
+    24 = 0
+    26 = 0
+    27 = 0
+    29 = 0
+    31 = 0
+    32 = 0
+    34 = 0
+    35 = 0
+    36 = 0
+    38 = 0
+    39 = 0
+    41 = 0
+    43 = 0
+    44 = 0
+    46 = 0
+    47 = 0
+    49 = 0
+    48 50:  0.500000
+    52 = 0
+    51 53:  0.500000
+    55 = 0
+    54 56:  0.500000
+    58 = 0
+    57 59:  0.500000
+    60 62:  0.500000
+    63 65:  0.500000
+    67 = 0
+    70 = 0
+    74 = 0
+    77 = 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.