From: Daniel Arndt Date: Tue, 19 Jul 2016 12:38:46 +0000 (+0200) Subject: Allow for non-primitive FiniteElements in MGTools::make_boundary_list X-Git-Tag: v8.5.0-rc1~855^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1212bf2cd6a52b58a28f32acc02225b7d47ef042;p=dealii.git Allow for non-primitive FiniteElements in MGTools::make_boundary_list --- diff --git a/include/deal.II/multigrid/mg_tools.h b/include/deal.II/multigrid/mg_tools.h index dae7f34d92..816ef7500c 100644 --- a/include/deal.II/multigrid/mg_tools.h +++ b/include/deal.II/multigrid/mg_tools.h @@ -197,10 +197,10 @@ namespace MGTools std::vector &boundary_indices, const ComponentMask &component_mask = ComponentMask()); - /** - * The same function as above, but return an IndexSet rather than a - * std::set on each level. - */ + /** + * The same function as above, but return an IndexSet rather than a + * std::set on each level. + */ template void make_boundary_list (const DoFHandler &mg_dof, diff --git a/source/multigrid/mg_tools.cc b/source/multigrid/mg_tools.cc index 183114aa0f..a65f2cb6e3 100644 --- a/source/multigrid/mg_tools.cc +++ b/source/multigrid/mg_tools.cc @@ -1201,21 +1201,91 @@ namespace MGTools std::vector > &boundary_indices, const ComponentMask &component_mask) { - // if for whatever reason we were - // passed an empty map, return - // immediately - if (function_map.size() == 0) - return; + Assert (boundary_indices.size() == dof.get_triangulation().n_global_levels(), + ExcDimensionMismatch (boundary_indices.size(), + dof.get_triangulation().n_global_levels())); - const unsigned int n_levels = dof.get_triangulation().n_global_levels(); + std::vector boundary_indicators; + typename std::map* >::const_iterator it + = function_map.begin(); + for (; it!=function_map.end(); ++it) + boundary_indicators.push_back(it->first); - (void)n_levels; + std::vector boundary_indexset; + make_boundary_list(dof, boundary_indicators, boundary_indexset, component_mask); + for (unsigned int i=0; i + void + make_boundary_list(const DoFHandler &dof, + const typename FunctionMap::type &function_map, + std::vector &boundary_indices, + const ComponentMask &component_mask) + { + Assert (boundary_indices.size() == dof.get_triangulation().n_global_levels(), + ExcDimensionMismatch (boundary_indices.size(), + dof.get_triangulation().n_global_levels())); + + std::vector boundary_indicators; + typename std::map* >::const_iterator it = function_map.begin(); + for (; it!=function_map.end(); ++it) + boundary_indicators.push_back(it->first); + + make_boundary_list (dof, boundary_indicators, boundary_indices, component_mask); + } + + + template <> + void + make_boundary_list( + const DoFHandler<1,1> &, + const std::vector &, + std::vector &, + const ComponentMask &) + { + Assert(false, ExcNotImplemented()); + } + + + + template <> + void + make_boundary_list( + const DoFHandler<1,2> &, + const std::vector &, + std::vector &, + const ComponentMask &) + { + Assert(false, ExcNotImplemented()); + } - AssertDimension (boundary_indices.size(), n_levels); + + template + void + make_boundary_list(const DoFHandler &dof, + const std::vector &boundary_indicators, + std::vector &boundary_indices, + const ComponentMask &component_mask) + { + boundary_indices.resize(dof.get_triangulation().n_global_levels()); + + // if for whatever reason we were passed an empty vector, return immediately + if (boundary_indicators.size() == 0) + return; + + for (unsigned int i=0; i local_dofs; local_dofs.reserve (DoFTools::max_dofs_per_face(dof)); @@ -1223,9 +1293,7 @@ namespace MGTools local_dofs.end (), DoFHandler::invalid_dof_index); - // First, deal with the simpler - // case when we have to identify - // all boundary dofs + // First, deal with the simpler case when we have to identif all boundary dofs if (component_mask.n_selected_components(n_components) == n_components) { typename DoFHandler::cell_iterator @@ -1247,13 +1315,11 @@ namespace MGTools const typename DoFHandler::face_iterator face = cell->face(face_no); const types::boundary_id bi = face->boundary_id(); - // Face is listed in - // boundary map - if (function_map.find(bi) != function_map.end()) + // Face is listed in boundary map + if (std::find(boundary_indicators.begin(), boundary_indicators.end(), bi) != boundary_indicators.end()) { face->get_mg_dof_indices(level, local_dofs); - for (unsigned int i=0; i &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 ComponentMask &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 DoFHandler::face_iterator face = cell->face(face_no); const types::boundary_id boundary_component = face->boundary_id(); - if (function_map.find(boundary_component) != function_map.end()) - // face is of the right component + if (std::find(boundary_indicators.begin(), boundary_indicators.end(), boundary_component) + != boundary_indicators.end()) + // we want to constrain this boundary { - // get indices, physical location and - // boundary values of dofs on this - // face + + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + { + const ComponentMask &nonzero_component_array + = cell->get_fe().get_nonzero_components (i); + // if we want to constrain one of the nonzero components, + // we have to constrain all of them + + bool selected = false; + for (unsigned int c=0; cget_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 + // 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; iget_fe().get_nonzero_components (i); + for (unsigned int c=0; c - void - make_boundary_list(const DoFHandler &dof, - const typename FunctionMap::type &function_map, - std::vector &boundary_indices, - const ComponentMask &component_mask) - { - Assert (boundary_indices.size() == dof.get_triangulation().n_global_levels(), - ExcDimensionMismatch (boundary_indices.size(), - dof.get_triangulation().n_global_levels())); - - std::vector > - my_boundary_indices (dof.get_triangulation().n_global_levels()); - make_boundary_list (dof, function_map, my_boundary_indices, component_mask); - for (unsigned int i=0; i - void - make_boundary_list(const DoFHandler &dof, - const std::vector &boundary_indicators, - std::vector &boundary_indices, - const ComponentMask &component_mask) - { - Assert (boundary_indices.size() == dof.get_triangulation().n_global_levels(), - ExcDimensionMismatch (boundary_indices.size(), - dof.get_triangulation().n_global_levels())); - - std::vector > - my_boundary_indices (dof.get_triangulation().n_global_levels()); - - //TODO get rid off this - typename FunctionMap::type homogen_bc; - ZeroFunction zero_function (dim); - for (unsigned int i=0;i void extract_non_interface_dofs (const DoFHandler &mg_dof_handler,