]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement extract_dofs for hp::DoFHandler. Make DoFRenumbering::component_wise actual...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 18 Dec 2009 12:52:37 +0000 (12:52 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 18 Dec 2009 12:52:37 +0000 (12:52 +0000)
git-svn-id: https://svn.dealii.org/trunk@20254 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/dofs/dof_renumbering.cc
deal.II/deal.II/source/dofs/dof_tools.cc

index 58bf0d7ba1df4a9bf781f8cd41c694845aaf35e8..7809a142c0e7dab8e1d4e3d40303f2c8cbd9b1fc 100644 (file)
@@ -538,7 +538,7 @@ namespace DoFRenumbering
 
   template <int dim>
   void
-  component_wise (hp::DoFHandler<dim>                 &dof_handler,
+  component_wise (hp::DoFHandler<dim>             &dof_handler,
                  const std::vector<unsigned int> &component_order_arg)
   {
     std::vector<unsigned int> renumbering (dof_handler.n_dofs(),
@@ -635,7 +635,7 @@ namespace DoFRenumbering
                          const ENDITERATOR& end,
                          const std::vector<unsigned int> &component_order_arg)
   {
-    const hp::FECollection<dim> fe_collection (start->get_fe ());
+    const hp::FECollection<dim> fe_collection (start->get_dof_handler().get_fe ());
 
                                     // do nothing if the FE has only
                                     // one component
@@ -645,9 +645,6 @@ namespace DoFRenumbering
        return 0;
       }
 
-    const FiniteElement<dim>& fe = fe_collection[0];
-    const unsigned int dofs_per_cell = fe.dofs_per_cell;
-
 //    Assert (new_indices.size() ==  start->dof_handler().n_dofs(),
 //       ExcDimensionMismatch(new_indices.size(),
 //                            start->dof_handler.n_dofs()));
@@ -660,19 +657,19 @@ namespace DoFRenumbering
                                     // store components in the order
                                     // found in the system.
     if (component_order.size() == 0)
-      for (unsigned int i=0; i<fe.n_components(); ++i)
+      for (unsigned int i=0; i<fe_collection.n_components(); ++i)
        component_order.push_back (i);
 
-    Assert (component_order.size() == fe.n_components(),
-           ExcDimensionMismatch(component_order.size(), fe.n_components()));
+    Assert (component_order.size() == fe_collection.n_components(),
+           ExcDimensionMismatch(component_order.size(), fe_collection.n_components()));
 
     for (unsigned int i=0; i< component_order.size(); ++i)
-      Assert(component_order[i] < fe.n_components(),
-            ExcIndexRange(component_order[i], 0, fe.n_components()));
+      Assert(component_order[i] < fe_collection.n_components(),
+            ExcIndexRange(component_order[i], 0, fe_collection.n_components()));
 
                                     // vector to hold the dof indices on
                                     // the cell we visit at a time
-    std::vector<unsigned int> local_dof_indices(dofs_per_cell);
+    std::vector<unsigned int> local_dof_indices;
 
                                     // prebuilt list to which component
                                     // a given dof on a cell
@@ -684,31 +681,29 @@ namespace DoFRenumbering
                                     // belongs. in this case, assign it
                                     // to the first vector component to
                                     // which it belongs
-    std::vector<unsigned int> component_list (dofs_per_cell);
-    for (unsigned int i=0; i<component_list.size(); ++i)
-      if (fe.is_primitive(i))
-       component_list[i]
-         = component_order[fe.system_to_component_index(i).first];
-      else
-       {
-         const unsigned int base_element =
-           fe.system_to_base_index(i).first.first;
-         const unsigned int base_multiplicity =
-           fe.system_to_base_index(i).first.second;
-                                          // sum up the number of
-                                          // components all the
-                                          // elements before that have
-         unsigned int c=0;
-         for (unsigned int b=0; b<base_element; ++b)
-           c += fe.base_element(b).n_components() *
-                fe.element_multiplicity(b);
-         for (unsigned int m=0; m<base_multiplicity; ++m)
-           c += fe.base_element(base_element).n_components();
+    std::vector<std::vector<unsigned int> > component_list (fe_collection.size());
+    for (unsigned int f=0; f<fe_collection.size(); ++f)
+      {
+       const FiniteElement<dim> & fe = fe_collection[f];
+       const unsigned int dofs_per_cell = fe.dofs_per_cell;
+       component_list[f].resize(dofs_per_cell);
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         if (fe.is_primitive(i))
+           component_list[f][i]
+             = component_order[fe.system_to_component_index(i).first];
+         else
+           {
+             const unsigned int comp = (std::find(fe.get_nonzero_components(i).begin(),
+                                                  fe.get_nonzero_components(i).end(),
+                                                  true) -
+                                        fe.get_nonzero_components(i).begin());
+
                                           // then associate this degree
                                           // of freedom with this
                                           // component
-         component_list[i] = component_order[c];
-       };
+             component_list[f][i] = component_order[comp];
+           };
+      }
 
                                     // set up a map where for each
                                     // component the respective degrees
@@ -722,15 +717,19 @@ namespace DoFRenumbering
                                     // multiply, so we will have to
                                     // take care of that
     std::vector<std::vector<unsigned int> >
-      component_to_dof_map (fe.n_components());
+      component_to_dof_map (fe_collection.n_components());
     for (;start!=end;++start)
       {
                                         // on each cell: get dof indices
                                         // and insert them into the global
                                         // list using their component
+       const unsigned int fe_index = start->active_fe_index();
+       const unsigned int dofs_per_cell =fe_collection[fe_index].dofs_per_cell;
+       local_dof_indices.resize (dofs_per_cell);
        start.get_dof_indices (local_dof_indices);
        for (unsigned int i=0; i<dofs_per_cell; ++i)
-         component_to_dof_map[component_list[i]].push_back (local_dof_indices[i]);
+         component_to_dof_map[component_list[fe_index][i]].
+           push_back (local_dof_indices[i]);
       };
 
                                     // now we've got all indices sorted
@@ -756,7 +755,8 @@ namespace DoFRenumbering
                                     // into the first one. The same
                                     // holds if several components were
                                     // joined into a single target.
-    for (unsigned int component=0; component<fe.n_components(); ++component)
+    for (unsigned int component=0; component<fe_collection.n_components();
+        ++component)
       {
        std::sort (component_to_dof_map[component].begin(),
                   component_to_dof_map[component].end());
@@ -770,7 +770,7 @@ namespace DoFRenumbering
                                     // components in the order the user
                                     // desired to see
     unsigned int next_free_index = 0;
-    for (unsigned int c=0; c<fe.n_components(); ++c)
+    for (unsigned int c=0; c<fe_collection.n_components(); ++c)
       {
        const unsigned int component = c;
 
index dfbac365bfa2a6d31bdf34cc7a4e01cd60905a5a..4cb761a01f0eba2cefbe848d3d19f77030af2189 100644 (file)
@@ -194,11 +194,8 @@ DoFTools::make_sparsity_pattern (const DH                &dof,
        ||
        (subdomain_id == cell->subdomain_id()))
       {
-       const unsigned int fe_index
-         = cell->active_fe_index();
-
-       const unsigned int dofs_per_cell
-         = fe_collection[fe_index].dofs_per_cell;
+       const unsigned int fe_index = cell->active_fe_index();
+       const unsigned int dofs_per_cell =fe_collection[fe_index].dofs_per_cell;
 
        dofs_on_this_cell.resize (dofs_per_cell);
        cell->get_dof_indices (dofs_on_this_cell);
@@ -2910,36 +2907,42 @@ namespace internal
                                // component can be assigned. Then, we sort
                                // them to the first selected component of the
                                // vector system.
-  template <int dim, int spacedim>
+  template <class DH>
   inline
   void
-  extract_dofs_by_component (const dealii::DoFHandler<dim,spacedim> &dof,
+  extract_dofs_by_component (const DH                   &dof,
                             const std::vector<bool>    &component_select,
                             const bool                  sort_by_blocks,
                             std::vector<unsigned char> &dofs_by_component)
   {
-    const FiniteElement<dim,spacedim> &fe = dof.get_fe();
-    Assert (fe.n_components() < 256, ExcNotImplemented());
+    const dealii::hp::FECollection<DH::dimension,DH::space_dimension> 
+      fe_collection (dof.get_fe());
+    Assert (fe_collection.n_components() < 256, ExcNotImplemented());
     Assert(dofs_by_component.size() == dof.n_dofs(),
           ExcDimensionMismatch(dofs_by_component.size(), dof.n_dofs()));
 
-
                                    // next set up a table for the
                                    // degrees of freedom on each of
                                    // the cells whether it is
                                    // something interesting or not
-    std::vector<unsigned char> local_component_association (fe.dofs_per_cell);
-    if (sort_by_blocks == true)
+    std::vector<std::vector<unsigned char> > local_component_association
+      (fe_collection.size());
+    for (unsigned int f=0; f<fe_collection.size(); ++f)
       {
-       for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-         local_component_association[i]
-           = fe.system_to_block_index(i).first;
-      }
-    else
-      for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-       if (fe.is_primitive(i))
-         local_component_association[i] = fe.system_to_component_index(i).first;
+       const FiniteElement<DH::dimension,DH::space_dimension> &fe = 
+         fe_collection[f];
+       local_component_association[f].resize(fe.dofs_per_cell);
+       if (sort_by_blocks == true)
+         {
+           for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+             local_component_association[f][i]
+               = fe.system_to_block_index(i).first;
+         }
        else
+         for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+           if (fe.is_primitive(i))
+             local_component_association[f][i] = fe.system_to_component_index(i).first;
+           else
                                         // if this shape function is
                                         // not primitive, then we have
                                         // to work harder. we have to
@@ -2951,57 +2954,50 @@ namespace internal
                                         // to do so, get the a list of
                                         // nonzero elements and see which are
                                         // actually active
-         {
-           const unsigned int first_comp = (std::find(fe.get_nonzero_components(i).begin(),
-                                                      fe.get_nonzero_components(i).end(),
-                                                      true) -
-                                            fe.get_nonzero_components(i).begin());
-           const unsigned int end_comp = (std::find(fe.get_nonzero_components(i).begin()+first_comp,
-                                                    fe.get_nonzero_components(i).end(),
-                                                    false)-
-                                          fe.get_nonzero_components(i).begin());
+             {
+               const unsigned int first_comp =
+                 (std::find(fe.get_nonzero_components(i).begin(),
+                            fe.get_nonzero_components(i).end(),
+                            true) -
+                  fe.get_nonzero_components(i).begin());
+               const unsigned int end_comp =
+                 (std::find(fe.get_nonzero_components(i).begin()+first_comp,
+                            fe.get_nonzero_components(i).end(),
+                            false)-
+                  fe.get_nonzero_components(i).begin());
 
                                           // now check whether any of
                                           // the components in between
                                           // is set
-           if (component_select.size() == 0 ||
-               (component_select[first_comp] == true ||
-                std::count(component_select.begin()+first_comp,
-                           component_select.begin()+end_comp, true) == 0))
-             local_component_association[i] = first_comp;
-           else
-             for (unsigned int c=first_comp; c<end_comp; ++c)
-               if (component_select[c] == true)
-                 {
-                   local_component_association[i] = c;
-                   break;
-                 }
-         }
+               if (component_select.size() == 0 ||
+                   (component_select[first_comp] == true ||
+                    std::count(component_select.begin()+first_comp,
+                               component_select.begin()+end_comp, true) == 0))
+                 local_component_association[f][i] = first_comp;
+               else
+                 for (unsigned int c=first_comp; c<end_comp; ++c)
+                   if (component_select[c] == true)
+                     {
+                       local_component_association[f][i] = c;
+                       break;
+                     }
+             }
+      }
 
                                    // then loop over all cells and do
                                    // the work
-    std::vector<unsigned int> indices(fe.dofs_per_cell);
-    typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator c;
-    for (c=dof.begin_active(); c!=dof.end(); ++ c)
+    std::vector<unsigned int> indices;
+    for (typename DH::active_cell_iterator c=dof.begin_active(); 
+        c!=dof.end(); ++ c)
       {
+       const unsigned int fe_index = c->active_fe_index();
+       const unsigned int dofs_per_cell = c->get_fe().dofs_per_cell;
+       indices.resize(dofs_per_cell);
        c->get_dof_indices(indices);
-       for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-         dofs_by_component[indices[i]] = local_component_association[i];
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         dofs_by_component[indices[i]] = local_component_association[fe_index][i];
       }
   }
-
-
-
-  template <int dim, int spacedim>
-  inline
-  void
-  extract_dofs_by_component (const dealii::hp::DoFHandler<dim,spacedim> &,
-                            const std::vector<bool>  &,
-                            const bool                ,
-                            std::vector<unsigned char>&)
-  {
-    Assert (false, ExcNotImplemented());
-  }
 }
 
 
@@ -3041,7 +3037,7 @@ void DoFTools::distribute_cell_to_dof_vector (
       std::vector<bool> component_mask (dof_handler.get_fe().n_components(),
                                        false);
       component_mask[component] = true;
-      internal::extract_dofs_by_component (dof_handler, component_mask, 
+      internal::extract_dofs_by_component (dof_handler, component_mask,
                                           false, component_dofs);
 
       for (unsigned int i=0; i<dof_data.size(); ++i)
@@ -3165,15 +3161,12 @@ DoFTools::extract_dofs (
 template <int dim, int spacedim>
 void
 DoFTools::extract_dofs (
-  const hp::DoFHandler<dim,spacedim>   &/*dof*/,
-  const std::vector<bool> &/*component_select*/,
-  std::vector<bool>       &/*selected_dofs*/,
-  const bool               /*count_by_blocks*/)
+  const hp::DoFHandler<dim,spacedim>   &dof,
+  const std::vector<bool> &component_select,
+  std::vector<bool>       &selected_dofs,
+  const bool               count_by_blocks)
 {
-//TODO[WB]: implement this function
-  Assert (false, ExcNotImplemented());
-/*
-  const FiniteElement<dim,spacedim> &fe = dof.get_fe();
+  const FiniteElement<dim,spacedim> &fe = dof.begin_active()->get_fe();
 
   if (count_by_blocks == true)
     {
@@ -3210,72 +3203,20 @@ DoFTools::extract_dofs (
                                   // preset all values by false
   std::fill_n (selected_dofs.begin(), dof.n_dofs(), false);
 
-                                   // next set up a table for the
-                                   // degrees of freedom on each of
-                                   // the cells whether it is
-                                   // something interesting or not
-  std::vector<bool> local_selected_dofs (fe.dofs_per_cell, false);
-  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-    if (count_by_blocks == true)
-      local_selected_dofs[i]
-        = component_select[fe.system_to_block_index(i).first];
-    else
-      if (fe.is_primitive(i))
-       local_selected_dofs[i]
-         = component_select[fe.system_to_component_index(i).first];
-      else
-                                        // if this shape function is
-                                        // not primitive, then we have
-                                        // to work harder. we have to
-                                        // find out whether _any_ of
-                                        // the vector components of
-                                        // this element is selected or
-                                        // not
-                                        //
-                                        // to do so, get the first and
-                                        // last vector components of
-                                        // the base element to which
-                                        // the local dof with index i
-                                        // belongs
-       {
-         unsigned int first_comp = 0;
-         const unsigned int this_base = fe.system_to_base_index(i).first.first;
-         const unsigned int this_multiplicity
-           = fe.system_to_base_index(i).first.second;
-
-         for (unsigned int b=0; b<this_base; ++b)
-           first_comp += fe.base_element(b).n_components() *
-                         fe.element_multiplicity(b);
-         for (unsigned int m=0; m<this_multiplicity; ++m)
-           first_comp += fe.base_element(this_base).n_components();
-         const unsigned int end_comp = first_comp +
-                                       fe.base_element(this_base).n_components();
-
-         Assert (first_comp < fe.n_components(), ExcInternalError());
-         Assert (end_comp <= fe.n_components(),  ExcInternalError());
-
-                                          // now check whether any of
-                                          // the components in between
-                                          // is set
-         for (unsigned int c=first_comp; c<end_comp; ++c)
-           if (component_select[c] == true)
-             {
-               local_selected_dofs[i] = true;
-               break;
-             }
-       }
+                               // if we count by blocks, we need to extract
+                               // the association of blocks with local dofs,
+                               // and then go through all the cells and set
+                               // the properties according to this
+                               // info. Otherwise, we let the function
+                               // extract_dofs_by_component function do the
+                               // job.
+  std::vector<unsigned char> dofs_by_component (dof.n_dofs());
+  internal::extract_dofs_by_component (dof, component_select, count_by_blocks,
+                                      dofs_by_component);
 
-                                   // then loop over all cells and do
-                                   // the work
-  std::vector<unsigned int> indices(fe.dofs_per_cell);
-  typename DoFHandler<dim,spacedim>::active_cell_iterator c;
-  for (c=dof.begin_active(); c!=dof.end(); ++ c)
-    {
-      c->get_dof_indices(indices);
-      for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-        selected_dofs[indices[i]] = local_selected_dofs[i];
-    }
-*/
+  for (unsigned int i=0; i<dof.n_dofs(); ++i)
+    if (component_select[dofs_by_component[i]] == true)
+      selected_dofs[i] = true;
 }
 
 

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.