[cig-commits] [commit] master: more fixes (a13708f)

[cig_noreply at geodynamics.org]

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

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

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

commit a13708f5cc51a0c1325f6f6c7b835430f7bdd50d
Author: Timo Heister <timo.heister at gmail.com>
Date:   Mon May 19 11:22:35 2014 -0400

    more fixes


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

a13708f5cc51a0c1325f6f6c7b835430f7bdd50d
 source/simulator/core.cc               | 6 +++---
 source/simulator/initial_conditions.cc | 9 ++++++---
 2 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index 7a83cf6..37420bf 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -1308,9 +1308,9 @@ namespace aspect
 
               current_linearization_point.block(introspection.block_indices.velocities)
                 = solution.block(introspection.block_indices.velocities);
-              current_linearization_point.block(introspection.block_indices.pressure)
-                = solution.block(introspection.block_indices.pressure);
-
+              if (introspection.block_indices.velocities != introspection.block_indices.pressure)
+                current_linearization_point.block(introspection.block_indices.pressure)
+                  = solution.block(introspection.block_indices.pressure);
 
               pcout << std::endl;
             }
diff --git a/source/simulator/initial_conditions.cc b/source/simulator/initial_conditions.cc
index ba2b699..05bfa60 100644
--- a/source/simulator/initial_conditions.cc
+++ b/source/simulator/initial_conditions.cc
@@ -180,6 +180,9 @@ namespace aspect
   template <int dim>
   void Simulator<dim>::compute_initial_pressure_field ()
   {
+    // Note that this code will overwrite the velocity solution with 0 if
+    // velocity and pressure are in the same block.
+
     // we'd like to interpolate the initial pressure onto the pressure
     // variable but that's a bit involved because the pressure may either
     // be an FE_Q (for which we can interpolate) or an FE_DGP (for which
@@ -207,7 +210,7 @@ namespace aspect
                                                                                std_cxx1x::cref (*adiabatic_conditions),
                                                                                std_cxx1x::_1),
                                                                                dim,
-                                                                               dim+2+parameters.n_compositional_fields),
+                                                                               introspection.n_components),
                                   system_tmp);
 
         // we may have hanging nodes, so apply constraints
@@ -281,9 +284,9 @@ namespace aspect
                     // for all other variables so that the whole thing remains
                     // invertible
                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                      if ((finite_element.system_to_component_index(i).first == dim)
+                      if ((finite_element.system_to_component_index(i).first == introspection.component_indices.pressure)
                           &&
-                          (finite_element.system_to_component_index(j).first == dim))
+                          (finite_element.system_to_component_index(j).first == introspection.component_indices.pressure))
                         local_mass_matrix(j,i) += (fe_values[introspection.extractors.pressure].value(i,point) *
                                                    fe_values[introspection.extractors.pressure].value(j,point) *
                                                    fe_values.JxW(point));