Implement MaterialModelInputs reinitialization function.

geodynamics · May 13, 2017 · b084299 · b084299
1 parent b26c788
commit b084299
Show file tree

Hide file tree

Showing 3 changed files with 32 additions and 47 deletions.
diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
@@ -435,14 +435,21 @@ namespace aspect
 
     std::vector<Tensor<1,dim> > velocity_values(n_q_points);
     std::vector<Tensor<1,dim> > fluid_velocity_values(n_q_points);
-    std::vector<std::vector<double> > composition_values (parameters.n_compositional_fields,std::vector<double> (n_q_points));
+    std::vector<std::vector<double> > composition_values (introspection.n_compositional_fields,std::vector<double> (n_q_points));
 
     double max_local_speed_over_meshsize = 0;
     double min_local_conduction_timestep = std::numeric_limits<double>::max();
 
     typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
+
+
+    MaterialModel::MaterialModelInputs<dim> in(n_q_points,
+                                               introspection.n_compositional_fields);
+    MaterialModel::MaterialModelOutputs<dim> out(n_q_points,
+                                                 introspection.n_compositional_fields);
+
     for (; cell!=endc; ++cell)
       if (cell->is_locally_owned())
         {
@@ -472,34 +479,14 @@ namespace aspect
 
           if (parameters.use_conduction_timestep)
             {
-              MaterialModel::MaterialModelInputs<dim> in(n_q_points, parameters.n_compositional_fields);
-              MaterialModel::MaterialModelOutputs<dim> out(n_q_points, parameters.n_compositional_fields);
-
-              in.strain_rate.resize(0); // we do not need the viscosity
-              in.position = fe_values.get_quadrature_points();
-              in.cell = &cell;
-
-              fe_values[introspection.extractors.pressure].get_function_values (solution,
-                                                                                in.pressure);
-              fe_values[introspection.extractors.temperature].get_function_values (solution,
-                                                                                   in.temperature);
-              fe_values[introspection.extractors.pressure].get_function_gradients (solution,
-                                                                                   in.pressure_gradient);
-
-              for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
-                fe_values[introspection.extractors.compositional_fields[c]].get_function_values (solution,
-                    composition_values[c]);
-
-              for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
-                {
-                  for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
-                    in.composition[q][c] = composition_values[c][q];
-
-                  in.velocity[q] = velocity_values[q];
-                }
+              in.reinit(fe_values,
+                        &cell,
+                        introspection,
+                        solution);
 
               material_model->evaluate(in, out);
 
+
               // Evaluate thermal diffusivity at each quadrature point and
               // calculate the corresponding conduction timestep, if applicable
               for (unsigned int q=0; q<n_q_points; ++q)

diff --git a/source/simulator/lateral_averaging.cc b/source/simulator/lateral_averaging.cc
@@ -60,6 +60,12 @@ 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)
@@ -70,27 +76,16 @@ namespace aspect
             {
 
               // 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());
+              in.reinit(fe_values,
+                        &cell,
+                        this->introspection(),
+                        this->get_solution());
 
               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)
             {

diff --git a/source/simulator/melt.cc b/source/simulator/melt.cc
@@ -915,6 +915,7 @@ namespace aspect
 
       MaterialModel::MaterialModelInputs<dim> in(quadrature.size(), this->n_compositional_fields());
       MaterialModel::MaterialModelOutputs<dim> out(quadrature.size(), this->n_compositional_fields());
+
       create_material_model_outputs(out);
 
       typename DoFHandler<dim>::active_cell_iterator cell = this->get_dof_handler().begin_active(),
@@ -934,11 +935,13 @@ namespace aspect
 
             fe_values[this->introspection().variable("fluid pressure").extractor_scalar()].get_function_gradients (
               solution, grad_p_f_values);
-            this->compute_material_model_input_values (solution,
-                                                       fe_values,
-                                                       cell,
-                                                       true,
-                                                       in);
+
+
+            in.reinit(fe_values,
+                      &cell,
+                      this->introspection(),
+                      this->get_solution());
+
 
             this->get_material_model().evaluate(in, out);