From f0ca7d4ec56420b1428da2d6ef9e3fff1d71ae80 Mon Sep 17 00:00:00 2001 From: kanschat Date: Wed, 17 Mar 2010 17:44:27 +0000 Subject: [PATCH] multigrid tests improved git-svn-id: https://svn.dealii.org/trunk@20841 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/multigrid/mg_constraints.h | 111 ++++++ .../include/multigrid/mg_dof_accessor.h | 5 +- deal.II/deal.II/source/multigrid/mg_tools.cc | 334 ++++++++++-------- deal.II/examples/step-39/step-39.cc | 2 +- 4 files changed, 303 insertions(+), 149 deletions(-) create mode 100644 deal.II/deal.II/include/multigrid/mg_constraints.h diff --git a/deal.II/deal.II/include/multigrid/mg_constraints.h b/deal.II/deal.II/include/multigrid/mg_constraints.h new file mode 100644 index 0000000000..dd009b95b6 --- /dev/null +++ b/deal.II/deal.II/include/multigrid/mg_constraints.h @@ -0,0 +1,111 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2010 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. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__mg_constraints_h +#define __deal2__mg_constraints_h + +#include +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +/** + * Collection of boundary constraints and refinement edge constraints + * for level vectors. + */ +class MGConstraints : public Subscriptor +{ + public: + /** + * Fill the internal data + * structures with values + * extracted from the dof + * handler. + * + * This function leaves + * #boundary_indices empty, since + * no boundary values are + * provided. + */ + template + void initialize(const MGDoFHandler& dof); + + /** + * Fill the internal data + * structures with values + * extracted from the dof + * handler, applying the boundary + * values provided. + */ + template + void initialize(const MGDoFHandler& dof, + const typename FunctionMap::type& function_map, + const std::vector& component_mask = std::vector()); + + /** + * Reset the data structures. + */ + void clear(); + + /** + * Determine whether a dof index + * is subject to a boundary + * constraint. + */ + bool at_boundary(unsigned int level, unsigned int index) const; + + /** + * Determine whether a dof index + * is at the refinement edge. + */ + bool at_refinement_edge(unsigned int level, unsigned int index) const; + + private: + /** + * The indices of boundary dofs + * for each level. + */ + std::vector > boundary_indices; + /** + * The degrees of freedom on the + * refinement edge between a + * level and coarser cells. + */ + std::vector > refinement_edge_indices; + + /** + * The degrees of freedom on the + * refinement edge between a + * level and coarser cells, which + * are also on the boundary. + * + * This is a subset of + * #refinement_edge_indices. + */ + std::vector > refinement_edge_boundary_indices; +}; + + +bool +MGConstraints::at_boundary(unsigned int level, unsigned int index) +{ + AssertIndexRange(level, boundary_indices.size()); + AssertIndexRange(level, boundary_indices.size()); +} + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/deal.II/include/multigrid/mg_dof_accessor.h b/deal.II/deal.II/include/multigrid/mg_dof_accessor.h index 948ffc2d28..01a42d0fa0 100644 --- a/deal.II/deal.II/include/multigrid/mg_dof_accessor.h +++ b/deal.II/deal.II/include/multigrid/mg_dof_accessor.h @@ -289,8 +289,11 @@ class MGDoFAccessor : public internal::MGDoFAccessor::Inheritance local_dofs; + + // First, deal with the simpler + // case when we have to identify + // all boundary dofs + if (component_mask_.size() == 0) + { + typename MGDoFHandler::cell_iterator + cell = dof.begin(), + endc = dof.end(); + for (; cell!=endc; ++cell) + { + const FiniteElement &fe = cell->get_fe(); + const unsigned int level = cell->level(); + local_dofs.resize(fe.dofs_per_face); + + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; + ++face_no) + { + const typename MGDoFHandler::face_iterator + face = cell->face(face_no); + face->get_mg_dof_indices(level, local_dofs); + const unsigned char bi = face->boundary_indicator(); + // Face is listed in + // boundary map + if (function_map.find(bi) != function_map.end()) + { + for (unsigned int i=0;i::type::const_iterator i=function_map.begin(); // i!=function_map.end(); ++i) // Assert (n_components == i->second->n_components, // ExcInvalidFE()); - - // set the component mask to either - // the original value or a vector - // of @p{true}s - const std::vector component_mask ((component_mask_.size() == 0) ? - std::vector (n_components, true) : - component_mask_); + + // set the component mask to either + // the original value or a vector + // of @p{true}s + const std::vector component_mask ((component_mask_.size() == 0) ? + std::vector (n_components, true) : + component_mask_); // Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0, // ExcComponentMismatch()); - - // field to store the indices - std::vector face_dofs; - face_dofs.reserve (DoFTools::max_dofs_per_face(dof)); - std::fill (face_dofs.begin (), face_dofs.end (), - DoFHandler::invalid_dof_index); - - typename MGDoFHandler::cell_iterator - cell = dof.begin(), - endc = dof.end(); - for (; cell!=endc; ++cell) - for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; - ++face_no) - { - const FiniteElement &fe = cell->get_fe(); - const unsigned int level = cell->level(); - - // we can presently deal only with - // primitive elements for boundary - // values. this does not preclude - // us using non-primitive elements - // in components that we aren't - // interested in, however. make - // sure that all shape functions - // that are non-zero for the - // components we are interested in, - // are in fact primitive - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) - { - const std::vector &nonzero_component_array - = cell->get_fe().get_nonzero_components (i); - for (unsigned int c=0; cget_fe().is_primitive (i), - ExcMessage ("This function can only deal with requested boundary " - "values that correspond to primitive (scalar) base " - "elements")); - } - - typename MGDoFHandler::face_iterator face = cell->face(face_no); - const unsigned char boundary_component = face->boundary_indicator(); - if (function_map.find(boundary_component) != function_map.end()) - // face is of the right component - { - // get indices, physical location and - // boundary values of dofs on this - // face - face_dofs.resize (fe.dofs_per_face); - face->get_mg_dof_indices (level, face_dofs); - if (fe_is_system) - { - // enter those dofs - // into the list that - // match the - // component - // signature. avoid - // the usual - // complication that - // we can't just use - // *_system_to_component_index - // for non-primitive - // FEs - for (unsigned int i=0; i local_dofs; + local_dofs.reserve (DoFTools::max_dofs_per_face(dof)); + std::fill (local_dofs.begin (), local_dofs.end (), + DoFHandler::invalid_dof_index); + + typename MGDoFHandler::cell_iterator + cell = dof.begin(), + endc = dof.end(); + for (; cell!=endc; ++cell) + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; + ++face_no) + { + const FiniteElement &fe = cell->get_fe(); + const unsigned int level = cell->level(); + + // we can presently deal only with + // primitive elements for boundary + // values. this does not preclude + // us using non-primitive elements + // in components that we aren't + // interested in, however. make + // sure that all shape functions + // that are non-zero for the + // components we are interested in, + // are in fact primitive + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + { + const std::vector &nonzero_component_array + = cell->get_fe().get_nonzero_components (i); + for (unsigned int c=0; cget_fe().is_primitive (i), + ExcMessage ("This function can only deal with requested boundary " + "values that correspond to primitive (scalar) base " + "elements")); + } + + typename MGDoFHandler::face_iterator face = cell->face(face_no); + const unsigned char boundary_component = face->boundary_indicator(); + if (function_map.find(boundary_component) != function_map.end()) + // face is of the right component + { + // get indices, physical location and + // boundary values of dofs on this + // face + local_dofs.resize (fe.dofs_per_face); + face->get_mg_dof_indices (level, local_dofs); + if (fe_is_system) + { + // enter those dofs + // into the list that + // match the + // component + // signature. avoid + // the usual + // complication that + // we can't just use + // *_system_to_component_index + // for non-primitive + // FEs + for (unsigned int i=0; i