]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Patch from Markus Buerg: Implement FEValuesViews::Vector::curl and friends.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 13 May 2010 13:22:15 +0000 (13:22 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 13 May 2010 13:22:15 +0000 (13:22 +0000)
git-svn-id: https://svn.dealii.org/trunk@21127 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/doc/news/changes.h

index 5d1b482199f9cdf5506621e063cc3ad306decc90..5eec7183926516e8c75e9b9b4bbe4ecdab41cbb5 100644 (file)
@@ -500,6 +500,19 @@ namespace FEValuesViews
                                        */
       typedef double                 divergence_type;
 
+                                      /**
+                                       * A typedef for the type of the curl
+                                       * of the view this class
+                                       * represents. Here, for a set of
+                                       * <code>dim</code> components of the
+                                       * finite element, the curl is of type
+                                       * <code>Tensor@<1, 1@></code> (i.e. a
+                                       * scalar) in 2d, and of type
+                                       * <code>Tensor@<1, 3@></code> (i.e. a
+                                       * vector) in 3d.
+                                       */
+      typedef Tensor<1, (spacedim == 3)? 3 : 1>     curl_type;
+
                                       /**
                                        * A typedef for the type of second
                                        * derivatives of the view this class
@@ -598,6 +611,17 @@ namespace FEValuesViews
       divergence (const unsigned int shape_function,
                  const unsigned int q_point) const;
 
+                                      /**
+                                       * Return the vector curl of
+                                       * the vector components selected by
+                                       * this view, for the shape function
+                                       * and quadrature point selected by the
+                                       * arguments.
+                                       */
+      curl_type
+      curl (const unsigned int shape_function,
+           const unsigned int q_point) const;
+
                                       /**
                                        * Return the Hessian (the tensor of
                                        * rank 2 of all second derivatives) of
@@ -689,7 +713,7 @@ namespace FEValuesViews
                                        *
                                        * There is no equivalent function such
                                        * as
-                                       * FEValuesBase::get_function_gradients
+                                       * FEValuesBase::get_function_divergences
                                        * in the FEValues classes but the
                                        * information can be obtained from
                                        * FEValuesBase::get_function_gradients,
@@ -699,6 +723,29 @@ namespace FEValuesViews
       void get_function_divergences (const InputVector& fe_function,
                                     std::vector<divergence_type>& divergences) const;
 
+                                      /**
+                                       * Return the curl of the selected
+                                       * vector components of the finite
+                                       * element function characterized by
+                                       * <tt>fe_function</tt> at the
+                                       * quadrature points of the cell, face
+                                       * or subface selected the last time
+                                       * the <tt>reinit</tt> function of the
+                                       * FEValues object was called, at the
+                                       * quadrature points.
+                                       *
+                                       * There is no equivalent function such
+                                       * as
+                                       * FEValuesBase::get_function_curls
+                                       * in the FEValues classes but the
+                                       * information can be obtained from
+                                       * FEValuesBase::get_function_gradients,
+                                       * of course.
+                                       */
+      template <class InputVector>
+      void get_function_curls (const InputVector& fe_function,
+                              std::vector<curl_type>& curls) const;
+
                                       /**
                                        * Return the Hessians of the selected
                                        * vector components of the finite
@@ -3730,6 +3777,120 @@ namespace FEValuesViews
 
   template <int dim, int spacedim>
   inline
+  typename Vector<dim,spacedim>::curl_type
+  Vector<dim,spacedim>::curl (const unsigned int shape_function,
+                             const unsigned int q_point) const
+  {
+                                    // this function works like in the case
+                                    // above
+     typedef FEValuesBase<dim,spacedim> FVB;
+
+     Assert (shape_function < fe_values.fe->dofs_per_cell,
+            ExcIndexRange (shape_function, 0, fe_values.fe->dofs_per_cell));
+     Assert (fe_values.update_flags & update_gradients,
+            typename FVB::ExcAccessToUninitializedField());
+                                     // same as for the scalar case except
+                                     // that we have one more index
+     const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+     if (snc == -2)
+        return curl_type ();
+
+     else
+        switch (dim) {
+           case 1:
+              return curl_type ();
+
+           case 2: {
+              if (snc != -1) {
+                 curl_type return_value;
+
+                 switch (shape_function_data[shape_function].single_nonzero_component_index) {
+                    case 0: {
+                       return_value[0] = -1.0 * fe_values.shape_gradients[snc][q_point][1];
+                       return return_value;
+                    }
+
+                    default: {
+                       return_value[0] = fe_values.shape_gradients[snc][q_point][2];
+                       return return_value;
+                    }
+                 }
+              }
+
+              else {
+                 curl_type return_value;
+
+                 return_value[0] = 0.0;
+
+                 if (shape_function_data[shape_function].is_nonzero_shape_function_component[0])
+                    return_value[0] -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][1];
+
+                 if (shape_function_data[shape_function].is_nonzero_shape_function_component[1])
+                    return_value[0] += fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][0];
+
+                 return return_value;
+              }
+           }
+
+           case 3: {
+              if (snc != -1) {
+                 curl_type return_value;
+
+                 switch (shape_function_data[shape_function].single_nonzero_component_index) {
+                    case 0: {
+                       return_value[0] = 0;
+                       return_value[1] = fe_values.shape_gradients[snc][q_point][2];
+                       return_value[2] = -1.0 * fe_values.shape_gradients[snc][q_point][1];
+                       return return_value;
+                    }
+
+                    case 1: {
+                       return_value[0] = -1.0 * fe_values.shape_gradients[snc][q_point][2];
+                       return_value[1] = 0;
+                       return_value[2] = fe_values.shape_gradients[snc][q_point][0];
+                       return return_value;
+                    }
+
+                    default: {
+                       return_value[0] = fe_values.shape_gradients[snc][q_point][1];
+                       return_value[1] = -1.0 * fe_values.shape_gradients[snc][q_point][0];
+                       return_value[2] = 0;
+                       return return_value;
+                    }
+                 }
+              }
+
+              else {
+                 curl_type return_value;
+
+                 for (unsigned int i = 0; i < dim; ++i)
+                    return_value[i] = 0.0;
+
+                 if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) {
+                    return_value[1] += fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][2];
+                    return_value[2] -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][1];
+                 }
+
+                 if (shape_function_data[shape_function].is_nonzero_shape_function_component[1]) {
+                    return_value[0] -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][2];
+                    return_value[2] += fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][0];
+                 }
+
+                 if (shape_function_data[shape_function].is_nonzero_shape_function_component[2]) {
+                    return_value[0] += fe_values.shape_gradients[shape_function_data[shape_function].row_index[2]][q_point][1];
+                    return_value[1] -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[2]][q_point][0];
+                 }
+
+                 return return_value;
+              }
+           }
+        }
+  }
+
+  
+  template <int dim, int spacedim>  
+  inline
   typename Vector<dim,spacedim>::hessian_type
   Vector<dim,spacedim>::hessian (const unsigned int shape_function,
                                 const unsigned int q_point) const
index c141fe27716537a907beed4d12ce8f4b42697b30..ff010357785250086f09956330f18a46554651cd 100644 (file)
@@ -755,9 +755,168 @@ namespace FEValuesViews
       }
   }
 
+  template <int dim, int spacedim>
+  template <class InputVector>
+  void
+  Vector<dim,spacedim>::get_function_curls (const InputVector &fe_function,
+                                           std::vector<curl_type> &curls) const {
+     typedef FEValuesBase<dim,spacedim> FVB;
+
+     Assert (fe_values.update_flags & update_gradients,
+            typename FVB::ExcAccessToUninitializedField());
+     Assert (curls.size() == fe_values.n_quadrature_points,
+            ExcDimensionMismatch (curls.size(), fe_values.n_quadrature_points));
+     Assert (fe_values.present_cell.get () != 0,
+            ExcMessage ("FEValues object is not reinit'ed to any cell"));
+     Assert (fe_function.size () == fe_values.present_cell->n_dofs_for_dof_handler (),
+            ExcDimensionMismatch (fe_function.size (), fe_values.present_cell->n_dofs_for_dof_handler ()));
+                                     // get function values of dofs on this
+                                     // cell
+     dealii::Vector<typename InputVector::value_type> dof_values (fe_values.dofs_per_cell);
+     fe_values.present_cell->get_interpolated_dof_values (fe_function, dof_values);
+
+     std::fill (curls.begin (), curls.end (), curl_type ());
+
+     switch (dim) {
+        case 2: {
+           for (unsigned int shape_function = 0; shape_function < fe_values.fe->dofs_per_cell; ++shape_function) {
+              const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+              if (snc == -2)
+                 // shape function is zero for the selected components
+                 continue;
+
+              const double value = dof_values (shape_function);
+
+              if (value == 0.)
+                 continue;
+
+              if (snc != -1) {
+                 const unsigned int comp = shape_function_data[shape_function].single_nonzero_component_index;
+                 const Tensor<1, spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[snc][0];
+
+                 switch (shape_function_data[shape_function].single_nonzero_component_index) {
+                    case 0: {
+                       for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+                          curls[q_point][0] = -1.0 * value * (*shape_gradient_ptr++)[1];
+
+                       break;
+                    }
+
+                    default:
+                       for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+                          curls[q_point][0] = value * (*shape_gradient_ptr)[0];
+                 }
+              }
+
+              else {
+                 for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+                    curls[q_point][0] = 0;
+
+                 if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) {
+                    const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][0];
+
+                    for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+                       curls[q_point][0] -= value * (*shape_gradient_ptr++)[1];
+                 }
+
+                 if (shape_function_data[shape_function].is_nonzero_shape_function_component[1]) {
+                    const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][0];
+
+                    for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+                       curls[q_point][0] += value * (*shape_gradient_ptr++)[0];
+                 }
+              }
+           }
+           break;
+        }
+
+        case 3: {
+           for (unsigned int shape_function = 0; shape_function < fe_values.fe->dofs_per_cell; ++shape_function) {
+              const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+              if (snc == -2)
+                 // shape function is zero for the selected components
+                 continue;
+
+              const double value = dof_values (shape_function);
+
+              if (value == 0.)
+                 continue;
+
+              if (snc != -1) {
+                 const unsigned int comp = shape_function_data[shape_function].single_nonzero_component_index;
+                 const Tensor<1, spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[snc][0];
+
+                 switch (shape_function_data[shape_function].single_nonzero_component_index) {
+                    case 0: {
+                       for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+                          curls[q_point][0] = 0;
+                          curls[q_point][1] = value * (*shape_gradient_ptr)[2];
+                          curls[q_point][2] = -1.0 * value * (*shape_gradient_ptr++)[1];
+                       }
+
+                       break;
+                    }
+
+                    case 1: {
+                       for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+                          curls[q_point][0] = -1.0 * value * (*shape_gradient_ptr)[2];
+                          curls[q_point][1] = 0;
+                          curls[q_point][2] = value * (*shape_gradient_ptr++)[0];
+                       }
+
+                       break;
+                    }
+
+                    default:
+                       for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+                          curls[q_point][0] = value * (*shape_gradient_ptr)[1];
+                          curls[q_point][1] = -1.0 * value * (*shape_gradient_ptr++)[0];
+                          curls[q_point][2] = 0;
+                       }
+                 }
+              }
+
+              else {
+                 for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+                    for (unsigned int d = 0; d < dim; ++d)
+                       curls[q_point][d] = 0;
+
+                 if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) {
+                    const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][0];
+
+                    for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+                       curls[q_point][1] += value * (*shape_gradient_ptr++)[2];
+                       curls[q_point][2] -= value * (*shape_gradient_ptr++)[1];
+                    }
+                 }
+
+                 if (shape_function_data[shape_function].is_nonzero_shape_function_component[1]) {
+                    const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][0];
+
+                    for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+                       curls[q_point][0] -= value * (*shape_gradient_ptr++)[2];
+                       curls[q_point][2] += value * (*shape_gradient_ptr++)[0];
+                    }
+                 }
+
+                 if (shape_function_data[shape_function].is_nonzero_shape_function_component[2]) {
+                    const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[2]][0];
+
+                    for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) {
+                       curls[q_point][0] += value * (*shape_gradient_ptr++)[1];
+                       curls[q_point][1] -= value * (*shape_gradient_ptr++)[0];
+                    }
+                 }
+              }
+           }
+        }
+     }
+  }
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim>  
   template <class InputVector>
   void
   Vector<dim,spacedim>::
index cae2c5aa8cf904b3c2d2e43185e733701062e8f7..e216d88a2589d57d4587f89e455d77acff15aab8 100644 (file)
@@ -636,6 +636,13 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li>
+  <p>New: The FEValuesViews::Vector class now has functions
+  FEValuesViews::Vector::curl and FEValuesViews::Vector::get_function_curls.
+  <br>
+  (Markus B&uuml;rg 2010/05/13)
+  </p></li>
+
   <li>
   <p>New: TriaAccessor::extent_in_direction() returns the length
   of an object in a given direction.

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.