]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
new fe_values::get_function* for nonprimitive elements
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 12 May 2005 07:13:50 +0000 (07:13 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 12 May 2005 07:13:50 +0000 (07:13 +0000)
git-svn-id: https://svn.dealii.org/trunk@10665 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_values.h
deal.II/deal.II/source/fe/fe_values.cc
deal.II/deal.II/source/fe/fe_values.instance.h

index c9fbc2505bf5b273fb028da556f4acef14f941ce..d40864b35236a36d78eae83debf854aa28c51d34 100644 (file)
@@ -706,6 +706,69 @@ class FEValuesBase : protected FEValuesData<dim>
                              const VectorSlice<const std::vector<unsigned int> >& indices,
                              std::vector<Vector<number> >& values) const;
 
+
+                                    /**
+                                     * Generate vector function
+                                     * values from an arbitrary
+                                     * vector.
+                                     *
+                                     * This function offers the
+                                     * possibility to extract
+                                     * function values in quadrature
+                                     * points from vectors not
+                                     * corresponding to a whole
+                                     * discretization.
+                                     *
+                                     * The length of the vector
+                                     * <tt>indices</tt> may even be a
+                                     * multiple of the number of dofs
+                                     * per cell. Then, the vectors in
+                                     * <tt>value</tt> should allow
+                                     * for the same multiple of the
+                                     * components of the finite
+                                     * element.
+                                     *
+                                     * Depending on the last
+                                     * argument, the outer vector of
+                                     * <tt>values</tt> has either the
+                                     * length of the quadrature rule
+                                     * (<tt>quadrature_points_fastest
+                                     * == false</tt>) or the length
+                                     * of components to be filled
+                                     * <tt>quadrature_points_fastest
+                                     * == true</tt>. If <tt>p</tt> is
+                                     * the xurrent quadrature point
+                                     * number and <tt>i</tt> is the
+                                     * vector component of the
+                                     * solution desired, the access
+                                     * to <tt>values</tt> is
+                                     * <tt>values[p][i]</tt> if
+                                     * <tt>quadrature_points_fastest
+                                     * == false</tt>, and
+                                     * <tt>values[i][p]</tt>
+                                     * otherwise.
+                                     *
+                                     * You may want to use this
+                                     * function, if you want to
+                                     * access just a single block
+                                     * from a BlockVector, if you
+                                     * have a multi-level vector or
+                                     * if you already have a local
+                                     * representation of your finite
+                                     * element data.
+                                     *
+                                     * Since this function allows for
+                                     * fairly general combinations of
+                                     * argument sizes, be aware that
+                                     * the checks on the arguments
+                                     * may not detect errors.
+                                     */
+    template <class InputVector, typename number>
+    void get_function_values (const InputVector& fe_function,
+                             const VectorSlice<const std::vector<unsigned int> >& indices,
+                             std::vector<std::vector<number> >& values,
+                             bool quadrature_points_fastest) const;
+
                                     /**
                                      * Compute the gradients of the finite
                                      * element function characterized
@@ -806,7 +869,8 @@ class FEValuesBase : protected FEValuesData<dim>
     template <class InputVector>
     void get_function_grads (const InputVector& fe_function,
                             const VectorSlice<const std::vector<unsigned int> >& indices,
-                            std::vector<std::vector<Tensor<1,dim> > >& gradients) const;
+                            std::vector<std::vector<Tensor<1,dim> > >& gradients,
+                            bool quadrature_points_fastest = false) const;
 
                                     /**
                                      * Compute the tensor of second
@@ -885,7 +949,8 @@ class FEValuesBase : protected FEValuesData<dim>
     template <class InputVector>
     void
     get_function_2nd_derivatives (const InputVector      &fe_function,
-                                  std::vector<std::vector<Tensor<2,dim> > > &second_derivatives) const;
+                                  std::vector<std::vector<Tensor<2,dim> > > &second_derivatives,
+                                 bool quadrature_points_fastest = false) const;
                                     //@}
     
                                     /**
index 3e9f1a6c98606365737686913a56ca38d9b2fa22..477e360bc38de05f7d5aef5956286ea862142a2e 100644 (file)
@@ -517,7 +517,7 @@ void FEValuesBase<dim>::get_function_values (
                                   // result may be a multiple of the
                                   // number of components of the
                                   // finite element
-  const unsigned int result_components = indices.size() / dofs_per_cell;
+  const unsigned int result_components = indices.size() * n_components / dofs_per_cell;
   
   for (unsigned i=0;i<values.size();++i)
     Assert (values[i].size() == result_components,
@@ -551,8 +551,102 @@ void FEValuesBase<dim>::get_function_values (
                * shape_value(shape_func, point));
        else
          for (unsigned int c=0; c<n_components; ++c)
-           values[point](c) += (fe_function(indices[shape_func]) *
-                                shape_value_component(shape_func, point, c));
+           values[point](c+mc*n_components)
+             += (fe_function(indices[shape_func])
+                 * shape_value_component(shape_func, point, c));
+}
+
+
+
+template <int dim>
+template <class InputVector, typename number>
+void FEValuesBase<dim>::get_function_values (
+  const InputVector& fe_function,
+  const VectorSlice<const std::vector<unsigned int> >& indices,
+  std::vector<std::vector<number> >& values,
+  bool quadrature_points_fastest) const
+{
+  const unsigned int n_components = fe->n_components();
+  
+                                  // Size of indices must be a
+                                  // multiple of dofs_per_cell such
+                                  // that an integer number of
+                                  // function values is generated in
+                                  // each point.
+  Assert (indices.size() % dofs_per_cell == 0,
+         ExcNotMultiple(indices.size(), dofs_per_cell));
+
+                                  // The number of components of the
+                                  // result may be a multiple of the
+                                  // number of components of the
+                                  // finite element
+  const unsigned int result_components = indices.size() * n_components / dofs_per_cell;
+
+                                  // Check if the value argument is
+                                  // initialized to the correct sizes
+  if (quadrature_points_fastest)
+    {
+      Assert (values.size() == result_components,
+             ExcDimensionMismatch(values.size(), result_components));
+      for (unsigned i=0;i<values.size();++i)
+       Assert (values[i].size() == n_quadrature_points,
+               ExcDimensionMismatch(values[i].size(), n_quadrature_points));
+    }
+  else
+    {
+      Assert(values.size() == n_quadrature_points,
+            ExcDimensionMismatch(values.size(), n_quadrature_points));
+      for (unsigned i=0;i<values.size();++i)
+       Assert (values[i].size() == result_components,
+               ExcDimensionMismatch(values[i].size(), result_components));
+    }
+  
+                                  // If the result has more
+                                  // components than the finite
+                                  // element, we need this number for
+                                  // loops filling all components
+  const unsigned int component_multiple = result_components / n_components;
+  
+  Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
+    
+                                  // initialize with zero
+  for (unsigned i=0;i<values.size();++i)
+    std::fill_n (values[i].begin(), values[i].size(), 0);
+
+                                  // add up contributions of trial
+                                  // functions. now check whether the
+                                  // shape function is primitive or
+                                  // not. if it is, then set its only
+                                  // non-zero component, otherwise
+                                  // loop over components
+  if (quadrature_points_fastest)
+    for (unsigned int mc = 0; mc < component_multiple; ++mc)
+      for (unsigned int point=0; point<n_quadrature_points; ++point)
+       for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+         if (fe->is_primitive(shape_func))
+           values[fe->system_to_component_index(shape_func).first
+                  +mc * n_components][point]
+             += (fe_function(indices[shape_func+mc*dofs_per_cell])
+                 * shape_value(shape_func, point));
+         else
+           for (unsigned int c=0; c<n_components; ++c)
+             values[c+mc*n_components][point]
+               += (fe_function(indices[shape_func])
+                   * shape_value_component(shape_func, point, c));
+  else
+    for (unsigned int mc = 0; mc < component_multiple; ++mc)
+      for (unsigned int point=0; point<n_quadrature_points; ++point)
+       for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+         if (fe->is_primitive(shape_func))
+           values[point][fe->system_to_component_index(shape_func).first
+                         +mc * n_components]
+             += (fe_function(indices[shape_func+mc*dofs_per_cell])
+                 * shape_value(shape_func, point));
+         else
+           for (unsigned int c=0; c<n_components; ++c)
+             values[point][c+mc*n_components]
+               += (fe_function(indices[shape_func])
+                   * shape_value_component(shape_func, point, c));
 }
 
 
@@ -675,7 +769,7 @@ get_function_grads (const InputVector                         &fe_function,
        {
          Tensor<1,dim> tmp = this->shape_grad(shape_func,point);
          tmp *= dof_values(shape_func);
-         gradients[point][fe->system_to_component_index(shape_func).first]
+         gradients[fe->system_to_component_index(shape_func).first][point]
            += tmp;
        }
       else
@@ -683,8 +777,8 @@ get_function_grads (const InputVector                         &fe_function,
          {
            Tensor<1,dim> tmp = this->shape_grad_component(shape_func,point,c);
            tmp *= dof_values(shape_func);
-           gradients[point][c] += tmp;
-         };
+           gradients[c][point] += tmp;
+         }
 }
 
 
@@ -694,12 +788,9 @@ template <class InputVector>
 void FEValuesBase<dim>::get_function_grads (
   const InputVector& fe_function,
   const VectorSlice<const std::vector<unsigned int> >& indices,
-  std::vector<std::vector<Tensor<1,dim> > >& values) const
+  std::vector<std::vector<Tensor<1,dim> > >& values,
+  bool quadrature_points_fastest) const
 {
-                                  // One value per quadrature point
-  Assert (n_quadrature_points == values.size(),
-         ExcDimensionMismatch(values.size(), n_quadrature_points));
-  
   const unsigned int n_components = fe->n_components();
   
                                   // Size of indices must be a
@@ -714,11 +805,26 @@ void FEValuesBase<dim>::get_function_grads (
                                   // result may be a multiple of the
                                   // number of components of the
                                   // finite element
-  const unsigned int result_components = indices.size() / dofs_per_cell;
+  const unsigned int result_components = indices.size() * n_components / dofs_per_cell;
   
-  for (unsigned i=0;i<values.size();++i)
-    Assert (values[i].size() == result_components,
-           ExcDimensionMismatch(values[i].size(), result_components));
+                                  // Check if the value argument is
+                                  // initialized to the correct sizes
+  if (quadrature_points_fastest)
+    {
+      Assert (values.size() == result_components,
+             ExcDimensionMismatch(values.size(), result_components));
+      for (unsigned i=0;i<values.size();++i)
+       Assert (values[i].size() == n_quadrature_points,
+               ExcDimensionMismatch(values[i].size(), n_quadrature_points));
+    }
+  else
+    {
+      Assert(values.size() == n_quadrature_points,
+            ExcDimensionMismatch(values.size(), n_quadrature_points));
+      for (unsigned i=0;i<values.size();++i)
+       Assert (values[i].size() == result_components,
+               ExcDimensionMismatch(values[i].size(), result_components));
+    }
 
                                   // If the result has more
                                   // components than the finite
@@ -738,18 +844,32 @@ void FEValuesBase<dim>::get_function_grads (
                                   // not. if it is, then set its only
                                   // non-zero component, otherwise
                                   // loop over components
-  for (unsigned int mc = 0; mc < component_multiple; ++mc)
-    for (unsigned int point=0; point<n_quadrature_points; ++point)
-      for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
-       if (fe->is_primitive(shape_func))
-         values[point][fe->system_to_component_index(shape_func).first
-                       +mc * n_components]
-           += fe_function(indices[shape_func+mc*dofs_per_cell])
-               * shape_grad(shape_func, point);
-       else
-         for (unsigned int c=0; c<n_components; ++c)
-           values[point][c] += (fe_function(indices[shape_func]) *
-                                shape_grad_component(shape_func, point, c));
+  if (quadrature_points_fastest)
+    for (unsigned int mc = 0; mc < component_multiple; ++mc)
+      for (unsigned int point=0; point<n_quadrature_points; ++point)
+       for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+         if (fe->is_primitive(shape_func))
+           values[point][fe->system_to_component_index(shape_func).first
+                         +mc * n_components]
+             += fe_function(indices[shape_func+mc*dofs_per_cell])
+             * shape_grad(shape_func, point);
+         else
+           for (unsigned int c=0; c<n_components; ++c)
+             values[point][c] += (fe_function(indices[shape_func]) *
+                                  shape_grad_component(shape_func, point, c));
+  else
+    for (unsigned int mc = 0; mc < component_multiple; ++mc)
+      for (unsigned int point=0; point<n_quadrature_points; ++point)
+       for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+         if (fe->is_primitive(shape_func))
+           values[fe->system_to_component_index(shape_func).first
+                  +mc * n_components][point]
+             += fe_function(indices[shape_func+mc*dofs_per_cell])
+             * shape_grad(shape_func, point);
+         else
+           for (unsigned int c=0; c<n_components; ++c)
+             values[c][point] += (fe_function(indices[shape_func]) *
+                                  shape_grad_component(shape_func, point, c));
 }
 
 
@@ -799,7 +919,8 @@ template <class InputVector>
 void
 FEValuesBase<dim>::
 get_function_2nd_derivatives (const InputVector                         &fe_function,
-                             std::vector<std::vector<Tensor<2,dim> > > &second_derivs) const
+                             std::vector<std::vector<Tensor<2,dim> > > &second_derivs,
+                             bool quadrature_points_fastest) const
 {
   Assert (n_quadrature_points == second_derivs.size(),
          ExcDimensionMismatch(second_derivs.size(), n_quadrature_points));
@@ -824,22 +945,40 @@ get_function_2nd_derivatives (const InputVector                         &fe_func
 
                                   // add up contributions of trial
                                   // functions
-  for (unsigned int point=0; point<n_quadrature_points; ++point)
-    for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
-      if (fe->is_primitive(shape_func))
-       {
-         Tensor<2,dim> tmp(shape_2nd_derivative(shape_func,point));
-         tmp *= dof_values(shape_func);
-         second_derivs[point][fe->system_to_component_index(shape_func).first]
-           += tmp;
-       }
-      else
-       for (unsigned int c=0; c<n_components; ++c)
+  if (quadrature_points_fastest)
+    for (unsigned int point=0; point<n_quadrature_points; ++point)
+      for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+       if (fe->is_primitive(shape_func))
+         {
+           Tensor<2,dim> tmp(shape_2nd_derivative(shape_func,point));
+           tmp *= dof_values(shape_func);
+           second_derivs[point][fe->system_to_component_index(shape_func).first]
+             += tmp;
+         }
+       else
+         for (unsigned int c=0; c<n_components; ++c)
+           {
+             Tensor<2,dim> tmp = this->shape_2nd_derivative_component(shape_func,point,c);
+             tmp *= dof_values(shape_func);
+             second_derivs[point][c] += tmp;
+           }
+  else
+    for (unsigned int point=0; point<n_quadrature_points; ++point)
+      for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
+       if (fe->is_primitive(shape_func))
          {
-           Tensor<2,dim> tmp = this->shape_2nd_derivative_component(shape_func,point,c);
+           Tensor<2,dim> tmp(shape_2nd_derivative(shape_func,point));
            tmp *= dof_values(shape_func);
-           second_derivs[point][c] += tmp;
-         };
+           second_derivs[fe->system_to_component_index(shape_func).first][point]
+             += tmp;
+         }
+       else
+         for (unsigned int c=0; c<n_components; ++c)
+           {
+             Tensor<2,dim> tmp = this->shape_2nd_derivative_component(shape_func,point,c);
+             tmp *= dof_values(shape_func);
+             second_derivs[c][point] += tmp;
+           }
 }
 
 
index ac9e17c272ed8ee104f8b3f86e46ea6a59909b98..f02bebe622631d89991f03cbf745f0efb60e16d5 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -45,6 +45,15 @@ void FEValuesBase<deal_II_dimension>::get_function_values<IN>
 (const IN&, const VectorSlice<const std::vector<unsigned int> >&,
  std::vector<Vector<float> > &) const;
 
+template
+void FEValuesBase<deal_II_dimension>::get_function_values<IN>
+(const IN&, const VectorSlice<const std::vector<unsigned int> >&,
+ std::vector<std::vector<double> > &, bool) const;
+template
+void FEValuesBase<deal_II_dimension>::get_function_values<IN>
+(const IN&, const VectorSlice<const std::vector<unsigned int> >&,
+ std::vector<std::vector<float> > &, bool) const;
+
 template
 void FEValuesBase<deal_II_dimension>::get_function_grads<IN>
 (const IN&, std::vector<Tensor<1,deal_II_dimension> > &) const;
@@ -59,7 +68,7 @@ void FEValuesBase<deal_II_dimension>::get_function_grads<IN>
 template
 void FEValuesBase<deal_II_dimension>::get_function_grads<IN>
 (const IN&, const VectorSlice<const std::vector<unsigned int> >&,
- std::vector<std::vector<Tensor<1,deal_II_dimension> > > &) const;
+ std::vector<std::vector<Tensor<1,deal_II_dimension> > > &, bool) const;
 
 template
 void FEValuesBase<deal_II_dimension>::get_function_2nd_derivatives<IN>
@@ -67,4 +76,4 @@ void FEValuesBase<deal_II_dimension>::get_function_2nd_derivatives<IN>
 
 template
 void FEValuesBase<deal_II_dimension>::get_function_2nd_derivatives<IN>
-(const IN&, std::vector<std::vector<Tensor<2,deal_II_dimension> > > &) const;
+(const IN&, std::vector<std::vector<Tensor<2,deal_II_dimension> > > &, bool) const;

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.