From 8c83cbdd66333113097fa1dae712503d803f270c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 23 Sep 2013 15:58:27 +0000 Subject: [PATCH] More editing. Return some values rather than having functions have several output arguments. git-svn-id: https://svn.dealii.org/trunk@30897 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/step-42.cc | 321 +++++++++++++--------------- 1 file changed, 153 insertions(+), 168 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index f1fd4d9db3..275c8c3623 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -22,34 +22,23 @@ */ // @sect3{Include files} -// We are using the the same -// include files as in step-41: -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +// The set of include files is not much of a surprise any more at this time: +#include +#include +#include +#include #include #include -#include -#include +#include + #include #include +#include #include #include -#include #include #include #include - #include #include #include @@ -57,28 +46,39 @@ #include #include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include + #include #include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include #include #include #include -#include -#include + #include #include -#include -#include +// This final include file provides the mkdir function +// that we will use to create a directory for output files, if necessary: #include -#include - namespace Step42 { using namespace dealii; @@ -87,8 +87,8 @@ namespace Step42 // This class has the the only purpose // to read in data from a picture file -// that has to be stored in pbm ascii -// format. This data will be bilinear +// stored in pbm ascii +// format. This data will be bilinearly // interpolated and provides in this way // a function which describes an obstacle. // @@ -116,9 +116,8 @@ namespace Step42 public: Input (const std::string &name) : - mpi_communicator(MPI_COMM_WORLD), pcout(std::cout, - (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), + (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)), obstacle_data(0), hx(0), hy(0), @@ -128,9 +127,8 @@ namespace Step42 read_obstacle(name); } - double - hv ( - int i, int j); + double hv (const int i, + const int j); double obstacle_function (const double x, @@ -140,7 +138,6 @@ namespace Step42 read_obstacle (const std::string name); private: - MPI_Comm mpi_communicator; ConditionalOStream pcout; std::vector obstacle_data; double hx, hy; @@ -187,8 +184,7 @@ namespace Step42 Vector X(4); Vector b(4); - double xx = 0.0; - double yy = 0.0; + double xx, yy; xx = ix * hx; yy = iy * hy; @@ -265,49 +261,47 @@ namespace Step42 // example we are using an elastoplastic // material behavior with linear, // isotropic hardening. -// For gamma = 0 we obtain perfect elastoplastic +// For $\gamma = 0$ we obtain perfect elastoplastic // behavior. template class ConstitutiveLaw { public: - ConstitutiveLaw ( - double _E, double _nu, double _sigma_0, double _gamma, - MPI_Comm _mpi_communicator, ConditionalOStream _pcout); + ConstitutiveLaw (const double _E, + const double _nu, + const double _sigma_0, + const double _gamma); + + bool + plast_linear_hardening (SymmetricTensor<4, dim> &stress_strain_tensor, + const SymmetricTensor<2, dim> &strain_tensor) const; void - plast_linear_hardening ( - SymmetricTensor<4, dim> &stress_strain_tensor, - const SymmetricTensor<2, dim> &strain_tensor, - unsigned int &elast_points, unsigned int &plast_points, - double &yield); - void - linearized_plast_linear_hardening ( - SymmetricTensor<4, dim> &stress_strain_tensor_linearized, + linearized_plast_linear_hardening (SymmetricTensor<4, dim> &stress_strain_tensor_linearized, SymmetricTensor<4, dim> &stress_strain_tensor, - const SymmetricTensor<2, dim> &strain_tensor); - inline SymmetricTensor<2, dim> - get_strain ( - const FEValues &fe_values, const unsigned int shape_func, + const SymmetricTensor<2, dim> &strain_tensor) const; + + SymmetricTensor<2, dim> + get_strain (const FEValues &fe_values, + const unsigned int shape_func, const unsigned int q_point) const; + void - set_sigma_0 ( - double sigma_hlp) + set_sigma_0 (double sigma_zero) { - sigma_0 = sigma_hlp; + sigma_0 = sigma_zero; } private: - SymmetricTensor<4, dim> stress_strain_tensor_mu; - SymmetricTensor<4, dim> stress_strain_tensor_kappa; - double E; - double nu; - double sigma_0; - double gamma; - double mu; - double kappa; - MPI_Comm mpi_communicator; - ConditionalOStream pcout; + const double E; + const double nu; + double sigma_0; + const double gamma; + const double mu; + const double kappa; + + const SymmetricTensor<4, dim> stress_strain_tensor_kappa; + const SymmetricTensor<4, dim> stress_strain_tensor_mu; }; // The constructor of the ConstitutiveLaw class sets the @@ -320,94 +314,84 @@ namespace Step42 // of the volumetric and deviator part. For further details // see the documentation above. template - ConstitutiveLaw::ConstitutiveLaw ( - double _E, double _nu, double _sigma_0, double _gamma, - MPI_Comm _mpi_communicator, ConditionalOStream _pcout) + ConstitutiveLaw::ConstitutiveLaw (double _E, + double _nu, + double _sigma_0, + double _gamma) : E(_E), nu(_nu), sigma_0(_sigma_0), gamma(_gamma), - mpi_communicator(_mpi_communicator), - pcout(_pcout) - { - mu = E / (2 * (1 + nu)); - kappa = E / (3 * (1 - 2 * nu)); - stress_strain_tensor_kappa = kappa + mu (E / (2 * (1 + nu))), + kappa (E / (3 * (1 - 2 * nu))), + stress_strain_tensor_kappa (kappa * outer_product(unit_symmetric_tensor(), - unit_symmetric_tensor()); - stress_strain_tensor_mu = 2 * mu + unit_symmetric_tensor())), + stress_strain_tensor_mu (2 * mu * (identity_tensor() - outer_product(unit_symmetric_tensor(), - unit_symmetric_tensor()) / 3.0); - } + unit_symmetric_tensor()) / 3.0)) + {} -// @sect3{ConstitutiveLaw::ConstitutiveLaw} +// @sect4{ConstitutiveLaw::ConstitutiveLaw} // Calculates the strain $\varepsilon(\varphi)=\dfrac{1}{2}\left(\nabla\varphi + \nabla\varphi^T$ // for the shape functions $\varphi$. template - inline SymmetricTensor<2, dim> - ConstitutiveLaw::get_strain ( - const FEValues &fe_values, const unsigned int shape_func, - const unsigned int q_point) const + SymmetricTensor<2, dim> + ConstitutiveLaw::get_strain (const FEValues &fe_values, + const unsigned int shape_func, + const unsigned int q_point) const { const FEValuesExtractors::Vector displacement(0); - SymmetricTensor<2, dim> tmp; - - tmp = fe_values[displacement].symmetric_gradient(shape_func, q_point); - - return tmp; + return fe_values[displacement].symmetric_gradient(shape_func, q_point); } -// @sect3{ConstitutiveLaw::plast_linear_hardening} +// @sect4{ConstitutiveLaw::plast_linear_hardening} // This is the implemented constitutive law. It projects the -// deviator part of the stresses in a quadrature point back to -// the yield stress plus the linear isotropic hardening. -// Also we sum up the elastic and the plastic quadrature -// points. We need this function to calculate the nonlinear +// deviatoric part of the stresses in a quadrature point back to +// the yield stress (i.e., the original yield stress $\sigma_0$ plus +// the term that describes linear isotropic hardening). +// We need this function to calculate the nonlinear // residual in // PlasticityContactProblem::residual_nl_system(TrilinosWrappers::MPI::Vector &u). +// +// The function returns whether the quadrature point is plastic to allow for +// some statistics downstream on how many of the quadrature points are +// plastic and how many are elastic. template - void - ConstitutiveLaw::plast_linear_hardening ( - SymmetricTensor<4, dim> &stress_strain_tensor, - const SymmetricTensor<2, dim> &strain_tensor, - unsigned int &elast_points, unsigned int &plast_points, double &yield) + bool + ConstitutiveLaw:: + plast_linear_hardening (SymmetricTensor<4, dim> &stress_strain_tensor, + const SymmetricTensor<2, dim> &strain_tensor) const { - if (dim == 3) - { - SymmetricTensor<2, dim> stress_tensor; - stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu) - * strain_tensor; - - SymmetricTensor<2, dim> deviator_stress_tensor = deviator( - stress_tensor); + Assert (dim == 3, ExcNotImplemented()); - double deviator_stress_tensor_norm = deviator_stress_tensor.norm(); + SymmetricTensor<2, dim> stress_tensor; + stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu) + * strain_tensor; - yield = 0; - stress_strain_tensor = stress_strain_tensor_mu; - double beta = 1.0; - if (deviator_stress_tensor_norm > sigma_0) - { - beta = sigma_0 / deviator_stress_tensor_norm; - stress_strain_tensor *= (gamma + (1 - gamma) * beta); - yield = 1; - plast_points += 1; - } - else - elast_points += 1; + const SymmetricTensor<2, dim> deviator_stress_tensor = deviator(stress_tensor); + const double deviator_stress_tensor_norm = deviator_stress_tensor.norm(); - stress_strain_tensor += stress_strain_tensor_kappa; + stress_strain_tensor = stress_strain_tensor_mu; + if (deviator_stress_tensor_norm > sigma_0) + { + const double beta = sigma_0 / deviator_stress_tensor_norm; + stress_strain_tensor *= (gamma + (1 - gamma) * beta); } + + stress_strain_tensor += stress_strain_tensor_kappa; + + return (deviator_stress_tensor_norm > sigma_0); } -// @sect3{ConstitutiveLaw::linearized_plast_linear_hardening} +// @sect4{ConstitutiveLaw::linearized_plast_linear_hardening} -// This function returns the linearized stress strain tensor -// in the solution $u^{i-1}$ of the previous Newton $i-1$ step. +// This function returns the linearized stress strain tensor, linearized +// around the solution $u^{i-1}$ of the previous Newton step $i-1$. // The parameter strain_tensor $\varepsilon(u^{i-1})$ is calculated // by $u^{i-1}$. It contains the derivative of the nonlinear // constitutive law. As the result this function returns @@ -418,39 +402,36 @@ namespace Step42 // where this function is used. template void - ConstitutiveLaw::linearized_plast_linear_hardening ( - SymmetricTensor<4, dim> &stress_strain_tensor_linearized, - SymmetricTensor<4, dim> &stress_strain_tensor, - const SymmetricTensor<2, dim> &strain_tensor) + ConstitutiveLaw:: + linearized_plast_linear_hardening (SymmetricTensor<4, dim> &stress_strain_tensor_linearized, + SymmetricTensor<4, dim> &stress_strain_tensor, + const SymmetricTensor<2, dim> &strain_tensor) const { - if (dim == 3) - { - SymmetricTensor<2, dim> stress_tensor; - stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu) - * strain_tensor; + Assert (dim == 3, ExcNotImplemented()); - SymmetricTensor<2, dim> deviator_stress_tensor = deviator( - stress_tensor); + SymmetricTensor<2, dim> stress_tensor; + stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu) + * strain_tensor; - double deviator_stress_tensor_norm = deviator_stress_tensor.norm(); + stress_strain_tensor = stress_strain_tensor_mu; + stress_strain_tensor_linearized = stress_strain_tensor_mu; - stress_strain_tensor = stress_strain_tensor_mu; - stress_strain_tensor_linearized = stress_strain_tensor_mu; - double beta = 1.0; - if (deviator_stress_tensor_norm > sigma_0) - { - beta = sigma_0 / deviator_stress_tensor_norm; - stress_strain_tensor *= (gamma + (1 - gamma) * beta); - stress_strain_tensor_linearized *= (gamma + (1 - gamma) * beta); - deviator_stress_tensor /= deviator_stress_tensor_norm; - stress_strain_tensor_linearized -= (1 - gamma) * beta * 2 * mu - * outer_product(deviator_stress_tensor, - deviator_stress_tensor); - } + SymmetricTensor<2, dim> deviator_stress_tensor = deviator(stress_tensor); + const double deviator_stress_tensor_norm = deviator_stress_tensor.norm(); - stress_strain_tensor += stress_strain_tensor_kappa; - stress_strain_tensor_linearized += stress_strain_tensor_kappa; + if (deviator_stress_tensor_norm > sigma_0) + { + const double beta = sigma_0 / deviator_stress_tensor_norm; + stress_strain_tensor *= (gamma + (1 - gamma) * beta); + stress_strain_tensor_linearized *= (gamma + (1 - gamma) * beta); + deviator_stress_tensor /= deviator_stress_tensor_norm; + stress_strain_tensor_linearized -= (1 - gamma) * beta * 2 * mu + * outer_product(deviator_stress_tensor, + deviator_stress_tensor); } + + stress_strain_tensor += stress_strain_tensor_kappa; + stress_strain_tensor_linearized += stress_strain_tensor_kappa; } namespace EquationData @@ -784,9 +765,7 @@ namespace Step42 TimerOutput::wall_times) { // double _E, double _nu, double _sigma_0, double _gamma - plast_lin_hard.reset( - new ConstitutiveLaw(e_modul, nu, sigma_0, gamma, - mpi_communicator, pcout)); + plast_lin_hard.reset(new ConstitutiveLaw(e_modul, nu, sigma_0, gamma)); degree = prm.get_integer("polynomial degree"); n_initial_refinements = prm.get_integer("number of initial refinements"); @@ -1189,10 +1168,17 @@ namespace Step42 SymmetricTensor<4, dim> stress_strain_tensor; SymmetricTensor<2, dim> stress_tensor; - plast_lin_hard->plast_linear_hardening(stress_strain_tensor, - strain_tensor[q_point], elast_points, plast_points, yield); + const bool q_point_is_plastic + = plast_lin_hard->plast_linear_hardening(stress_strain_tensor, + strain_tensor[q_point]); + if (q_point_is_plastic) + { + ++plast_points; + ++cell_constitution(cell_number); + } + else + ++elast_points; - cell_constitution(cell_number) += yield; for (unsigned int i = 0; i < dofs_per_cell; ++i) { cell_rhs(i) -= (strain_tensor[q_point] @@ -1205,8 +1191,8 @@ namespace Step42 rhs_values = 0; cell_rhs(i) += ((fe_values[displacement].value(i, q_point) * rhs_values) * fe_values.JxW(q_point)); - }; - }; + } + } for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) @@ -1216,9 +1202,8 @@ namespace Step42 { fe_values_face.reinit(cell, face); - right_hand_side.vector_value_list( - fe_values_face.get_quadrature_points(), - right_hand_side_values_face); + right_hand_side.vector_value_list(fe_values_face.get_quadrature_points(), + right_hand_side_values_face); for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) @@ -1246,7 +1231,7 @@ namespace Step42 { cell_constitution(cell_number) = 0; cell_number += 1; - }; + } cell_constitution /= n_q_points; cell_constitution.compress(VectorOperation::add); @@ -1255,9 +1240,9 @@ namespace Step42 // constraints_hanging_nodes.condense(system_rhs_lambda); - unsigned int sum_elast_points = Utilities::MPI::sum(elast_points, + const unsigned int sum_elast_points = Utilities::MPI::sum(elast_points, mpi_communicator); - unsigned int sum_plast_points = Utilities::MPI::sum(plast_points, + const unsigned int sum_plast_points = Utilities::MPI::sum(plast_points, mpi_communicator); pcout << " Number of elastic quadrature points: " << sum_elast_points << " and plastic quadrature points: " << sum_plast_points -- 2.39.5