]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Avoid some unnecessary work in FE_Q setup. Give compiler some hint in FEPolyTensor.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 28 May 2009 09:57:02 +0000 (09:57 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 28 May 2009 09:57:02 +0000 (09:57 +0000)
git-svn-id: https://svn.dealii.org/trunk@18885 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/fe_poly_tensor.cc
deal.II/deal.II/source/fe/fe_q.cc

index a6cc1f226b9135bfbb14dc72d699613f456b0757..9abae0427ac3f379ab9cb55263ac1976fe66196d 100644 (file)
@@ -330,8 +330,12 @@ FE_PolyTensor<POLY,dim,spacedim>::get_data (
                data->shape_values[i][k] = values[i];
            else
              for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-               for (unsigned int j=0; j<this->dofs_per_cell; ++j)
-                 data->shape_values[i][k] += inverse_node_matrix(j,i) * values[j];
+               {
+                 Tensor<1,dim> add_values;
+                 for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+                   add_values += inverse_node_matrix(j,i) * values[j];
+                 data->shape_values[i][k] = add_values;
+               }
          }
        
        if (flags & update_gradients)
@@ -341,8 +345,12 @@ FE_PolyTensor<POLY,dim,spacedim>::get_data (
                data->shape_grads[i][k] = grads[i];
            else
              for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-               for (unsigned int j=0; j<this->dofs_per_cell; ++j)
-                 data->shape_grads[i][k] += inverse_node_matrix(j,i) * grads[j];
+               {
+                 Tensor<2,dim> add_grads;
+                 for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+                   add_grads += inverse_node_matrix(j,i) * grads[j];
+                 data->shape_grads[i][k] = add_grads;
+               }
          }
        
       }
index e3fa057a36ee5723ebd03b0a2808ef86ceac5457..8b27f0e2f0802d2376f799b29b466b863e25391b 100644 (file)
@@ -1530,6 +1530,60 @@ FE_Q<dim,spacedim>::initialize_embedding ()
   const std::vector<unsigned int> &index_map=
     this->poly_space.get_numbering();
 
+  const double zero_threshold = 2e-13*this->degree*this->degree*dim;
+
+                                  // precompute subcell interpolation
+                                  // matrix
+  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+    for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+      {
+       const Point<dim> p_subcell
+         = FE_Q_Helper::generate_unit_point (index_map[j], this->dofs_per_cell,
+                                             dealii::internal::int2type<dim>());
+       const double
+         subcell_value = this->poly_space.compute_value (i, p_subcell);
+
+                                                // cut off values that are
+                                                // too small. note that we
+                                                // have here Lagrange
+                                                // interpolation functions,
+                                                // so they should be zero
+                                                // at almost all points,
+                                                // and one at the others,
+                                                // at least on the
+                                                // subcells. so set them to
+                                                // their exact values
+                                                //
+                                                // the actual cut-off value
+                                                // is somewhat fuzzy, but
+                                                // it works for
+                                                // 2e-13*degree^2*dim (see
+                                                // above), which is kind of
+                                                // reasonable given that we
+                                                // compute the values of
+                                                // the polynomials via an
+                                                // degree-step recursion
+                                                // and then multiply the
+                                                // 1d-values. this gives us
+                                                // a linear growth in
+                                                // degree*dim, times a
+                                                // small constant.
+       if (std::fabs(subcell_value) < zero_threshold)
+         subcell_interpolation(j, i) = 0.;
+       else if (std::fabs(subcell_value-1) < zero_threshold)
+         subcell_interpolation(j, i) = 1.;
+       else                    
+                                                    // we have put our
+                                                    // evaluation
+                                                    // points onto the
+                                                    // interpolation
+                                                    // points, so we
+                                                    // should either
+                                                    // get zeros or
+                                                    // ones!
+         Assert (false, ExcInternalError());
+      }
+
   for (unsigned int ref=0; ref<RefinementCase<dim>::isotropic_refinement; ++ref)
     for (unsigned int child=0; child<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref+1)); ++child)
       {
@@ -1543,66 +1597,43 @@ FE_Q<dim,spacedim>::initialize_embedding ()
              = FE_Q_Helper::generate_unit_point (index_map[j], this->dofs_per_cell,
                                                  dealii::internal::int2type<dim>());
            const Point<dim> p_cell =
-             GeometryInfo<dim>::child_to_cell_coordinates (p_subcell, child, RefinementCase<dim>(ref+1));
+             GeometryInfo<dim>::child_to_cell_coordinates (p_subcell, child, 
+                                                           RefinementCase<dim>(ref+1));
 
            for (unsigned int i=0; i<this->dofs_per_cell; ++i)
              {
                const double
-                 cell_value    = this->poly_space.compute_value (i, p_cell),
-                 subcell_value = this->poly_space.compute_value (i, p_subcell);
-
-                                                // cut off values that
-                                                // are too small. note
-                                                // that we have here
-                                                // Lagrange
-                                                // interpolation
-                                                // functions, so they
-                                                // should be zero at
-                                                // almost all points,
-                                                // and one at the
-                                                // others, at least on
-                                                // the subcells. so set
-                                                // them to their exact
-                                                // values
+                 cell_value    = this->poly_space.compute_value (i, p_cell);
+
+                                                // cut off values that are
+                                                // too small. note that we
+                                                // have here Lagrange
+                                                // interpolation functions,
+                                                // so they should be zero
+                                                // at almost all points,
+                                                // and one at the others,
+                                                // at least on the
+                                                // subcells. so set them to
+                                                // their exact values
                                                 //
-                                                // the actual cut-off
-                                                // value is somewhat
-                                                // fuzzy, but it works
-                                                // for
-                                                // 1e-14*degree*dim,
-                                                // which is kind of
-                                                // reasonable given
-                                                // that we compute the
-                                                // values of the
-                                                // polynomials via an
-                                                // degree-step
-                                                // recursion and then
-                                                // multiply the
-                                                // 1d-values. this
-                                                // gives us a linear
-                                                // growth in
+                                                // the actual cut-off value
+                                                // is somewhat fuzzy, but
+                                                // it works for
+                                                // 2e-13*degree^2*dim (see
+                                                // above), which is kind of
+                                                // reasonable given that we
+                                                // compute the values of
+                                                // the polynomials via an
+                                                // degree-step recursion
+                                                // and then multiply the
+                                                // 1d-values. this gives us
+                                                // a linear growth in
                                                 // degree*dim, times a
                                                 // small constant.
-               if (std::fabs(cell_value) < 2e-13*this->degree*this->degree*dim)
+               if (std::fabs(cell_value) < zero_threshold)
                  cell_interpolation(j, i) = 0.;
                else
                  cell_interpolation(j, i) = cell_value;
-
-               if (std::fabs(subcell_value) < 2e-13*this->degree*this->degree*dim)
-                 subcell_interpolation(j, i) = 0.;
-               else
-                 if (std::fabs(subcell_value-1) < 2e-13*this->degree*this->degree*dim)
-                   subcell_interpolation(j, i) = 1.;
-                 else                  
-                                                    // we have put our
-                                                    // evaluation
-                                                    // points onto the
-                                                    // interpolation
-                                                    // points, so we
-                                                    // should either
-                                                    // get zeros or
-                                                    // ones!
-                   Assert (false, ExcInternalError());
              }
          }
 
@@ -1631,7 +1662,7 @@ FE_Q<dim,spacedim>::initialize_embedding ()
                                         // here
        for (unsigned int i=0; i<this->dofs_per_cell; ++i)
          for (unsigned int j=0; j<this->dofs_per_cell; ++j)
-           if (std::fabs(this->prolongation[ref][child](i,j)) < 2e-13*this->degree*dim)
+           if (std::fabs(this->prolongation[ref][child](i,j)) < zero_threshold)
              this->prolongation[ref][child](i,j) = 0.;
 
                                         // and make sure that the row
@@ -1643,7 +1674,7 @@ FE_Q<dim,spacedim>::initialize_embedding ()
            double sum = 0;
            for (unsigned int col=0; col<this->dofs_per_cell; ++col)
              sum += this->prolongation[ref][child](row,col);
-           Assert (std::fabs(sum-1.) < 2e-13*this->degree*this->degree*dim,
+           Assert (std::fabs(sum-1.) < zero_threshold,
                    ExcInternalError());
          }
       }

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.