[cig-commits] [commit] master: implement normalize_pressure (c61d9c7)

[cig_noreply at geodynamics.org]

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

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

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

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

    implement normalize_pressure


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

c61d9c77d214b3c50ffd4f58fc2800acd95c7ac0
 source/simulator/helper_functions.cc | 46 ++++++++++++++++++++++++++++++------
 1 file changed, 39 insertions(+), 7 deletions(-)

diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
index 488e298..3bfb383 100644
--- a/source/simulator/helper_functions.cc
+++ b/source/simulator/helper_functions.cc
@@ -535,11 +535,6 @@ 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());
-
     double my_pressure = 0.0;
     double my_area = 0.0;
     if (parameters.pressure_normalization == "surface")
@@ -636,7 +631,42 @@ namespace aspect
     distributed_vector = vector;
 
     if (parameters.use_locally_conservative_discretization == false)
-      distributed_vector.block(introspection.block_indices.pressure).add(pressure_adjustment);
+      {
+        if (introspection.block_indices.velocities != introspection.block_indices.pressure)
+          distributed_vector.block(introspection.block_indices.pressure).add(pressure_adjustment);
+        else
+          {
+            // velocity and pressure are in the same block, so we have to modify the values manually
+            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())
+                {
+                  // identify the first pressure dof
+                  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)
+                  {
+                      unsigned int support_point_index
+                          = finite_element.component_to_system_index(introspection.component_indices.pressure,
+                                                                                             /*dof index within component=*/ j);
+
+                  // make sure that this DoF is really owned by the current processor
+                  // and that it is in fact a pressure dof
+                  Assert (dof_handler.locally_owned_dofs().is_element(local_dof_indices[support_point_index]),
+                          ExcInternalError());
+
+                  Assert (introspection.block_indices.velocities == introspection.block_indices.pressure
+                      || local_dof_indices[support_point_index] >= vector.block(0).size(),
+                          ExcInternalError());
+
+                  // then adjust its value
+                  distributed_vector(local_dof_indices[support_point_index]) += pressure_adjustment;
+                  }
+                }
+          }
+      }
     else
       {
         // this case is a bit more complicated: if the condition above is false
@@ -668,7 +698,9 @@ namespace aspect
               // and that it is in fact a pressure dof
               Assert (dof_handler.locally_owned_dofs().is_element(local_dof_indices[first_pressure_dof]),
                       ExcInternalError());
-              Assert (local_dof_indices[first_pressure_dof] >= vector.block(0).size(),
+
+              Assert (introspection.block_indices.velocities == introspection.block_indices.pressure
+                  || local_dof_indices[first_pressure_dof] >= vector.block(0).size(),
                       ExcInternalError());
 
               // then adjust its value