]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fixed bug in FEValuesViews::SymmetricTensor< 2, dim, spacedim >::get_function_divergences
authormcbride <mcbride@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 9 Sep 2010 08:55:23 +0000 (08:55 +0000)
committermcbride <mcbride@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 9 Sep 2010 08:55:23 +0000 (08:55 +0000)
git-svn-id: https://svn.dealii.org/trunk@21902 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/fe_values.cc

index 3c62232cce0e0c33b1058825479265d71931f953..1dd877aeeb4b294d76e63a838b8c360ce5fd6b8c 100644 (file)
@@ -1121,79 +1121,90 @@ namespace FEValuesViews
   void
   SymmetricTensor<2, dim, spacedim>::
   get_function_divergences(const InputVector &fe_function,
-                          std::vector<divergence_type> &divergences) const
+                           std::vector<divergence_type> &divergences) const
   {
     typedef FEValuesBase<dim, spacedim> FVB;
     Assert(fe_values.update_flags & update_gradients,
-          typename FVB::ExcAccessToUninitializedField());
+           typename FVB::ExcAccessToUninitializedField());
     Assert(divergences.size() == fe_values.n_quadrature_points,
-          ExcDimensionMismatch(divergences.size(), fe_values.n_quadrature_points));
+           ExcDimensionMismatch(divergences.size(), fe_values.n_quadrature_points));
     Assert(fe_values.present_cell.get() != 0,
-          ExcMessage("FEValues object is not reinit'ed to any cell"));
+           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()));
+           ExcDimensionMismatch(fe_function.size(),
+                                fe_values.present_cell->n_dofs_for_dof_handler()));
 
-                                    // get function values of dofs
-                                    // on this cell
+                                     // 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(divergences.begin(), divergences.end(), divergence_type());
 
     for (unsigned int shape_function = 0;
-        shape_function < fe_values.fe->dofs_per_cell; ++shape_function)
+         shape_function < fe_values.fe->dofs_per_cell; ++shape_function)
       {
-       const int snc = shape_function_data[shape_function].single_nonzero_component;
+        const int snc = shape_function_data[shape_function].single_nonzero_component;
 
-       if (snc == -2)
-                                          // shape function is zero for the
-                                          // selected components
-         continue;
+        if (snc == -2)
+                                           // shape function is zero for the
+                                           // selected components
+          continue;
 
-       const double value = dof_values(shape_function);
-       if (value == 0.)
-         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;
+        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];
-           for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
-                ++q_point, ++shape_gradient_ptr) {
-             for (unsigned int j = 0; j < dim; ++j)
-               divergences[q_point][value_type::unrolled_to_component_indices(comp)[0]]
-                 += value * (*shape_gradient_ptr)[j];
-           }
+            const Tensor < 1, spacedim> * shape_gradient_ptr =
+              &fe_values.shape_gradients[snc][0];
+
+            const unsigned int ii = value_type::unrolled_to_component_indices(comp)[0];
+            const unsigned int jj = value_type::unrolled_to_component_indices(comp)[1];
+
+            for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
+                 ++q_point, ++shape_gradient_ptr) {
+
+                divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
+
+                if (ii != jj)
+                    divergences[q_point][jj] += value * (*shape_gradient_ptr)[ii];
+            }
          }
-       else
+        else
          {
-           for (unsigned int d = 0; d < value_type::n_independent_components; ++d)
-             if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
-               {
+            for (unsigned int d = 0; d < value_type::n_independent_components; ++d)
+              if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
+                {
                  Assert (false, ExcNotImplemented());
 
-// the following implementation needs to be looked over -- I think it
-// can't be right, because we are in a case where there is no single
-// nonzero component
+                    // the following implementation needs to be looked over -- I think it
+                    // can't be right, because we are in a case where there is no single
+                    // nonzero component
+                    //
+                    // the following is not implemented! we need to consider the interplay between
+                    // mutliple non-zero entries in shape function and the representation
+                    // as a symmetric second-order tensor
+
                  const unsigned int comp =
-                   shape_function_data[shape_function].single_nonzero_component_index;
+                    shape_function_data[shape_function].single_nonzero_component_index;
 
                  const Tensor < 1, spacedim> * shape_gradient_ptr =
-                   &fe_values.shape_gradients[shape_function_data[shape_function].
-                                              row_index[d]][0];
+                    &fe_values.shape_gradients[shape_function_data[shape_function].
+                                               row_index[d]][0];
                  for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
                       ++q_point, ++shape_gradient_ptr) {
-                   for (unsigned int j = 0; j < dim; ++j)
+                    for (unsigned int j = 0; j < dim; ++j)
                      {
-                       const unsigned int vector_component = value_type::component_to_unrolled_index (TableIndices<2>(comp,j));
-                       divergences[q_point][vector_component] += value * (*shape_gradient_ptr++)[j];
+                        const unsigned int vector_component = value_type::component_to_unrolled_index (TableIndices<2>(comp,j));
+                        divergences[q_point][vector_component] += value * (*shape_gradient_ptr++)[j];
                      }
                  }
-               }
+                }
          }
       }
   }

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.