[cig-commits] [commit] master: implement denormalize_pressure() for direct solver (28bdda2)

[cig_noreply at geodynamics.org]

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

On branch  : master
Link       : https://github.com/geodynamics/aspect/compare/75a1eead49c5165ff3ecee1573c3f4577ef534bb...62599834f9944b32d7d25d117a9a394c565642ea

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

commit 28bdda282bce666d5e65b77e93c085e3def01dbf
Author: Timo Heister <timo.heister at gmail.com>
Date:   Wed Jun 11 11:14:23 2014 -0400

    implement denormalize_pressure() for direct solver


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

28bdda282bce666d5e65b77e93c085e3def01dbf
 source/simulator/helper_functions.cc | 42 ++++++++++++++++++++++++++++++------
 1 file changed, 36 insertions(+), 6 deletions(-)

diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
index fa77a95..27fdf89 100644
--- a/source/simulator/helper_functions.cc
+++ b/source/simulator/helper_functions.cc
@@ -717,13 +717,43 @@ namespace aspect
     if (parameters.pressure_normalization == "no")
       return;
 
-    // TODO: pressure normalization currently does not work if velocity and
-    // pressure are in the same block.
-    Assert(introspection.block_indices.velocities != introspection.block_indices.pressure,
-           ExcNotImplemented());
-
     if (parameters.use_locally_conservative_discretization == false)
-      vector.block (introspection.block_indices.pressure).add (-1.0 * pressure_adjustment);
+      {
+        if (introspection.block_indices.velocities != introspection.block_indices.pressure)
+          vector.block(introspection.block_indices.pressure).add(-1.0 * pressure_adjustment);
+        else
+          {
+            // velocity and pressure are in the same block, so we have to modify the values manually
+
+            LinearAlgebra::BlockVector distributed_vector (introspection.index_sets.stokes_partitioning,
+                                                           mpi_communicator);
+            distributed_vector = vector;
+
+            std::vector<types::global_dof_index> local_dof_indices (finite_element.dofs_per_cell);
+            typename DoFHandler<dim>::active_cell_iterator
+            cell = dof_handler.begin_active(),
+            endc = dof_handler.end();
+            for (; cell != endc; ++cell)
+              if (cell->is_locally_owned())
+                {
+                  cell->get_dof_indices (local_dof_indices);
+                  for (unsigned int j=0; j<finite_element.base_element(introspection.base_elements.pressure).dofs_per_cell; ++j)
+                    {
+                      const unsigned int local_dof_index
+                        = finite_element.component_to_system_index(introspection.component_indices.pressure,
+                                                                   /*dof index within component=*/ j);
+
+                      // then adjust its value. Note that because we end up touching
+                      // entries more than once, we are not simply incrementing
+                      // distributed_vector but copy from the unchanged vector.
+                      distributed_vector(local_dof_indices[support_point_index])
+                        = vector(local_dof_indices[local_dof_index]) - pressure_adjustment;
+                    }
+                }
+            distributed_vector.compress(VectorOperation::insert);
+            vector = distributed_vector;
+          }
+      }
     else
       {
         // this case is a bit more complicated: if the condition above is false