]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
compute_row_length_vector for systems
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 24 Aug 2005 07:47:39 +0000 (07:47 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 24 Aug 2005 07:47:39 +0000 (07:47 +0000)
git-svn-id: https://svn.dealii.org/trunk@11312 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_tools.h
deal.II/deal.II/source/dofs/dof_tools.cc

index 65b8abc891926197986d7b73c1d5a788073aea54..7b5a2bc4e02e9f798aa6092edaf7004855207069 100644 (file)
@@ -225,8 +225,15 @@ class DoFTools
                                      * <b>not</b> coupling in the
                                      * differential equation can be
                                      * eliminated by the two optional
-                                     * tables.
+                                     * tables, which may reduce the
+                                     * amount of pre-allocated
+                                     * memory dramatically.
                                      *
+                                     * @todo This function is only
+                                     * implemented for primitive
+                                     * elements. Implementation in 1D
+                                     * is missing completely.
+                                     * 
                                      * @param dofs The DoFHandler
                                      * @param row_lengths The vector
                                      * containing the resulting row
@@ -254,8 +261,8 @@ class DoFTools
     void compute_row_length_vector(
       const DoFHandler<dim>& dofs,
       std::vector<unsigned int>& row_lengths,
-      Table<2,Coupling> couplings /*= typename Table<2,Coupling>()*/,
-      Table<2,Coupling> flux_couplings /*= Table<2,Coupling>()*/);
+      const Table<2,Coupling>& couplings /*= typename Table<2,Coupling>()*/,
+      const Table<2,Coupling>& flux_couplings /*= Table<2,Coupling>()*/);
     
                                     /**
                                      * Locate non-zero entries of the
@@ -330,7 +337,7 @@ class DoFTools
 
                                     /**
                                      * Locate non-zero entries for
-                                     * ector valued finite elements.
+                                     * vector valued finite elements.
                                      * This function does mostly the
                                      * same as the previous
                                      * @p make_sparsity_pattern, but
@@ -372,16 +379,16 @@ class DoFTools
                                      * contained in there.
                                      *
                                      * This function is designed to
-                                     * accept a mask, like the one
+                                     * accept a coupling pattern, like the one
                                      * shown above, through the
-                                     * @p mask parameter, which
-                                     * contains boolean values. It
+                                     * @p couplings parameter, which
+                                     * contains values of type #Coupling. It
                                      * builds the matrix structure
                                      * just like the previous
                                      * function, but does not create
                                      * matrix elements if not
-                                     * specified by the mask. If the
-                                     * mask is symmetric, then so
+                                     * specified by the coupling pattern. If the
+                                     * couplings are symmetric, then so
                                      * will be the resulting sparsity
                                      * pattern.
                                      *
@@ -402,11 +409,24 @@ class DoFTools
                                      * more than one component (in
                                      * deal.II speak: they are
                                      * non-primitive). In this case,
-                                     * the mask element correspoding
-                                     * to the first non-zero
-                                     * component is taken, and the
-                                     * mask elements of all following
-                                     * components are ignored.
+                                     * the coupling element
+                                     * correspoding to the first
+                                     * non-zero component is taken
+                                     * and additional ones for this
+                                     * component are ignored.
+                                     */
+    template <int dim, class SparsityPattern>
+    static
+    void
+    make_sparsity_pattern (const DoFHandler<dim>&    dof,
+                          const Table<2, Coupling>& coupling,
+                          SparsityPattern&          sparsity_pattern);
+                                    /**
+                                     * @deprecated This is the old
+                                     * form of the previous
+                                     * function. It generates a table
+                                     * of #Coupling values and calls
+                                     * it.
                                      */
     template <int dim, class SparsityPattern>
     static
@@ -1554,12 +1574,39 @@ class DoFTools
 
 // ---------------------- inline and template functions --------------------
 
+template <int dim, class SparsityPattern>
+inline
+void
+DoFTools::make_sparsity_pattern (
+  const DoFHandler<dim>                 &dof,
+  const std::vector<std::vector<bool> > &mask,
+  SparsityPattern                       &sparsity_pattern)
+{
+  const unsigned int ncomp = dof.get_fe().n_components();
+  
+  Assert (mask.size() == ncomp,
+         ExcDimensionMismatch(mask.size(), ncomp));
+  for (unsigned int i=0; i<mask.size(); ++i)
+    Assert (mask[i].size() == ncomp,
+           ExcDimensionMismatch(mask[i].size(), ncomp));
+                                  // Create a coupling table out of the mask
+  Table<2, Coupling> couplings(ncomp, ncomp);
+  for (unsigned int i=0;i<ncomp;++i)
+    for (unsigned int j=0;j<ncomp;++j)
+      if (mask[i][j])
+       couplings(i,j) = always;
+                                  // Call the new function
+  make_sparsity_pattern(dof, couplings, sparsity_pattern);
+}
+
+
 template <int dim, class Comp>
+inline
 void
-DoFTools::
-map_support_points_to_dofs (const Mapping<dim>       &mapping,
-                           const DoFHandler<dim>    &dof_handler,
-                           std::map<Point<dim>, unsigned int, Comp> &point_to_index_map)
+DoFTools::map_support_points_to_dofs (
+  const Mapping<dim>       &mapping,
+  const DoFHandler<dim>    &dof_handler,
+  std::map<Point<dim>, unsigned int, Comp> &point_to_index_map)
 {
                                   // let the checking of arguments be
                                   // done by the function first
index 74a3a82a58a737ede55701b24bc5ca4c8a4bc66c..92fb934be2d026641e5f919188d94b64a60c3ab6 100644 (file)
 
 #if  deal_II_dimension == 1
 
-
+// Specialization for 1D
 template <int dim>
 void
 DoFTools::
 compute_row_length_vector(const DoFHandler<dim>     &dofs,
                           std::vector<unsigned int> &row_lengths,
-                          const Coupling)
+                          const Coupling flux_coupling)
 {
   const FiniteElement<dim>& fe = dofs.get_fe();
   Assert (row_lengths.size() == dofs.n_dofs(),
@@ -68,12 +68,87 @@ compute_row_length_vector(const DoFHandler<dim>     &dofs,
                                        // dof on this cell
       for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
        row_lengths[cell_indices[i]] += fe.dofs_per_cell;
+      
+                                      // If fluxes couple, add
+                                      // coupling to neighbor cells
+      if (flux_coupling != none)
+       for (unsigned int face=0;face<GeometryInfo<dim>::faces_per_cell;++face)
+         {
+           if (cell->at_boundary(face)) continue;
+           for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+             row_lengths[cell_indices[i]] += fe.dofs_per_cell;
+       }
+    }
+}
+
+
+// Specialization for 1D
+//TODO:[GK] Implement this.
+template <int dim>
+void
+DoFTools::compute_row_length_vector(
+  const DoFHandler<dim>& dofs,
+  std::vector<unsigned int>& row_lengths,
+  const Table<2,Coupling>& couplings,
+  const Table<2,Coupling>& flux_couplings)
+{
+  const FiniteElement<dim>& fe = dofs.get_fe();
+  Assert (fe.is_primitive(), FiniteElementBase<dim>::ExcFENotPrimitive());
+  Assert (row_lengths.size() == dofs.n_dofs(),
+         ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs()));
+  Assert (couplings.n_rows()==fe.n_components(),
+         ExcDimensionMismatch(couplings.n_rows(), fe.n_components()));
+  Assert (couplings.n_cols()==fe.n_components(),
+         ExcDimensionMismatch(couplings.n_cols(), fe.n_components()));
+  Assert (flux_couplings.n_rows()==fe.n_components(),
+         ExcDimensionMismatch(flux_couplings.n_rows(), fe.n_components()));
+  Assert (flux_couplings.n_cols()==fe.n_components(),
+         ExcDimensionMismatch(flux_couplings.n_cols(), fe.n_components()));
+  
+                                  // Function starts here by
+                                  // resetting the counters.
+  std::fill(row_lengths.begin(), row_lengths.end(), 0);
+  
+                                  // Two vectors containing the
+                                  // number of dofs of a cell, split
+                                  // by components.
+  std::vector<unsigned int> dofs_cell(fe.n_components());
+  
+  for (unsigned int comp=0;comp<fe.n_components();++comp)
+      dofs_cell[comp] = fe.base_element(fe.component_to_base(comp).first).dofs_per_cell;
+  
+  const typename DoFHandler<dim>::cell_iterator end = dofs.end();
+  typename DoFHandler<dim>::active_cell_iterator cell;
+  std::vector<unsigned int> cell_indices(fe.dofs_per_cell);
+  
+  for (cell = dofs.begin_active(); cell != end; ++cell)
+    {
+      cell->get_dof_indices(cell_indices);
+      
+                                       // each dof can couple with each other
+                                       // dof on this cell
+      for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+       for (unsigned int comp=0;comp<fe.n_components();++comp)
+         if (couplings(fe.system_to_component_index(i).first,comp) != none)
+           row_lengths[cell_indices[i]] += dofs_cell[comp];
+      
+                                      // If fluxes couple, add
+                                      // coupling to neighbor cells
+      for (unsigned int face=0;face<GeometryInfo<dim>::faces_per_cell;++face)
+       {
+         if (cell->at_boundary(face)) continue;
+         for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+           for (unsigned int comp=0;comp<fe.n_components();++comp)
+             if (flux_couplings(fe.system_to_component_index(i).first,comp) != none)
+               row_lengths[cell_indices[i]] += dofs_cell[comp];
+       }
     }
 }
 
 
 #else
 
+// Template for 2D and 3D. For 1D see specialization above
 template <int dim>
 void
 DoFTools::compute_row_length_vector(
@@ -243,21 +318,236 @@ DoFTools::compute_row_length_vector(
   user_flags_triangulation.load_user_flags(old_flags);
 }
 
-#endif
-
 
-//TODO:[GK] This function is not finished yet!!!
+//TODO:[GK] Extend this to non-primitive elements
+// This is the template for 2D and 3D. See version for 1D above
 template <int dim>
 void
 DoFTools::compute_row_length_vector(
   const DoFHandler<dim>& dofs,
   std::vector<unsigned int>& row_lengths,
-  Table<2,Coupling> couplings,
-  Table<2,Coupling> flux_couplings)
+  const Table<2,Coupling>& couplings,
+  const Table<2,Coupling>& flux_couplings)
 {
-  Assert(false, ExcNotImplemented());
+  const FiniteElement<dim>& fe = dofs.get_fe();
+  Assert (fe.is_primitive(), FiniteElementBase<dim>::ExcFENotPrimitive());
+  Assert (row_lengths.size() == dofs.n_dofs(),
+         ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs()));
+  Assert (couplings.n_rows()==fe.n_components(),
+         ExcDimensionMismatch(couplings.n_rows(), fe.n_components()));
+  Assert (couplings.n_cols()==fe.n_components(),
+         ExcDimensionMismatch(couplings.n_cols(), fe.n_components()));
+  Assert (flux_couplings.n_rows()==fe.n_components(),
+         ExcDimensionMismatch(flux_couplings.n_rows(), fe.n_components()));
+  Assert (flux_couplings.n_cols()==fe.n_components(),
+         ExcDimensionMismatch(flux_couplings.n_cols(), fe.n_components()));
+  
+                                  // Function starts here by
+                                  // resetting the counters.
+  std::fill(row_lengths.begin(), row_lengths.end(), 0);
+                                  // We need the user flags, so we
+                                  // save them for later restoration
+  std::vector<bool> old_flags;
+                                  // We need a non-constant
+                                  // triangulation for the user
+                                  // flags. Since we restore them in
+                                  // the end, this cast is safe.
+  Triangulation<dim>& user_flags_triangulation =
+    const_cast<Triangulation<dim>&> (dofs.get_tria());
+  user_flags_triangulation.save_user_flags(old_flags);
+  user_flags_triangulation.clear_user_flags();
+  
+  const typename DoFHandler<dim>::cell_iterator end = dofs.end();
+  typename DoFHandler<dim>::active_cell_iterator cell;
+  std::vector<unsigned int> cell_indices(fe.dofs_per_cell);
+  std::vector<unsigned int> neighbor_indices(fe.dofs_per_cell);
+
+                                  // Two vectors containing the
+                                  // number of dofs of a cell and of
+                                  // a face, split by components.
+  std::vector<unsigned int> dofs_cell(fe.n_components());
+  std::vector<unsigned int> dofs_face(fe.n_components());
+
+  for (unsigned int comp=0;comp<fe.n_components();++comp)
+    {
+      dofs_cell[comp] = fe.base_element(fe.component_to_base(comp).first).dofs_per_cell;
+      dofs_face[comp] = fe.base_element(fe.component_to_base(comp).first).dofs_per_face;
+    }
+  
+                                  // We loop over cells and go from
+                                  // cells to lower dimensional
+                                  // objects. This is the only way to
+                                  // cope withthe fact, that an
+                                  // unknown number of cells may
+                                  // share an object of dimension
+                                  // smaller than dim-1.
+  for (cell = dofs.begin_active(); cell != end; ++cell)
+    {
+      cell->get_dof_indices(cell_indices);
+      unsigned int i = 0;
+                                      // First, dofs on
+                                      // vertices. We assume that
+                                      // each vertex dof couples
+                                      // with all dofs on
+                                      // adjacent grid cells.
+
+                                      // Adding all dofs of the cells
+                                      // will add dofs of the faces
+                                      // of the cell adjacent to the
+                                      // vertex twice. Therefore, we
+                                      // subtract these here and add
+                                      // them in a loop over the
+                                      // faces below.
+
+                                      // in 1d, faces and vertices
+                                      // are identical. Nevertheless,
+                                      // this will only work if
+                                      // dofs_per_face is zero and
+                                      // dofs_per_vertex is
+                                      // arbitrary, not the other way
+                                      // round.
+      unsigned int increment;
+      while (i < fe.first_line_index)
+       {
+         for (unsigned int comp=0;comp<fe.n_components();++comp)
+           if (couplings(fe.system_to_component_index(i).first,comp) != none)
+             {
+               increment = dofs_cell[comp] - dim * dofs_face[comp];
+               row_lengths[cell_indices[i]] += increment;
+             }
+         ++i;
+       }
+                                      // From now on, if an object is
+                                      // a cell, its dofs only couple
+                                      // inside the cell. Since the
+                                      // faces are handled below, we
+                                      // have to subtract ALL faces
+                                      // in this case.
+      
+                                      // In all other cases we
+                                      // subtract adjacent faces to be
+                                      // added in the loop below.
+      while (i < fe.first_quad_index)
+       {
+         for (unsigned int comp=0;comp<fe.n_components();++comp)
+           if (couplings(fe.system_to_component_index(i).first,comp) != none)
+             {
+               increment = (dim>1)
+                           ? dofs_cell[comp] - (dim-1) * dofs_face[comp]
+                           : dofs_cell[comp] - GeometryInfo<dim>::faces_per_cell * dofs_face[comp];
+               row_lengths[cell_indices[i]] += increment;
+             }
+         ++i;
+       }
+      
+                                      // Now quads in 2D and 3D
+      while (i < fe.first_hex_index)
+       {
+         for (unsigned int comp=0;comp<fe.n_components();++comp)
+           if (couplings(fe.system_to_component_index(i).first,comp) != none)
+             {
+               increment = (dim>2)
+                           ? dofs_cell[comp] - (dim-2) * dofs_face[comp]
+                           : dofs_cell[comp] - GeometryInfo<dim>::faces_per_cell * dofs_face[comp];
+               row_lengths[cell_indices[i]] += increment;
+             }
+         ++i;
+       }
+      
+                                      // Finally, cells in 3D
+      while (i < fe.dofs_per_cell)
+       {
+         for (unsigned int comp=0;comp<fe.n_components();++comp)
+           if (couplings(fe.system_to_component_index(i).first,comp) != none)
+             {
+               increment = dofs_cell[comp] - GeometryInfo<dim>::faces_per_cell * dofs_face[comp];
+               row_lengths[cell_indices[i]] += increment;
+             }
+         ++i;
+       }
+      
+                                  // At this point, we have
+                                  // counted all dofs
+                                  // contributiong from cells
+                                  // coupled topologically to the
+                                  // adjacent cells, but we
+                                  // subtracted some faces.
+  
+                                  // Now, let's go by the faces
+                                  // and add the missing
+                                  // contribution as well as the
+                                  // flux contributions.
+      for (unsigned int iface=0;iface<GeometryInfo<dim>::faces_per_cell;++iface)
+       {
+         if (cell->at_boundary(iface))
+           {
+             for (unsigned int i=0;i<fe.dofs_per_cell;++i)
+               row_lengths[cell_indices[i]] += fe.dofs_per_face;
+             continue;
+           }
+         
+         const typename DoFHandler<dim>::cell_iterator neighbor = cell->neighbor(iface);
+//       const bool neighbor_same_level = (neighbor->level() == cell->level());
+//       const unsigned int
+//         nface = neighbor_same_level
+//         ? cell->neighbor_of_neighbor(iface)
+//         : cell->neighbor_of_coarser_neighbor(iface).first;
+         typename DoFHandler<dim>::face_iterator face = cell->face(iface);
+         
+                                          // Flux couplings are
+                                          // computed from both sides
+                                          // for simplicity.
+         
+                                          // The dofs on the common face
+                                          // will be handled below,
+                                          // therefore, we subtract them
+                                          // here.
+         for (unsigned int comp=0;comp<fe.n_components();++comp)
+           for (unsigned int i=0;i<fe.dofs_per_cell;++i)
+             if (flux_couplings(fe.system_to_component_index(i).first,comp) != none)
+               {
+                 unsigned int increment = dofs_cell[comp] - dofs_face[comp];
+                 row_lengths[cell_indices[i]] += increment;
+               }
+         
+                                          // Do this only once per
+                                          // face and not on the
+                                          // hanging faces.
+         if (face->user_flag_set() || neighbor->has_children())
+           continue;
+         face->set_user_flag();
+                                          // At this point, we assume
+                                          // that each cell added its
+                                          // dofs minus the face to
+                                          // the couplings of the
+                                          // face dofs. Since we
+                                          // subtracted two faces, we
+                                          // have to re-add one.
+         
+                                          // If one side of the face
+                                          // is refined, all the fine
+                                          // face dofs couple with
+                                          // the coarse one.
+         
+                                          // Wolfgang, do they couple
+                                          // with each other by
+                                          // constraints?
+         neighbor->get_dof_indices(neighbor_indices);
+         for (unsigned int comp=0;comp<fe.n_components();++comp)
+           for (unsigned int i=0;i<fe.dofs_per_cell;++i)
+             if (couplings(fe.system_to_component_index(i).first,comp) != none)
+               {
+                 row_lengths[cell_indices[i]] += dofs_face[comp];
+                 row_lengths[neighbor_indices[i]] += dofs_face[comp];
+               }
+       }
+    }  
+  user_flags_triangulation.load_user_flags(old_flags);
 }
 
+#endif
+
+
 
 template <int dim, class SparsityPattern>
 void
@@ -292,9 +582,9 @@ DoFTools::make_sparsity_pattern (
 template <int dim, class SparsityPattern>
 void
 DoFTools::make_sparsity_pattern (
-  const DoFHandler<dim>                 &dof,
-  const std::vector<std::vector<bool> > &mask,
-  SparsityPattern                       &sparsity)
+  const DoFHandler<dim>&   dof,
+  const Table<2,Coupling>& couplings,
+  SparsityPattern&         sparsity)
 {
   const unsigned int n_dofs = dof.n_dofs();
   const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
@@ -303,11 +593,10 @@ DoFTools::make_sparsity_pattern (
          ExcDimensionMismatch (sparsity.n_rows(), n_dofs));
   Assert (sparsity.n_cols() == n_dofs,
          ExcDimensionMismatch (sparsity.n_cols(), n_dofs));
-  Assert (mask.size() == dof.get_fe().n_components(),
-         ExcDimensionMismatch(mask.size(), dof.get_fe().n_components()));
-  for (unsigned int i=0; i<mask.size(); ++i)
-    Assert (mask[i].size() == dof.get_fe().n_components(),
-           ExcDimensionMismatch(mask[i].size(), dof.get_fe().n_components()));
+  Assert (couplings.n_rows() == dof.get_fe().n_components(),
+         ExcDimensionMismatch(couplings.n_rows(), dof.get_fe().n_components()));
+  Assert (couplings.n_cols() == dof.get_fe().n_components(),
+         ExcDimensionMismatch(couplings.n_cols(), dof.get_fe().n_components()));
 
                                   // first build a mask for each dof,
                                   // not like the one given which
@@ -322,9 +611,8 @@ DoFTools::make_sparsity_pattern (
     for (unsigned int j=0; j<dofs_per_cell; ++j)
       if (dof.get_fe().is_primitive(i) &&
           dof.get_fe().is_primitive(j))
-        dof_mask[i][j] = mask
-                         [dof.get_fe().system_to_component_index(i).first]
-                         [dof.get_fe().system_to_component_index(j).first];
+        dof_mask[i][j] = (couplings(dof.get_fe().system_to_component_index(i).first,
+                                   dof.get_fe().system_to_component_index(j).first) != none);
       else
         {
           const unsigned int first_nonzero_comp_i
@@ -344,7 +632,7 @@ DoFTools::make_sparsity_pattern (
           Assert (first_nonzero_comp_j < dof.get_fe().n_components(),
                   ExcInternalError());          
           
-          dof_mask[i][j] = mask[first_nonzero_comp_i][first_nonzero_comp_j];
+          dof_mask[i][j] = (couplings(first_nonzero_comp_i,first_nonzero_comp_j) != none);
         }
   
 
@@ -3105,6 +3393,10 @@ DoFTools::compute_row_length_vector(
   const DoFHandler<deal_II_dimension>& dofs, std::vector<unsigned int>& row_lengths,
   const Coupling flux_coupling);
 
+template void
+DoFTools::compute_row_length_vector(
+  const DoFHandler<deal_II_dimension>& dofs, std::vector<unsigned int>& row_lengths,
+  const Table<2,Coupling>& couplings, const Table<2,Coupling>& flux_couplings);
 
 template void
 DoFTools::make_sparsity_pattern<deal_II_dimension,SparsityPattern>
@@ -3128,27 +3420,19 @@ DoFTools::make_sparsity_pattern<deal_II_dimension,CompressedBlockSparsityPattern
 
 template void 
 DoFTools::make_sparsity_pattern<deal_II_dimension,SparsityPattern>
-(const DoFHandler<deal_II_dimension> &dof,
- const std::vector<std::vector<bool> > &mask,
- SparsityPattern             &sparsity);
+(const DoFHandler<deal_II_dimension>&, const Table<2,Coupling>&, SparsityPattern&);
 
 template void 
 DoFTools::make_sparsity_pattern<deal_II_dimension,CompressedSparsityPattern>
-(const DoFHandler<deal_II_dimension> &dof,
- const std::vector<std::vector<bool> > &mask,
- CompressedSparsityPattern             &sparsity);
+(const DoFHandler<deal_II_dimension>&, const Table<2,Coupling>&, CompressedSparsityPattern&);
 
 template void 
 DoFTools::make_sparsity_pattern<deal_II_dimension,BlockSparsityPattern>
-(const DoFHandler<deal_II_dimension> &dof,
- const std::vector<std::vector<bool> > &mask,
- BlockSparsityPattern        &sparsity);
+(const DoFHandler<deal_II_dimension>&, const Table<2,Coupling>&, BlockSparsityPattern&);
 
 template void 
 DoFTools::make_sparsity_pattern<deal_II_dimension,CompressedBlockSparsityPattern>
-(const DoFHandler<deal_II_dimension> &dof,
- const std::vector<std::vector<bool> > &mask,
- CompressedBlockSparsityPattern        &sparsity);
+(const DoFHandler<deal_II_dimension>&, const Table<2,Coupling>&, CompressedBlockSparsityPattern&);
 
 #if deal_II_dimension > 1
 template void

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.