From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 17 Feb 2003 15:59:20 +0000 (+0000)
Subject: Fix another bug with non-primitive elements.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e2edf401c7abddaf50292bea833b782d0fd123df;p=dealii-svn.git

Fix another bug with non-primitive elements.


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

diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h
index 19d93ac6fc..cdba28b0af 100644
--- a/deal.II/deal.II/include/numerics/vectors.h
+++ b/deal.II/deal.II/include/numerics/vectors.h
@@ -620,6 +620,20 @@ class VectorTools
 				      * of the finite element used by
 				      * @p{dof}.
 				      *
+				      * If the finite element used has
+				      * shape functions that are
+				      * non-zero in more than one
+				      * component (in deal.II speak:
+				      * they are non-primitive), then
+				      * these components can presently
+				      * not be used for interpolating
+				      * boundary values. Thus, the
+				      * elements in the component mask
+				      * corresponding to the
+				      * components of these
+				      * non-primitive shape functions
+				      * must be @p{false}.
+				      *
 				      * See the general doc for more
 				      * information.
 				      */
@@ -922,6 +936,10 @@ class VectorTools
 				     /**
 				      * Exception
 				      */
+    DeclException0 (ExcFENotPrimitive);
+				     /**
+				      * Exception
+				      */
     DeclException0 (ExcInvalidBoundaryIndicator);
 				     /**
 				      * Exception
diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc
index 116ba292f2..1bf6f6a19d 100644
--- a/deal.II/deal.II/source/numerics/vectors.cc
+++ b/deal.II/deal.II/source/numerics/vectors.cc
@@ -865,14 +865,120 @@ interpolate_boundary_values (const Mapping<dim>            &mapping,
 	    
 	    if (fe_is_system)
 	      {
-		function_map.find(boundary_component)->second->vector_value_list (dof_locations,
-										  dof_values_system);
+		function_map.find(boundary_component)->second
+                  ->vector_value_list (dof_locations, dof_values_system);
 		
-						 // enter into list
+						 // 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<face_dofs.size(); ++i)
-		  if (component_mask[fe.face_system_to_component_index(i).first])
-		    boundary_values[face_dofs[i]]
-		      = dof_values_system[i](fe.face_system_to_component_index(i).first);
+                  {
+                    unsigned int component;
+                    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))
+                               :
+                               static_cast<unsigned int>(-1))));
+                        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 = (std::find (fe.get_nonzero_components(cell_i).begin(),
+                                                fe.get_nonzero_components(cell_i).end(),
+                                                true)
+                                     -
+                                     fe.get_nonzero_components(cell_i).begin());
+                      }
+                    
+                    if (component_mask[component] == true)
+                      boundary_values[face_dofs[i]] = dof_values_system[i](component);
+                  } 
 	      }
 	    else
 					       // fe has only one component,
diff --git a/tests/bits/Makefile b/tests/bits/Makefile
index 41d99e11ee..f774ac799b 100644
--- a/tests/bits/Makefile
+++ b/tests/bits/Makefile
@@ -27,6 +27,7 @@ anna_2.exe              : anna_2.g.$(OBJEXT)              $(libraries)
 anna_3.exe              : anna_3.g.$(OBJEXT)              $(libraries)
 anna_4.exe              : anna_4.g.$(OBJEXT)              $(libraries)
 anna_5.exe              : anna_5.g.$(OBJEXT)              $(libraries)
+anna_6.exe              : anna_6.g.$(OBJEXT)              $(libraries)
 dof_tools_01a.exe       : dof_tools_01a.g.$(OBJEXT)       $(libraries)
 dof_tools_01b.exe       : dof_tools_01b.g.$(OBJEXT)       $(libraries)
 dof_tools_01c.exe       : dof_tools_01c.g.$(OBJEXT)       $(libraries)
@@ -52,7 +53,7 @@ denis_1.exe             : denis_1.g.$(OBJEXT)             $(libraries)
 unit_support_points.exe : unit_support_points.g.$(OBJEXT) $(libraries)
 parameter_handler_1.exe : parameter_handler_1.g.$(OBJEXT) $(libraries)
 
-tests = anna_1 anna_2 anna_3 anna_4 anna_5 \
+tests = anna_1 anna_2 anna_3 anna_4 anna_5 anna_6 \
         dof_tools_01a dof_tools_01b dof_tools_01c dof_tools_01d \
         dof_tools_02a dof_tools_02b dof_tools_02c dof_tools_02d \
 	dof_tools_03 dof_tools_04 dof_tools_05 dof_tools_06 \
diff --git a/tests/bits/anna_6.cc b/tests/bits/anna_6.cc
new file mode 100644
index 0000000000..07b2aaaa1b
--- /dev/null
+++ b/tests/bits/anna_6.cc
@@ -0,0 +1,326 @@
+//----------------------------  anna_6.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2002, 2003 by the deal.II authors and Anna Schneebeli
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  anna_6.cc  ---------------------------
+
+/*	Test code to help fixing 
+
+DoFTools::extract_boundary_dofs
+and
+	
+VectorTools::interpolate_boundary_values
+	   
+	   
+	   
+The underlying FE is a FESystem with one Nedelec
+component and one scalar component.
+	
+We would like to use
+DoFTools::extract_boundary_dofs
+with a component_mask disabling the
+Nedelec and the scalar component respectively.
+	
+Furthermore, we would like to use
+VectorTools::interpolate_boundary_values
+with a component_mask enabling only
+the scalar component.
+	
+author: Anna Schneebeli, February 2003
+  
+*/
+
+			
+			
+
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+
+#include <grid/tria.h>
+
+#include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+
+				
+#include <fe/fe_system.h>	
+#include <fe/fe_nedelec.h>
+#include <fe/fe_q.h>
+
+#include <fstream>
+#include <iostream>
+
+
+template <int dim>
+class ImposeBC
+{
+  public:
+    ImposeBC();
+    ~ImposeBC();
+    void run ();
+    
+  private:
+   
+    void get_ready ();
+    void test_extract_boundary_DoFs ();
+    void test_interpolate_BC ();
+	
+    Triangulation<dim>   triangulation;
+	
+                                     // We use a FE-System with 2 components:
+                                     // a vector-valued one for the
+                                     // field-variable u and a scalar one for
+                                     // the pressure-variable p
+    FESystem<dim>        fe;
+    DoFHandler<dim>      dof_handler;
+    
+    unsigned int 		n_u_dofs;
+    unsigned int 		n_p_dofs;
+	
+};			
+
+
+// some boundary function for the scalar component
+template <int dim>
+class BoundaryFunction : public Function<dim>
+{
+  public:
+    BoundaryFunction ();
+    
+    virtual void vector_value (const Point<dim>   &p,
+                               Vector<double> &values) const;
+	
+};
+
+
+template <int dim>
+BoundaryFunction<dim>::BoundaryFunction () : 
+                Function<dim> (dim+1) {};
+
+
+
+template <int dim>
+inline
+void BoundaryFunction<dim>::vector_value (const Point<dim>   &p,
+                                          Vector<double> &values) const
+{
+	
+  Assert (values.size() == dim+1, 
+	  ExcDimensionMismatch (values.size(), dim+1));
+				  
+  values.clear ();
+  values(dim) = 1.;
+};
+
+
+
+
+
+                                 // Construct FE with first component: Nedelec-Element,
+                                 // second component: Q1_Element
+template <int dim>
+ImposeBC<dim>::ImposeBC() :
+                fe (FE_Nedelec<dim>(1), 1, FE_Q<dim>(1), 1),
+                dof_handler (triangulation)
+			
+				 
+{};
+
+
+
+template <int dim>
+ImposeBC<dim>::~ImposeBC() 
+{
+  dof_handler.clear ();
+};
+
+
+template <int dim>
+void ImposeBC<dim>::get_ready () 
+{  
+  dof_handler.distribute_dofs (fe);
+  std::vector<unsigned int>	dofs_per_comp;
+  DoFTools::count_dofs_per_component(dof_handler, dofs_per_comp);
+	
+                                   // For an FESystem with Nedelec-elements as
+                                   // first component and bilinear elements as
+                                   // component we have:
+                                   // dofs_per_comp[0] = dofs_per_comp[1] = # Ned-DoFs
+                                   // dofs_per_comp[2] = # Q1-DoFs
+  n_u_dofs = dofs_per_comp[0];
+  n_p_dofs = dofs_per_comp[2];
+	
+};
+
+
+template <int dim>
+void ImposeBC<dim>::test_extract_boundary_DoFs () 
+{  	
+	
+  std::map<unsigned int, double> boundary_values;	
+  std::vector<bool> 	bc_component_select(dim + 1);
+  	 
+                                   // extract boundary DoFs for the Nedelec-component
+                                   // and impose zero boundary condition						
+  bc_component_select[0] = true;
+  bc_component_select[1] = true;
+  bc_component_select[2] = false;
+  			
+  std::vector<bool> ned_boundary_dofs (dof_handler.n_dofs());
+  std::set<unsigned char> boundary_indicators;
+  boundary_indicators.insert (0);
+  DoFTools::extract_boundary_dofs (dof_handler, 
+                                   bc_component_select,
+                                   ned_boundary_dofs,
+                                   boundary_indicators);
+								   
+  
+  for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
+    if (ned_boundary_dofs[i] == true)
+      boundary_values[i] = 0.;
+	  
+
+};
+
+
+				
+template <int dim>
+void ImposeBC<dim>::test_interpolate_BC () 
+{  	
+	
+  std::map<unsigned int, double> boundary_values;	
+  std::vector<bool> 	bc_component_select(dim + 1, false);
+  
+	  
+                                   // impose inhomogeneous boundary condition
+                                   // on the scalar variable
+  bc_component_select.back() = true;
+  
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                            0,
+                                            BoundaryFunction<dim>(),
+                                            boundary_values,
+                                            bc_component_select);
+						 
+						 
+	
+                                   // check 
+                                   // (the pressure is assumed to be set to 1
+                                   // on the boundary)					 
+  std::vector<bool> p_boundary_dofs (dof_handler.n_dofs());
+  std::set<unsigned char> boundary_indicators;
+  boundary_indicators.insert (0);
+  DoFTools::extract_boundary_dofs (dof_handler, 
+                                   bc_component_select,
+                                   p_boundary_dofs,
+                                   boundary_indicators);
+  for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
+    {
+                                       // error: pressure boundary DoF
+                                       // i has not been set to the
+                                       // correct value
+                                       //
+                                       // or:
+                                       //
+                                       // nedelec boundary DoF i has
+                                       // wrongly been set to some
+                                       // value
+      Assert ((p_boundary_dofs[i] && boundary_values[i] == 1.)
+              ||
+              (!(p_boundary_dofs[i])  && boundary_values[i] != 1.),
+              ExcInternalError());
+
+      deallog << boundary_values[i] << ' ';
+    }
+  deallog << std::endl;
+};
+
+
+
+
+template <int dim>
+void ImposeBC<dim>::run () 
+{ 
+  GridGenerator::hyper_cube(triangulation, -1,1);
+  triangulation.refine_global (1);
+  triangulation.begin_active()->set_refine_flag ();
+  triangulation.execute_coarsening_and_refinement ();
+  
+  deallog << "   Number of active cells:       "
+          << triangulation.n_active_cells()
+          << std::endl;
+	
+  get_ready ();
+ 
+  deallog << "  Total number of degrees of freedom: "
+          << dof_handler.n_dofs()
+          << std::endl
+          << "   Number of degrees of freedom for the field variable U: "
+          << n_u_dofs
+          << std::endl
+          << "   Number of degrees of freedom for the pressure variable p: "
+          << n_p_dofs
+          << std::endl;
+
+  test_extract_boundary_DoFs ();
+  test_interpolate_BC ();
+	  
+};
+
+
+int main () 
+{
+  try
+    {
+      std::ofstream logfile("anna_6.output");
+      logfile.precision (2);      
+      deallog.attach(logfile);
+      deallog.depth_console(0);
+
+      ImposeBC<2>().run ();
+      ImposeBC<3>().run ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+		<< exc.what() << std::endl
+		<< "Aborting!" << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+      
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+		<< "Aborting!" << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+      return 1;
+    };
+
+  return 0;
+};
diff --git a/tests/bits/anna_6.checked b/tests/bits/anna_6.checked
new file mode 100644
index 0000000000..e91fa2436c
--- /dev/null
+++ b/tests/bits/anna_6.checked
@@ -0,0 +1,11 @@
+
+DEAL::   Number of active cells:       7
+DEAL::  Total number of degrees of freedom: 36
+DEAL::   Number of degrees of freedom for the field variable U: 22
+DEAL::   Number of degrees of freedom for the pressure variable p: 14
+DEAL::1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 1.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
+DEAL::   Number of active cells:       15
+DEAL::  Total number of degrees of freedom: 151
+DEAL::   Number of degrees of freedom for the field variable U: 105
+DEAL::   Number of degrees of freedom for the pressure variable p: 105
+DEAL::1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0