[cig-commits] [commit] master: Updating location and naming in DT postprocessor (ac54b57)

[cig_noreply at geodynamics.org]

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

On branch  : master
Link       : https://github.com/geodynamics/aspect/compare/eef6299b61e1952be22819a80f6a6e4ee344b677...09cd1e0322ded4719e699729a269b602229b86b3

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

commit ac54b5710c54a09a558373a4464b440c084de73d
Author: Jacky Austermann <jaustermann at fas.harvard.edu>
Date:   Fri Jan 16 17:01:54 2015 -0500

    Updating location and naming in DT postprocessor


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

ac54b5710c54a09a558373a4464b440c084de73d
 source/postprocess/dynamic_topography.cc               | 16 ++++++++++------
 source/postprocess/visualization/dynamic_topography.cc |  9 ++++++---
 2 files changed, 16 insertions(+), 9 deletions(-)

diff --git a/source/postprocess/dynamic_topography.cc b/source/postprocess/dynamic_topography.cc
index e481627..dda7e1b 100644
--- a/source/postprocess/dynamic_topography.cc
+++ b/source/postprocess/dynamic_topography.cc
@@ -120,11 +120,13 @@ namespace aspect
               
               double dynamic_topography_x_volume = 0;
               double volume = 0;
-              Point<dim> location;
 
+              // Compute the integral of the dynamic topography function 
+              // over the entire cell, by looping over all quadrature points 
+              // (currently, there is only one, but the code is generic).
               for (unsigned int q=0; q<quadrature_formula.size(); ++q)
                 {
-                  location = fe_values.quadrature_point(q);
+                  Point<dim> location = fe_values.quadrature_point(q);
                   const double viscosity = out.viscosities[q];
                   const double density   = out.densities[q];
 
@@ -145,20 +147,22 @@ namespace aspect
                  volume += fe_values.JxW(q);  
                }
                
-               const double dynamic_topography_total = dynamic_topography_x_volume / volume;
+               const double dynamic_topography_cell_average = dynamic_topography_x_volume / volume;
                // Compute the associated surface area to later compute the surfaces weighted integral
-               double surface = 0;  
+               double surface;
+               Point<dim> midpoint_at_surface;  
                for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
                  if (cell->at_boundary(f) && this->get_geometry_model().depth (cell->face(f)->center()) < cell->face(f)->minimum_vertex_distance()/3)
                      {
                      fe_face_values.reinit(cell,f);
                      surface = fe_face_values.JxW(0);
+                     midpoint_at_surface = cell->face(f)->center(); 
                      }
 
-               integrated_topography += dynamic_topography_total*surface;
+               integrated_topography += dynamic_topography_cell_average*surface;
                integrated_surface_area += surface;
 
-               stored_values.push_back (std::make_pair(location, dynamic_topography_total));
+               stored_values.push_back (std::make_pair(midpoint_at_surface, dynamic_topography_cell_average));
             }
 
       // Calculate surface weighted average dynamic topography
diff --git a/source/postprocess/visualization/dynamic_topography.cc b/source/postprocess/visualization/dynamic_topography.cc
index e3d1d5c..9872f57 100644
--- a/source/postprocess/visualization/dynamic_topography.cc
+++ b/source/postprocess/visualization/dynamic_topography.cc
@@ -127,6 +127,9 @@ namespace aspect
                 double dynamic_topography_x_volume = 0;
                 double volume = 0;
 
+                // Compute the integral of the dynamic topography function 
+                // over the entire cell, by looping over all quadrature points 
+                // (currently, there is only one, but the code is generic).
                 for (unsigned int q=0; q<quadrature_formula.size(); ++q)
                   {
                     const Point<dim> location = fe_values.quadrature_point(q);
@@ -150,7 +153,7 @@ namespace aspect
                     volume += fe_values.JxW(q); 
                  }
        
-                 const double dynamic_topography_total = dynamic_topography_x_volume / volume;
+                 const double dynamic_topography_cell_average = dynamic_topography_x_volume / volume;
                  // Compute the associated surface area to later compute the surfaces weighted integral
                  double surface = 0;
                  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
@@ -160,10 +163,10 @@ namespace aspect
                        surface = fe_face_values.JxW(0);
                        }
 
-                 integrated_topography += dynamic_topography_total*surface;
+                 integrated_topography += dynamic_topography_cell_average*surface;
                  integrated_surface_area += surface;
        
-                 (*return_value.second)(cell_index) = dynamic_topography_total;
+                 (*return_value.second)(cell_index) = dynamic_topography_cell_average;
                }
 
         // Calculate surface weighted average dynamic topography