From 2f5d03576d224bc04da4718cc730cd2a9e3ffb34 Mon Sep 17 00:00:00 2001
From: guido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 17 Dec 2001 22:44:40 +0000
Subject: [PATCH] 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
---
 deal.II/deal.II/source/dofs/dof_tools.cc | 149 ++++++++++++++++++-----
 deal.II/deal.II/source/fe/fe_dgq.cc      | 106 +++++++---------
 2 files changed, 162 insertions(+), 93 deletions(-)

diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc
index 074c521712..5477f0f375 100644
--- 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 f27602159a..12970c699f 100644
--- 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;
 }
-- 
2.39.5