Skip to content

Commit

Permalink
Merge pull request #1689 from ian-r-rose/use_material_model_reinit
Browse files Browse the repository at this point in the history 
Use MaterialModelInputs<dim>::reinit() where appropriate.
  • Loading branch information
@bangerth
bangerth committed May 16, 2017
2 parents 04524c1 + 0f689ca commit 32ae628
Show file tree
Hide file tree
Showing 13 changed files with 38 additions and 323 deletions.
6 changes: 4 additions & 2 deletions include/aspect/material_model/interface.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -194,7 +194,8 @@ namespace aspect
MaterialModelInputs(const FEValuesBase<dim,dim> &fe_values,
const typename DoFHandler<dim>::active_cell_iterator *cell,
const Introspection<dim> &introspection,
const LinearAlgebra::BlockVector &solution_vector);
const LinearAlgebra::BlockVector &solution_vector,
const bool use_strain_rates = true);

/**
* Function to re-initialize and populate the pre-existing arrays
Expand All @@ -203,7 +204,8 @@ namespace aspect
void reinit(const FEValuesBase<dim,dim> &fe_values,
const typename DoFHandler<dim>::active_cell_iterator *cell,
const Introspection<dim> &introspection,
const LinearAlgebra::BlockVector &solution_vector);
const LinearAlgebra::BlockVector &solution_vector,
const bool use_strain_rates = true);


/**
Expand Down
13 changes: 9 additions & 4 deletions source/material_model/interface.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,7 +245,8 @@ namespace aspect
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)
const LinearAlgebra::BlockVector &solution_vector,
const bool use_strain_rate)
:
position(fe_values.get_quadrature_points()),
temperature(fe_values.n_quadrature_points, numbers::signaling_nan<double>()),
Expand All @@ -257,7 +258,7 @@ namespace aspect
cell(cell_x)
{
// Call the function reinit to populate the new arrays.
this->reinit(fe_values, cell, introspection, solution_vector);
this->reinit(fe_values, cell, introspection, solution_vector, use_strain_rate);
}


Expand All @@ -266,14 +267,18 @@ namespace aspect
MaterialModelInputs<dim>::reinit(const FEValuesBase<dim,dim> &fe_values,
const typename DoFHandler<dim>::active_cell_iterator *cell_x,
const Introspection<dim> &introspection,
const LinearAlgebra::BlockVector &solution_vector)
const LinearAlgebra::BlockVector &solution_vector,
const bool use_strain_rate)
{
// 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);
if (use_strain_rate)
fe_values[introspection.extractors.velocities].get_function_symmetric_gradients (solution_vector,this->strain_rate);
else
this->strain_rate.resize(0);

// 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));
Expand Down
24 changes: 2 additions & 22 deletions source/mesh_refinement/density.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,28 +76,8 @@ namespace aspect
if (cell->is_locally_owned())
{
fe_values.reinit(cell);

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.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(),
prelim_composition_values[c]);

in.position = fe_values.get_quadrature_points();
in.strain_rate.resize(0);// we are not reading the viscosity
for (unsigned int i=0; i<quadrature.size(); ++i)
{
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
in.composition[i][c] = prelim_composition_values[c][i];
}
in.cell = &cell;
// Set use_strain_rates to false since we don't need viscosity
in.reinit(fe_values, &cell, this->introspection(), this->get_solution(), false);

this->get_material_model().evaluate(in, out);

Expand Down
24 changes: 3 additions & 21 deletions source/mesh_refinement/thermal_energy_density.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,29 +70,11 @@ namespace aspect
if (cell->is_locally_owned())
{
fe_values.reinit(cell);
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.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(),
prelim_composition_values[c]);
// Set use_strain_rates to false since we don't need viscosity
in.reinit(fe_values, &cell, this->introspection(), this->get_solution(), false);
this->get_material_model().evaluate(in, out);

cell->get_dof_indices (local_dof_indices);
in.position = fe_values.get_quadrature_points();
in.strain_rate.resize(0);// we are not reading the viscosity
for (unsigned int i=0; i<quadrature.size(); ++i)
{
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
in.composition[i][c] = prelim_composition_values[c][i];
}
in.cell = &cell;

this->get_material_model().evaluate(in, out);

// for each temperature dof, write into the output
// vector the density. note that quadrature points and
Expand Down
25 changes: 1 addition & 24 deletions source/mesh_refinement/viscosity.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,30 +76,7 @@ namespace aspect
if (cell->is_locally_owned())
{
fe_values.reinit(cell);

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_symmetric_gradients (this->get_solution(),
in.strain_rate);
fe_values[this->introspection().extractors.velocities].get_function_values (this->get_solution(),
in.velocity);
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(),
prelim_composition_values[c]);

in.position = fe_values.get_quadrature_points();
for (unsigned int i=0; i<quadrature.size(); ++i)
{
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
in.composition[i][c] = prelim_composition_values[c][i];
}
in.cell = &cell;

in.reinit(fe_values, &cell, this->introspection(), this->get_solution());
this->get_material_model().evaluate(in, out);

cell->get_dof_indices (local_dof_indices);
Expand Down
78 changes: 4 additions & 74 deletions source/postprocess/geoid.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,31 +117,8 @@ namespace aspect
if (cell->is_locally_owned())
{
fe_values.reinit (cell);
// get the various components of the solution, then
// evaluate the material properties there
fe_values[this->introspection().extractors.temperature]
.get_function_values (this->get_solution(), in.temperature);
fe_values[this->introspection().extractors.pressure]
.get_function_values (this->get_solution(), in.pressure);
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);

in.position = fe_values.get_quadrature_points();

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<fe_values.n_quadrature_points; ++i)
{
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
in.composition[i][c] = composition_values[c][i];
}
in.cell = &cell;
// Set use_strain_rates to false since we don't need viscosity
in.reinit(fe_values, &cell, this->introspection(), this->get_solution(), false);

this->get_material_model().evaluate(in, out);

Expand Down Expand Up @@ -263,31 +240,7 @@ namespace aspect

// focus on the boundary cell
fe_values.reinit (cell);
// get the various components of the solution, then
// evaluate the material properties there
fe_values[this->introspection().extractors.temperature]
.get_function_values (this->get_solution(), in.temperature);
fe_values[this->introspection().extractors.pressure]
.get_function_values (this->get_solution(), in.pressure);
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);

in.position = fe_values.get_quadrature_points();

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<fe_values.n_quadrature_points; ++i)
{
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
in.composition[i][c] = composition_values[c][i];
}
in.cell = &cell;
in.reinit(fe_values, &cell, this->introspection(), this->get_solution());

this->get_material_model().evaluate(in, out);

Expand Down Expand Up @@ -349,30 +302,7 @@ namespace aspect

// focus on the boundary cell's upper face if on the top boundary and lower face if on the bottom boundary
fe_face_values.reinit(cell, face_idx);
// get the various components of the solution, then
// evaluate the material properties there
fe_face_values[this->introspection().extractors.temperature]
.get_function_values (this->get_solution(), in_face.temperature);
fe_face_values[this->introspection().extractors.pressure]
.get_function_values (this->get_solution(), in_face.pressure);
fe_face_values[this->introspection().extractors.velocities]
.get_function_values (this->get_solution(), in_face.velocity);
fe_face_values[this->introspection().extractors.velocities]
.get_function_symmetric_gradients (this->get_solution(), in_face.strain_rate);
fe_face_values[this->introspection().extractors.pressure]
.get_function_gradients (this->get_solution(), in_face.pressure_gradient);

in_face.position = fe_face_values.get_quadrature_points();

for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
fe_face_values[this->introspection().extractors.compositional_fields[c]]
.get_function_values(this->get_solution(),
face_composition_values[c]);
for (unsigned int i=0; i<fe_face_values.n_quadrature_points; ++i)
{
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
in_face.composition[i][c] = face_composition_values[c][i];
}
in_face.reinit(fe_face_values, &cell, this->introspection(), this->get_solution(), false);

this->get_material_model().evaluate(in_face, out_face);

Expand Down
33 changes: 4 additions & 29 deletions source/postprocess/heat_flux_densities.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,38 +71,13 @@ namespace aspect
if (cell->at_boundary(f))
{
fe_face_values.reinit (cell, f);
fe_face_values[this->introspection().extractors.temperature].get_function_gradients (this->get_solution(),
temperature_gradients);
fe_face_values[this->introspection().extractors.temperature].get_function_values (this->get_solution(),
in.temperature);
fe_face_values[this->introspection().extractors.pressure].get_function_values (this->get_solution(),
in.pressure);
fe_face_values[this->introspection().extractors.velocities].get_function_values (this->get_solution(),
in.velocity);
fe_face_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_face_values[this->introspection().extractors.compositional_fields[c]].get_function_values(this->get_solution(),
composition_values[c]);

in.position = fe_face_values.get_quadrature_points();

// since we are not reading the viscosity and the viscosity
// is the only coefficient that depends on the strain rate,
// we need not compute the strain rate. set the corresponding
// array to empty, to prevent accidental use and skip the
// evaluation of the strain rate in evaluate().
in.strain_rate.resize(0);

for (unsigned int i=0; i<fe_face_values.n_quadrature_points; ++i)
{
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
in.composition[i][c] = composition_values[c][i];
}
in.cell = &cell;
// Set use_strain_rates to false since we don't need viscosity
in.reinit(fe_face_values, &cell, this->introspection(), this->get_solution(), false);

this->get_material_model().evaluate(in, out);

//Get the temperature gradients from the solution.
fe_face_values[this->introspection().extractors.temperature].get_function_gradients (this->get_solution(), temperature_gradients);

double local_normal_flux = 0;
double local_area = 0;
Expand Down

0 comments on commit 32ae628

Please sign in to comment.