]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Let shape_value return a reference to the data, and not only the value. This allows...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 5 Jun 2009 11:10:59 +0000 (11:10 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 5 Jun 2009 11:10:59 +0000 (11:10 +0000)
git-svn-id: https://svn.dealii.org/trunk@18914 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_values.h
deal.II/deal.II/source/fe/fe_tools.cc
deal.II/deal.II/source/numerics/matrices.cc

index 5efbd524025a18d59098d3f502b5858a4a4d4d69..ee712fc875c93c2a4119516a55e43093a9d9f93d 100644 (file)
@@ -1298,8 +1298,8 @@ class FEValuesBase : protected FEValuesData<dim,spacedim>,
                                      * the quadrature point at which
                                      * function is to be computed
                                      */
-    double shape_value (const unsigned int function_no,
-                       const unsigned int point_no) const;
+    const double & shape_value (const unsigned int function_no,
+                               const unsigned int point_no) const;
 
                                     /**
                                      * Compute one vector component of
@@ -3789,7 +3789,7 @@ operator[] (const FEValuesExtractors::Vector &vector) const
 
 template <int dim, int spacedim>  
 inline
-double
+const double &
 FEValuesBase<dim,spacedim>::shape_value (const unsigned int i,
                                         const unsigned int j) const
 {
index 5df260a7ce2e5d5cf6bdd976e92b1c07521c15ce..9e2c36dbd28bce862999bc012ab8878f68a20143 100644 (file)
@@ -637,9 +637,9 @@ FETools::compute_embedding_matrices(const FiniteElement<dim,spacedim>& fe,
                                       // children.
       fine.reinit(tria.begin_active());
       FullMatrix<number> A(nq*nd, n);
-      for (unsigned int d=0;d<nd;++d)
-       for (unsigned int k=0;k<nq;++k)
-         for (unsigned int j=0;j<n;++j)
+      for (unsigned int j=0;j<n;++j)
+       for (unsigned int d=0;d<nd;++d)
+         for (unsigned int k=0;k<nq;++k)
            A(k*nd+d,j) = fine.shape_value_component(j,k,d);
 
       Householder<double> H(A);
@@ -664,6 +664,7 @@ FETools::compute_embedding_matrices(const FiniteElement<dim,spacedim>& fe,
          coarse.reinit(tria.begin(0));
 
          FullMatrix<double> &this_matrix = matrices[ref_case-1][cell_number];
+         v_coarse = 0;
       
                                           // Compute this once for each
                                           // coarse grid basis function
@@ -675,9 +676,17 @@ FETools::compute_embedding_matrices(const FiniteElement<dim,spacedim>& fe,
                                               // function values of the
                                               // coarse grid function in
                                               // each quadrature point.
-             for (unsigned int d=0;d<nd;++d)
-               for (unsigned int k=0;k<nq;++k)
-                 v_coarse(k*nd+d) = coarse.shape_value_component(i,k,d);
+             if (fe.is_primitive())
+               {
+                 const unsigned int d = fe.system_to_component_index(i).first;
+                 const double * phi_i = &coarse.shape_value (i,0);
+                 for (unsigned int k=0;k<nq;++k)
+                   v_coarse(k*nd+d) = phi_i[k];
+               }
+             else
+               for (unsigned int d=0;d<nd;++d)
+                 for (unsigned int k=0;k<nq;++k)
+                   v_coarse(k*nd+d) = coarse.shape_value_component(i,k,d);
 
                                               // solve the least squares
                                               // problem.
@@ -810,9 +819,9 @@ FETools::compute_face_embedding_matrices(const FiniteElement<dim,spacedim>& fe,
                                   // children.
   fine.reinit(tria.begin_active());
   FullMatrix<number> A(nq*nd, n);
-  for (unsigned int d=0;d<nd;++d)
-    for (unsigned int k=0;k<nq;++k)
-      for (unsigned int j=0;j<n;++j)
+  for (unsigned int j=0;j<n;++j)
+    for (unsigned int d=0;d<nd;++d)
+      for (unsigned int k=0;k<nq;++k)
        A(k*nd+d,j) = fine.shape_value_component(face_f_dofs[j],k,d);
 
   Householder<double> H(A);
@@ -849,7 +858,7 @@ FETools::compute_face_embedding_matrices(const FiniteElement<dim,spacedim>& fe,
                                           // each quadrature point.
          for (unsigned int d=0;d<nd;++d)
            for (unsigned int k=0;k<nq;++k)
-             v_coarse(k*nd+d) = coarse.shape_value_component(face_c_dofs[i],k,d);
+             v_coarse(k*nd+d) = coarse.shape_value_component (face_c_dofs[i],k,d);
 
                                           // solve the least squares
                                           // problem.
@@ -904,9 +913,55 @@ FETools::compute_projection_matrices(const FiniteElement<dim,spacedim>& fe,
   const unsigned int nd = fe.n_components();
   const unsigned int degree = fe.degree;
 
+                                  // prepare FEValues, quadrature etc on
+                                  // coarse cell
+  MappingCartesian<dim> mapping;
+  QGauss<dim> q_fine(degree+1);
+  const unsigned int nq = q_fine.size();
+
+                                  // create mass matrix on coarse cell.
+  FullMatrix<number> mass(n, n);
+  {
+                                  // set up a triangulation for coarse cell
+    Triangulation<dim,spacedim> tr;
+    GridGenerator::hyper_cube (tr, 0, 1);
+
+    FEValues<dim> coarse (mapping, fe, q_fine,
+                         update_JxW_values | update_values);
+
+    typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
+      = tr.begin(0);
+    coarse.reinit (coarse_cell);
+
+    const std::vector<double> & JxW = coarse.get_JxW_values();
+    for (unsigned int i=0;i<n;++i)
+      for (unsigned int j=0;j<n;++j)
+       if (fe.is_primitive())
+         {
+           const double * coarse_i = &coarse.shape_value(i,0);
+           const double * coarse_j = &coarse.shape_value(j,0);
+           double mass_ij = 0;
+           for (unsigned int k=0;k<nq;++k)
+             mass_ij += JxW[k] * coarse_i[k] * coarse_j[k];
+           mass(i,j) = mass_ij;
+         }
+       else
+         {
+           double mass_ij = 0;
+           for (unsigned int d=0;d<nd;++d)
+             for (unsigned int k=0;k<nq;++k)
+               mass_ij += JxW[k] * coarse.shape_value_component(i,k,d)
+                                 * coarse.shape_value_component(j,k,d);
+           mass(i,j) = mass_ij;
+         }
+
+                                  // invert mass matrix
+    mass.gauss_jordan();
+  }
+
                                   // loop over all possible refinement cases
   for (unsigned int ref_case = RefinementCase<dim>::cut_x;
-       ref_case<RefinementCase<dim>::isotropic_refinement+1; ++ref_case)
+       ref_case < RefinementCase<dim>::isotropic_refinement+1; ++ref_case)
     {
       const unsigned int
        nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
@@ -918,45 +973,20 @@ FETools::compute_projection_matrices(const FiniteElement<dim,spacedim>& fe,
          Assert(matrices[ref_case-1][i].m() == n,
                 ExcDimensionMismatch(matrices[ref_case-1][i].m(),n));
        }
-  
+
+                                  // create a respective refinement on the
+                                  // triangulation
       Triangulation<dim,spacedim> tr;
       GridGenerator::hyper_cube (tr, 0, 1);
       tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
       tr.execute_coarsening_and_refinement();
-  
-      MappingCartesian<dim> mapping;
-      QGauss<dim> q_fine(degree+1);
-      const unsigned int nq = q_fine.size();
-  
-      FEValues<dim> coarse (mapping, fe, q_fine,
-                           update_quadrature_points |
-                           update_JxW_values |
-                           update_values);
+
       FEValues<dim> fine (mapping, fe, q_fine,
-                         update_quadrature_points |
-                         update_JxW_values |
+                         update_quadrature_points | update_JxW_values |
                          update_values);
-  
+
       typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
        = tr.begin(0);
-                                      // Compute the coarse level mass
-                                      // matrix
-      coarse.reinit(coarse_cell);
-      FullMatrix<number> A(n, n);
-      for (unsigned int k=0;k<nq;++k)
-       for (unsigned int i=0;i<n;++i)
-         for (unsigned int j=0;j<n;++j)
-           if (fe.is_primitive())
-             A(i,j) += coarse.JxW(k)
-                       * coarse.shape_value(i,k)
-                       * coarse.shape_value(j,k);
-           else
-             for (unsigned int d=0;d<nd;++d)
-               A(i,j) = coarse.JxW(k)
-                        * coarse.shape_value_component(i,k,d)
-                        * coarse.shape_value_component(j,k,d);
-  
-      Householder<double> H(A);
 
       Vector<number> v_coarse(n);
       Vector<number> v_fine(n);
@@ -977,37 +1007,39 @@ FETools::compute_projection_matrices(const FiniteElement<dim,spacedim>& fe,
       
                                           // Build RHS
 
+         const std::vector<double> & JxW = fine.get_JxW_values();
+
                                           // Outer loop over all fine
                                           // grid shape functions phi_j
          for (unsigned int j=0;j<fe.dofs_per_cell;++j)
            {
-             v_fine = 0.;
-                                              // Loop over all quadrature points
-             for (unsigned int k=0;k<fine.n_quadrature_points;++k)
+             for (unsigned int i=0; i<fe.dofs_per_cell;++i)
                {
-                                                  // integrate the scalar
-                                                  // product
-                                                  // (phi_i,phi_j) for
-                                                  // all coarse shape
-                                                  // functions to get the
-                                                  // right hand side
-                 for (unsigned int i=0;i<fe.dofs_per_cell;++i)
+                 if (fe.is_primitive())
                    {
-                     if (fe.is_primitive())
-                       v_fine(i) += fine.JxW(k)
-                                    * coarse.shape_value(i,k)
-                                    * fine.shape_value(j,k);
-                     else
-                       for (unsigned int d=0;d<nd;++d)
-                         v_fine(i) += fine.JxW(k)
-                                      * coarse.shape_value_component(i,k,d)
-                                      * fine.shape_value_component(j,k,d);
+                     const double * coarse_i = &coarse.shape_value(i,0);
+                     const double * fine_j = &fine.shape_value(j,0);
+
+                     double update = 0;
+                     for (unsigned int k=0; k<nq; ++k)
+                       update += JxW[k] * coarse_i[k] * fine_j[k];
+                     v_fine(i) = update;
+                   }
+                 else
+                   {
+                     double update = 0;
+                     for (unsigned int d=0; d<nd; ++d)
+                       for (unsigned int k=0; k<nq; ++k)
+                         update += JxW[k] * coarse.shape_value_component(i,k,d) 
+                                          * fine.shape_value_component(j,k,d);
+                     v_fine(i) = update;
                    }
                }
+
                                               // RHS ready. Solve system
                                               // and enter row into
                                               // matrix
-             H.least_squares(v_coarse, v_fine);
+             mass.vmult (v_coarse, v_fine);
              for (unsigned int i=0;i<fe.dofs_per_cell;++i)
                this_matrix(i,j) = v_coarse(i);
            }
index cc47a687d437563c0d4be81b78ee771f378745fc..58d8142dc2573ef2e3bd60ea34670f8e822f2a61 100644 (file)
 #include <base/thread_management.h>
 #include <base/work_stream.h>
 #include <base/multithread_info.h>
-#include <grid/tria_iterator.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
+#include <base/quadrature.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
-#include <grid/tria_iterator.h>
-#include <base/geometry_info.h>
-#include <base/quadrature.h>
 #include <fe/fe.h>
 #include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+#include <grid/tria_iterator.h>
+#include <grid/geometry_info.h>
 #include <hp/fe_values.h>
 #include <hp/mapping_collection.h>
 #include <numerics/matrices.h>
@@ -34,7 +34,6 @@
 #include <lac/block_vector.h>
 #include <lac/sparse_matrix.h>
 #include <lac/block_sparse_matrix.h>
-#include <fe/mapping_q1.h>
 
 #include <algorithm>
 #include <set>
@@ -56,10 +55,10 @@ namespace internal
       {
          Scratch (const FiniteElement<dim,spacedim> &fe,
                   const UpdateFlags         update_flags,
-                  const Function<spacedim>      *coefficient,
-                  const Function<spacedim>      *rhs_function,
+                  const Function<spacedim> *coefficient,
+                  const Function<spacedim> *rhs_function,
                   const Quadrature<dim>    &quadrature,
-                  const Mapping<dim,spacedim>       &mapping)
+                  const Mapping<dim,spacedim> &mapping)
                          :
                          fe_collection (fe),
                          quadrature_collection (quadrature),
@@ -81,10 +80,10 @@ namespace internal
 
          Scratch (const ::dealii::hp::FECollection<dim,spacedim> &fe,
                   const UpdateFlags         update_flags,
-                  const Function<spacedim>      *coefficient,
-                  const Function<spacedim>      *rhs_function,
-                  const ::dealii::hp::QCollection<dim>    &quadrature,
-                  const ::dealii::hp::MappingCollection<dim,spacedim>       &mapping)
+                  const Function<spacedim> *coefficient,
+                  const Function<spacedim> *rhs_function,
+                  const ::dealii::hp::QCollection<dim> &quadrature,
+                  const ::dealii::hp::MappingCollection<dim,spacedim> &mapping)
                          :
                          fe_collection (fe),
                          quadrature_collection (quadrature),
@@ -122,8 +121,8 @@ namespace internal
                          update_flags (data.update_flags)
            {}
        
-         const ::dealii::hp::FECollection<dim,spacedim> fe_collection;
-         const ::dealii::hp::QCollection<dim> quadrature_collection;
+         const ::dealii::hp::FECollection<dim,spacedim>      fe_collection;
+         const ::dealii::hp::QCollection<dim>                quadrature_collection;
          const ::dealii::hp::MappingCollection<dim,spacedim> mapping_collection;
        
          ::dealii::hp::FEValues<dim,spacedim> x_fe_values;
@@ -133,6 +132,8 @@ namespace internal
          std::vector<double>                  rhs_values;
          std::vector<dealii::Vector<double> > rhs_vector_values;
 
+         std::vector<double> old_JxW;
+
          const Function<spacedim>   *coefficient;
          const Function<spacedim>   *rhs_function;
 
@@ -164,7 +165,7 @@ namespace internal
       
       const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
                         n_q_points    = fe_values.n_quadrature_points;
-      const FiniteElement<dim,spacedim>    &fe  = fe_values.get_fe();
+      const FiniteElement<dim,spacedim> &fe  = fe_values.get_fe();
       const unsigned int n_components  = fe.n_components();
 
       Assert(data.rhs_function == 0 ||
@@ -177,15 +178,13 @@ namespace internal
             ::dealii::MatrixCreator::ExcComponentMismatch());
 
       copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
-      copy_data.cell_matrix = 0;
-
       copy_data.cell_rhs.reinit (dofs_per_cell);
-      copy_data.cell_rhs = 0;
       
       copy_data.dof_indices.resize (dofs_per_cell);
       cell->get_dof_indices (copy_data.dof_indices);
 
-      if (data.rhs_function != 0)
+      const bool use_rhs_function = data.rhs_function != 0;
+      if (use_rhs_function)
        {
          if (data.rhs_function->n_components==1)
            {
@@ -202,7 +201,8 @@ namespace internal
            }
        }
 
-      if (data.coefficient != 0)
+      const bool use_coefficient = data.coefficient != 0;
+      if (use_coefficient)
        {
          if (data.coefficient->n_components==1)
            {
@@ -220,175 +220,116 @@ namespace internal
        }
 
 
-      if (data.coefficient != 0)
-       {
-         if (data.coefficient->n_components==1)
-           {
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
+      double add_data;
+      const std::vector<double> & JxW = fe_values.get_JxW_values();
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       if (fe.is_primitive ())
+         {
+           const unsigned int component_i =
+             fe.system_to_component_index(i).first;
+           const double * phi_i = &fe_values.shape_value(i,0);
+           add_data = 0;
+
+                                  // use symmetry in the mass matrix here:
+                                  // just need to calculate the diagonal
+                                  // and half of the elements above the
+                                  // diagonal
+           for (unsigned int j=i; j<dofs_per_cell; ++j)
+             if ((n_components==1) ||
+                 (fe.system_to_component_index(j).first ==
+                  component_i))
                {
-                 const unsigned int component_i =
-                   fe.system_to_component_index(i).first;
-                 for (unsigned int j=0; j<dofs_per_cell; ++j)
-                   if ((n_components==1) ||
-                       (fe.system_to_component_index(j).first ==
-                        component_i))
-                     for (unsigned int point=0; point<n_q_points; ++point)
-                       copy_data.cell_matrix(i,j)
-                         += (fe_values.shape_value(i,point) * 
-                             fe_values.shape_value(j,point) * 
-                             fe_values.JxW(point) *
-                             data.coefficient_values[point]);
-
-                 if (data.rhs_function != 0)
+                 const double * phi_j = &fe_values.shape_value(j,0);
+                 add_data = 0;
+                 if (use_coefficient)
                    {
-                     if (data.rhs_function->n_components==1)
+                     if (data.coefficient->n_components==1)
                        for (unsigned int point=0; point<n_q_points; ++point)
-                         copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
-                           data.rhs_values[point] * fe_values.JxW(point);
+                         add_data += (phi_i[point] * phi_j[point] * JxW[point] * 
+                                      data.coefficient_values[point]);
                      else
                        for (unsigned int point=0; point<n_q_points; ++point)
-                         copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
-                           data.rhs_vector_values[point](component_i) * 
-                           fe_values.JxW(point);
+                         add_data += (phi_i[point] * phi_j[point] * JxW[point] * 
+                                      data.coefficient_vector_values[point](component_i));
                    }
+                 else
+                   for (unsigned int point=0; point<n_q_points; ++point)
+                     add_data += phi_i[point] * phi_j[point] * JxW[point]; 
+
+                                  // this is even ok for i==j, since then
+                                  // we just write the same value twice.
+                 copy_data.cell_matrix(i,j) = add_data;
+                 copy_data.cell_matrix(j,i) = add_data;
                }
-           }
-         else
-           {
-             if (fe.is_primitive ())
-               {
-                 for (unsigned int i=0; i<dofs_per_cell; ++i)
+
+           if (use_rhs_function)
+             {
+               add_data = 0;
+               if (data.rhs_function->n_components==1)
+                 for (unsigned int point=0; point<n_q_points; ++point)
+                   add_data += phi_i[point] * JxW[point] *
+                               data.rhs_values[point];
+               else
+                 for (unsigned int point=0; point<n_q_points; ++point)
+                   add_data += phi_i[point] * JxW[point] *
+                               data.rhs_vector_values[point](component_i);
+               copy_data.cell_rhs(i) = add_data;
+             }
+         }
+       else
+         {
+                                  // non-primitive vector-valued FE, using
+                                  // symmetry again
+           for (unsigned int j=i; j<dofs_per_cell; ++j)
+             {
+               add_data = 0;
+               for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+                 if (fe.get_nonzero_components(i)[comp_i] &&
+                     fe.get_nonzero_components(j)[comp_i])
                    {
-                     const unsigned int component_i =
-                       fe.system_to_component_index(i).first;
-                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       if ((n_components==1) ||
-                           (fe.system_to_component_index(j).first ==
-                            component_i))
-                         for (unsigned int point=0; point<n_q_points; ++point)
-                           copy_data.cell_matrix(i,j) +=
-                             (fe_values.shape_value(i,point) * 
-                              fe_values.shape_value(j,point) * 
-                              fe_values.JxW(point) *
-                              data.coefficient_vector_values[point](component_i));
-
-                     if (data.rhs_function != 0)
+                     if (use_coefficient)
                        {
-                         if (data.rhs_function->n_components==1)
+                         if (data.coefficient->n_components==1)
                            for (unsigned int point=0; point<n_q_points; ++point)
-                             copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
-                               data.rhs_values[point] * fe_values.JxW(point);
+                             add_data += (fe_values.shape_value_component(i,point,comp_i) *
+                                          fe_values.shape_value_component(j,point,comp_i) *
+                                          JxW[point] * 
+                                          data.coefficient_values[point]);
                          else
                            for (unsigned int point=0; point<n_q_points; ++point)
-                             copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
-                               data.rhs_vector_values[point](component_i) * 
-                               fe_values.JxW(point);
+                             add_data += (fe_values.shape_value_component(i,point,comp_i) *
+                                          fe_values.shape_value_component(j,point,comp_i) *
+                                          JxW[point] * 
+                                          data.coefficient_vector_values[point](comp_i));
                        }
+                     else
+                       for (unsigned int point=0; point<n_q_points; ++point)
+                         add_data += fe_values.shape_value_component(i,point,comp_i) *
+                           fe_values.shape_value_component(j,point,comp_i) * JxW[point];
                    }
-               }
-             else
-                                                // non-primitive vector-valued FE
-               {
-                 for (unsigned int i=0; i<dofs_per_cell; ++i)
-                   for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
-                     if (fe.get_nonzero_components(i)[comp_i])
-                       {
-                         for (unsigned int j=0; j<dofs_per_cell; ++j)
-                           for (unsigned int comp_j = 0; comp_j < n_components; ++comp_j)
-                             if (fe.get_nonzero_components(j)[comp_j])
-                               if (comp_i == comp_j)
-                                 for (unsigned int point=0; point<n_q_points; ++point)
-                                   copy_data.cell_matrix(i,j) += 
-                                     (fe_values.shape_value_component(i,point,comp_i) * 
-                                      fe_values.shape_value_component(j,point,comp_j) * 
-                                      fe_values.JxW(point) *
-                                      data.coefficient_vector_values[point](comp_i));
-                             
-                         if (data.rhs_function != 0)
-                           {
-                             if (data.rhs_function->n_components==1)
-                               for (unsigned int point=0; point<n_q_points; ++point)
-                                 copy_data.cell_rhs(i) += 
-                                   fe_values.shape_value_component(i,point,comp_i) *
-                                   data.rhs_values[point] * fe_values.JxW(point);
-                             else
-                               for (unsigned int point=0; point<n_q_points; ++point)
-                                 copy_data.cell_rhs(i) += 
-                                   fe_values.shape_value_component(i,point,comp_i) *
-                                   data.rhs_vector_values[point](comp_i) * 
-                                   fe_values.JxW(point);
-                           }
-                       }
-               }
-           }
-       }
-      else
-                                        // no coefficient
-       {
-         if (fe.is_primitive ())
-           {
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               {
-                 const unsigned int component_i =
-                   fe.system_to_component_index(i).first;
-                 for (unsigned int j=0; j<dofs_per_cell; ++j)
-                   if ((n_components==1) ||
-                       (fe.system_to_component_index(j).first ==
-                        component_i))
-                     for (unsigned int point=0; point<n_q_points; ++point)
-                       copy_data.cell_matrix(i,j) +=
-                         (fe_values.shape_value(i,point) * 
-                          fe_values.shape_value(j,point) * 
-                          fe_values.JxW(point));
-
-                 if (data.rhs_function != 0)
+
+               copy_data.cell_matrix(i,j) = add_data;
+               copy_data.cell_matrix(j,i) = add_data;
+             }
+             
+           if (use_rhs_function)
+             {
+               add_data = 0;
+               for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) 
+                 if (fe.get_nonzero_components(i)[comp_i])
                    {
                      if (data.rhs_function->n_components==1)
                        for (unsigned int point=0; point<n_q_points; ++point)
-                         copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
-                           data.rhs_values[point] * fe_values.JxW(point);
+                         add_data += fe_values.shape_value_component(i,point,comp_i) * 
+                                     JxW[point] * data.rhs_values[point];
                      else
                        for (unsigned int point=0; point<n_q_points; ++point)
-                         copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
-                           data.rhs_vector_values[point](component_i) * 
-                           fe_values.JxW(point);
-                   }
-               }
-           }
-         else
-                                            // non-primitive FE, no coefficient
-           {
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
-                 if (fe.get_nonzero_components(i)[comp_i])
-                   {
-                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       for (unsigned int comp_j = 0; comp_j < n_components; ++comp_j)
-                         if (fe.get_nonzero_components(j)[comp_j])
-                           if (comp_i == comp_j)
-                             for (unsigned int point=0; point<n_q_points; ++point)
-                               copy_data.cell_matrix(i,j) += 
-                                 (fe_values.shape_value_component(i,point,comp_i) * 
-                                  fe_values.shape_value_component(j,point,comp_j) * 
-                                  fe_values.JxW(point));
-
-                     if (data.rhs_function != 0)
-                       {
-                         if (data.rhs_function->n_components==1)
-                           for (unsigned int point=0; point<n_q_points; ++point)
-                             copy_data.cell_rhs(i) += 
-                               fe_values.shape_value_component(i,point,comp_i) *
-                               data.rhs_values[point] * fe_values.JxW(point);
-                         else
-                           for (unsigned int point=0; point<n_q_points; ++point)
-                             copy_data.cell_rhs(i) += 
-                               fe_values.shape_value_component(i,point,comp_i) *
-                               data.rhs_vector_values[point](comp_i) * 
-                               fe_values.JxW(point);
-                       }
+                         add_data += fe_values.shape_value_component(i,point,comp_i) * 
+                                     JxW[point] * data.rhs_vector_values[point](comp_i);
                    }
-           }
-       }
+               copy_data.cell_rhs(i) = add_data;
+             }
+         }
     }
 
 
@@ -418,16 +359,13 @@ namespace internal
             ::dealii::MatrixCreator::ExcComponentMismatch());
 
       copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
-      copy_data.cell_matrix = 0;
-
       copy_data.cell_rhs.reinit (dofs_per_cell);
-      copy_data.cell_rhs = 0;
-      
       copy_data.dof_indices.resize (dofs_per_cell);
       cell->get_dof_indices (copy_data.dof_indices);
 
 
-      if (data.rhs_function != 0)
+      const bool use_rhs_function = data.rhs_function != 0;
+      if (use_rhs_function)
        {
          if (data.rhs_function->n_components==1)
            {
@@ -444,7 +382,8 @@ namespace internal
            }
        }
 
-      if (data.coefficient != 0)
+      const bool use_coefficient = data.coefficient != 0;
+      if (use_coefficient)
        {
          if (data.coefficient->n_components==1)
            {
@@ -462,179 +401,120 @@ namespace internal
        }
 
       
-      if (data.coefficient != 0)
-       {
-         if (data.coefficient->n_components==1)
-           {
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
+      const std::vector<double> & JxW = fe_values.get_JxW_values();
+      double add_data;
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       if (fe.is_primitive ())
+         {
+           const unsigned int component_i =
+             fe.system_to_component_index(i).first;
+           const Tensor<1,spacedim> * grad_phi_i = 
+             &fe_values.shape_grad(i,0);
+
+                                  // can use symmetry
+           for (unsigned int j=i; j<dofs_per_cell; ++j)
+             if ((n_components==1) ||
+                 (fe.system_to_component_index(j).first ==
+                  component_i))
                {
-                 const unsigned int component_i =
-                   fe.system_to_component_index(i).first;
-                 for (unsigned int j=0; j<dofs_per_cell; ++j)
-                   if ((n_components==1) ||
-                       (fe.system_to_component_index(j).first ==
-                        component_i))
-                     for (unsigned int point=0; point<n_q_points; ++point)
-                       copy_data.cell_matrix(i,j)
-                         += (fe_values.shape_grad(i,point) * 
-                             fe_values.shape_grad(j,point) * 
-                             fe_values.JxW(point) *
-                             data.coefficient_values[point]);
-
-                 if (data.rhs_function != 0)
+                 const Tensor<1,spacedim> * grad_phi_j = 
+                   & fe_values.shape_grad(j,0);
+                 add_data = 0;
+                 if (use_coefficient)
                    {
-                     if (data.rhs_function->n_components==1)
+                     if (data.coefficient->n_components==1)
                        for (unsigned int point=0; point<n_q_points; ++point)
-                         copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
-                           data.rhs_values[point] * fe_values.JxW(point);
+                         add_data += ((grad_phi_i[point]*grad_phi_j[point]) * 
+                                      JxW[point] * 
+                                      data.coefficient_values[point]);
                      else
                        for (unsigned int point=0; point<n_q_points; ++point)
-                         copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
-                           data.rhs_vector_values[point](component_i) * 
-                           fe_values.JxW(point);
+                         add_data += ((grad_phi_i[point]*grad_phi_j[point]) * 
+                                      JxW[point] * 
+                                      data.coefficient_vector_values[point](component_i));
                    }
+                 else
+                   for (unsigned int point=0; point<n_q_points; ++point)
+                     add_data += (grad_phi_i[point]*grad_phi_j[point]) * 
+                                 JxW[point]; 
+
+                 copy_data.cell_matrix(i,j) = add_data;
+                 copy_data.cell_matrix(j,i) = add_data;
                }
-           }
-         else
-           {
-             if (fe.is_primitive ())
-               {
-                 for (unsigned int i=0; i<dofs_per_cell; ++i)
+
+           if (use_rhs_function)
+             {
+               const double * phi_i = &fe_values.shape_value(i,0);
+               add_data = 0;
+               if (data.rhs_function->n_components==1)
+                 for (unsigned int point=0; point<n_q_points; ++point)
+                   add_data += phi_i[point] * JxW[point] *
+                               data.rhs_values[point];
+               else
+                 for (unsigned int point=0; point<n_q_points; ++point)
+                   add_data += phi_i[point] * JxW[point] *
+                               data.rhs_vector_values[point](component_i);
+               copy_data.cell_rhs(i) = add_data;
+             }
+         }
+       else
+         {
+                                  // non-primitive vector-valued FE
+           for (unsigned int j=i; j<dofs_per_cell; ++j)
+             {
+               add_data = 0;
+               for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+                 if (fe.get_nonzero_components(i)[comp_i] &&
+                     fe.get_nonzero_components(j)[comp_i])
                    {
-                     const unsigned int component_i =
-                       fe.system_to_component_index(i).first;
-                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       if ((n_components==1) ||
-                           (fe.system_to_component_index(j).first ==
-                            component_i))
-                         for (unsigned int point=0; point<n_q_points; ++point)
-                           copy_data.cell_matrix(i,j) +=
-                             (fe_values.shape_grad(i,point) * 
-                              fe_values.shape_grad(j,point) * 
-                              fe_values.JxW(point) *
-                              data.coefficient_vector_values[point](component_i));
-
-                     if (data.rhs_function != 0)
+                     if (use_coefficient)
                        {
-                         if (data.rhs_function->n_components==1)
+                         if (data.coefficient->n_components==1)
                            for (unsigned int point=0; point<n_q_points; ++point)
-                             copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
-                               data.rhs_values[point] * fe_values.JxW(point);
+                             add_data += ((fe_values.shape_grad_component(i,point,comp_i) * 
+                                           fe_values.shape_grad_component(j,point,comp_i)) * 
+                                          JxW[point] * 
+                                          data.coefficient_values[point]);
                          else
                            for (unsigned int point=0; point<n_q_points; ++point)
-                             copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
-                               data.rhs_vector_values[point](component_i) * 
-                               fe_values.JxW(point);
+                             add_data += ((fe_values.shape_grad_component(i,point,comp_i) * 
+                                           fe_values.shape_grad_component(j,point,comp_i)) * 
+                                          JxW[point] * 
+                                          data.coefficient_vector_values[point](comp_i));
                        }
+                     else
+                       for (unsigned int point=0; point<n_q_points; ++point)
+                         add_data += (fe_values.shape_grad_component(i,point,comp_i) * 
+                                      fe_values.shape_grad_component(j,point,comp_i)) * 
+                                     JxW[point];
                    }
-               }
-             else
-                                                // non-primitive vector-valued FE
-               {
-                 for (unsigned int i=0; i<dofs_per_cell; ++i)
-                   for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
-                     if (fe.get_nonzero_components(i)[comp_i])
-                       {
-                         for (unsigned int j=0; j<dofs_per_cell; ++j)
-                           for (unsigned int comp_j = 0; comp_j < n_components; ++comp_j)
-                             if (fe.get_nonzero_components(j)[comp_j])
-                               if (comp_i == comp_j)
-                                 for (unsigned int point=0; point<n_q_points; ++point)
-                                   copy_data.cell_matrix(i,j) += 
-                                     (fe_values.shape_grad_component(i,point,comp_i) * 
-                                      fe_values.shape_grad_component(j,point,comp_j) * 
-                                      fe_values.JxW(point) *
-                                      data.coefficient_vector_values[point](comp_i));
-                             
-                         if (data.rhs_function != 0)
-                           {
-                             if (data.rhs_function->n_components==1)
-                               for (unsigned int point=0; point<n_q_points; ++point)
-                                 copy_data.cell_rhs(i) += 
-                                   fe_values.shape_value_component(i,point,comp_i) *
-                                   data.rhs_values[point] * fe_values.JxW(point);
-                             else
-                               for (unsigned int point=0; point<n_q_points; ++point)
-                                 copy_data.cell_rhs(i) += 
-                                   fe_values.shape_value_component(i,point,comp_i) *
-                                   data.rhs_vector_values[point](comp_i) * 
-                                   fe_values.JxW(point);
-                           }
-                       }
-               }
-           }
-       }
-      else
-                                        // no coefficient
-       {
-         if (fe.is_primitive ())
-           {
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               {
-                 const unsigned int component_i =
-                   fe.system_to_component_index(i).first;
-                 for (unsigned int j=0; j<dofs_per_cell; ++j)
-                   if ((n_components==1) ||
-                       (fe.system_to_component_index(j).first ==
-                        component_i))
-                     for (unsigned int point=0; point<n_q_points; ++point)
-                       copy_data.cell_matrix(i,j) +=
-                         (fe_values.shape_grad(i,point) * 
-                          fe_values.shape_grad(j,point) * 
-                          fe_values.JxW(point));
-
-                 if (data.rhs_function != 0)
+
+               copy_data.cell_matrix(i,j) = add_data;
+               copy_data.cell_matrix(j,i) = add_data;
+             }
+             
+           if (use_rhs_function)
+             {
+               add_data = 0;
+               for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i) 
+                 if (fe.get_nonzero_components(i)[comp_i])
                    {
                      if (data.rhs_function->n_components==1)
                        for (unsigned int point=0; point<n_q_points; ++point)
-                         copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
-                           data.rhs_values[point] * fe_values.JxW(point);
+                         add_data += fe_values.shape_value_component(i,point,comp_i) * 
+                                     JxW[point] * data.rhs_values[point];
                      else
                        for (unsigned int point=0; point<n_q_points; ++point)
-                         copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
-                           data.rhs_vector_values[point](component_i) * 
-                           fe_values.JxW(point);
-                   }
-               }
-           }
-         else
-                                            // non-primitive FE, no coefficient
-           {
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
-                 if (fe.get_nonzero_components(i)[comp_i])
-                   {
-                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       for (unsigned int comp_j = 0; comp_j < n_components; ++comp_j)
-                         if (fe.get_nonzero_components(j)[comp_j])
-                           if (comp_i == comp_j)
-                             for (unsigned int point=0; point<n_q_points; ++point)
-                               copy_data.cell_matrix(i,j) += 
-                                 (fe_values.shape_grad_component(i,point,comp_i) * 
-                                  fe_values.shape_grad_component(j,point,comp_j) * 
-                                  fe_values.JxW(point));
-
-                     if (data.rhs_function != 0)
-                       {
-                         if (data.rhs_function->n_components==1)
-                           for (unsigned int point=0; point<n_q_points; ++point)
-                             copy_data.cell_rhs(i) += 
-                               fe_values.shape_value_component(i,point,comp_i) *
-                               data.rhs_values[point] * fe_values.JxW(point);
-                         else
-                           for (unsigned int point=0; point<n_q_points; ++point)
-                             copy_data.cell_rhs(i) += 
-                               fe_values.shape_value_component(i,point,comp_i) *
-                               data.rhs_vector_values[point](comp_i) * 
-                               fe_values.JxW(point);
-                       }
+                         add_data += fe_values.shape_value_component(i,point,comp_i) * 
+                                     JxW[point] * data.rhs_vector_values[point](comp_i);
                    }
-           }
-       }
+               copy_data.cell_rhs(i) = add_data;
+             }
+         }
     }
-    
 
-    
+
+
     template <typename MatrixType,
              typename VectorType>
     void copy_local_to_global (const AssemblerData::CopyData &data,
@@ -652,8 +532,7 @@ namespace internal
              (data.cell_rhs.size() == dofs_per_cell),
              ExcInternalError());
 
-      matrix->add (data.dof_indices, data.cell_matrix);
-
+      matrix->add(data.dof_indices, data.cell_matrix);
       if (right_hand_side != 0)
        for (unsigned int i=0; i<dofs_per_cell; ++i)
          (*right_hand_side)(data.dof_indices[i]) += data.cell_rhs(i);

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.