[DO NOT MERGE] use current_constraints for SparsityPattern

doc/modules/changes/20170801_tjhei

@@ -0,0 +1,5 @@
+New: Sparsity patterns of the system matrix now uses the current_constraints
+instead of the static constraints. This allows new kind of constraints in the
+system and reduces the number of nonzero entries in the system matrix.
+<br>
+(Timo Heister, 2017/08/01)

include/aspect/simulator.h

@@ -1430,6 +1430,7 @@ namespace aspect
 //TODO: use n_compositional_field separate preconditioners
       std_cxx11::shared_ptr<LinearAlgebra::PreconditionILU>     C_preconditioner;
 
+      bool                                                      rebuild_sparsity_and_matrices;
       bool                                                      rebuild_stokes_matrix;
       bool                                                      assemble_newton_stokes_matrix;
       bool                                                      assemble_newton_stokes_system;

include/aspect/simulator_access.h

@@ -657,6 +657,13 @@ namespace aspect
        */
       bool simulator_is_initialized () const;
 
+      /**
+       * Return the value used for rescaling the pressure in the linear
+       * solver.
+       */
+      double
+      get_pressure_scaling () const;
+
       /**
        * A convenience function that copies the values of the compositional
        * fields at the quadrature point q given as input parameter to the

source/simulator/core.cc

@@ -807,6 +807,16 @@ namespace aspect
     // into current_constraints and then add to current_constraints
     compute_current_constraints ();
 
+    // If needed, construct sparsity patterns and matrices with the current
+    // constraints. Of course we need to force assembly too.
+    if (rebuild_sparsity_and_matrices)
+      {
+        rebuild_sparsity_and_matrices = false;
+        setup_system_matrix (introspection.index_sets.system_partitioning);
+        setup_system_preconditioner (introspection.index_sets.system_partitioning);
+        rebuild_stokes_matrix = rebuild_stokes_preconditioner = true;
+      }
+
 
     // TODO: do this in a more efficient way (TH)? we really only need
     // to make sure that the time dependent velocity boundary conditions
@@ -1398,10 +1408,9 @@ namespace aspect
     constraints.close();
     signals.post_compute_no_normal_flux_constraints(triangulation);
 
-    // finally initialize vectors, matrices, etc.
-
-    setup_system_matrix (introspection.index_sets.system_partitioning);
-    setup_system_preconditioner (introspection.index_sets.system_partitioning);
+    // Finally initialize vectors. We delay construction of the sparsity
+    // patterns and matrices until we have current_constraints.
+    rebuild_sparsity_and_matrices = true;
 
     system_rhs.reinit(introspection.index_sets.system_partitioning, mpi_communicator);
     solution.reinit(introspection.index_sets.system_partitioning, introspection.index_sets.system_relevant_partitioning, mpi_communicator);

source/simulator/simulator_access.cc

@@ -596,6 +596,12 @@ namespace aspect
     return (simulator != NULL);
   }
 
+  template <int dim>
+  double
+  SimulatorAccess<dim>::get_pressure_scaling () const
+  {
+    return (simulator->pressure_scaling);
+  }
 }