From c3ccd5b31ea02ed5ea39e0e8365272cbd178243b Mon Sep 17 00:00:00 2001 From: j507 Date: Mon, 16 Sep 2024 17:10:09 +0200 Subject: [PATCH] (M) Postprocessing output --- applications/simple_shear.cc | 5 +- include/gCP/postprocessing.h | 34 +- .../gradient_crystal_plasticity_solver.cc | 2 + source/postprocessing.cc | 428 +++++++----------- 4 files changed, 186 insertions(+), 283 deletions(-) diff --git a/applications/simple_shear.cc b/applications/simple_shear.cc index 38d9e6f..545399b 100644 --- a/applications/simple_shear.cc +++ b/applications/simple_shear.cc @@ -163,7 +163,10 @@ homogenization( mapping), postprocessor( fe_field, - crystals_data), + crystals_data, + parameters.solver_parameters.dimensionless_form_parameters, + false, + parameters.output.flag_output_dimensionless_quantities), simple_shear( fe_field, mapping, diff --git a/include/gCP/postprocessing.h b/include/gCP/postprocessing.h index 6b70a3f..d5e667c 100644 --- a/include/gCP/postprocessing.h +++ b/include/gCP/postprocessing.h @@ -32,10 +32,12 @@ class Postprocessor : public dealii::DataPostprocessor { public: Postprocessor( - std::shared_ptr> &fe_field, - std::shared_ptr> &crystals_data, - const bool flag_light_output = false, - const bool flag_output_fluctuations = false); + std::shared_ptr> &fe_field, + std::shared_ptr> &crystals_data, + const RunTimeParameters::DimensionlessForm ¶meters, + const bool flag_light_output = false, + const bool flag_output_dimensionless_quantities = false, + const bool flag_output_fluctuations = false); virtual void evaluate_vector_field( const dealii::DataPostprocessorInputs::Vector &inputs, @@ -57,25 +59,29 @@ class Postprocessor : public dealii::DataPostprocessor const dealii::SymmetricTensor<2,dim> macroscopic_strain); private: - std::shared_ptr> fe_field; + std::shared_ptr> fe_field; - std::shared_ptr> crystals_data; + std::shared_ptr> crystals_data; - std::shared_ptr> hooke_law; + std::shared_ptr> hooke_law; - std::vector> voigt_indices; + std::vector> voigt_indices; - dealii::SymmetricTensor<2,dim> macroscopic_strain; + dealii::SymmetricTensor<2,dim> macroscopic_strain; - const dealii::SymmetricTensor<4,dim> deviatoric_projector; + const dealii::SymmetricTensor<4,dim> deviatoric_projector; - const dealii::SymmetricTensor<4,3> deviatoric_projector_3d; + const dealii::SymmetricTensor<4,3> deviatoric_projector_3d; - bool flag_light_output; + RunTimeParameters::DimensionlessForm parameters; - bool flag_output_fluctuations; + bool flag_light_output; - bool flag_init_was_called; + bool flag_output_dimensionless_quantities; + + bool flag_output_fluctuations; + + bool flag_init_was_called; dealii::SymmetricTensor<2,3> convert_2d_to_3d( dealii::SymmetricTensor<2,dim> symmetric_tensor) const; diff --git a/source/gradient_crystal_plasticity/gradient_crystal_plasticity_solver.cc b/source/gradient_crystal_plasticity/gradient_crystal_plasticity_solver.cc index 704a3cf..c58fc46 100644 --- a/source/gradient_crystal_plasticity/gradient_crystal_plasticity_solver.cc +++ b/source/gradient_crystal_plasticity/gradient_crystal_plasticity_solver.cc @@ -66,6 +66,8 @@ nonlinear_solver_logger( postprocessor( fe_field, crystals_data, + parameters.dimensionless_form_parameters, + true, true), flag_init_was_called(false) { diff --git a/source/postprocessing.cc b/source/postprocessing.cc index 13cdba7..4d94770 100644 --- a/source/postprocessing.cc +++ b/source/postprocessing.cc @@ -17,10 +17,12 @@ namespace Postprocessing template Postprocessor::Postprocessor( - std::shared_ptr> &fe_field, - std::shared_ptr> &crystals_data, - const bool flag_light_output, - const bool flag_output_fluctuations) + std::shared_ptr> &fe_field, + std::shared_ptr> &crystals_data, + const RunTimeParameters::DimensionlessForm ¶meters, + const bool flag_light_output, + const bool flag_output_dimensionless_quantities, + const bool flag_output_fluctuations) : fe_field(fe_field), crystals_data(crystals_data), @@ -37,7 +39,10 @@ deviatoric_projector_3d( dealii::outer_product( dealii::unit_symmetric_tensor<3>(), dealii::unit_symmetric_tensor<3>())), +parameters(parameters), flag_light_output(flag_light_output), +flag_output_dimensionless_quantities( + flag_output_dimensionless_quantities), flag_output_fluctuations(flag_output_fluctuations) { macroscopic_strain = 0.; @@ -49,6 +54,12 @@ flag_output_fluctuations(flag_output_fluctuations) voigt_indices[4] = std::make_pair(0,2); voigt_indices[5] = std::make_pair(0,1); + if (flag_output_dimensionless_quantities) + { + this->parameters.characteristic_quantities = + RunTimeParameters::CharacteristicQuantities(); + } + if (flag_light_output) { flag_init_was_called = true; @@ -247,50 +258,59 @@ void Postprocessor::evaluate_vector_field( { AssertThrow(flag_init_was_called, dealii::ExcMessage("The Postprocessor instance has" - " not been initialized.")); + " not been initialized.")); + // Get data const typename dealii::DoFHandler::cell_iterator current_cell = inputs.template get_cell(); - - const unsigned int material_id = current_cell->material_id(); - - const unsigned int n_q_points = inputs.solution_values.size(); - + const unsigned int material_id = current_cell->material_id(); + const unsigned int n_quadrature_points = + inputs.solution_values.size(); const unsigned int n_components = fe_field->get_n_components(); - - const unsigned int n_slips = crystals_data->get_n_slips(); - - const unsigned int n_crystals = crystals_data->get_n_crystals(); + const unsigned int n_slips = crystals_data->get_n_slips(); + const unsigned int n_crystals = crystals_data->get_n_crystals(); + const double &dimensionless_number = + parameters.dimensionless_numbers[0]; (void)n_components; - Assert(inputs.solution_gradients.size() == n_q_points, - dealii::ExcInternalError()); - - Assert(computed_quantities.size() == n_q_points, - dealii::ExcInternalError()); - + Assert(inputs.solution_gradients.size() == n_quadrature_points, + dealii::ExcInternalError()); + Assert(computed_quantities.size() == n_quadrature_points, + dealii::ExcInternalError()); Assert(inputs.solution_values[0].size() == n_components, - dealii::ExcInternalError()); + dealii::ExcInternalError()); // Reset - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) - for (unsigned int d = 0; d < computed_quantities[0].size(); ++d) - computed_quantities[q_point](d) = 0.0; - - - dealii::Tensor<2,dim> displacement_gradient; - dealii::SymmetricTensor<2,dim> strain_tensor; - dealii::SymmetricTensor<2,3> strain_tensor_3d; - dealii::SymmetricTensor<2,dim> plastic_strain_tensor; - dealii::SymmetricTensor<2,3> plastic_strain_tensor_3d; - dealii::SymmetricTensor<2,dim> elastic_strain_tensor; - dealii::SymmetricTensor<2,3> elastic_strain_tensor_3d; - dealii::SymmetricTensor<2,3> stress_tensor; - double equivalent_edge_dislocation_density; - double equivalent_screw_dislocation_density; + for (unsigned int quadrature_point_id = 0; + quadrature_point_id < n_quadrature_points; + ++quadrature_point_id) + { + for (unsigned int index = 0; index < computed_quantities[0].size(); + ++index) + { + computed_quantities[quadrature_point_id](index) = 0.0; + } + } - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + // Local instances + dealii::Tensor<2,dim> displacement_gradient; + dealii::SymmetricTensor<2,dim> strain_tensor; + dealii::SymmetricTensor<2,dim> elastic_strain_tensor; + dealii::SymmetricTensor<2,dim> plastic_strain_tensor; + dealii::SymmetricTensor<2,3> strain_tensor_3d; + dealii::SymmetricTensor<2,3> plastic_strain_tensor_3d; + dealii::SymmetricTensor<2,3> elastic_strain_tensor_3d; + dealii::SymmetricTensor<2,3> stress_tensor; + double equivalent_edge_dislocation_density; + double equivalent_screw_dislocation_density; + + const unsigned int index_offset = + fe_field->is_decohesion_allowed() ? dim * n_crystals : dim; + + for (unsigned int quadrature_point_index = 0; + quadrature_point_index < n_quadrature_points; + ++quadrature_point_index) { // Reset displacement_gradient = 0.0; @@ -304,305 +324,183 @@ void Postprocessor::evaluate_vector_field( equivalent_edge_dislocation_density = 0.0; equivalent_screw_dislocation_density = 0.0; - /*! - * @note This if-else can probably be done in a way more elegant - * manner - */ - if (fe_field->is_decohesion_allowed()) + // Displacement + for (unsigned int index = 0; index < dim; ++index) { - // Displacement - for (unsigned int d = 0; d < dim; ++d) + if (fe_field->is_decohesion_allowed()) { for (unsigned int crystal_id = 0; crystal_id < n_crystals; ++crystal_id) { - computed_quantities[q_point](d) += - inputs.solution_values[q_point](d + dim * crystal_id); + computed_quantities[quadrature_point_index](index) += + parameters.characteristic_quantities.displacement * + inputs.solution_values[quadrature_point_index]( + index + dim * crystal_id); if (!flag_light_output) { - displacement_gradient[d] += - inputs.solution_gradients[q_point][d + dim * crystal_id]; + displacement_gradient[index] += + inputs.solution_gradients[quadrature_point_index][ + index + dim * crystal_id]; } } } - - for (unsigned int slip_id = 0; - slip_id < n_slips; ++slip_id) + else { - for (unsigned int crystal_id = 0; - crystal_id < n_crystals; ++crystal_id) + computed_quantities[quadrature_point_index](index) = + parameters.characteristic_quantities.displacement * + inputs.solution_values[quadrature_point_index](index); + + if (!flag_light_output) { + displacement_gradient[index] = + inputs.solution_gradients[quadrature_point_index][index]; + } + } + } + + // Slips + for (unsigned int slip_id = 0; + slip_id < n_slips; ++slip_id) + { + for (unsigned int crystal_id = 0; + crystal_id < n_crystals; ++crystal_id) + { // Slips - computed_quantities[q_point](dim + slip_id) += - inputs.solution_values[q_point]( - dim * n_crystals + slip_id + n_slips * crystal_id); + computed_quantities[quadrature_point_index](dim + slip_id) += + inputs.solution_values[quadrature_point_index]( + index_offset + slip_id + n_slips * crystal_id); if (!flag_light_output) { // Equivalent plastic strain - computed_quantities[q_point](dim + n_slips) += - inputs.solution_values[q_point]( - dim * n_crystals + slip_id + n_slips * crystal_id); + computed_quantities[quadrature_point_index](dim + n_slips) += + inputs.solution_values[quadrature_point_index]( + index_offset + slip_id + n_slips * crystal_id); // Equivalent absolute plastic strain - computed_quantities[q_point](dim + n_slips + 1) += - std::abs(inputs.solution_values[q_point]( - dim * n_crystals + slip_id + n_slips * crystal_id)); + computed_quantities[quadrature_point_index]( + dim + n_slips + 1) += + std::abs(inputs.solution_values[quadrature_point_index]( + index_offset + slip_id + n_slips * crystal_id)); // Equivalent edge dislocation density equivalent_edge_dislocation_density += - std::pow(inputs.solution_gradients[q_point][ - dim * n_crystals + slip_id + n_slips * crystal_id] * + std::pow(inputs.solution_gradients[quadrature_point_index][ + index_offset + slip_id + n_slips * crystal_id] * crystals_data->get_slip_direction(crystal_id, slip_id), 2); // Equivalent screw dislocation density equivalent_screw_dislocation_density += - std::pow(inputs.solution_gradients[q_point][ - dim * n_crystals + slip_id + n_slips * crystal_id] * + std::pow(inputs.solution_gradients[quadrature_point_index][ + index_offset + slip_id + n_slips * crystal_id] * crystals_data->get_slip_orthogonal(crystal_id, slip_id), 2); // Plastic strain tensor plastic_strain_tensor += - inputs.solution_values[q_point]( - dim * n_crystals + slip_id + n_slips * crystal_id) * + inputs.solution_values[quadrature_point_index]( + index_offset + slip_id + n_slips * crystal_id) * crystals_data->get_symmetrized_schmid_tensor( crystal_id, slip_id); } - } } + } - if (!flag_light_output) - { - strain_tensor = dealii::symmetrize(displacement_gradient) + - macroscopic_strain; + if (!flag_light_output) + { + computed_quantities[quadrature_point_index](dim + n_slips + 2) = + parameters.characteristic_quantities.dislocation_density * + std::sqrt(equivalent_edge_dislocation_density); - strain_tensor_3d = convert_2d_to_3d(strain_tensor); + computed_quantities[quadrature_point_index](dim + n_slips + 3) = + parameters.characteristic_quantities.dislocation_density * + std::sqrt(equivalent_screw_dislocation_density); - plastic_strain_tensor_3d = - convert_2d_to_3d(plastic_strain_tensor); + strain_tensor = + dealii::symmetrize(displacement_gradient) + + macroscopic_strain; - computed_quantities[q_point](dim + n_slips + 2) = - std::sqrt(equivalent_edge_dislocation_density); + elastic_strain_tensor = + strain_tensor - dimensionless_number * plastic_strain_tensor; - computed_quantities[q_point](dim + n_slips + 3) = - std::sqrt(equivalent_screw_dislocation_density); + strain_tensor_3d = convert_2d_to_3d(strain_tensor); - elastic_strain_tensor = - strain_tensor - plastic_strain_tensor; + plastic_strain_tensor_3d = + convert_2d_to_3d(plastic_strain_tensor); - elastic_strain_tensor_3d = - convert_2d_to_3d(elastic_strain_tensor); + elastic_strain_tensor_3d = + convert_2d_to_3d(elastic_strain_tensor); - stress_tensor = - hooke_law->get_stiffness_tetrad_3d(material_id) * - convert_2d_to_3d(elastic_strain_tensor); + stress_tensor = + hooke_law->get_stiffness_tetrad_3d(material_id) * + convert_2d_to_3d(elastic_strain_tensor); - // Von-Mises stress - computed_quantities[q_point](dim + n_slips + 4) = - get_von_mises_stress(stress_tensor); + // Von-Mises stress + computed_quantities[quadrature_point_index](dim + n_slips + 4) = + parameters.characteristic_quantities.stress * + get_von_mises_stress(stress_tensor); - // Von-Mises plastic strain - computed_quantities[q_point](dim + n_slips + 5) = - get_von_mises_plastic_strain(plastic_strain_tensor); + // Von-Mises plastic strain + computed_quantities[quadrature_point_index](dim + n_slips + 5) = + get_von_mises_plastic_strain(plastic_strain_tensor); + for (unsigned int i = 0; i < voigt_indices.size(); ++i) + { // Stress components - for (unsigned int i = 0; i < voigt_indices.size(); ++i) - computed_quantities[q_point](dim + n_slips + 6 + i) = + computed_quantities[quadrature_point_index](dim + n_slips + 6 + i) = + parameters.characteristic_quantities.stress * stress_tensor[voigt_indices[i].first][voigt_indices[i].second]; // Strain components - for (unsigned int i = 0; i < voigt_indices.size(); ++i) - computed_quantities[q_point](dim + n_slips + 12 + i) = - (i < 3 ? 1.0 : 2.0) * - strain_tensor_3d[voigt_indices[i].first][voigt_indices[i].second]; + computed_quantities[quadrature_point_index](dim + n_slips + 12 + i) = + (i < 3 ? 1.0 : 2.0) * + parameters.characteristic_quantities.strain * + strain_tensor_3d[voigt_indices[i].first][voigt_indices[i].second]; // Elastic strain components - for (unsigned int i = 0; i < voigt_indices.size(); ++i) - computed_quantities[q_point](dim + n_slips + 18 + i) = - (i < 3 ? 1.0 : 2.0) * - elastic_strain_tensor_3d[voigt_indices[i].first][voigt_indices[i].second]; + computed_quantities[quadrature_point_index](dim + n_slips + 18 + i) = + (i < 3 ? 1.0 : 2.0) * + parameters.characteristic_quantities.strain * + elastic_strain_tensor_3d[voigt_indices[i].first][voigt_indices[i].second]; // Plastic strain components - for (unsigned int i = 0; i < voigt_indices.size(); ++i) - computed_quantities[q_point](dim + n_slips + 24 + i) = - (i < 3 ? 1.0 : 2.0) * - plastic_strain_tensor_3d[voigt_indices[i].first][voigt_indices[i].second]; - - if (flag_output_fluctuations) - { - strain_tensor = dealii::symmetrize(displacement_gradient); - - strain_tensor_3d = convert_2d_to_3d(strain_tensor); - - elastic_strain_tensor = - strain_tensor - plastic_strain_tensor; - - stress_tensor = - hooke_law->get_stiffness_tetrad_3d(material_id) * - convert_2d_to_3d(elastic_strain_tensor); - - // Von-Mises stress - computed_quantities[q_point](dim + n_slips + 30) = - get_von_mises_stress(stress_tensor); - - // Stress components - for (unsigned int i = 0; i < voigt_indices.size(); ++i) - computed_quantities[q_point](dim + n_slips + 31 + i) = - stress_tensor[voigt_indices[i].first][voigt_indices[i].second]; - - // Strain components - for (unsigned int i = 0; i < voigt_indices.size(); ++i) - computed_quantities[q_point](dim + n_slips + 36 + i) = - (i < 3 ? 1.0 : 2.0) * - strain_tensor_3d[voigt_indices[i].first][voigt_indices[i].second]; - } + computed_quantities[quadrature_point_index](dim + n_slips + 24 + i) = + (i < 3 ? 1.0 : 2.0) * + plastic_strain_tensor_3d[voigt_indices[i].first][voigt_indices[i].second]; } - } - else - { - // Displacement - for (unsigned int d = 0; d < dim; ++d) - { - computed_quantities[q_point](d) = - inputs.solution_values[q_point](d); - if (!flag_light_output) - { - displacement_gradient[d] = - inputs.solution_gradients[q_point][d]; - } - } - - for (unsigned int slip_id = 0; - slip_id < n_slips; ++slip_id) + if (flag_output_fluctuations) { - for (unsigned int crystal_id = 0; - crystal_id < n_crystals; ++crystal_id) - { - // Slips - computed_quantities[q_point](dim + slip_id) += - inputs.solution_values[q_point]( - dim + slip_id + n_slips * crystal_id); - - if (!flag_light_output) - { - // Equivalent plastic strain - computed_quantities[q_point](dim + n_slips) += - inputs.solution_values[q_point]( - dim + slip_id + n_slips * crystal_id); - - // Equivalent absolute plastic strain - computed_quantities[q_point](dim + n_slips + 1) += - std::abs(inputs.solution_values[q_point]( - dim + slip_id + n_slips * crystal_id)); - - // Equivalent edge dislocation density - equivalent_edge_dislocation_density += - std::pow(inputs.solution_gradients[q_point][ - dim + slip_id + n_slips * crystal_id] * - crystals_data->get_slip_direction(crystal_id, slip_id), 2); - - // Equivalent screw dislocation density - equivalent_screw_dislocation_density += - std::pow(inputs.solution_gradients[q_point][ - dim + slip_id + n_slips * crystal_id] * - crystals_data->get_slip_orthogonal(crystal_id, slip_id), 2); - - // Plastic strain tensor - plastic_strain_tensor += - inputs.solution_values[q_point]( - dim + slip_id + n_slips * crystal_id) * - crystals_data->get_symmetrized_schmid_tensor( - crystal_id, slip_id); - } - } - } - - if (!flag_light_output) - { - strain_tensor = dealii::symmetrize(displacement_gradient) + - macroscopic_strain; + strain_tensor = dealii::symmetrize(displacement_gradient); strain_tensor_3d = convert_2d_to_3d(strain_tensor); - plastic_strain_tensor_3d = - convert_2d_to_3d(plastic_strain_tensor); - - computed_quantities[q_point](dim + n_slips + 2) = - std::sqrt(equivalent_edge_dislocation_density); - - computed_quantities[q_point](dim + n_slips + 3) = - std::sqrt(equivalent_screw_dislocation_density); - elastic_strain_tensor = - strain_tensor - plastic_strain_tensor; - - elastic_strain_tensor_3d = - convert_2d_to_3d(elastic_strain_tensor); + strain_tensor - dimensionless_number * + plastic_strain_tensor; stress_tensor = + parameters.characteristic_quantities.stress * hooke_law->get_stiffness_tetrad_3d(material_id) * convert_2d_to_3d(elastic_strain_tensor); // Von-Mises stress - computed_quantities[q_point](dim + n_slips + 4) = + computed_quantities[quadrature_point_index](dim + n_slips + 30) = + parameters.characteristic_quantities.stress * get_von_mises_stress(stress_tensor); - // Von-Mises plastic strain - computed_quantities[q_point](dim + n_slips + 5) = - get_von_mises_plastic_strain(plastic_strain_tensor); - - // Stress components - for (unsigned int i = 0; i < voigt_indices.size(); ++i) - computed_quantities[q_point](dim + n_slips + 6 + i) = - stress_tensor[voigt_indices[i].first][voigt_indices[i].second]; - - // Strain components - for (unsigned int i = 0; i < voigt_indices.size(); ++i) - computed_quantities[q_point](dim + n_slips + 12 + i) = - (i < 3 ? 1.0 : 2.0) * - strain_tensor_3d[voigt_indices[i].first][voigt_indices[i].second]; - - // Elastic strain components for (unsigned int i = 0; i < voigt_indices.size(); ++i) - computed_quantities[q_point](dim + n_slips + 18 + i) = - (i < 3 ? 1.0 : 2.0) * - elastic_strain_tensor_3d[voigt_indices[i].first][voigt_indices[i].second]; - - // Plastic strain components - for (unsigned int i = 0; i < voigt_indices.size(); ++i) - computed_quantities[q_point](dim + n_slips + 24 + i) = - (i < 3 ? 1.0 : 2.0) * - plastic_strain_tensor_3d[voigt_indices[i].first][voigt_indices[i].second]; - - if (flag_output_fluctuations) { - strain_tensor = dealii::symmetrize(displacement_gradient); - - strain_tensor_3d = convert_2d_to_3d(strain_tensor); - - elastic_strain_tensor = - strain_tensor - plastic_strain_tensor; - - stress_tensor = - hooke_law->get_stiffness_tetrad_3d(material_id) * - convert_2d_to_3d(elastic_strain_tensor); - - // Von-Mises stress - computed_quantities[q_point](dim + n_slips + 30) = - get_von_mises_stress(stress_tensor); - // Stress components - for (unsigned int i = 0; i < voigt_indices.size(); ++i) - computed_quantities[q_point](dim + n_slips + 31 + i) = - stress_tensor[voigt_indices[i].first][voigt_indices[i].second]; + computed_quantities[quadrature_point_index](dim + n_slips + 31 + i) = + parameters.characteristic_quantities.stress * + stress_tensor[voigt_indices[i].first][voigt_indices[i].second]; // Strain components - for (unsigned int i = 0; i < voigt_indices.size(); ++i) - computed_quantities[q_point](dim + n_slips + 37 + i) = - (i < 3 ? 1.0 : 2.0) * - strain_tensor_3d[voigt_indices[i].first][voigt_indices[i].second]; + computed_quantities[quadrature_point_index](dim + n_slips + 37 + i) = + (i < 3 ? 1.0 : 2.0) * + parameters.characteristic_quantities.strain * + strain_tensor_3d[voigt_indices[i].first][voigt_indices[i].second]; } } } @@ -1503,12 +1401,6 @@ void Homogenization::compute_macroscopic_stiffness_tetrad() template class gCP::Postprocessing::Postprocessor<2>; template class gCP::Postprocessing::Postprocessor<3>; -template class gCP::Postprocessing::ResidualPostprocessor<2>; -template class gCP::Postprocessing::ResidualPostprocessor<3>; - -template class gCP::Postprocessing::RatePostprocessor<2>; -template class gCP::Postprocessing::RatePostprocessor<3>; - template class gCP::Postprocessing::TrialstressPostprocessor<2>; template class gCP::Postprocessing::TrialstressPostprocessor<3>;