[cig-commits] [commit] master: Use introspection in several more places. Also rename AdvectionField object 'advf' rather than 'torc' for clarity. (86cf52d)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu May 22 06:05:09 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/5b2b586a4f29499284711e77f2cb8d5a0bc64afe...80f462a3e073ac705f736fd43243a28f9833c866
>---------------------------------------------------------------
commit 86cf52de759746c1f788e90e603111eb36762327
Author: Jonathan Perry-Houts <jperryh2 at uoregon.edu>
Date: Wed May 21 16:30:08 2014 -0700
Use introspection in several more places. Also rename AdvectionField object 'advf' rather than 'torc' for clarity.
>---------------------------------------------------------------
86cf52de759746c1f788e90e603111eb36762327
source/mesh_refinement/density.cc | 4 ++--
source/mesh_refinement/nonadiabatic_temperature.cc | 4 ++--
source/mesh_refinement/thermal_energy_density.cc | 4 ++--
source/mesh_refinement/viscosity.cc | 4 ++--
source/postprocess/visualization/density.cc | 8 ++++----
.../postprocess/visualization/friction_heating.cc | 10 +++++-----
source/postprocess/visualization/melt_fraction.cc | 8 ++++----
.../visualization/nonadiabatic_pressure.cc | 2 +-
.../visualization/nonadiabatic_temperature.cc | 4 ++--
source/postprocess/visualization/seismic_vp.cc | 8 ++++----
source/postprocess/visualization/seismic_vs.cc | 8 ++++----
source/postprocess/visualization/specific_heat.cc | 8 ++++----
source/postprocess/visualization/strain_rate.cc | 4 ++--
.../visualization/thermal_expansivity.cc | 8 ++++----
.../visualization/thermodynamic_phase.cc | 8 ++++----
source/postprocess/visualization/viscosity.cc | 10 +++++-----
.../postprocess/visualization/viscosity_ratio.cc | 10 +++++-----
source/simulator/helper_functions.cc | 2 +-
source/simulator/initial_conditions.cc | 23 ++++++++++++----------
source/simulator/nullspace.cc | 6 +++---
20 files changed, 73 insertions(+), 70 deletions(-)
diff --git a/source/mesh_refinement/density.cc b/source/mesh_refinement/density.cc
index e90f1aa..ab52dd0 100644
--- a/source/mesh_refinement/density.cc
+++ b/source/mesh_refinement/density.cc
@@ -102,7 +102,7 @@ namespace aspect
for (unsigned int i=0; i<this->get_fe().base_element(this->introspection().base_elements.temperature).dofs_per_cell; ++i)
{
const unsigned int system_local_dof
- = this->get_fe().component_to_system_index(/*temperature component=*/dim+1,
+ = this->get_fe().component_to_system_index(this->introspection().component_indices.temperature,
/*dof index within component=*/i);
vec_distributed(local_dof_indices[system_local_dof])
@@ -123,7 +123,7 @@ namespace aspect
vec,
indicators,
//TODO: replace by the appropriate component mask
- dim+1);
+ this->introspection().component_indices.temperature);
// Scale gradient in each cell with the correct power of h. Otherwise,
// error indicators do not reduce when refined if there is a density
diff --git a/source/mesh_refinement/nonadiabatic_temperature.cc b/source/mesh_refinement/nonadiabatic_temperature.cc
index c2bf06c..69a0d40 100644
--- a/source/mesh_refinement/nonadiabatic_temperature.cc
+++ b/source/mesh_refinement/nonadiabatic_temperature.cc
@@ -80,7 +80,7 @@ namespace aspect
for (unsigned int i=0; i<this->get_fe().base_element(this->introspection().base_elements.temperature).dofs_per_cell; ++i)
{
const unsigned int system_local_dof
- = this->get_fe().component_to_system_index(/*temperature component=*/dim+1,
+ = this->get_fe().component_to_system_index(this->introspection().component_indices.temperature,
/*dof index within component=*/i);
vec_distributed(local_dof_indices[system_local_dof])
@@ -99,7 +99,7 @@ namespace aspect
vec,
indicators,
//TODO: replace by the appropriate component mask
- dim+1);
+ this->introspection().component_indices.temperature);
// Scale gradient in each cell with the correct power of h. Otherwise,
// error indicators do not reduce when refined if there is a density
diff --git a/source/mesh_refinement/thermal_energy_density.cc b/source/mesh_refinement/thermal_energy_density.cc
index b33a151..0ba63fc 100644
--- a/source/mesh_refinement/thermal_energy_density.cc
+++ b/source/mesh_refinement/thermal_energy_density.cc
@@ -100,7 +100,7 @@ namespace aspect
for (unsigned int i=0; i<this->get_fe().base_element(this->introspection().base_elements.temperature).dofs_per_cell; ++i)
{
const unsigned int system_local_dof
- = this->get_fe().component_to_system_index(/*temperature component=*/dim+1,
+ = this->get_fe().component_to_system_index(this->introspection().component_indices.temperature,
/*dof index within component=*/i);
vec_distributed(local_dof_indices[system_local_dof])
@@ -121,7 +121,7 @@ namespace aspect
vec,
indicators,
//TODO: get this from introspection
- dim+1);
+ this->introspection().component_indices.temperature);
// Scale gradient in each cell with the correct power of h. Otherwise,
// error indicators do not reduce when refined if there is a density
diff --git a/source/mesh_refinement/viscosity.cc b/source/mesh_refinement/viscosity.cc
index 3c7078c..c61c20f 100644
--- a/source/mesh_refinement/viscosity.cc
+++ b/source/mesh_refinement/viscosity.cc
@@ -102,7 +102,7 @@ namespace aspect
for (unsigned int i=0; i<this->get_fe().base_element(this->introspection().base_elements.temperature).dofs_per_cell; ++i)
{
const unsigned int system_local_dof
- = this->get_fe().component_to_system_index(/*temperature component=*/dim+1,
+ = this->get_fe().component_to_system_index(this->introspection().component_indices.temperature,
/*dof index within component=*/i);
vec_distributed(local_dof_indices[system_local_dof])
@@ -121,7 +121,7 @@ namespace aspect
vec,
indicators,
//TODO: replace by the appropriate component mask
- dim+1);
+ this->introspection().component_indices.temperature);
// Scale gradient in each cell with the correct power of h. Otherwise,
// error indicators do not reduce when refined if there is a viscosity
diff --git a/source/postprocess/visualization/density.cc b/source/postprocess/visualization/density.cc
index ffe8d3e..45c425f 100644
--- a/source/postprocess/visualization/density.cc
+++ b/source/postprocess/visualization/density.cc
@@ -54,7 +54,7 @@ namespace aspect
const unsigned int n_quadrature_points = uh.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+ Assert (uh[0].size() == this->introspection().n_components, ExcInternalError());
typename MaterialModel::Interface<dim>::MaterialModelInputs in(n_quadrature_points,
this->n_compositional_fields());
@@ -65,11 +65,11 @@ namespace aspect
in.strain_rate.resize(0); // we do not need the viscosity
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
- in.pressure[q]=uh[q][dim];
- in.temperature[q]=uh[q][dim+1];
+ in.pressure[q]=uh[q][this->introspection().component_indices.pressure];
+ in.temperature[q]=uh[q][this->introspection().component_indices.temperature];
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
- in.composition[q][c] = uh[q][dim+2+c];
+ in.composition[q][c] = uh[q][this->introspection().component_indices.compositional_fields[0]+c];
}
this->get_material_model().evaluate(in, out);
diff --git a/source/postprocess/visualization/friction_heating.cc b/source/postprocess/visualization/friction_heating.cc
index 888ea2a..798ace8 100644
--- a/source/postprocess/visualization/friction_heating.cc
+++ b/source/postprocess/visualization/friction_heating.cc
@@ -54,8 +54,8 @@ namespace aspect
const unsigned int n_quadrature_points = uh.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
- Assert (duh[0].size() == dim+2+this->n_compositional_fields(),ExcInternalError());
+ Assert (uh[0].size() == this->introspection().n_components, ExcInternalError());
+ Assert (duh[0].size() == this->introspection().n_components, ExcInternalError());
typename MaterialModel::Interface<dim>::MaterialModelInputs in(n_quadrature_points,
this->n_compositional_fields());
@@ -70,11 +70,11 @@ namespace aspect
grad_u[d] = duh[q][d];
in.strain_rate[q] = symmetrize (grad_u);
- in.pressure[q]=uh[q][dim];
- in.temperature[q]=uh[q][dim+1];
+ in.pressure[q]=uh[q][this->introspection().component_indices.pressure];
+ in.temperature[q]=uh[q][this->introspection().component_indices.temperature];
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
- in.composition[q][c] = uh[q][dim+2+c];
+ in.composition[q][c] = uh[q][this->introspection().component_indices.compositional_fields[0]+c];
}
diff --git a/source/postprocess/visualization/melt_fraction.cc b/source/postprocess/visualization/melt_fraction.cc
index 2aacbb2..a1b9533 100644
--- a/source/postprocess/visualization/melt_fraction.cc
+++ b/source/postprocess/visualization/melt_fraction.cc
@@ -55,16 +55,16 @@ namespace aspect
const unsigned int n_quadrature_points = uh.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+ Assert (uh[0].size() == this->introspection().n_components, ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
- const double pressure = uh[q][dim];
- const double temperature = uh[q][dim+1];
+ const double pressure = uh[q][this->introspection().component_indices.pressure];
+ const double temperature = uh[q][this->introspection().component_indices.temperature];
std::vector<double> composition(this->n_compositional_fields());
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
- composition[c] = uh[q][dim+2+c];
+ composition[c] = uh[q][this->introspection().component_indices.compositional_fields[0]+c];
// anhydrous melting of peridotite after Katz, 2003
const double T_solidus = A1 + 273.15
diff --git a/source/postprocess/visualization/nonadiabatic_pressure.cc b/source/postprocess/visualization/nonadiabatic_pressure.cc
index 003231a..d1a45eb 100644
--- a/source/postprocess/visualization/nonadiabatic_pressure.cc
+++ b/source/postprocess/visualization/nonadiabatic_pressure.cc
@@ -54,7 +54,7 @@ namespace aspect
const unsigned int n_quadrature_points = uh.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+ Assert (uh[0].size() == this->introspection().n_components, ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
diff --git a/source/postprocess/visualization/nonadiabatic_temperature.cc b/source/postprocess/visualization/nonadiabatic_temperature.cc
index 3df7641..181373d 100644
--- a/source/postprocess/visualization/nonadiabatic_temperature.cc
+++ b/source/postprocess/visualization/nonadiabatic_temperature.cc
@@ -54,11 +54,11 @@ namespace aspect
const unsigned int n_quadrature_points = uh.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+ Assert (uh[0].size() == this->introspection().n_components, ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
- const double temperature = uh[q][dim+1];
+ const double temperature = uh[q][this->introspection().component_indices.temperature];
computed_quantities[q](0) = temperature - this->get_adiabatic_conditions().temperature(evaluation_points[q]);
}
diff --git a/source/postprocess/visualization/seismic_vp.cc b/source/postprocess/visualization/seismic_vp.cc
index 3f45b16..fb1ff83 100644
--- a/source/postprocess/visualization/seismic_vp.cc
+++ b/source/postprocess/visualization/seismic_vp.cc
@@ -54,15 +54,15 @@ namespace aspect
const unsigned int n_quadrature_points = uh.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+ Assert (uh[0].size() == this->introspection().n_components, ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
- const double pressure = uh[q][dim];
- const double temperature = uh[q][dim+1];
+ const double pressure = uh[q][this->introspection().component_indices.pressure];
+ const double temperature = uh[q][this->introspection().component_indices.temperature];
std::vector<double> composition(this->n_compositional_fields());
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
- composition[c] = uh[q][dim+2+c];
+ composition[c] = uh[q][this->introspection().component_indices.compositional_fields[0]+c];
computed_quantities[q](0) = this->get_material_model().seismic_Vp(temperature,
pressure,
diff --git a/source/postprocess/visualization/seismic_vs.cc b/source/postprocess/visualization/seismic_vs.cc
index ce536b9..183a709 100644
--- a/source/postprocess/visualization/seismic_vs.cc
+++ b/source/postprocess/visualization/seismic_vs.cc
@@ -54,15 +54,15 @@ namespace aspect
const unsigned int n_quadrature_points = uh.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+ Assert (uh[0].size() == this->introspection().n_components, ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
- const double pressure = uh[q][dim];
- const double temperature = uh[q][dim+1];
+ const double pressure = uh[q][this->introspection().component_indices.pressure];
+ const double temperature = uh[q][this->introspection().component_indices.temperature];
std::vector<double> composition(this->n_compositional_fields());
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
- composition[c] = uh[q][dim+2+c];
+ composition[c] = uh[q][this->introspection().component_indices.compositional_fields[0]+c];
computed_quantities[q](0) = this->get_material_model().seismic_Vs(temperature,
pressure,
diff --git a/source/postprocess/visualization/specific_heat.cc b/source/postprocess/visualization/specific_heat.cc
index c82dda7..118b07f 100644
--- a/source/postprocess/visualization/specific_heat.cc
+++ b/source/postprocess/visualization/specific_heat.cc
@@ -54,7 +54,7 @@ namespace aspect
const unsigned int n_quadrature_points = uh.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+ Assert (uh[0].size() == this->introspection().n_components, ExcInternalError());
typename MaterialModel::Interface<dim>::MaterialModelInputs in(n_quadrature_points,
this->n_compositional_fields());
@@ -65,11 +65,11 @@ namespace aspect
in.strain_rate.resize(0); // we do not need the viscosity
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
- in.pressure[q]=uh[q][dim];
- in.temperature[q]=uh[q][dim+1];
+ in.pressure[q]=uh[q][this->introspection().component_indices.pressure];
+ in.temperature[q]=uh[q][this->introspection().component_indices.temperature];
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
- in.composition[q][c] = uh[q][dim+2+c];
+ in.composition[q][c] = uh[q][this->introspection().component_indices.compositional_fields[0]+c];
}
this->get_material_model().evaluate(in, out);
diff --git a/source/postprocess/visualization/strain_rate.cc b/source/postprocess/visualization/strain_rate.cc
index 304df9a..2e8bdc3 100644
--- a/source/postprocess/visualization/strain_rate.cc
+++ b/source/postprocess/visualization/strain_rate.cc
@@ -54,8 +54,8 @@ namespace aspect
const unsigned int n_quadrature_points = uh.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
- Assert (duh[0].size() == dim+2+this->n_compositional_fields(),ExcInternalError());
+ Assert (uh[0].size() == this->introspection().n_components, ExcInternalError());
+ Assert (duh[0].size() == this->introspection().n_components, ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
diff --git a/source/postprocess/visualization/thermal_expansivity.cc b/source/postprocess/visualization/thermal_expansivity.cc
index d3364ed..f670db6 100644
--- a/source/postprocess/visualization/thermal_expansivity.cc
+++ b/source/postprocess/visualization/thermal_expansivity.cc
@@ -54,7 +54,7 @@ namespace aspect
const unsigned int n_quadrature_points = uh.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+ Assert (uh[0].size() == this->introspection().n_components, ExcInternalError());
typename MaterialModel::Interface<dim>::MaterialModelInputs in(n_quadrature_points,
this->n_compositional_fields());
@@ -66,11 +66,11 @@ namespace aspect
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
//in.strain_rate[q] =
- in.pressure[q]=uh[q][dim];
- in.temperature[q]=uh[q][dim+1];
+ in.pressure[q]=uh[q][this->introspection().component_indices.pressure];
+ in.temperature[q]=uh[q][this->introspection().component_indices.temperature];
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
- in.composition[q][c] = uh[q][dim+2+c];
+ in.composition[q][c] = uh[q][this->introspection().component_indices.compositional_fields[0]+c];
}
diff --git a/source/postprocess/visualization/thermodynamic_phase.cc b/source/postprocess/visualization/thermodynamic_phase.cc
index 3a330fd..f7463f6 100644
--- a/source/postprocess/visualization/thermodynamic_phase.cc
+++ b/source/postprocess/visualization/thermodynamic_phase.cc
@@ -54,15 +54,15 @@ namespace aspect
const unsigned int n_quadrature_points = uh.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+ Assert (uh[0].size() == this->introspection().n_components, ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
- const double pressure = uh[q][dim];
- const double temperature = uh[q][dim+1];
+ const double pressure = uh[q][this->introspection().component_indices.pressure];
+ const double temperature = uh[q][this->introspection().component_indices.temperature];
std::vector<double> composition(this->n_compositional_fields());
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
- composition[c] = uh[q][dim+2+c];
+ composition[c] = uh[q][this->introspection().component_indices.compositional_fields[0]+c];
computed_quantities[q](0) = this->get_material_model().thermodynamic_phase(temperature,
pressure,
diff --git a/source/postprocess/visualization/viscosity.cc b/source/postprocess/visualization/viscosity.cc
index b1036fd..2392eaa 100644
--- a/source/postprocess/visualization/viscosity.cc
+++ b/source/postprocess/visualization/viscosity.cc
@@ -54,8 +54,8 @@ namespace aspect
const unsigned int n_quadrature_points = uh.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
- Assert (duh[0].size() == dim+2+this->n_compositional_fields(),ExcInternalError());
+ Assert (uh[0].size() == this->introspection().n_components, ExcInternalError());
+ Assert (duh[0].size() == this->introspection().n_components, ExcInternalError());
typename MaterialModel::Interface<dim>::MaterialModelInputs in(n_quadrature_points,
this->n_compositional_fields());
@@ -70,11 +70,11 @@ namespace aspect
grad_u[d] = duh[q][d];
in.strain_rate[q] = symmetrize (grad_u);
- in.pressure[q]=uh[q][dim];
- in.temperature[q]=uh[q][dim+1];
+ in.pressure[q]=uh[q][this->introspection().component_indices.pressure];
+ in.temperature[q]=uh[q][this->introspection().component_indices.temperature];
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
- in.composition[q][c] = uh[q][dim+2+c];
+ in.composition[q][c] = uh[q][this->introspection().component_indices.compositional_fields[0]+c];
}
this->get_material_model().evaluate(in, out);
diff --git a/source/postprocess/visualization/viscosity_ratio.cc b/source/postprocess/visualization/viscosity_ratio.cc
index 24adbfd..e3ed171 100644
--- a/source/postprocess/visualization/viscosity_ratio.cc
+++ b/source/postprocess/visualization/viscosity_ratio.cc
@@ -54,8 +54,8 @@ namespace aspect
const unsigned int n_quadrature_points = uh.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
- Assert (duh[0].size() == dim+2+this->n_compositional_fields(),ExcInternalError());
+ Assert (uh[0].size() == this->introspection().n_components, ExcInternalError());
+ Assert (duh[0].size() == this->introspection().n_components, ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
@@ -64,8 +64,8 @@ namespace aspect
for (unsigned int d=0; d<dim; ++d)
grad_u[d] = duh[q][d];
- const double pressure = uh[q][dim];
- const double temperature = uh[q][dim+1];
+ const double pressure = uh[q][this->introspection().component_indices.pressure];
+ const double temperature = uh[q][this->introspection().component_indices.temperature];
const SymmetricTensor<2,dim> strain_rate = symmetrize (grad_u);
const SymmetricTensor<2,dim> compressible_strain_rate
@@ -77,7 +77,7 @@ namespace aspect
std::vector<double> composition(this->n_compositional_fields());
for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
- composition[c] = uh[q][dim+2+c];
+ composition[c] = uh[q][this->introspection().component_indices.compositional_fields[0]+c];
computed_quantities[q](0) = this->get_material_model().viscosity_ratio(temperature,
pressure,
diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
index 8fa78da..4eeaf0c 100644
--- a/source/simulator/helper_functions.cc
+++ b/source/simulator/helper_functions.cc
@@ -730,7 +730,7 @@ namespace aspect
ExcNotImplemented());
if (parameters.use_locally_conservative_discretization == false)
- vector.block (1).add (-1.0 * pressure_adjustment);
+ vector.block (introspection.block_indices.pressure).add (-1.0 * pressure_adjustment);
else
{
// this case is a bit more complicated: if the condition above is false
diff --git a/source/simulator/initial_conditions.cc b/source/simulator/initial_conditions.cc
index 5c9d4c3..20b8c61 100644
--- a/source/simulator/initial_conditions.cc
+++ b/source/simulator/initial_conditions.cc
@@ -72,12 +72,15 @@ namespace aspect
// should probably ask AdvectionField how many fields there actually are
for (unsigned int n=0; n<1+parameters.n_compositional_fields; ++n)
{
- AdvectionField torc = (n==0) ? AdvectionField::temperature()
- : AdvectionField::composition(n-1);
- initial_solution.reinit(system_rhs, false);
+ AdvectionField advf = ((n == 0)
+ ?
+ AdvectionField::temperature()
+ :
+ AdvectionField::composition(n-1));
- const unsigned int base_element = torc.base_element(introspection);
+ initial_solution.reinit(system_rhs, false);
+ const unsigned int base_element = advf.base_element(introspection);
// get the temperature/composition support points
const std::vector<Point<dim> > support_points
@@ -104,18 +107,18 @@ namespace aspect
for (unsigned int i=0; i<finite_element.base_element(base_element).dofs_per_cell; ++i)
{
const unsigned int system_local_dof
- = finite_element.component_to_system_index(torc.component_index(introspection),
- /*dof index within component=*/i);
+ = finite_element.component_to_system_index(advf.component_index(introspection),
+ /*dof index within component=*/i);
const double value =
- (torc.is_temperature()
+ (advf.is_temperature()
?
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 (!torc.is_temperature())
+ if (!advf.is_temperature())
Assert (value >= 0,
ExcMessage("Invalid initial conditions: Composition is negative"));
@@ -141,7 +144,7 @@ namespace aspect
// we should not have written at all into any of the blocks with
// the exception of the current temperature or composition block
for (unsigned int b=0; b<initial_solution.n_blocks(); ++b)
- if (b != torc.block_index(introspection))
+ if (b != advf.block_index(introspection))
Assert (initial_solution.block(b).l2_norm() == 0,
ExcInternalError());
@@ -168,7 +171,7 @@ namespace aspect
constraints.distribute(initial_solution);
// copy temperature/composition block only
- const unsigned int blockidx = torc.block_index(introspection);
+ const unsigned int blockidx = advf.block_index(introspection);
solution.block(blockidx) = initial_solution.block(blockidx);
old_solution.block(blockidx) = initial_solution.block(blockidx);
old_old_solution.block(blockidx) = initial_solution.block(blockidx);
diff --git a/source/simulator/nullspace.cc b/source/simulator/nullspace.cc
index 9124454..35824b8 100644
--- a/source/simulator/nullspace.cc
+++ b/source/simulator/nullspace.cc
@@ -229,10 +229,10 @@ namespace aspect
parameters.n_compositional_fields);
for (unsigned int i=0; i< q_points.size(); i++)
{
- in.pressure[i] = fe_vals[i][dim];
- in.temperature[i] = fe_vals[i][dim+1];
+ in.pressure[i] = fe_vals[i][introspection.component_indices.pressure];
+ in.temperature[i] = fe_vals[i][introspection.component_indices.temperature];
for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
- in.composition[i][c] = fe_vals[i][dim+2+c];
+ in.composition[i][c] = fe_vals[i][introspection.component_indices.compositional_fields[0]+c];
in.position[i] = q_points[i];
}
More information about the CIG-COMMITS
mailing list