diff --git a/include/aspect/material_model/interface.h b/include/aspect/material_model/interface.h index 9e9974585a3..6a86a999c9b 100644 --- a/include/aspect/material_model/interface.h +++ b/include/aspect/material_model/interface.h @@ -1247,10 +1247,10 @@ namespace aspect std::vector dilation_lhs_term; /** - * A scalar value per evaluation point corresponding to the RHS term + * A scalar value per direction per evaluation point corresponding to the RHS term * due to plastic dilation. */ - std::vector dilation_rhs_term; + std::vector> dilation_rhs_term; }; diff --git a/source/material_model/interface.cc b/source/material_model/interface.cc index 7c7cd1b4455..4ef6eb0bbb1 100644 --- a/source/material_model/interface.cc +++ b/source/material_model/interface.cc @@ -999,11 +999,14 @@ namespace aspect namespace { - std::vector make_prescribed_dilation_outputs_names() + std::vector make_prescribed_dilation_outputs_names(int dim) { std::vector names; names.emplace_back("dilation_lhs_term"); names.emplace_back("dilation_rhs_term"); + names.emplace_back("dilation_rhs_term"); + if (dim == 3) + names.emplace_back("dilation_rhs_term"); return names; } } @@ -1012,9 +1015,9 @@ namespace aspect template PrescribedPlasticDilation::PrescribedPlasticDilation (const unsigned int n_points) - : NamedAdditionalMaterialOutputs(make_prescribed_dilation_outputs_names()), + : NamedAdditionalMaterialOutputs(make_prescribed_dilation_outputs_names(dim)), dilation_lhs_term(n_points, numbers::signaling_nan()), - dilation_rhs_term(n_points, numbers::signaling_nan()) + dilation_rhs_term(dim, std::vector(n_points, numbers::signaling_nan())) {} @@ -1022,14 +1025,20 @@ namespace aspect template std::vector PrescribedPlasticDilation::get_nth_output(const unsigned int idx) const { - AssertIndexRange (idx, 2); + AssertIndexRange (idx, 4); switch (idx) { case 0: return dilation_lhs_term; case 1: - return dilation_rhs_term; + return dilation_rhs_term[0]; + + case 2: + return dilation_rhs_term[1]; + + case 3: + return dilation_rhs_term[2]; default: AssertThrow(false, ExcInternalError()); diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc index 0e54fad8dee..fd267df8ef7 100644 --- a/source/material_model/visco_plastic.cc +++ b/source/material_model/visco_plastic.cc @@ -299,9 +299,15 @@ namespace aspect // should cancel out (RHS = LHS * p) so as to satisfy the loading-unloading conditions. plastic_dilation->dilation_lhs_term[i] = dilation_lhs_term; if (dilation_rhs_term - dilation_lhs_term * in.pressure[i] > 0) - plastic_dilation->dilation_rhs_term[i] = dilation_rhs_term; + { + for (unsigned int dim_i = 0; dim_i < dim; ++dim_i) + plastic_dilation->dilation_rhs_term[dim_i][i] = dilation_rhs_term; + } else - plastic_dilation->dilation_rhs_term[i] = dilation_lhs_term * in.pressure[i]; + { + for (unsigned int dim_i = 0; dim_i < dim; ++dim_i) + plastic_dilation->dilation_rhs_term[dim_i][i] = dilation_lhs_term * in.pressure[i]; + } } } diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc index 516116489fd..489e1e91a27 100644 --- a/source/simulator/assemblers/newton_stokes.cc +++ b/source/simulator/assemblers/newton_stokes.cc @@ -461,13 +461,17 @@ namespace aspect // dilation term, but also the LHS dilation term should be included in the // system residual if (enable_prescribed_dilation) - data.local_rhs(i) += ( - - pressure_scaling - * (prescribed_dilation->dilation_rhs_term[q] - - prescribed_dilation->dilation_lhs_term[q] * - scratch.material_model_inputs.pressure[q]) - * scratch.phi_p[i] - ) * JxW; + { + const unsigned int index_direction=fe.system_to_component_index(i).first; + if (introspection.is_stokes_component(index_direction+1)) + data.local_rhs(i) += ( + - pressure_scaling + * (prescribed_dilation->dilation_rhs_term[index_direction][q] - + prescribed_dilation->dilation_lhs_term[q] * + scratch.material_model_inputs.pressure[q]) + * scratch.phi_p[i] + ) * JxW; + } } // and then the matrix, if necessary @@ -665,7 +669,10 @@ namespace aspect Assert(!this->get_parameters().enable_prescribed_dilation || (outputs.template get_additional_output_object>()->dilation_lhs_term.size() == n_points && - outputs.template get_additional_output_object>()->dilation_rhs_term.size() == n_points), + outputs.template get_additional_output_object>()->dilation_rhs_term.size() == dim && + outputs.template get_additional_output_object>()->dilation_rhs_term[0].size() == n_points && + outputs.template get_additional_output_object>()->dilation_rhs_term[1].size() == n_points && + (dim == 2 || (dim == 3 && outputs.template get_additional_output_object>()->dilation_rhs_term[2].size() == n_points))), ExcInternalError()); if (this->get_newton_handler().parameters.newton_derivative_scaling_factor != 0) diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc index 7fd884b0394..fd75f07db6c 100644 --- a/source/simulator/assemblers/stokes.cc +++ b/source/simulator/assemblers/stokes.cc @@ -451,11 +451,20 @@ namespace aspect * JxW; if (prescribed_dilation != nullptr) - data.local_rhs(i) += ( - - pressure_scaling - * prescribed_dilation->dilation_rhs_term[q] - * scratch.phi_p[i] - ) * JxW; + { + const unsigned int index_direction=fe.system_to_component_index(i).first; + // Only want the velocity components and not the pressure one (which is the last one), so add 1 + if (introspection.is_stokes_component(index_direction+1)) + { + data.local_rhs(i) += ( + - pressure_scaling + * (prescribed_dilation->dilation_rhs_term[index_direction][q] - + prescribed_dilation->dilation_lhs_term[q] * + scratch.material_model_inputs.pressure[q]) + * scratch.phi_p[i] + ) * JxW; + } + } if (scratch.rebuild_stokes_matrix) for (unsigned int j=0; jget_parameters().enable_prescribed_dilation || (outputs.template get_additional_output_object>()->dilation_lhs_term.size() == n_points && - outputs.template get_additional_output_object>()->dilation_rhs_term.size() == n_points), + outputs.template get_additional_output_object>()->dilation_rhs_term.size() == dim && + outputs.template get_additional_output_object>()->dilation_rhs_term[0].size() == n_points && + outputs.template get_additional_output_object>()->dilation_rhs_term[1].size() == n_points && + (dim == 2 || (dim == 3 && outputs.template get_additional_output_object>()->dilation_rhs_term[2].size() == n_points))), ExcInternalError()); // Elasticity: diff --git a/tests/prescribed_dilation.cc b/tests/prescribed_dilation.cc index cb4afcf89b8..465feede6ae 100644 --- a/tests/prescribed_dilation.cc +++ b/tests/prescribed_dilation.cc @@ -207,8 +207,11 @@ namespace aspect } if (prescribed_dilation) { - prescribed_dilation->dilation_rhs_term[i] = x; - prescribed_dilation->dilation_lhs_term[i] = 0.; + prescribed_dilation->dilation_rhs_term[0][i] = x; + prescribed_dilation->dilation_rhs_term[1][i] = x; + if (dim == 3) + prescribed_dilation->dilation_rhs_term[2][i] = x; + prescribed_dilation->dilation_lhs_term[i] = 0; } }