<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.
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],
// do
if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
return;
-
+
// this is only
// implemented, if the
// FE is a Nedelec
// 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])))
{
// do
if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
return;
-
+
// this is only
// implemented, if the
// FE is a Nedelec
// 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];
// do
if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
return;
-
+
// This is only
// implemented, if the
// FE is a Nedelec
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])))
{
// do
if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
return;
-
+
// This is only
// implemented, if the
// FE is a Nedelec
}
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];
// do
if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
return;
-
+
// This is only
// implemented, if the
// FE is a Raviart-Thomas
// 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);
// 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)
{
// it in
// the
// current
- // form to
- // make
+ // form
+ // (with
+ // [dim-2])
+ // to make
// sure
// that
// compilers
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);
}
--- /dev/null
+//----------------------------------------------------------------------
+// $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>();
+}