[cig-commits] [commit] master: Some more generalizations now possible with the extended Introspection class. (946dc1a)

[cig_noreply at geodynamics.org]

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

On branch  : master
Link       : https://github.com/geodynamics/aspect/compare/555391b308f972e57b2488cb024ed6bda549ebf0...485a82a477d4e359644ce65c4155a55cdde612d7

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

commit 946dc1a814a657df0a26f943fc5ac02dae5b4142
Author: Wolfgang Bangerth <bangerth at math.tamu.edu>
Date:   Sun May 11 21:27:31 2014 -0500

    Some more generalizations now possible with the extended Introspection
    class.


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

946dc1a814a657df0a26f943fc5ac02dae5b4142
 source/simulator/initial_conditions.cc | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/source/simulator/initial_conditions.cc b/source/simulator/initial_conditions.cc
index de8e0c2..857bda4 100644
--- a/source/simulator/initial_conditions.cc
+++ b/source/simulator/initial_conditions.cc
@@ -67,14 +67,18 @@ namespace aspect
     // we need to track whether we need to normalize the totality of fields
     bool normalize_composition = false;
 
+//TODO: The code here is confusing. We should be using something
+// like the TemperatureOrComposition class instead of just a single
+// integer 'n'
     for (unsigned int n=0; n<1+parameters.n_compositional_fields; ++n)
       {
         initial_solution.reinit(system_rhs, false);
 
         // base element in the finite element is 2 for temperature (n=0) and 3 for
         // compositional fields (n>0)
-//TODO: can we use introspection here, instead of the hard coded numbers?
-        const unsigned int base_element = (n==0 ? 2 : 3);
+        const unsigned int base_element = (n==0 ?
+                                           introspection.base_elements.temperature :
+                                           introspection.base_elements.compositional_fields);
 
         // get the temperature/composition support points
         const std::vector<Point<dim> > support_points
@@ -100,19 +104,20 @@ namespace aspect
               cell->get_dof_indices (local_dof_indices);
               for (unsigned int i=0; i<finite_element.base_element(base_element).dofs_per_cell; ++i)
                 {
+//TODO: Use introspection here
                   const unsigned int system_local_dof
                     = finite_element.component_to_system_index(/*temperature/composition component=*/dim+1+n,
                         /*dof index within component=*/i);
 
                   const double value =
-                    (base_element == 2
+                    (n == 0
                      ?
                      initial_conditions->initial_temperature(fe_values.quadrature_point(i))
                      :
                      compositional_initial_conditions->initial_composition(fe_values.quadrature_point(i),n-1));
                   initial_solution(local_dof_indices[system_local_dof]) = value;
 
-                  if (base_element != 2)
+                  if (n > 0)
                     Assert (value >= 0,
                             ExcMessage("Invalid initial conditions: Composition is negative"));
 
@@ -123,7 +128,7 @@ namespace aspect
                       double sum = 0;
                       for (unsigned int m=0; m<parameters.normalized_fields.size(); ++m)
                         sum += compositional_initial_conditions->initial_composition(fe_values.quadrature_point(i),parameters.normalized_fields[m]);
-                      if (abs(sum) > 1.0+1e-6)
+                      if (std::abs(sum) > 1.0+1e-6)
                         {
                           max_sum_comp = std::max(sum, max_sum_comp);
                           normalize_composition = true;