]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Allow for non-primitive FiniteElements in MGTools::make_boundary_list
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Tue, 19 Jul 2016 12:38:46 +0000 (14:38 +0200)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Wed, 20 Jul 2016 00:22:49 +0000 (02:22 +0200)
include/deal.II/multigrid/mg_tools.h
source/multigrid/mg_tools.cc

index dae7f34d926739e9e2100b7cb8cccb2952b45e2b..816ef7500cdd56ff53ad99db7e875bce56e567f3 100644 (file)
@@ -197,10 +197,10 @@ namespace MGTools
                       std::vector<IndexSet>                 &boundary_indices,
                       const ComponentMask               &component_mask = ComponentMask());
 
-   /**
-   * The same function as above, but return an IndexSet rather than a
-   * std::set<unsigned int> on each level.
-   */
+  /**
+  * The same function as above, but return an IndexSet rather than a
+  * std::set<unsigned int> on each level.
+  */
   template <int dim, int spacedim>
   void
   make_boundary_list (const DoFHandler<dim,spacedim>      &mg_dof,
index 183114aa0fc099028175c5c6601b0653a828334b..a65f2cb6e3f0a2aa4d9cd48c23d94bf40b9d628d 100644 (file)
@@ -1201,21 +1201,91 @@ namespace MGTools
     std::vector<std::set<types::global_dof_index> > &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<types::boundary_id> boundary_indicators;
+    typename std::map<types::boundary_id, const Function<dim>* >::const_iterator it
+      = function_map.begin();
+    for (; it!=function_map.end(); ++it)
+      boundary_indicators.push_back(it->first);
 
-    (void)n_levels;
+    std::vector<IndexSet>  boundary_indexset;
+    make_boundary_list(dof, boundary_indicators, boundary_indexset, component_mask);
+    for (unsigned int i=0; i<dof.get_triangulation().n_global_levels(); ++i)
+      {
+        IndexSet::ElementIterator it_element = boundary_indexset[i].begin();
+        for (; it_element != boundary_indexset[i].end(); ++it_element)
+          boundary_indices[i].insert(*it_element);
+      }
+  }
 
 
-    const unsigned int n_components = DoFTools::n_components(dof);
-    const bool          fe_is_system = (n_components != 1);
+  template <int dim, int spacedim>
+  void
+  make_boundary_list(const DoFHandler<dim,spacedim> &dof,
+                     const typename FunctionMap<dim>::type &function_map,
+                     std::vector<IndexSet> &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<types::boundary_id> boundary_indicators;
+    typename std::map<types::boundary_id, const Function<dim>* >::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<types::boundary_id> &,
+    std::vector<IndexSet> &,
+    const ComponentMask &)
+  {
+    Assert(false, ExcNotImplemented());
+  }
+
+
+
+  template <>
+  void
+  make_boundary_list(
+    const DoFHandler<1,2> &,
+    const std::vector<types::boundary_id> &,
+    std::vector<IndexSet> &,
+    const ComponentMask &)
+  {
+    Assert(false, ExcNotImplemented());
+  }
 
-    AssertDimension (boundary_indices.size(), n_levels);
+
+  template <int dim, int spacedim>
+  void
+  make_boundary_list(const DoFHandler<dim,spacedim> &dof,
+                     const std::vector<types::boundary_id> &boundary_indicators,
+                     std::vector<IndexSet> &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<dof.get_triangulation().n_global_levels(); ++i)
+      if (boundary_indices[i].size() == 0)
+        boundary_indices[i] = IndexSet (dof.n_dofs(i));
+
+    const unsigned int n_components = DoFTools::n_components(dof);
+    const bool         fe_is_system = (n_components != 1);
 
     std::vector<types::global_dof_index> local_dofs;
     local_dofs.reserve (DoFTools::max_dofs_per_face(dof));
@@ -1223,9 +1293,7 @@ namespace MGTools
                local_dofs.end (),
                DoFHandler<dim,spacedim>::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<dim,spacedim>::cell_iterator
@@ -1247,13 +1315,11 @@ namespace MGTools
                   const typename DoFHandler<dim,spacedim>::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.dofs_per_face; ++i)
-                        boundary_indices[level].insert(local_dofs[i]);
+                      boundary_indices[level].add_indices(local_dofs.begin(), local_dofs.end());
                     }
                 }
           }
@@ -1278,191 +1344,76 @@ namespace MGTools
                 const FiniteElement<dim> &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; i<cell->get_fe().dofs_per_cell; ++i)
-                  {
-                    const ComponentMask &nonzero_component_array
-                      = cell->get_fe().get_nonzero_components (i);
-                    for (unsigned int c=0; c<n_components; ++c)
-                      if ((nonzero_component_array[c] == true)
-                          &&
-                          (component_mask[c] == true))
-                        Assert (cell->get_fe().is_primitive (i),
-                                ExcMessage ("This function can only deal with requested boundary "
-                                            "values that correspond to primitive (scalar) base "
-                                            "elements"));
-                  }
-
                 typename DoFHandler<dim,spacedim>::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; i<cell->get_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; c<n_components; ++c)
+                          if (nonzero_component_array[c] == true
+                              && component_mask[c]== true)
+                            {
+                              selected = true;
+                              break;
+                            }
+                        if (selected)
+                          for (unsigned int c=0; c<n_components; ++c)
+                            Assert (nonzero_component_array[c] == false || component_mask[c] == true,
+                                    ExcMessage ("You are using a non-primitive FiniteElement "
+                                                "and try to constrain just some of its components!"));
+                      }
+
+                    // 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
+                        // 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.size(); ++i)
                           {
-                            unsigned int component;
+                            unsigned int component = numbers::invalid_unsigned_int;
                             if (fe.is_primitive())
                               component = fe.face_system_to_component_index(i).first;
                             else
                               {
-                                // non-primitive
-                                // case. make
-                                // sure that
-                                // this
-                                // particular
-                                // shape
-                                // function
-                                // _is_
-                                // primitive,
-                                // and get at
-                                // it's
-                                // component. use
-                                // usual
-                                // trick to
-                                // transfer
-                                // face dof
-                                // index to
-                                // cell dof
-                                // index
-                                const unsigned int cell_i
-                                  = (dim == 1 ?
-                                     i
-                                     :
-                                     (dim == 2 ?
-                                      (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
-                                      :
-                                      (dim == 3 ?
-                                       (i<4*fe.dofs_per_vertex ?
-                                        i
-                                        :
-                                        (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
-                                         i+4*fe.dofs_per_vertex
-                                         :
-                                         i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
-                                       :
-                                       numbers::invalid_unsigned_int)));
-                                Assert (cell_i < fe.dofs_per_cell, ExcInternalError());
-
-                                // make sure
-                                // that if
-                                // this is
-                                // not a
-                                // primitive
-                                // shape function,
-                                // then all
-                                // the
-                                // corresponding
-                                // components
-                                // in the
-                                // mask are
-                                // not set
-//                         if (!fe.is_primitive(cell_i))
-//                           for (unsigned int c=0; c<n_components; ++c)
-//                             if (fe.get_nonzero_components(cell_i)[c])
-//                               Assert (component_mask[c] == false,
-//                                       ExcFENotPrimitive());
-
-// let's pick the first of possibly more than one non-zero
-// components. if shape function is non-primitive, then we will ignore
-// the result in the following anyway, otherwise there's only one
-// non-zero component which we will use
-                                component = fe.get_nonzero_components(cell_i).first_selected_component();
+                                // Just pick the first of the components
+                                // We already know that either all or none
+                                // of the components are selected
+                                const ComponentMask &nonzero_component_array
+                                  = cell->get_fe().get_nonzero_components (i);
+                                for (unsigned int c=0; c<n_components; ++c)
+                                  if (nonzero_component_array[c] == true)
+                                    {
+                                      component = c;
+                                      break;
+                                    }
                               }
-
+                            Assert(component!=numbers::invalid_unsigned_int, ExcInternalError());
                             if (component_mask[component] == true)
-                              boundary_indices[level].insert(local_dofs[i]);
+                              boundary_indices[level].add_index(local_dofs[i]);
                           }
                       }
                     else
                       for (unsigned int i=0; i<local_dofs.size(); ++i)
-                        boundary_indices[level].insert(local_dofs[i]);
+                        boundary_indices[level].add_indices(local_dofs.begin(), local_dofs.end());
                   }
               }
       }
   }
 
 
-  template <int dim, int spacedim>
-  void
-  make_boundary_list(const DoFHandler<dim,spacedim> &dof,
-                     const typename FunctionMap<dim>::type &function_map,
-                     std::vector<IndexSet> &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<std::set<types::global_dof_index> >
-    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<dof.get_triangulation().n_global_levels(); ++i)
-      {
-        boundary_indices[i] = IndexSet (dof.n_dofs(i));
-        boundary_indices[i].add_indices (my_boundary_indices[i].begin(),
-                                         my_boundary_indices[i].end());
-      }
-  }
-
-
-  template <int dim, int spacedim>
-  void
-  make_boundary_list(const DoFHandler<dim,spacedim> &dof,
-                     const std::vector<types::boundary_id> &boundary_indicators,
-                     std::vector<IndexSet> &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<std::set<types::global_dof_index> >
-    my_boundary_indices (dof.get_triangulation().n_global_levels());
-
-    //TODO get rid off this
-    typename FunctionMap<dim>::type homogen_bc;
-    ZeroFunction<dim> zero_function (dim);
-    for (unsigned int i=0;i<boundary_indicators.size();++i)
-       homogen_bc[i] = &zero_function;
-
-    make_boundary_list (dof, homogen_bc, my_boundary_indices, component_mask);
-    for (unsigned int i=0; i<dof.get_triangulation().n_global_levels(); ++i)
-      {
-        if (boundary_indices[i].size()!=dof.n_dofs(i))
-          boundary_indices[i] = IndexSet (dof.n_dofs(i));
-        boundary_indices[i].add_indices (my_boundary_indices[i].begin(),
-                                         my_boundary_indices[i].end());
-      }
-  }
-
-
   template <int dim, int spacedim>
   void
   extract_non_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.