]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix another bug with non-primitive elements.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 17 Feb 2003 15:59:20 +0000 (15:59 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 17 Feb 2003 15:59:20 +0000 (15:59 +0000)
git-svn-id: https://svn.dealii.org/trunk@7125 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/vectors.h
deal.II/deal.II/source/numerics/vectors.cc
tests/bits/Makefile
tests/bits/anna_6.cc [new file with mode: 0644]
tests/bits/anna_6.checked [new file with mode: 0644]

index 19d93ac6fc82d4d7fe504803be55b98158c27cfd..cdba28b0afff8edd4efdb137c82d05b4e5679831 100644 (file)
@@ -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
index 116ba292f24096941b05588af15547166ae18fa9..1bf6f6a19d7e82c26f3a52e042f30f1592048b31 100644 (file)
@@ -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,
index 41d99e11ee5314c127409c34a3ba062b65a3f4d1..f774ac799b75b82a4dff491e07b059835ee287aa 100644 (file)
@@ -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 (file)
index 0000000..07b2aaa
--- /dev/null
@@ -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 (file)
index 0000000..e91fa24
--- /dev/null
@@ -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 

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.