From: bangerth Date: Mon, 26 Nov 2012 04:57:29 +0000 (+0000) Subject: Add a version of DoFTools::make_zero_boundary_constraints that works on only part... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5f95d45fca4f01bec66f0d13d3c1287efa31ea97;p=dealii-svn.git Add a version of DoFTools::make_zero_boundary_constraints that works on only part of the boundary. git-svn-id: https://svn.dealii.org/trunk@27683 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 3a5b562c21..ec306a6bdc 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -154,6 +154,11 @@ DoFHandler, in particular removal of specializations.

Specific improvements

    +
  1. New: There is now a version of DoFTools::make_zero_boundary_constraints() +that accepts a boundary indicator as argument. +
    +(Wolfgang Bangerth, 2012/11/25) +
  2. 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 diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h index 7fb3d61df6..c39b590c47 100644 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@ -2655,39 +2655,67 @@ namespace DoFTools /** * 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 class DH> void make_zero_boundary_constraints (const DH &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 class DH> + void + make_zero_boundary_constraints (const DH &dof, + ConstraintMatrix &zero_boundary_constraints, + const ComponentMask &component_mask = ComponentMask()); + + /** * Map a coupling table from the * user friendly organization by diff --git a/deal.II/source/dofs/dof_tools.cc b/deal.II/source/dofs/dof_tools.cc index 8a08604065..b9da248100 100644 --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@ -6167,6 +6167,7 @@ namespace DoFTools template class DH> void make_zero_boundary_constraints (const DH &dof, + const types::boundary_id boundary_indicator, ConstraintMatrix &zero_boundary_constraints, const ComponentMask &component_mask) { @@ -6194,10 +6195,15 @@ namespace DoFTools { const FiniteElement &fe = cell->get_fe(); - typename DH::face_iterator face = cell->face(face_no); + const typename DH::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 @@ -6231,6 +6237,19 @@ namespace DoFTools } } + + + template class DH> + void + make_zero_boundary_constraints (const DH &dof, + ConstraintMatrix &zero_boundary_constraints, + const ComponentMask &component_mask) + { + make_zero_boundary_constraints(dof, numbers::invalid_boundary_id, + zero_boundary_constraints, component_mask); + } + + template void make_cell_patches( Sparsity &block_list, diff --git a/deal.II/source/dofs/dof_tools.inst.in b/deal.II/source/dofs/dof_tools.inst.in index b0f889f3f6..1b9a510043 100644 --- a/deal.II/source/dofs/dof_tools.inst.in +++ b/deal.II/source/dofs/dof_tools.inst.in @@ -1010,6 +1010,14 @@ DoFTools::make_zero_boundary_constraints ConstraintMatrix &, const ComponentMask &); +template +void +DoFTools::make_zero_boundary_constraints +(const DoFHandler &, + const types::boundary_id , + ConstraintMatrix &, + const ComponentMask &); + template void DoFTools::make_zero_boundary_constraints @@ -1017,6 +1025,15 @@ DoFTools::make_zero_boundary_constraints ConstraintMatrix &, const ComponentMask &); +template +void +DoFTools::make_zero_boundary_constraints +(const hp::DoFHandler &, + const types::boundary_id , + ConstraintMatrix &, + const ComponentMask &); + + #if deal_II_dimension < 3 template diff --git a/tests/bits/make_boundary_constraints_02.cc b/tests/bits/make_boundary_constraints_02.cc new file mode 100644 index 0000000000..4e5ab2e5ad --- /dev/null +++ b/tests/bits/make_boundary_constraints_02.cc @@ -0,0 +1,188 @@ +//---------------------------- 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +void test () +{ + deallog << dim << "D" << std::endl; + + Triangulation triangulation; + const Point p2 = (dim == 1 ? + Point(1.) : + (dim == 2 ? + Point(1.,1.) : + Point(1.,1.,1.))); + GridGenerator::hyper_rectangle (triangulation, + Point(), 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::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 fe(1); + DoFHandler 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 A(sparsity), B(sparsity); + Vector a1 (dof_handler.n_dofs()); + + // initialize object denoting zero + // boundary values and boundary + // constraints + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 1, + ConstantFunction(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 local_dofs (fe.dofs_per_cell); + FullMatrix local_matrix (fe.dofs_per_cell, fe.dofs_per_cell); + Vector local_vector (fe.dofs_per_cell); + for (typename DoFHandler::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; idistribute_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; + }; +} diff --git a/tests/bits/make_boundary_constraints_02/cmp/generic b/tests/bits/make_boundary_constraints_02/cmp/generic new file mode 100644 index 0000000000..e6ee237e6b --- /dev/null +++ b/tests/bits/make_boundary_constraints_02/cmp/generic @@ -0,0 +1,16 @@ + +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