Merge pull request #1537 from pmbremner/material_constructor

Create a new MaterialModel constructor.
geodynamics · May 10, 2017 · ca3e368 · ca3e368
2 parents 5873817 + cdcb02c
commit ca3e368
Show file tree

Hide file tree

Showing 3 changed files with 71 additions and 27 deletions.
diff --git a/include/aspect/material_model/interface.h b/include/aspect/material_model/interface.h
@@ -33,6 +33,10 @@
 
 namespace aspect
 {
+  template <int dim>
+  struct Introspection;
+
+
   /**
    * A namespace in which we define everything that has to do with modeling
    * convecting material, including descriptions of material parameters such
@@ -180,6 +184,18 @@ namespace aspect
       MaterialModelInputs(const unsigned int n_points,
                           const unsigned int n_comp);
 
+
+      /**
+       * Constructor. Both initialize and populate the various arrays of this
+       * structure with the FEValues and introspection objects and
+       * the solution_vector.
+       */
+      MaterialModelInputs(const FEValuesBase<dim,dim> &fe_values,
+                          const typename DoFHandler<dim>::active_cell_iterator *cell,
+                          const Introspection<dim> &introspection,
+                          const LinearAlgebra::BlockVector &solution_vector);
+
+
       /**
        * Vector with global positions where the material has to be evaluated
        * in evaluate().

diff --git a/source/material_model/interface.cc b/source/material_model/interface.cc
@@ -255,6 +255,42 @@ namespace aspect
     {}
 
 
+    template <int dim>
+    MaterialModelInputs<dim>::MaterialModelInputs(const FEValuesBase<dim,dim> &fe_values,
+                                                  const typename DoFHandler<dim>::active_cell_iterator *cell_x,
+                                                  const Introspection<dim> &introspection,
+                                                  const LinearAlgebra::BlockVector &solution_vector)
+      :
+      position(fe_values.get_quadrature_points()),
+      temperature(fe_values.n_quadrature_points, numbers::signaling_nan<double>()),
+      pressure(fe_values.n_quadrature_points, numbers::signaling_nan<double>()),
+      pressure_gradient(fe_values.n_quadrature_points, numbers::signaling_nan<Tensor<1,dim> >()),
+      velocity(fe_values.n_quadrature_points, numbers::signaling_nan<Tensor<1,dim> >()),
+      composition(fe_values.n_quadrature_points, std::vector<double>(introspection.n_compositional_fields, numbers::signaling_nan<double>())),
+      strain_rate(fe_values.n_quadrature_points, numbers::signaling_nan<SymmetricTensor<2,dim> >()),
+      cell(cell_x)
+    {
+      // Populate the newly allocated arrays
+      fe_values[introspection.extractors.temperature].get_function_values (solution_vector, this->temperature);
+      fe_values[introspection.extractors.velocities].get_function_values (solution_vector, this->velocity);
+      fe_values[introspection.extractors.pressure].get_function_values (solution_vector, this->pressure);
+      fe_values[introspection.extractors.pressure].get_function_gradients (solution_vector, this->pressure_gradient);
+      fe_values[introspection.extractors.velocities].get_function_symmetric_gradients (solution_vector,this->strain_rate);
+
+      // Vectors for evaluating the the compositional field parts of the finite element solution
+      std::vector<std::vector<double> > composition_values (introspection.n_compositional_fields, std::vector<double> (fe_values.n_quadrature_points));
+      for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
+        {
+          fe_values[introspection.extractors.compositional_fields[c]].get_function_values(solution_vector,composition_values[c]);
+        }
+
+      for (unsigned int i=0; i<fe_values.n_quadrature_points; ++i)
+        {
+          this->position[i] = fe_values.quadrature_point(i);
+          for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
+            this->composition[i][c] = composition_values[c][i];
+        }
+    }
 
     template <int dim>
     MaterialModelOutputs<dim>::MaterialModelOutputs(const unsigned int n_points,

diff --git a/source/simulator/lateral_averaging.cc b/source/simulator/lateral_averaging.cc
@@ -60,11 +60,6 @@ namespace aspect
     cell = this->get_dof_handler().begin_active(),
     endc = this->get_dof_handler().end();
 
-    MaterialModel::MaterialModelInputs<dim> in(n_q_points,
-                                               this->n_compositional_fields());
-    MaterialModel::MaterialModelOutputs<dim> out(n_q_points,
-                                                 this->n_compositional_fields());
-
     fctr.setup(quadrature_formula.size());
 
     for (; cell!=endc; ++cell)
@@ -73,32 +68,29 @@ namespace aspect
           fe_values.reinit (cell);
           if (fctr.need_material_properties())
             {
-              fe_values[this->introspection().extractors.pressure].get_function_values (this->get_solution(),
-                                                                                        in.pressure);
-              fe_values[this->introspection().extractors.temperature].get_function_values (this->get_solution(),
-                                                                                           in.temperature);
-              fe_values[this->introspection().extractors.velocities].get_function_values (this->get_solution(),
-                                                                                          in.velocity);
-              fe_values[this->introspection().extractors.velocities].get_function_symmetric_gradients (this->get_solution(),
-                  in.strain_rate);
-              fe_values[this->introspection().extractors.pressure].get_function_gradients (this->get_solution(),
-                                                                                           in.pressure_gradient);
-              for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
-                fe_values[this->introspection().extractors.compositional_fields[c]].get_function_values(this->get_solution(),
-                    composition_values[c]);
-
-              for (unsigned int i=0; i<n_q_points; ++i)
-                {
-                  in.position[i] = fe_values.quadrature_point(i);
-                  for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
-                    in.composition[i][c] = composition_values[c][i];
-                }
-              in.cell = &cell;
+
+              // get the density at each quadrature point if necessary
+              MaterialModel::MaterialModelInputs<dim> in(fe_values,
+                                                         &cell,
+                                                         this->introspection(),
+                                                         this->get_solution());
+              MaterialModel::MaterialModelOutputs<dim> out(n_q_points,
+                                                           this->n_compositional_fields());
 
               this->get_material_model().evaluate(in, out);
+              fctr(in, out, fe_values, this->get_solution(), output_values);
+
+            }
+          else
+            {
+              MaterialModel::MaterialModelInputs<dim> in(n_q_points,
+                                                         this->n_compositional_fields());
+              MaterialModel::MaterialModelOutputs<dim> out(n_q_points,
+                                                           this->n_compositional_fields());
+
+              fctr(in, out, fe_values, this->get_solution(), output_values);
             }
 
-          fctr(in, out, fe_values, this->get_solution(), output_values);
 
           for (unsigned int q = 0; q < n_q_points; ++q)
             {