[cig-commits] [commit] master: Addressed comments. (5c9362d)

[cig_noreply at geodynamics.org]

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

On branch  : master
Link       : https://github.com/geodynamics/aspect/compare/78abb106b73747fb364d8481c9df27e32aed4597...ea45ef11705f505f2d2d63c77aa04f2d7be9a4c1

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

commit 5c9362d13f5c744a79cdec6eda792929859b8e1b
Author: Rene Gassmoeller <R.Gassmoeller at mailbox.org>
Date:   Mon Aug 4 15:25:52 2014 +0200

    Addressed comments.


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

5c9362d13f5c744a79cdec6eda792929859b8e1b
 source/simulator/helper_functions.cc | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
index c88d483..3679ce4 100644
--- a/source/simulator/helper_functions.cc
+++ b/source/simulator/helper_functions.cc
@@ -814,8 +814,13 @@ namespace aspect
     Assert(introspection.block_indices.velocities != introspection.block_indices.pressure,
            ExcNotImplemented());
 
-    const QGauss<dim-1> quadrature_formula (parameters.stokes_velocity_degree+1);
+    // In the following we integrate the normal velocity over every surface
+    // of the model. This integral is part of the correction term that needs
+    // to be added to the pressure right hand side. To calculate the normal
+    // velocity we need the positions and normals of every quadrature point on
+    // the surface.
 
+    const QGauss<dim-1> quadrature_formula (parameters.stokes_velocity_degree+1);
     FEFaceValues<dim> fe_face_values (mapping,
                                       finite_element,
                                       quadrature_formula,
@@ -829,20 +834,21 @@ namespace aspect
       cell = dof_handler.begin_active(),
       endc = dof_handler.end();
 
-    // for every surface face that is part of a geometry boundary
-    // and that is owned by this processor,
+    // for every surface face that is part of a geometry boundary with
+    // prescribed velocity and that is owned by this processor,
     // integrate the normal velocity magnitude.
     for (; cell!=endc; ++cell)
       if (cell->is_locally_owned())
         for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
           if (cell->face(f)->at_boundary() &&
-              parameters.prescribed_velocity_boundary_indicators.find(cell->face(f)->boundary_indicator())!=
-              parameters.prescribed_velocity_boundary_indicators.end())
+              (parameters.prescribed_velocity_boundary_indicators.find(cell->face(f)->boundary_indicator())!=
+              parameters.prescribed_velocity_boundary_indicators.end()))
             {
               fe_face_values.reinit (cell, f);
 
               // for each of the quadrature points, evaluate the
-              // normal velocity and add it to the integral
+              // normal velocity by calling the boundary conditions and add
+              // it to the integral
               for (unsigned int q=0; q<quadrature_formula.size(); ++q)
                 {
                   const Tensor<1,dim> velocity =
@@ -857,7 +863,7 @@ namespace aspect
         Utilities::MPI::sum (local_normal_velocity_integral,mpi_communicator);
 
     const double mean       = vector.block(introspection.block_indices.pressure).mean_value();
-    const double correction = (local_normal_velocity_integral - mean * vector.block(introspection.block_indices.pressure).size()) / global_volume;
+    const double correction = (global_normal_velocity_integral - mean * vector.block(introspection.block_indices.pressure).size()) / global_volume;
 
     vector.block(introspection.block_indices.pressure).add(correction, pressure_shape_function_integrals.block(introspection.block_indices.pressure));
   }