[cig-commits] commit 1931 by heister to /var/svn/dealii/aspect

[dealii.demon at gmail.com]

Revision 1931

improve matrix assembly performance

U   trunk/aspect/doc/modules/changes.h
U   trunk/aspect/source/simulator/assembly.cc


http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1931&peg=1931

Diff:
Modified: trunk/aspect/doc/modules/changes.h
===================================================================
--- trunk/aspect/doc/modules/changes.h	2013-10-01 17:20:34 UTC (rev 1930)
+++ trunk/aspect/doc/modules/changes.h	2013-10-01 21:24:41 UTC (rev 1931)
@@ -8,6 +8,11 @@
 </p>
 
 <ol>
+  <li>Fixed: performance of matrix assembly has been improved significantly,
+  especially in 3d.
+  <br>
+  (Timo Heister, 2013/10/01)
+
   <li>New: HDF5/XDMF will only output mesh data when the mesh changes,
   reducing total data output significantly. XDMF serialization is also
   properly implemented.

Modified: trunk/aspect/source/simulator/assembly.cc
===================================================================
--- trunk/aspect/source/simulator/assembly.cc	2013-10-01 17:20:34 UTC (rev 1930)
+++ trunk/aspect/source/simulator/assembly.cc	2013-10-01 21:24:41 UTC (rev 1931)
@@ -1151,6 +1151,14 @@
     const unsigned int dofs_per_cell = scratch.finite_element_values.get_fe().dofs_per_cell;
     const unsigned int n_q_points    = scratch.finite_element_values.n_quadrature_points;
 
+    const unsigned int solution_component
+      = (temperature_or_composition.is_temperature()
+        ?
+        introspection.component_indices.temperature
+        :
+        introspection.component_indices.compositional_fields[temperature_or_composition.compositional_variable]
+       );
+
     const FEValuesExtractors::Scalar solution_field
       = (temperature_or_composition.is_temperature()
          ?
@@ -1254,6 +1262,10 @@
     for (unsigned int q=0; q<n_q_points; ++q)
       {
         for (unsigned int k=0; k<dofs_per_cell; ++k)
+          // We only need to look up values of shape functions if they
+          // belong to 'our' component. They are zero otherwise anyway.
+          // Note that we later only look at the values that we do set here.
+          if (cell->get_fe().system_to_component_index(k).first == solution_component)
           {
             scratch.grad_phi_field[k] = scratch.finite_element_values[solution_field].gradient (k,q);
             scratch.phi_field[k]      = scratch.finite_element_values[solution_field].value (k, q);
@@ -1297,6 +1309,7 @@
                                                   (time_step + old_time_step)) : 1.0;
 
         for (unsigned int i=0; i<dofs_per_cell; ++i)
+          if (cell->get_fe().system_to_component_index(i).first == solution_component)
           {
             data.local_rhs(i) += (field_term_for_rhs * scratch.phi_field[i]
                                   + time_step *
@@ -1306,6 +1319,7 @@
                                  scratch.finite_element_values.JxW(q);
 
             for (unsigned int j=0; j<dofs_per_cell; ++j)
+              if (cell->get_fe().system_to_component_index(j).first == solution_component)
               {
                 data.local_matrix(i,j)
                 += (