From 3fd37f27535ad816b010a9238c67fd0903cb6d28 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 10 Aug 2011 22:53:44 +0000 Subject: [PATCH] Fix a bug in compute_no_normal_flux_constraints. git-svn-id: https://svn.dealii.org/trunk@24044 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 6 ++ .../deal.II/numerics/vectors.templates.h | 57 ++++++++---- tests/deal.II/no_flux_07.cc | 89 +++++++++++++++++++ tests/deal.II/no_flux_07/cmp/generic | 55 ++++++++++++ 4 files changed, 191 insertions(+), 16 deletions(-) create mode 100644 tests/deal.II/no_flux_07.cc create mode 100644 tests/deal.II/no_flux_07/cmp/generic diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index f27eaee2dc..6ef4cab8f3 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -252,6 +252,12 @@ and DoF handlers embedded in higher dimensional space have been added.

Specific improvements

    +
  1. 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. +
    +(Wolfgang Bangerth, Thomas Geenen, Timo Heister, 2011/08/10)
  2. New: The function FullMatrix::triple_product() adds triple products like Schur complements to existing matrices. diff --git a/deal.II/include/deal.II/numerics/vectors.templates.h b/deal.II/include/deal.II/numerics/vectors.templates.h index 8c06e64637..6c12bd2f89 100644 --- a/deal.II/include/deal.II/numerics/vectors.templates.h +++ b/deal.II/include/deal.II/numerics/vectors.templates.h @@ -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::epsilon()) constraints.add_entry (dof_indices.dof_indices[1], @@ -3177,7 +3177,7 @@ project_boundary_values_curl_conforming (const DoFHandler& dof_handler, // do if (dynamic_cast*> (&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& 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& dof_handler, // do if (dynamic_cast*> (&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& 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& dof_handler, // do if (dynamic_cast*> (&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& 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& dof_handler, // do if (dynamic_cast*> (&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& 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& dof_ // do if (dynamic_cast*> (&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 &dof_handler, // then construct constraints // from this: const internal::VectorTools::VectorDoFTuple & - 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 &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 &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 &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 index 0000000000..9c9263ffb6 --- /dev/null +++ b/tests/deal.II/no_flux_07.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +template +void test (const Triangulation& tr, + const FiniteElement& fe) +{ + DoFHandler dof(tr); + dof.distribute_dofs(fe); + + std::set 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 +void test_hyper_cube() +{ + Triangulation tr; + + // create a hypercube, then shear + // its top surface to make the + // result non-square + GridGenerator::hyper_cube(tr); + Point shift; + for (unsigned int i=GeometryInfo::vertices_per_cell/2; + i < GeometryInfo::vertices_per_cell; ++i) + tr.begin_active()->vertex(i)[0] += 0.5; + + FESystem fe (FE_Q(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 index 0000000000..1109ca2ea3 --- /dev/null +++ b/tests/deal.II/no_flux_07/cmp/generic @@ -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 -- 2.39.5