make_flux_sparsity_pattern more sophisticated; need check in 3D and cleanup

git-svn-id: https://svn.dealii.org/trunk@5339 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc
index 074c521712a84fe4d8bb5c9d212baeae9310308f..5477f0f3756bd96688c56d4171942588d823c6c6 100644 (file)
--- a/deal.II/deal.II/source/dofs/dof_tools.cc
+++ b/deal.II/deal.II/source/dofs/dof_tools.cc
@@ -444,7 +444,8 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
                                      const FullMatrix<double>& flux_mask)
 {
   const unsigned int n_dofs = dof.n_dofs();
-  const unsigned int n_comp = dof.get_fe().n_components();
+  const FiniteElement<dim>& fe = dof.get_fe();
+  const unsigned int n_comp = fe.n_components();
   
   Assert (sparsity.n_rows() == n_dofs,
          ExcDimensionMismatch (sparsity.n_rows(), n_dofs));
@@ -459,28 +460,35 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
   Assert (flux_mask.n() == n_comp,
          ExcDimensionMismatch (flux_mask.n(), n_comp));
   
-  const unsigned int total_dofs = dof.get_fe().dofs_per_cell;
+  const unsigned int total_dofs = fe.dofs_per_cell;
   std::vector<unsigned int> dofs_on_this_cell(total_dofs);
   std::vector<unsigned int> dofs_on_other_cell(total_dofs);
+  vector2d<bool> support_on_face(total_dofs, GeometryInfo<dim>::faces_per_cell);
+  
   typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
                                                 endc = dof.end();
 
 
-  std::vector<std::vector<bool> > int_dof_mask(total_dofs,
-                                              std::vector<bool>(total_dofs, false));
-  std::vector<std::vector<bool> > flux_dof_mask(total_dofs,
-                                               std::vector<bool>(total_dofs, false));
+  vector2d<char> int_dof_mask(total_dofs, total_dofs);
+  vector2d<char> flux_dof_mask(total_dofs, total_dofs);
+
   for (unsigned int i=0; i<total_dofs; ++i)
-    for (unsigned int j=0; j<total_dofs; ++j)
-      {
-       unsigned int ii = dof.get_fe().system_to_component_index(i).first;
-       unsigned int jj = dof.get_fe().system_to_component_index(j).first;
-       
-       if (int_mask(ii,jj) != 0)
-         int_dof_mask[i][j] = true;
-       if (flux_mask(ii,jj) != 0)
-         flux_dof_mask[i][j] = true;
-      }
+    {
+      for (unsigned int j=0; j<total_dofs; ++j)
+       {
+         unsigned int ii = fe.system_to_component_index(i).first;
+         unsigned int jj = fe.system_to_component_index(j).first;
+         
+         if (int_mask(ii,jj) != 0)
+           int_dof_mask(i,j) = 'f';
+         if (flux_mask(ii,jj) == 1.)
+           flux_dof_mask(i,j) = 'a';
+         if (flux_mask(ii,jj) == 2.)
+           flux_dof_mask(i,j) = 'f';
+       }
+      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell;++f)
+       support_on_face(i,f) = fe.has_support_on_face(i,f);
+    }
   
   (const_cast<Triangulation<dim>& > (dof.get_tria())).clear_user_flags();
   
@@ -490,7 +498,7 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
                                       // make sparsity pattern for this cell
       for (unsigned int i=0; i<total_dofs; ++i)
        for (unsigned int j=0; j<total_dofs; ++j)
-         if (int_dof_mask[i][j])
+         if (int_dof_mask(i,j))
            sparsity.add (dofs_on_this_cell[i],
                          dofs_on_this_cell[j]);
 
@@ -503,7 +511,30 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
          if (cell_face->user_flag_set ())
            continue;
 
-         if (! cell->at_boundary (face) )
+         if (cell->at_boundary (face) )
+           {
+             for (unsigned int i=0; i<total_dofs; ++i)
+               {
+                 const bool i_non_zero_i = support_on_face (i, face);
+                 for (unsigned int j=0; j<total_dofs; ++j)
+                   {
+                     const bool j_non_zero_i = support_on_face (j, face);
+                     
+                     if (flux_dof_mask(i,j) == 'f')
+                       {
+                         sparsity.add (dofs_on_this_cell[i],
+                                       dofs_on_this_cell[j]);
+                       }
+                     if (flux_dof_mask(i,j) == 'a')
+                       {
+                         if (i_non_zero_i && j_non_zero_i)
+                           sparsity.add (dofs_on_this_cell[i],
+                                         dofs_on_this_cell[j]);
+                       }
+                   }
+               }
+           }
+         else
            {
              typename DoFHandler<dim>::cell_iterator
                neighbor = cell->neighbor(face);
@@ -527,14 +558,22 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
                      sub_neighbor->get_dof_indices (dofs_on_other_cell);
                      for (unsigned int i=0; i<total_dofs; ++i)
                        {
+                         const bool i_non_zero_i = support_on_face (i, face);
+                         const bool i_non_zero_e = support_on_face (i, neighbor_face);
                          for (unsigned int j=0; j<total_dofs; ++j)
                            {
-                             if (flux_dof_mask[i][j])
+                             const bool j_non_zero_i = support_on_face (j, face);
+                             const bool j_non_zero_e  =support_on_face (j, neighbor_face);
+                             if (flux_dof_mask(i,j) == 'f'
+                                 || (flux_dof_mask(i,j) == 'a'
+                                      && i_non_zero_i && j_non_zero_e))
                                sparsity.add (dofs_on_this_cell[i],
                                              dofs_on_other_cell[j]);
-                             if (flux_dof_mask[j][i])
-                               sparsity.add (dofs_on_other_cell[j],
-                                             dofs_on_this_cell[i]);
+                             if (flux_dof_mask(j,i) == 'f'
+                                 || (flux_dof_mask(j,i) == 'a'
+                                     && i_non_zero_e && j_non_zero_i))
+                               sparsity.add (dofs_on_other_cell[i],
+                                             dofs_on_this_cell[j]);
                            }
                        }
                      sub_neighbor->face(neighbor_face)->set_user_flag ();
@@ -543,14 +582,68 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
                  neighbor->get_dof_indices (dofs_on_other_cell);
                  for (unsigned int i=0; i<total_dofs; ++i)
                    {
+                     const bool i_non_zero_i = support_on_face (i, face);
+                     const bool i_non_zero_e = support_on_face (i, neighbor_face);
                      for (unsigned int j=0; j<total_dofs; ++j)
                        {
-                         if (flux_dof_mask[i][j])
-                           sparsity.add (dofs_on_this_cell[i],
-                                         dofs_on_other_cell[j]);
-                         if (flux_dof_mask[j][i])
-                           sparsity.add (dofs_on_other_cell[j],
-                                         dofs_on_this_cell[i]);
+                         const bool j_non_zero_i = support_on_face (j, face);
+                         const bool j_non_zero_e = support_on_face (j, neighbor_face);
+                         if (flux_dof_mask(i,j)==0)
+                           flux_dof_mask(i,j) = 'n';
+                         
+                         if (flux_dof_mask(i,j) == 'f')
+                           {
+                             sparsity.add (dofs_on_this_cell[i],
+                                           dofs_on_other_cell[j]);
+                             sparsity.add (dofs_on_other_cell[i],
+                                           dofs_on_this_cell[j]);
+                             sparsity.add (dofs_on_this_cell[i],
+                                           dofs_on_this_cell[j]);
+                             sparsity.add (dofs_on_other_cell[i],
+                                           dofs_on_other_cell[j]);
+                           }
+                         if (flux_dof_mask(i,j) == 'a')
+                           {
+                             if (i_non_zero_i && j_non_zero_e)
+                               sparsity.add (dofs_on_this_cell[i],
+                                             dofs_on_other_cell[j]);
+                             if (i_non_zero_e && j_non_zero_i)
+                               sparsity.add (dofs_on_other_cell[i],
+                                             dofs_on_this_cell[j]);
+                             if (i_non_zero_i && j_non_zero_i)
+                               sparsity.add (dofs_on_this_cell[i],
+                                             dofs_on_this_cell[j]);
+                             if (i_non_zero_e && j_non_zero_e)
+                               sparsity.add (dofs_on_other_cell[i],
+                                             dofs_on_other_cell[j]); 
+                           }
+
+                         if (flux_dof_mask(j,i) == 'f')
+                           {
+                             sparsity.add (dofs_on_this_cell[j],
+                                           dofs_on_other_cell[i]);
+                             sparsity.add (dofs_on_other_cell[j],
+                                           dofs_on_this_cell[i]);
+                             sparsity.add (dofs_on_this_cell[j],
+                                           dofs_on_this_cell[i]);
+                             sparsity.add (dofs_on_other_cell[j],
+                                           dofs_on_other_cell[i]);
+                           }
+                         if (flux_dof_mask(j,i) == 'a')
+                           {
+                             if (j_non_zero_i && i_non_zero_e)
+                               sparsity.add (dofs_on_this_cell[j],
+                                             dofs_on_other_cell[i]);
+                             if (j_non_zero_e && i_non_zero_i)
+                               sparsity.add (dofs_on_other_cell[j],
+                                             dofs_on_this_cell[i]);
+                             if (j_non_zero_i && i_non_zero_i)
+                               sparsity.add (dofs_on_this_cell[j],
+                                             dofs_on_this_cell[i]);
+                             if (j_non_zero_e && i_non_zero_e)
+                               sparsity.add (dofs_on_other_cell[j],
+                                             dofs_on_other_cell[i]); 
+                           }
                        }
                    }
                  neighbor->face(neighbor_face)->set_user_flag (); 

diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc
index f27602159a2da41abc2a87e17bcfff5c1ec976ab..12970c699f8f0ef46cf15bafa4792c8aa1b4eec0 100644 (file)
--- a/deal.II/deal.II/source/fe/fe_dgq.cc
+++ b/deal.II/deal.II/source/fe/fe_dgq.cc
@@ -588,75 +588,51 @@ bool
 FE_DGQ<dim>::has_support_on_face (unsigned int shape_index,
                                  const unsigned int face_index) const
 {
-  if (dim==1)
-    return true;
-  const unsigned int cell_start = (dim==2)
-                                 ? first_quad_index
-                                 : first_hex_index;
-  const unsigned int face_start = (dim==2)
-                                 ? first_line_index
-                                 : first_quad_index;
+
+  unsigned int n = degree+1;
+  unsigned int n2 = n*n;
   
-                                  // Interior degrees of
-                                  // freedom correspond to
-                                  // shape functions with
-                                  // support inside the cell.
-  if (shape_index >= cell_start)
-    return false;
-                                  // Shape functions are sorted
-                                  // by face. If we dived by
-                                  // the number of shapes per
-                                  // face, the result must be
-                                  // equal to the face index.
-  if (shape_index >= face_start)
-    {
-      shape_index -= first_line_index;
-      shape_index /=  dofs_per_face;
-      return (shape_index == face_index);
-    }
-                                  // Only degrees of freedom on
-                                  // a vertex are left.
-  shape_index /=  dofs_per_vertex;
-                                  // Use a table to figure out
-                                  // which face is neighbor to
-                                  // which vertex.
-  switch (100*dim+10*face_index+shape_index)
+  switch (dim)
     {
-      case 200:
-      case 230:
-      case 201:
-      case 211:
-      case 212:
-      case 222:
-      case 223:
-      case 233:
-      case 300:
-      case 320:
-      case 350:
-      case 301:
-      case 321:
-      case 331:
-      case 302:
-      case 332:
-      case 342:
-      case 303:
-      case 343:
-      case 353:
-      case 314:
-      case 324:
-      case 354:
-      case 315:
-      case 325:
-      case 335:
-      case 316:
-      case 336:
-      case 346:
-      case 317:
-      case 347:
-      case 357:
+      case 1:
+                                        // This is not correct, but it
+                                        // should not matter in 1D.
        return true;
-      default:
+      case 2:
+       if (face_index==0 && shape_index < n)
+         return true;
+       if (face_index==1 && (shape_index % n) == degree)
+         return true;
+       if (face_index==2 && shape_index >= dofs_per_cell-n)
+         return true;
+       if (face_index==3 && (shape_index % n) == 0)
+         return true;
        return false;
+      case 3:
+       if (true)
+         {
+           const unsigned int in2 = shape_index % n2;
+           
+                                            // y=0
+           if (face_index==0 && in2 < n )
+             return true;
+                                            // y=1
+           if (face_index==1 && in2 >= n2-n)
+             return true;
+                                            // z=0
+           if (face_index==2 && shape_index < n2)
+             return true;
+                                            // x=1
+           if (face_index==3 && (shape_index % n) == n-1)
+             return true;
+                                            // z=1
+           if (face_index==4 && shape_index >= dofs_per_cell - n2)
+             return true;
+                                            // x=0
+           if (face_index==5 && (shape_index % n) == 0)
+             return true;
+           return false;
+         }
     }
   return true;
 }