<h3>Specific improvements</h3>
<ol>
+<li> New: There is now a version of DoFTools::make_zero_boundary_constraints()
+that accepts a boundary indicator as argument.
+<br>
+(Wolfgang Bangerth, 2012/11/25)
+
<li> Fixed: The DoFTools::make_flux_sparsity_pattern() function
had a bug that triggered in 1d whenever there were neighboring
cells that differed in refinement level by more than one. This
/**
* Make a constraint matrix for the
* constraints that result from zero
- * boundary values.
+ * boundary values on the given boundary indicator.
*
* This function constrains all
- * degrees of freedom on the
- * boundary. Optionally, you can
- * add a component mask, which
- * restricts this functionality
- * to a subset of an FESystem.
- *
- * For non-@ref GlossPrimitive "primitive"
- * shape functions, any degree of freedom
- * is affected that belongs to a
- * shape function where at least
- * one of its nonzero components
- * is affected.
- *
- * The last argument indicates which
- * components of the solution
- * vector should be constrained to zero
- * (see @ref GlossComponentMask).
- *
- * This function is used
- * in step-36, for
- * example.
+ * degrees of freedom on the given part of the
+ * boundary.
+ *
+ * A variant of this function with different arguments is used
+ * in step-36.
+ *
+ * @param dof The DoFHandler to work on.
+ * @param boundary_indicator The indicator of that part of the boundary
+ * for which constraints should be computed. If this number equals
+ * numbers::invalid_boundary_id then all boundaries of the domain
+ * will be treated.
+ * @param zero_boundary_constraints The constraint object into which the
+ * constraints will be written. The new constraints due to zero boundary
+ * values will simply be added, preserving any other constraints
+ * previously present. However, this will only work if the previous
+ * content of that object consists of constraints on degrees of freedom
+ * that are not located on the boundary treated here. If there are
+ * previously existing constraints for degrees of freedom located on the
+ * boundary, then this would constitute a conflict. See the @ref constraints
+ * module for handling the case where there are conflicting constraints
+ * on individual degrees of freedom.
+ * @param component_mask An optional component mask that
+ * restricts the functionality of this function to a subset of an FESystem.
+ * For non-@ref GlossPrimitive "primitive"
+ * shape functions, any degree of freedom
+ * is affected that belongs to a
+ * shape function where at least
+ * one of its nonzero components
+ * is affected by the component mask (see @ref GlossComponentMask). If
+ * this argument is omitted, all components of the finite element with
+ * degrees of freedom at the boundary will be considered.
*
* @ingroup constraints
*/
template <int dim, int spacedim, template <int, int> class DH>
void
make_zero_boundary_constraints (const DH<dim,spacedim> &dof,
+ const types::boundary_id boundary_indicator,
ConstraintMatrix &zero_boundary_constraints,
const ComponentMask &component_mask = ComponentMask());
+ /**
+ * Do the same as the previous function, except do it for all
+ * parts of the boundary, not just those with a particular boundary
+ * indicator. This function is then equivalent to calling the previous
+ * one with numbers::invalid_boundary_id as second argument.
+ *
+ * This function is used in step-36, for example.
+ *
+ * @ingroup constraints
+ */
+ template <int dim, int spacedim, template <int, int> class DH>
+ void
+ make_zero_boundary_constraints (const DH<dim,spacedim> &dof,
+ ConstraintMatrix &zero_boundary_constraints,
+ const ComponentMask &component_mask = ComponentMask());
+
+
/**
* Map a coupling table from the
* user friendly organization by
template <int dim, int spacedim, template <int,int> class DH>
void
make_zero_boundary_constraints (const DH<dim, spacedim> &dof,
+ const types::boundary_id boundary_indicator,
ConstraintMatrix &zero_boundary_constraints,
const ComponentMask &component_mask)
{
{
const FiniteElement<dim,spacedim> &fe = cell->get_fe();
- typename DH<dim,spacedim>::face_iterator face = cell->face(face_no);
+ const typename DH<dim,spacedim>::face_iterator face = cell->face(face_no);
- // if face is on the boundary
- if (face->at_boundary ())
+ // if face is on the boundary and satisfies the correct
+ // boundary id property
+ if (face->at_boundary ()
+ &&
+ ((boundary_indicator == numbers::invalid_boundary_id)
+ ||
+ (face->boundary_indicator() == boundary_indicator)))
{
// get indices and physical
// location on this face
}
}
+
+
+ template <int dim, int spacedim, template <int,int> class DH>
+ void
+ make_zero_boundary_constraints (const DH<dim, spacedim> &dof,
+ ConstraintMatrix &zero_boundary_constraints,
+ const ComponentMask &component_mask)
+ {
+ make_zero_boundary_constraints(dof, numbers::invalid_boundary_id,
+ zero_boundary_constraints, component_mask);
+ }
+
+
template <class DH, class Sparsity>
void make_cell_patches(
Sparsity &block_list,
ConstraintMatrix &,
const ComponentMask &);
+template
+void
+DoFTools::make_zero_boundary_constraints
+(const DoFHandler<deal_II_dimension> &,
+ const types::boundary_id ,
+ ConstraintMatrix &,
+ const ComponentMask &);
+
template
void
DoFTools::make_zero_boundary_constraints
ConstraintMatrix &,
const ComponentMask &);
+template
+void
+DoFTools::make_zero_boundary_constraints
+(const hp::DoFHandler<deal_II_dimension> &,
+ const types::boundary_id ,
+ ConstraintMatrix &,
+ const ComponentMask &);
+
+
#if deal_II_dimension < 3
template
--- /dev/null
+//---------------------------- make_boundary_constraints_02.cc ------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2009, 2012 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.
+//
+//---------------------------- make_boundary_constraints_02.cc ------------------
+
+
+// check DoFTools::make_zero_boundary_constraints by comparing
+// apply_boundary_values with a zero function to
+// make_zero_boundary_constraints
+//
+// this is the version of the test with only one boundary
+// indicator selected
+
+#include "../tests.h"
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <fstream>
+#include <iomanip>
+
+
+template <int dim>
+void test ()
+{
+ deallog << dim << "D" << std::endl;
+
+ Triangulation<dim> triangulation;
+ const Point<dim> p2 = (dim == 1 ?
+ Point<dim>(1.) :
+ (dim == 2 ?
+ Point<dim>(1.,1.) :
+ Point<dim>(1.,1.,1.)));
+ GridGenerator::hyper_rectangle (triangulation,
+ Point<dim>(), p2, true);
+
+ // refine the mesh in a random way
+ triangulation.refine_global (4-dim);
+ for (unsigned int i=0; i<11-2*dim; ++i)
+ {
+ typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active();
+ for (unsigned int index=0; cell != triangulation.end(); ++cell, ++index)
+ if (index % (3*dim) == 0)
+ cell->set_refine_flag();
+ triangulation.execute_coarsening_and_refinement ();
+ }
+ deallog << "Number of cells: "
+ << triangulation.n_active_cells() << std::endl;
+
+ // assign quadrature, set up a
+ // DoFHandler, and distribute dofs
+ FE_Q<dim> fe(1);
+ DoFHandler<dim> dof_handler (triangulation);
+ dof_handler.distribute_dofs (fe);
+ deallog << "Number of dofs : "
+ << dof_handler.n_dofs() << std::endl;
+
+ // then set up a sparsity pattern
+ // and two matrices and vectors on
+ // top of it.
+ SparsityPattern sparsity (dof_handler.n_dofs(),
+ dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, sparsity);
+ sparsity.compress ();
+ SparseMatrix<double> A(sparsity), B(sparsity);
+ Vector<double> a1 (dof_handler.n_dofs());
+
+ // initialize object denoting zero
+ // boundary values and boundary
+ // constraints
+ std::map<unsigned int, double> boundary_values;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 1,
+ ConstantFunction<dim>(1.),
+ boundary_values);
+ ConstraintMatrix constraints;
+ constraints.clear();
+ DoFTools::make_zero_boundary_constraints (dof_handler,
+ 1,
+ constraints);
+ constraints.close();
+
+ // then fill two matrices by
+ // setting up bogus matrix entries
+ // and for A applying constraints
+ // right away and for B applying
+ // constraints later on
+ std::vector<unsigned int> local_dofs (fe.dofs_per_cell);
+ FullMatrix<double> local_matrix (fe.dofs_per_cell, fe.dofs_per_cell);
+ Vector<double> local_vector (fe.dofs_per_cell);
+ for (typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active();
+ cell != dof_handler.end(); ++cell)
+ {
+ cell->get_dof_indices (local_dofs);
+ local_matrix = 0;
+
+ // create local matrices
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ local_matrix(i,j) = (i+1.)*(j+1.)*(local_dofs[i]+1.)*(local_dofs[j]+1.);
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ local_vector(i) = (i+1.)*(local_dofs[i]+1.);
+
+ // for matrix A and vector a1 apply boundary
+ // values
+ MatrixTools::local_apply_boundary_values (boundary_values, local_dofs,
+ local_matrix, local_vector,
+ true);
+ cell->distribute_local_to_global (local_matrix, A);
+ cell->distribute_local_to_global (local_vector, a1);
+
+ // for matrix B cast constraint
+ // matrix
+ constraints.distribute_local_to_global (local_matrix,
+ local_dofs,
+ B);
+ }
+
+ // here comes the check: compare
+ // the l1_norm of matrices A and B,
+ // their difference should be zero
+ deallog << "|A| = " << A.l1_norm() << std::endl;
+ deallog << "|B| = " << B.l1_norm() << std::endl;
+ Assert (A.l1_norm() - B.l1_norm() == 0,
+ ExcInternalError());
+}
+
+int main ()
+{
+ std::ofstream logfile("make_boundary_constraints_02/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ test<1> ();
+ test<2> ();
+ test<3> ();
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::1D
+DEAL::Number of cells: 115
+DEAL::Number of dofs : 116
+DEAL::|A| = 118784.
+DEAL::|B| = 118784.
+DEAL::2D
+DEAL::Number of cells: 1186
+DEAL::Number of dofs : 1537
+DEAL::|A| = 2.33581e+08
+DEAL::|B| = 2.33581e+08
+DEAL::3D
+DEAL::Number of cells: 736
+DEAL::Number of dofs : 1311
+DEAL::|A| = 2.10521e+09
+DEAL::|B| = 2.10521e+09
projects / dealii-svn.git / commitdiff