[cig-commits] [commit] master: fix block indices in several places (95c9804)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/aspect

On branch  : master
Link       : https://github.com/geodynamics/aspect/compare/76d4275352ef2cae5de9a073acd1c03a92c2670c...4f3d06fd1f3754419813db37ec9ef7f0f6f3cb15

>---------------------------------------------------------------

commit 95c9804ddcbe1e30cd3b8dae43ae96e75f2ebae7
Author: Timo Heister <timo.heister at gmail.com>
Date:   Mon May 12 17:03:05 2014 -0400

    fix block indices in several places


>---------------------------------------------------------------

95c9804ddcbe1e30cd3b8dae43ae96e75f2ebae7
 source/simulator/core.cc          | 30 +++++++++++++++++++++---------
 source/simulator/introspection.cc | 34 ++++++++++++++++++++++++----------
 2 files changed, 45 insertions(+), 19 deletions(-)

diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index 298a075..d7e30fb 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -474,14 +474,16 @@ namespace aspect
     statistics.add_value("Number of mesh cells",
                          triangulation.n_global_active_cells());
 
-    statistics.add_value("Number of Stokes degrees of freedom",
-                         introspection.system_dofs_per_block[0] +
-                         introspection.system_dofs_per_block[1]);
+    unsigned int n_stokes_dofs = introspection.system_dofs_per_block[0];
+    if (introspection.block_indices.velocities != introspection.block_indices.pressure)
+      n_stokes_dofs += introspection.system_dofs_per_block[introspection.block_indices.pressure];
+
+    statistics.add_value("Number of Stokes degrees of freedom", n_stokes_dofs);
     statistics.add_value("Number of temperature degrees of freedom",
-                         introspection.system_dofs_per_block[2]);
+                         introspection.system_dofs_per_block[introspection.block_indices.temperature]);
     if (parameters.n_compositional_fields > 0)
       statistics.add_value("Number of composition degrees of freedom",
-                           introspection.system_dofs_per_block[3]);
+                           introspection.system_dofs_per_block[introspection.block_indices.compositional_fields[0]]);
 
 
     // then interpolate the current boundary velocities. this adds to
@@ -887,10 +889,19 @@ namespace aspect
                                     introspection.system_dofs_per_block,
                                     introspection.components_to_blocks);
     {
-      const types::global_dof_index n_u = introspection.system_dofs_per_block[0],
-                                    n_p = introspection.system_dofs_per_block[1],
-                                    n_T = introspection.system_dofs_per_block[2];
-      std::vector<types::global_dof_index> n_C (parameters.n_compositional_fields+1);
+      types::global_dof_index n_u = introspection.system_dofs_per_block[0],
+          n_p = introspection.system_dofs_per_block[introspection.block_indices.pressure],
+          n_T = introspection.system_dofs_per_block[introspection.block_indices.temperature];
+
+      Assert(!parameters.direct_stokes_solver ||
+          (introspection.block_indices.velocities == introspection.block_indices.pressure),
+          ExcInternalError());
+
+      // only count pressure once if they are in the same block
+      if (parameters.direct_stokes_solver)
+        n_p = 0;
+
+      std::vector<types::global_dof_index> n_C (parameters.n_compositional_fields);
       for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
         n_C[c] = introspection.system_dofs_per_block
                  [introspection.block_indices.compositional_fields[c]];
@@ -908,6 +919,7 @@ namespace aspect
           introspection.index_sets.system_partitioning.push_back(system_index_set.get_view(n_u,n_u+n_p));
         }
       introspection.index_sets.system_partitioning.push_back(system_index_set.get_view(n_u+n_p,n_u+n_p+n_T));
+
       introspection.index_sets.stokes_partitioning.clear ();
       if (parameters.direct_stokes_solver)
         {
diff --git a/source/simulator/introspection.cc b/source/simulator/introspection.cc
index b9aeaad..6d52362 100644
--- a/source/simulator/introspection.cc
+++ b/source/simulator/introspection.cc
@@ -58,16 +58,28 @@ namespace aspect
   {
     template <int dim>
     std::vector<unsigned int>
-    component_to_block_mapping (const unsigned int n_components)
+    component_to_block_mapping (const unsigned int n_components,
+        const bool split_vel_pressure)
     {
       // set up a mapping between vector components to the blocks they
-      // correspond to. each variable has its own block except
-      // for the velocities which are all mapped into block 0
+      // correspond to.
       std::vector<unsigned int> components_to_blocks (n_components, 0U);
-      components_to_blocks[dim] = 1;
-      components_to_blocks[dim+1] = 2;
-      for (unsigned int i=dim+2; i<n_components; ++i)
-        components_to_blocks[i] = i-dim+1;
+      if (split_vel_pressure)
+        {
+          // each variable has its own block except
+          // for the velocities which are all mapped into block 0
+          components_to_blocks[dim] = 1;
+          components_to_blocks[dim+1] = 2;
+          for (unsigned int i=dim+2; i<n_components; ++i)
+            components_to_blocks[i] = i-dim+1;
+        }
+      else
+        {
+          // here velocity and pressure is block 0:
+          components_to_blocks[dim+1] = 1;
+          for (unsigned int i=dim+2; i<n_components; ++i)
+            components_to_blocks[i] = i-dim;
+        }
 
       return components_to_blocks;
     }
@@ -79,12 +91,12 @@ namespace aspect
       const bool split_vel_pressure)
     :
     n_components (dim+2+n_compositional_fields),
-    n_blocks (3+n_compositional_fields),
+    n_blocks (((split_vel_pressure)?3:2)+n_compositional_fields),
     extractors (n_compositional_fields),
     component_indices (n_compositional_fields),
     block_indices (n_compositional_fields, split_vel_pressure),
     base_elements (n_compositional_fields),
-    components_to_blocks (component_to_block_mapping<dim>(n_components)),
+    components_to_blocks (component_to_block_mapping<dim>(n_components, split_vel_pressure)),
     system_dofs_per_block (n_blocks)
   {}
 
@@ -119,7 +131,9 @@ namespace aspect
     velocities(0),
     pressure (split_vel_pressure?1:0),
     temperature (split_vel_pressure?2:1),
-    compositional_fields (half_open_sequence(3, 3+n_compositional_fields))
+    compositional_fields (half_open_sequence(
+        (split_vel_pressure?3:2),
+        (split_vel_pressure?3:2)+n_compositional_fields))
   {}