// the class PlasticityContactProblem.
template <int dim>
- class Input
+ class Input
+ {
+ public:
+ Input (const std::string &name)
+ :
+ mpi_communicator(MPI_COMM_WORLD),
+ pcout(std::cout,
+ (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
+ obstacle_data(0),
+ hx(0),
+ hy(0),
+ nx(0),
+ ny(0)
{
- public:
- Input (const std::string &name)
- :
- mpi_communicator(MPI_COMM_WORLD),
- pcout(std::cout,
- (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
- obstacle_data(0),
- hx(0),
- hy(0),
- nx(0),
- ny(0)
- {
- read_obstacle(name);
- }
+ read_obstacle(name);
+ }
- double
- hv (
- int i, int j);
+ double
+ hv (
+ int i, int j);
- double
- obstacle_function (const double x,
- const double y);
+ double
+ obstacle_function (const double x,
+ const double y);
- void
- read_obstacle (const std::string name);
+ void
+ read_obstacle (const std::string name);
- private:
- MPI_Comm mpi_communicator;
- ConditionalOStream pcout;
- std::vector<double> obstacle_data;
- double hx, hy;
- int nx, ny;
- };
+ private:
+ MPI_Comm mpi_communicator;
+ ConditionalOStream pcout;
+ std::vector<double> obstacle_data;
+ double hx, hy;
+ int nx, ny;
+ };
// This function is used in obstacle_function ()
// to provide the proper value of the obstacle.
template <int dim>
- double
- Input<dim>::hv (const int i,
- const int j)
- {
- assert(i >= 0 && i < nx);
- assert(j >= 0 && j < ny);
- return obstacle_data[nx * (ny - 1 - j) + i]; // i indiziert x-werte, j indiziert y-werte
- }
+ double
+ Input<dim>::hv (const int i,
+ const int j)
+ {
+ assert(i >= 0 && i < nx);
+ assert(j >= 0 && j < ny);
+ return obstacle_data[nx * (ny - 1 - j) + i]; // i indiziert x-werte, j indiziert y-werte
+ }
// obstacle_function () calculates the bilinear interpolated
// value in the point (x,y).
template <int dim>
- double
- Input<dim>::obstacle_function (const double x,
- const double y)
- {
- int ix = (int) (x / hx);
- int iy = (int) (y / hy);
-
- if (ix < 0)
- ix = 0;
+ double
+ Input<dim>::obstacle_function (const double x,
+ const double y)
+ {
+ int ix = (int) (x / hx);
+ int iy = (int) (y / hy);
- if (iy < 0)
- iy = 0;
+ if (ix < 0)
+ ix = 0;
- if (ix >= nx - 1)
- ix = nx - 2;
+ if (iy < 0)
+ iy = 0;
- if (iy >= ny - 1)
- iy = ny - 2;
+ if (ix >= nx - 1)
+ ix = nx - 2;
- double val = 0.0;
- {
- FullMatrix<double> H(4, 4);
- Vector<double> X(4);
- Vector<double> b(4);
-
- double xx = 0.0;
- double yy = 0.0;
-
- xx = ix * hx;
- yy = iy * hy;
- H(0, 0) = xx;
- H(0, 1) = yy;
- H(0, 2) = xx * yy;
- H(0, 3) = 1.0;
- b(0) = hv(ix, iy);
-
- xx = (ix + 1) * hx;
- yy = iy * hy;
- H(1, 0) = xx;
- H(1, 1) = yy;
- H(1, 2) = xx * yy;
- H(1, 3) = 1.0;
- b(1) = hv(ix + 1, iy);
-
- xx = (ix + 1) * hx;
- yy = (iy + 1) * hy;
- H(2, 0) = xx;
- H(2, 1) = yy;
- H(2, 2) = xx * yy;
- H(2, 3) = 1.0;
- b(2) = hv(ix + 1, iy + 1);
-
- xx = ix * hx;
- yy = (iy + 1) * hy;
- H(3, 0) = xx;
- H(3, 1) = yy;
- H(3, 2) = xx * yy;
- H(3, 3) = 1.0;
- b(3) = hv(ix, iy + 1);
-
- H.gauss_jordan();
- H.vmult(X, b);
-
- val = X(0) * x + X(1) * y + X(2) * x * y + X(3);
- }
+ if (iy >= ny - 1)
+ iy = ny - 2;
- return val;
+ double val = 0.0;
+ {
+ FullMatrix<double> H(4, 4);
+ Vector<double> X(4);
+ Vector<double> b(4);
+
+ double xx = 0.0;
+ double yy = 0.0;
+
+ xx = ix * hx;
+ yy = iy * hy;
+ H(0, 0) = xx;
+ H(0, 1) = yy;
+ H(0, 2) = xx * yy;
+ H(0, 3) = 1.0;
+ b(0) = hv(ix, iy);
+
+ xx = (ix + 1) * hx;
+ yy = iy * hy;
+ H(1, 0) = xx;
+ H(1, 1) = yy;
+ H(1, 2) = xx * yy;
+ H(1, 3) = 1.0;
+ b(1) = hv(ix + 1, iy);
+
+ xx = (ix + 1) * hx;
+ yy = (iy + 1) * hy;
+ H(2, 0) = xx;
+ H(2, 1) = yy;
+ H(2, 2) = xx * yy;
+ H(2, 3) = 1.0;
+ b(2) = hv(ix + 1, iy + 1);
+
+ xx = ix * hx;
+ yy = (iy + 1) * hy;
+ H(3, 0) = xx;
+ H(3, 1) = yy;
+ H(3, 2) = xx * yy;
+ H(3, 3) = 1.0;
+ b(3) = hv(ix, iy + 1);
+
+ H.gauss_jordan();
+ H.vmult(X, b);
+
+ val = X(0) * x + X(1) * y + X(2) * x * y + X(3);
}
+ return val;
+ }
+
// As mentioned above this function reads in the
// obstacle data and stores them in the std::vector
// obstacle_data. It will be used only in run ().
template <int dim>
- void
- Input<dim>::read_obstacle (const std::string name)
- {
- std::ifstream f(name.c_str());
+ void
+ Input<dim>::read_obstacle (const std::string name)
+ {
+ std::ifstream f(name.c_str());
- std::string temp;
- f >> temp >> nx >> ny;
- assert(nx > 0 && ny > 0);
+ std::string temp;
+ f >> temp >> nx >> ny;
+ assert(nx > 0 && ny > 0);
- for (int k = 0; k < nx * ny; k++)
- {
- double val;
- f >> val;
- obstacle_data.push_back(val);
- }
+ for (int k = 0; k < nx * ny; k++)
+ {
+ double val;
+ f >> val;
+ obstacle_data.push_back(val);
+ }
- hx = 1.0 / (nx - 1);
- hy = 1.0 / (ny - 1);
+ hx = 1.0 / (nx - 1);
+ hy = 1.0 / (ny - 1);
- pcout << "Resolution of the scanned obstacle picture: " << nx << " x "
+ pcout << "Resolution of the scanned obstacle picture: " << nx << " x "
<< ny << std::endl;
- }
+ }
// @sect3{The <code>ConstitutiveLaw</code> class template}
// For gamma = 0 we obtain perfect elastoplastic
// behavior.
template <int dim>
- class ConstitutiveLaw
+ class ConstitutiveLaw
+ {
+ public:
+ ConstitutiveLaw (
+ double _E, double _nu, double _sigma_0, double _gamma,
+ MPI_Comm _mpi_communicator, ConditionalOStream _pcout);
+
+ 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,
+ SymmetricTensor<4, dim> &stress_strain_tensor,
+ const SymmetricTensor<2, dim> &strain_tensor);
+ inline SymmetricTensor<2, dim>
+ get_strain (
+ const FEValues<dim> &fe_values, const unsigned int shape_func,
+ const unsigned int q_point) const;
+ void
+ set_sigma_0 (
+ double sigma_hlp)
{
- public:
- ConstitutiveLaw (
- double _E, double _nu, double _sigma_0, double _gamma,
- MPI_Comm _mpi_communicator, ConditionalOStream _pcout);
-
- 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,
- SymmetricTensor<4, dim> &stress_strain_tensor,
- const SymmetricTensor<2, dim> &strain_tensor);
- inline SymmetricTensor<2, dim>
- get_strain (
- const FEValues<dim> &fe_values, const unsigned int shape_func,
- const unsigned int q_point) const;
- void
- set_sigma_0 (
- double sigma_hlp)
- {
- sigma_0 = sigma_hlp;
- }
+ sigma_0 = sigma_hlp;
+ }
- 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;
- };
+ 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;
+ };
// The constructor of the ConstitutiveLaw class sets the
// required material parameter for our deformable body:
// of the volumetric and deviator part. For further details
// see the documentation above.
template <int dim>
- ConstitutiveLaw<dim>::ConstitutiveLaw (
- double _E, double _nu, double _sigma_0, double _gamma,
- MPI_Comm _mpi_communicator, ConditionalOStream _pcout)
- :
- 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
- * outer_product(unit_symmetric_tensor<dim>(),
- unit_symmetric_tensor<dim>());
- stress_strain_tensor_mu = 2 * mu
- * (identity_tensor<dim>()
- - outer_product(unit_symmetric_tensor<dim>(),
- unit_symmetric_tensor<dim>()) / 3.0);
- }
+ ConstitutiveLaw<dim>::ConstitutiveLaw (
+ double _E, double _nu, double _sigma_0, double _gamma,
+ MPI_Comm _mpi_communicator, ConditionalOStream _pcout)
+ :
+ 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
+ * outer_product(unit_symmetric_tensor<dim>(),
+ unit_symmetric_tensor<dim>());
+ stress_strain_tensor_mu = 2 * mu
+ * (identity_tensor<dim>()
+ - outer_product(unit_symmetric_tensor<dim>(),
+ unit_symmetric_tensor<dim>()) / 3.0);
+ }
// @sect3{ConstitutiveLaw::ConstitutiveLaw}
// Calculates the strain $\varepsilon(\varphi)=\dfrac{1}{2}\left(\nabla\varphi + \nabla\varphi^T$
// for the shape functions $\varphi$.
template <int dim>
- inline SymmetricTensor<2, dim>
- ConstitutiveLaw<dim>::get_strain (
- const FEValues<dim> &fe_values, const unsigned int shape_func,
- const unsigned int q_point) const
- {
- const FEValuesExtractors::Vector displacement(0);
- SymmetricTensor<2, dim> tmp;
+ inline SymmetricTensor<2, dim>
+ ConstitutiveLaw<dim>::get_strain (
+ const FEValues<dim> &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);
+ tmp = fe_values[displacement].symmetric_gradient(shape_func, q_point);
- return tmp;
- }
+ return tmp;
+ }
// @sect3{ConstitutiveLaw::plast_linear_hardening}
// residual in
// PlasticityContactProblem::residual_nl_system(TrilinosWrappers::MPI::Vector &u).
template <int dim>
- void
- ConstitutiveLaw<dim>::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)
- {
- if (dim == 3)
- {
- SymmetricTensor<2, dim> stress_tensor;
- stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
- * strain_tensor;
+ void
+ ConstitutiveLaw<dim>::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)
+ {
+ 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);
+ SymmetricTensor<2, dim> deviator_stress_tensor = deviator(
+ stress_tensor);
- double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
+ double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
- 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;
+ 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;
- stress_strain_tensor += stress_strain_tensor_kappa;
- }
- }
+ stress_strain_tensor += stress_strain_tensor_kappa;
+ }
+ }
// @sect3{ConstitutiveLaw::linearized_plast_linear_hardening}
// PlasticityContactProblem::assemble_nl_system(TrilinosWrappers::MPI::Vector &u)
// where this function is used.
template <int dim>
- void
- ConstitutiveLaw<dim>::linearized_plast_linear_hardening (
- SymmetricTensor<4, dim> &stress_strain_tensor_linearized,
- SymmetricTensor<4, dim> &stress_strain_tensor,
- const SymmetricTensor<2, dim> &strain_tensor)
- {
- if (dim == 3)
- {
- SymmetricTensor<2, dim> stress_tensor;
- stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
- * strain_tensor;
+ void
+ ConstitutiveLaw<dim>::linearized_plast_linear_hardening (
+ SymmetricTensor<4, dim> &stress_strain_tensor_linearized,
+ SymmetricTensor<4, dim> &stress_strain_tensor,
+ const SymmetricTensor<2, dim> &strain_tensor)
+ {
+ 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);
+ SymmetricTensor<2, dim> deviator_stress_tensor = deviator(
+ stress_tensor);
- double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
+ double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
- 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);
- }
+ 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);
+ }
- stress_strain_tensor += stress_strain_tensor_kappa;
- stress_strain_tensor_linearized += stress_strain_tensor_kappa;
- }
- }
+ stress_strain_tensor += stress_strain_tensor_kappa;
+ stress_strain_tensor_linearized += stress_strain_tensor_kappa;
+ }
+ }
namespace EquationData
{
// It possible to apply an additional body force
// but in here it is set to zero.
template <int dim>
- class RightHandSide : public Function<dim>
+ class RightHandSide : public Function<dim>
+ {
+ public:
+ RightHandSide ()
+ :
+ Function<dim>(dim)
{
- public:
- RightHandSide ()
- :
- Function<dim>(dim)
- {
- }
+ }
- virtual double
- value (
- const Point<dim> &p, const unsigned int component = 0) const;
+ virtual double
+ value (
+ const Point<dim> &p, const unsigned int component = 0) const;
- virtual void
- vector_value (
- const Point<dim> &p, Vector<double> &values) const;
- };
+ virtual void
+ vector_value (
+ const Point<dim> &p, Vector<double> &values) const;
+ };
template <int dim>
- double
- RightHandSide<dim>::value (
- const Point<dim> &p, const unsigned int component) const
- {
- double return_value = 0.0;
+ double
+ RightHandSide<dim>::value (
+ const Point<dim> &p, const unsigned int component) const
+ {
+ double return_value = 0.0;
- if (component == 0)
- return_value = 0.0;
- if (component == 1)
- return_value = 0.0;
- if (component == 2)
- return_value = 0.0;
+ if (component == 0)
+ return_value = 0.0;
+ if (component == 1)
+ return_value = 0.0;
+ if (component == 2)
+ return_value = 0.0;
- return return_value;
- }
+ return return_value;
+ }
template <int dim>
- void
- RightHandSide<dim>::vector_value (
- const Point<dim> &p, Vector<double> &values) const
- {
- for (unsigned int c = 0; c < this->n_components; ++c)
- values(c) = RightHandSide<dim>::value(p, c);
- }
+ void
+ RightHandSide<dim>::vector_value (
+ const Point<dim> &p, Vector<double> &values) const
+ {
+ for (unsigned int c = 0; c < this->n_components; ++c)
+ values(c) = RightHandSide<dim>::value(p, c);
+ }
// This function class is used to describe the prescribed displacements
// at the boundary. But again we set this to zero.
template <int dim>
- class BoundaryValues : public Function<dim>
+ class BoundaryValues : public Function<dim>
+ {
+ public:
+ BoundaryValues ()
+ :
+ Function<dim>(dim)
{
- public:
- BoundaryValues ()
- :
- Function<dim>(dim)
- {
- }
- ;
+ }
+ ;
- virtual double
- value (
- const Point<dim> &p, const unsigned int component = 0) const;
+ virtual double
+ value (
+ const Point<dim> &p, const unsigned int component = 0) const;
- virtual void
- vector_value (
- const Point<dim> &p, Vector<double> &values) const;
- };
+ virtual void
+ vector_value (
+ const Point<dim> &p, Vector<double> &values) const;
+ };
template <int dim>
- double
- BoundaryValues<dim>::value (
- const Point<dim> &p, const unsigned int component) const
- {
- double return_value = 0;
+ double
+ BoundaryValues<dim>::value (
+ const Point<dim> &p, const unsigned int component) const
+ {
+ double return_value = 0;
- if (component == 0)
- return_value = 0.0;
- if (component == 1)
- return_value = 0.0;
- if (component == 2)
- return_value = 0.0;
+ if (component == 0)
+ return_value = 0.0;
+ if (component == 1)
+ return_value = 0.0;
+ if (component == 2)
+ return_value = 0.0;
- return return_value;
- }
+ return return_value;
+ }
template <int dim>
- void
- BoundaryValues<dim>::vector_value (
- const Point<dim> &p, Vector<double> &values) const
- {
- for (unsigned int c = 0; c < this->n_components; ++c)
- values(c) = BoundaryValues<dim>::value(p, c);
- }
+ void
+ BoundaryValues<dim>::vector_value (
+ const Point<dim> &p, Vector<double> &values) const
+ {
+ for (unsigned int c = 0; c < this->n_components; ++c)
+ values(c) = BoundaryValues<dim>::value(p, c);
+ }
// This function is obviously implemented to
// define the obstacle that penetrates our deformable
// body. You can choose between two ways to define
// your obstacle: to read it from a file or to use
// a function (here a ball).
-// z_max_domain is the z value of the surface of the work piece
+// z_max_domain is the z value of the surface of the work piece
template <int dim>
- class Obstacle : public Function<dim>
+ class Obstacle : public Function<dim>
+ {
+ public:
+ Obstacle (
+ std_cxx1x::shared_ptr<Input<dim> > const &_input,
+ bool _use_read_obstacle, double z_max_domain)
+ :
+ Function<dim>(dim),
+ input_obstacle_copy(_input),
+ use_read_obstacle(_use_read_obstacle),
+ z_max_domain(z_max_domain)
{
- public:
- Obstacle (
- std_cxx1x::shared_ptr<Input<dim> > const &_input,
- bool _use_read_obstacle, double z_max_domain)
- :
- Function<dim>(dim),
- input_obstacle_copy(_input),
- use_read_obstacle(_use_read_obstacle),
- z_max_domain(z_max_domain)
- {
- }
+ }
- virtual double
- value (
- const Point<dim> &p, const unsigned int component = 0) const;
+ virtual double
+ value (
+ const Point<dim> &p, const unsigned int component = 0) const;
- virtual void
- vector_value (
- const Point<dim> &p, Vector<double> &values) const;
+ virtual void
+ vector_value (
+ const Point<dim> &p, Vector<double> &values) const;
- private:
- std_cxx1x::shared_ptr<Input<dim> > const &input_obstacle_copy;
- bool use_read_obstacle;
- double z_max_domain;
- };
+ private:
+ std_cxx1x::shared_ptr<Input<dim> > const &input_obstacle_copy;
+ bool use_read_obstacle;
+ double z_max_domain;
+ };
template <int dim>
- double
- Obstacle<dim>::value (
- const Point<dim> &p, const unsigned int component) const
- {
- if (component == 0)
- return p(0);
- if (component == 1)
- return p(1);
+ double
+ Obstacle<dim>::value (
+ const Point<dim> &p, const unsigned int component) const
+ {
+ if (component == 0)
+ return p(0);
+ if (component == 1)
+ return p(1);
- //component==2:
- if (use_read_obstacle)
- {
- if (p(0) >= 0.0 && p(0) <= 1.0 && p(1) >= 0.0 && p(1) <= 1.0)
- return z_max_domain + 0.999
- - input_obstacle_copy->obstacle_function(p(0), p(1));
- else
- return 10000.0;
- }
- else
- {
- //sphere:
- return -std::sqrt(
- 0.36 - (p(0) - 0.5) * (p(0) - 0.5)
- - (p(1) - 0.5) * (p(1) - 0.5)) + z_max_domain + 0.59;
- }
- }
+ //component==2:
+ if (use_read_obstacle)
+ {
+ if (p(0) >= 0.0 && p(0) <= 1.0 && p(1) >= 0.0 && p(1) <= 1.0)
+ return z_max_domain + 0.999
+ - input_obstacle_copy->obstacle_function(p(0), p(1));
+ else
+ return 10000.0;
+ }
+ else
+ {
+ //sphere:
+ return -std::sqrt(
+ 0.36 - (p(0) - 0.5) * (p(0) - 0.5)
+ - (p(1) - 0.5) * (p(1) - 0.5)) + z_max_domain + 0.59;
+ }
+ }
template <int dim>
- void
- Obstacle<dim>::vector_value (
- const Point<dim> &p, Vector<double> &values) const
- {
- for (unsigned int c = 0; c < this->n_components; ++c)
- values(c) = Obstacle<dim>::value(p, c);
- }
+ void
+ Obstacle<dim>::vector_value (
+ const Point<dim> &p, Vector<double> &values) const
+ {
+ for (unsigned int c = 0; c < this->n_components; ++c)
+ values(c) = Obstacle<dim>::value(p, c);
+ }
}
// @sect3{The <code>PlasticityContactProblem</code> class template}
// operator for the constitutive law.
template <int dim>
- class PlasticityContactProblem
+ class PlasticityContactProblem
+ {
+ public:
+ PlasticityContactProblem (
+ const ParameterHandler &prm);
+ void
+ run ();
+
+ static void
+ declare (
+ ParameterHandler &prm);
+
+ private:
+ void
+ make_grid ();
+ void
+ setup_system ();
+ void
+ assemble_nl_system (
+ TrilinosWrappers::MPI::Vector &u);
+ void
+ residual_nl_system (
+ TrilinosWrappers::MPI::Vector &u);
+ void
+ assemble_mass_matrix_diagonal (
+ TrilinosWrappers::SparseMatrix &mass_matrix);
+ void
+ update_solution_and_constraints ();
+ void
+ dirichlet_constraints ();
+ void
+ solve ();
+ void
+ solve_newton ();
+ void
+ refine_grid ();
+ void
+ move_mesh (
+ const TrilinosWrappers::MPI::Vector &_complete_displacement) const;
+ void
+ output_results (
+ const std::string &title);
+ void
+ output_contact_force (
+ const unsigned int cycle);
+
+ double to_refine_factor;
+ double to_coarsen_factor;
+ unsigned int cycle;
+
+ MPI_Comm mpi_communicator;
+
+ parallel::distributed::Triangulation<dim> triangulation;
+
+ FE_Q<dim> u;
+ FESystem<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+ // We are using the SolutionTransfer class to interpolate the
+ // solution on the new refined mesh. It appears in th refine_grid()
+ // and the run() function.
+ std_cxx1x::shared_ptr<
+ parallel::distributed::SolutionTransfer<dim,
+ TrilinosWrappers::MPI::Vector> > soltrans;
+
+ IndexSet locally_owned_dofs;
+ IndexSet locally_relevant_dofs;
+
+ unsigned int number_iterations;
+
+ ConstraintMatrix constraints;
+ ConstraintMatrix constraints_hanging_nodes;
+ ConstraintMatrix constraints_dirichlet_hanging_nodes;
+
+ TrilinosWrappers::SparseMatrix system_matrix_newton;
+
+ TrilinosWrappers::MPI::Vector solution;
+ TrilinosWrappers::MPI::Vector system_rhs_newton;
+ TrilinosWrappers::MPI::Vector system_rhs_lambda;
+ TrilinosWrappers::MPI::Vector resid_vector;
+ TrilinosWrappers::MPI::Vector diag_mass_matrix_vector;
+ Vector<float> cell_constitution;
+ IndexSet active_set;
+
+ ConditionalOStream pcout;
+
+ TrilinosWrappers::PreconditionAMG::AdditionalData additional_data;
+ TrilinosWrappers::PreconditionAMG preconditioner_u;
+
+ std_cxx1x::shared_ptr<Input<dim> > input_obstacle;
+ std_cxx1x::shared_ptr<ConstitutiveLaw<dim> > plast_lin_hard;
+
+ double sigma_0; // Yield stress
+ double gamma; // Parameter for the linear isotropic hardening
+ double e_modul; // E-Modul
+ double nu; // Poisson ratio
+
+ TimerOutput computing_timer;
+
+ unsigned int degree;
+ unsigned int n_initial_refinements;
+ struct RefinementStrategy
{
- public:
- PlasticityContactProblem (
- const ParameterHandler &prm);
- void
- run ();
-
- static void
- declare (
- ParameterHandler &prm);
-
- private:
- void
- make_grid ();
- void
- setup_system ();
- void
- assemble_nl_system (
- TrilinosWrappers::MPI::Vector &u);
- void
- residual_nl_system (
- TrilinosWrappers::MPI::Vector &u);
- void
- assemble_mass_matrix_diagonal (
- TrilinosWrappers::SparseMatrix &mass_matrix);
- void
- update_solution_and_constraints ();
- void
- dirichlet_constraints ();
- void
- solve ();
- void
- solve_newton ();
- void
- refine_grid ();
- void
- move_mesh (
- const TrilinosWrappers::MPI::Vector &_complete_displacement) const;
- void
- output_results (
- const std::string &title);
- void
- output_contact_force (
- const unsigned int cycle);
-
- double to_refine_factor;
- double to_coarsen_factor;
- unsigned int cycle;
-
- MPI_Comm mpi_communicator;
-
- parallel::distributed::Triangulation<dim> triangulation;
-
- FE_Q<dim> u;
- FESystem<dim> fe;
- DoFHandler<dim> dof_handler;
-
- // We are using the SolutionTransfer class to interpolate the
- // solution on the new refined mesh. It appears in th refine_grid()
- // and the run() function.
- std_cxx1x::shared_ptr<
- parallel::distributed::SolutionTransfer<dim,
- TrilinosWrappers::MPI::Vector> > soltrans;
-
- IndexSet locally_owned_dofs;
- IndexSet locally_relevant_dofs;
-
- unsigned int number_iterations;
-
- ConstraintMatrix constraints;
- ConstraintMatrix constraints_hanging_nodes;
- ConstraintMatrix constraints_dirichlet_hanging_nodes;
-
- TrilinosWrappers::SparseMatrix system_matrix_newton;
-
- TrilinosWrappers::MPI::Vector solution;
- TrilinosWrappers::MPI::Vector system_rhs_newton;
- TrilinosWrappers::MPI::Vector system_rhs_lambda;
- TrilinosWrappers::MPI::Vector resid_vector;
- TrilinosWrappers::MPI::Vector diag_mass_matrix_vector;
- Vector<float> cell_constitution;
- IndexSet active_set;
-
- ConditionalOStream pcout;
-
- TrilinosWrappers::PreconditionAMG::AdditionalData additional_data;
- TrilinosWrappers::PreconditionAMG preconditioner_u;
-
- std_cxx1x::shared_ptr<Input<dim> > input_obstacle;
- std_cxx1x::shared_ptr<ConstitutiveLaw<dim> > plast_lin_hard;
-
- double sigma_0; // Yield stress
- double gamma; // Parameter for the linear isotropic hardening
- double e_modul; // E-Modul
- double nu; // Poisson ratio
-
- TimerOutput computing_timer;
-
- unsigned int degree;
- unsigned int n_initial_refinements;
- struct RefinementStrategy
- {
- enum value
- {
- refine_global,
- refine_percentage,
- refine_fix_dofs
- };
- };
- typename RefinementStrategy::value refinement_strategy;
- unsigned int n_cycles;
- std::string obstacle_filename;
- std::string output_dir;
- bool transfer_solution;
- std::string base_mesh;
+ enum value
+ {
+ refine_global,
+ refine_percentage,
+ refine_fix_dofs
+ };
};
+ typename RefinementStrategy::value refinement_strategy;
+ unsigned int n_cycles;
+ std::string obstacle_filename;
+ std::string output_dir;
+ bool transfer_solution;
+ std::string base_mesh;
+ };
// @sect3{Implementation of the <code>PlasticityContactProblem</code> class}
// above. As before, we will write everything
template <int dim>
- PlasticityContactProblem<dim>::PlasticityContactProblem (
- const ParameterHandler &prm)
- :
- mpi_communicator(MPI_COMM_WORLD),
- triangulation(mpi_communicator),
- u(QGaussLobatto<1>(prm.get_integer("polynomial degree") + 1)),
- fe(u, dim),
- dof_handler(triangulation),
- pcout(std::cout,
- (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
- sigma_0(400.0),
- gamma(0.01),
- e_modul(2.0e+5),
- nu(0.3),
- computing_timer(MPI_COMM_WORLD, pcout, TimerOutput::never,
- TimerOutput::wall_times)
- {
- // double _E, double _nu, double _sigma_0, double _gamma
- plast_lin_hard.reset(
- new ConstitutiveLaw<dim>(e_modul, nu, sigma_0, gamma,
- mpi_communicator, pcout));
-
- degree = prm.get_integer("polynomial degree");
- n_initial_refinements = prm.get_integer("number of initial refinements");
- std::string strat = prm.get("refinement strategy");
- if (strat == "global")
- refinement_strategy = RefinementStrategy::refine_global;
- else if (strat == "percentage")
- refinement_strategy = RefinementStrategy::refine_percentage;
- else
- throw ExcNotImplemented();
-
- n_cycles = prm.get_integer("number of cycles");
- obstacle_filename = prm.get("obstacle filename");
- output_dir = prm.get("output directory");
- if (output_dir != "" && *(output_dir.rbegin()) != '/')
- output_dir += "/";
- mkdir(output_dir.c_str(), 0777);
-
- transfer_solution = prm.get_bool("transfer solution");
- base_mesh = prm.get("base mesh");
-
- pcout << " Using output directory '" << output_dir << "'" << std::endl;
- pcout << " FE degree " << degree << std::endl;
- pcout << " Obstacle '" << obstacle_filename << "'" << std::endl;
- pcout << " transfer solution "
+ PlasticityContactProblem<dim>::PlasticityContactProblem (
+ const ParameterHandler &prm)
+ :
+ mpi_communicator(MPI_COMM_WORLD),
+ triangulation(mpi_communicator),
+ u(QGaussLobatto<1>(prm.get_integer("polynomial degree") + 1)),
+ fe(u, dim),
+ dof_handler(triangulation),
+ pcout(std::cout,
+ (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
+ sigma_0(400.0),
+ gamma(0.01),
+ e_modul(2.0e+5),
+ nu(0.3),
+ computing_timer(MPI_COMM_WORLD, pcout, TimerOutput::never,
+ TimerOutput::wall_times)
+ {
+ // double _E, double _nu, double _sigma_0, double _gamma
+ plast_lin_hard.reset(
+ new ConstitutiveLaw<dim>(e_modul, nu, sigma_0, gamma,
+ mpi_communicator, pcout));
+
+ degree = prm.get_integer("polynomial degree");
+ n_initial_refinements = prm.get_integer("number of initial refinements");
+ std::string strat = prm.get("refinement strategy");
+ if (strat == "global")
+ refinement_strategy = RefinementStrategy::refine_global;
+ else if (strat == "percentage")
+ refinement_strategy = RefinementStrategy::refine_percentage;
+ else
+ throw ExcNotImplemented();
+
+ n_cycles = prm.get_integer("number of cycles");
+ obstacle_filename = prm.get("obstacle filename");
+ output_dir = prm.get("output directory");
+ if (output_dir != "" && *(output_dir.rbegin()) != '/')
+ output_dir += "/";
+ mkdir(output_dir.c_str(), 0777);
+
+ transfer_solution = prm.get_bool("transfer solution");
+ base_mesh = prm.get("base mesh");
+
+ pcout << " Using output directory '" << output_dir << "'" << std::endl;
+ pcout << " FE degree " << degree << std::endl;
+ pcout << " Obstacle '" << obstacle_filename << "'" << std::endl;
+ pcout << " transfer solution "
<< (transfer_solution ? "true" : "false") << std::endl;
- }
+ }
// @sect4{PlasticityContactProblem::declare}
template <int dim>
- void
- PlasticityContactProblem<dim>::declare (
- ParameterHandler &prm)
- {
- prm.declare_entry("polynomial degree", "1", Patterns::Integer(),
- "polynomial degree of the FE_Q finite element space, typically 1 or 2");
- prm.declare_entry("number of initial refinements", "2",
- Patterns::Integer(),
- "number of initial global refinements before the first computation");
- prm.declare_entry("refinement strategy", "percentage",
- Patterns::Selection("global|percentage|fix dofs"),
- "refinement strategy for each cycle:\n"
- " global: one global refinement\n"
- "percentage: fixed percentage gets refined using kelly\n"
- " fix dofs: tries to achieve 2^initial_refinement*300 dofs after cycle 1 (only use 2 cycles!). Changes the coarse mesh!");
- prm.declare_entry("number of cycles", "5", Patterns::Integer(),
- "number of adaptive cycles to run");
- prm.declare_entry("obstacle filename", "", Patterns::Anything(),
- "obstacle file to read, use 'obstacle_file.pbm' or leave empty to use a sphere");
- prm.declare_entry("output directory", "", Patterns::Anything(),
- "directory to put output files (graphical output and benchmark statistics), leave empty to put into current directory");
- prm.declare_entry("transfer solution", "false", Patterns::Bool(),
- "decide if the solution should be used as a starting guess for the finer mesh, use 0 otherwise.");
- prm.declare_entry("base mesh", "box",
- Patterns::Selection("box|half sphere"),
- "select the shape of the work piece: 'box' or 'half sphere'");
+ void
+ PlasticityContactProblem<dim>::declare (
+ ParameterHandler &prm)
+ {
+ prm.declare_entry("polynomial degree", "1", Patterns::Integer(),
+ "polynomial degree of the FE_Q finite element space, typically 1 or 2");
+ prm.declare_entry("number of initial refinements", "2",
+ Patterns::Integer(),
+ "number of initial global refinements before the first computation");
+ prm.declare_entry("refinement strategy", "percentage",
+ Patterns::Selection("global|percentage|fix dofs"),
+ "refinement strategy for each cycle:\n"
+ " global: one global refinement\n"
+ "percentage: fixed percentage gets refined using kelly\n"
+ " fix dofs: tries to achieve 2^initial_refinement*300 dofs after cycle 1 (only use 2 cycles!). Changes the coarse mesh!");
+ prm.declare_entry("number of cycles", "5", Patterns::Integer(),
+ "number of adaptive cycles to run");
+ prm.declare_entry("obstacle filename", "", Patterns::Anything(),
+ "obstacle file to read, use 'obstacle_file.pbm' or leave empty to use a sphere");
+ prm.declare_entry("output directory", "", Patterns::Anything(),
+ "directory to put output files (graphical output and benchmark statistics), leave empty to put into current directory");
+ prm.declare_entry("transfer solution", "false", Patterns::Bool(),
+ "decide if the solution should be used as a starting guess for the finer mesh, use 0 otherwise.");
+ prm.declare_entry("base mesh", "box",
+ Patterns::Selection("box|half sphere"),
+ "select the shape of the work piece: 'box' or 'half sphere'");
- }
+ }
Point<3>
rotate_half_sphere (
- const Point<3> &in)
+ const Point<3> &in)
{
return Point<3>(in(2), in(1), -in(0));
}
// @sect4{PlasticityContactProblem::make_grid}
template <int dim>
- void
- PlasticityContactProblem<dim>::make_grid ()
- {
-
- if (base_mesh == "half sphere")
- {
- Point<dim> center(0, 0, 0);
- double radius = 0.8;
- GridGenerator::half_hyper_ball(triangulation, center, radius);
- GridTools::transform(&rotate_half_sphere, triangulation);
- Point<dim> shift(0.5, 0.5, 0.5);
- GridTools::shift(shift, triangulation);
- static HyperBallBoundary<dim> boundary_description(
- Point<dim>(0.5, 0.5, 0.5), radius);
- triangulation.set_boundary(0, boundary_description);
-
- triangulation.refine_global(n_initial_refinements);
-
- to_refine_factor = 0.3;
- to_coarsen_factor = 0.03;
- return;
- }
-
- Point<dim> p1(0, 0, 0);
- Point<dim> p2(1.0, 1.0, 1.0);
+ void
+ PlasticityContactProblem<dim>::make_grid ()
+ {
- GridGenerator::hyper_rectangle(triangulation, p1, p2);
- to_refine_factor = 0.3;
- to_coarsen_factor = 0.03;
+ if (base_mesh == "half sphere")
+ {
+ Point<dim> center(0, 0, 0);
+ double radius = 0.8;
+ GridGenerator::half_hyper_ball(triangulation, center, radius);
+ GridTools::transform(&rotate_half_sphere, triangulation);
+ Point<dim> shift(0.5, 0.5, 0.5);
+ GridTools::shift(shift, triangulation);
+ static HyperBallBoundary<dim> boundary_description(
+ Point<dim>(0.5, 0.5, 0.5), radius);
+ triangulation.set_boundary(0, boundary_description);
+
+ triangulation.refine_global(n_initial_refinements);
+
+ to_refine_factor = 0.3;
+ to_coarsen_factor = 0.03;
+ return;
+ }
- Triangulation<3>::active_cell_iterator cell =
- triangulation.begin_active(), endc = triangulation.end();
+ Point<dim> p1(0, 0, 0);
+ Point<dim> p2(1.0, 1.0, 1.0);
- /* boundary_indicators:
- _______
- / 1 /|
- /______ / |
- | | 8|
- | 8 | /
- |_______|/
- 6
+ GridGenerator::hyper_rectangle(triangulation, p1, p2);
+ to_refine_factor = 0.3;
+ to_coarsen_factor = 0.03;
- The boundary indicators of the sides of the cube are 8.
- The boundary indicator of the bottom is indicated with 6
- and the top with 1.
- */
+ Triangulation<3>::active_cell_iterator cell =
+ triangulation.begin_active(), endc = triangulation.end();
- for (; cell != endc; ++cell)
- for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
- ++face)
- {
- if (cell->face(face)->center()[2] == p2(2))
- cell->face(face)->set_boundary_indicator(1);
- if (cell->face(face)->center()[0] == p1(0)
- || cell->face(face)->center()[0] == p2(0)
- || cell->face(face)->center()[1] == p1(1)
- || cell->face(face)->center()[1] == p2(1))
- cell->face(face)->set_boundary_indicator(8);
- if (cell->face(face)->center()[2] == p1(2))
- cell->face(face)->set_boundary_indicator(6);
- }
+ /* boundary_indicators:
+ _______
+ / 1 /|
+ /______ / |
+ | | 8|
+ | 8 | /
+ |_______|/
+ 6
+
+ The boundary indicators of the sides of the cube are 8.
+ The boundary indicator of the bottom is indicated with 6
+ and the top with 1.
+ */
+
+ for (; cell != endc; ++cell)
+ for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+ ++face)
+ {
+ if (cell->face(face)->center()[2] == p2(2))
+ cell->face(face)->set_boundary_indicator(1);
+ if (cell->face(face)->center()[0] == p1(0)
+ || cell->face(face)->center()[0] == p2(0)
+ || cell->face(face)->center()[1] == p1(1)
+ || cell->face(face)->center()[1] == p2(1))
+ cell->face(face)->set_boundary_indicator(8);
+ if (cell->face(face)->center()[2] == p1(2))
+ cell->face(face)->set_boundary_indicator(6);
+ }
- triangulation.refine_global(n_initial_refinements);
- }
+ triangulation.refine_global(n_initial_refinements);
+ }
template <int dim>
- void
- PlasticityContactProblem<dim>::setup_system ()
+ void
+ PlasticityContactProblem<dim>::setup_system ()
+ {
+ // setup dofs
{
- // setup dofs
- {
- TimerOutput::Scope t(computing_timer, "Setup: distribute DoFs");
- dof_handler.distribute_dofs(fe);
+ TimerOutput::Scope t(computing_timer, "Setup: distribute DoFs");
+ dof_handler.distribute_dofs(fe);
- locally_owned_dofs = dof_handler.locally_owned_dofs();
- locally_relevant_dofs.clear();
- DoFTools::extract_locally_relevant_dofs(dof_handler,
- locally_relevant_dofs);
- }
+ locally_owned_dofs = dof_handler.locally_owned_dofs();
+ locally_relevant_dofs.clear();
+ DoFTools::extract_locally_relevant_dofs(dof_handler,
+ locally_relevant_dofs);
+ }
- // setup hanging nodes and Dirichlet constraints
- {
- TimerOutput::Scope t(computing_timer, "Setup: constraints");
- constraints_hanging_nodes.reinit(locally_relevant_dofs);
- DoFTools::make_hanging_node_constraints(dof_handler,
- constraints_hanging_nodes);
- constraints_hanging_nodes.close();
-
- pcout << " Number of active cells: "
- << triangulation.n_global_active_cells() << std::endl
- << " Number of degrees of freedom: " << dof_handler.n_dofs()
- << std::endl;
-
- dirichlet_constraints();
- }
+ // setup hanging nodes and Dirichlet constraints
+ {
+ TimerOutput::Scope t(computing_timer, "Setup: constraints");
+ constraints_hanging_nodes.reinit(locally_relevant_dofs);
+ DoFTools::make_hanging_node_constraints(dof_handler,
+ constraints_hanging_nodes);
+ constraints_hanging_nodes.close();
+
+ pcout << " Number of active cells: "
+ << triangulation.n_global_active_cells() << std::endl
+ << " Number of degrees of freedom: " << dof_handler.n_dofs()
+ << std::endl;
+
+ dirichlet_constraints();
+ }
- // Initialization for matrices and vectors
- {
- TimerOutput::Scope t(computing_timer, "Setup: vectors");
- solution.reinit(locally_relevant_dofs, mpi_communicator);
- system_rhs_newton.reinit(locally_owned_dofs, mpi_communicator);
- system_rhs_lambda.reinit(system_rhs_newton);
- resid_vector.reinit(system_rhs_newton);
- diag_mass_matrix_vector.reinit(system_rhs_newton);
- cell_constitution.reinit(triangulation.n_active_cells());
- active_set.clear();
- active_set.set_size(locally_relevant_dofs.size());
- }
+ // Initialization for matrices and vectors
+ {
+ TimerOutput::Scope t(computing_timer, "Setup: vectors");
+ solution.reinit(locally_relevant_dofs, mpi_communicator);
+ system_rhs_newton.reinit(locally_owned_dofs, mpi_communicator);
+ system_rhs_lambda.reinit(system_rhs_newton);
+ resid_vector.reinit(system_rhs_newton);
+ diag_mass_matrix_vector.reinit(system_rhs_newton);
+ cell_constitution.reinit(triangulation.n_active_cells());
+ active_set.clear();
+ active_set.set_size(locally_relevant_dofs.size());
+ }
- // setup sparsity pattern
- {
- TimerOutput::Scope t(computing_timer, "Setup: matrix");
- TrilinosWrappers::SparsityPattern sp(locally_owned_dofs,
- mpi_communicator);
+ // setup sparsity pattern
+ {
+ TimerOutput::Scope t(computing_timer, "Setup: matrix");
+ TrilinosWrappers::SparsityPattern sp(locally_owned_dofs,
+ mpi_communicator);
- DoFTools::make_sparsity_pattern(dof_handler, sp,
- constraints_dirichlet_hanging_nodes, false,
- Utilities::MPI::this_mpi_process(mpi_communicator));
+ DoFTools::make_sparsity_pattern(dof_handler, sp,
+ constraints_dirichlet_hanging_nodes, false,
+ Utilities::MPI::this_mpi_process(mpi_communicator));
- sp.compress();
+ sp.compress();
- system_matrix_newton.reinit(sp);
+ system_matrix_newton.reinit(sp);
- // we are going to reuse the system
- // matrix for assembling the diagonal
- // of the mass matrix so that we do not
- // need to allocate two sparse matrices
- // at the same time:
- TrilinosWrappers::SparseMatrix & mass_matrix = system_matrix_newton;
- assemble_mass_matrix_diagonal(mass_matrix);
- const unsigned int start = (system_rhs_newton.local_range().first),
- end = (system_rhs_newton.local_range().second);
- for (unsigned int j = start; j < end; j++)
- diag_mass_matrix_vector(j) = mass_matrix.diag_element(j);
+ // we are going to reuse the system
+ // matrix for assembling the diagonal
+ // of the mass matrix so that we do not
+ // need to allocate two sparse matrices
+ // at the same time:
+ TrilinosWrappers::SparseMatrix &mass_matrix = system_matrix_newton;
+ assemble_mass_matrix_diagonal(mass_matrix);
+ const unsigned int start = (system_rhs_newton.local_range().first),
+ end = (system_rhs_newton.local_range().second);
+ for (unsigned int j = start; j < end; j++)
+ diag_mass_matrix_vector(j) = mass_matrix.diag_element(j);
- number_iterations = 0;
+ number_iterations = 0;
- diag_mass_matrix_vector.compress(VectorOperation::insert);
+ diag_mass_matrix_vector.compress(VectorOperation::insert);
- // remove the mass matrix entries from the matrix:
- mass_matrix = 0;
- }
+ // remove the mass matrix entries from the matrix:
+ mass_matrix = 0;
}
+ }
template <int dim>
- void
- PlasticityContactProblem<dim>::assemble_nl_system (
- TrilinosWrappers::MPI::Vector &u)
- {
- TimerOutput::Scope t(computing_timer, "Assembling");
+ void
+ PlasticityContactProblem<dim>::assemble_nl_system (
+ TrilinosWrappers::MPI::Vector &u)
+ {
+ TimerOutput::Scope t(computing_timer, "Assembling");
- QGauss<dim> quadrature_formula(fe.degree + 1);
- QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
+ QGauss<dim> quadrature_formula(fe.degree + 1);
+ QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
- FEValues<dim> fe_values(fe, quadrature_formula,
- UpdateFlags(
- update_values | update_gradients | update_q_points
- | update_JxW_values));
+ FEValues<dim> fe_values(fe, quadrature_formula,
+ UpdateFlags(
+ update_values | update_gradients | update_q_points
+ | update_JxW_values));
- FEFaceValues<dim> fe_values_face(fe, face_quadrature_formula,
- update_values | update_quadrature_points | update_JxW_values);
+ FEFaceValues<dim> fe_values_face(fe, face_quadrature_formula,
+ update_values | update_quadrature_points | update_JxW_values);
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
- const unsigned int n_q_points = quadrature_formula.size();
- const unsigned int n_face_q_points = face_quadrature_formula.size();
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
+ const unsigned int n_face_q_points = face_quadrature_formula.size();
- const EquationData::RightHandSide<dim> right_hand_side;
- std::vector<Vector<double> > right_hand_side_values(n_q_points,
- Vector<double>(dim));
- std::vector<Vector<double> > right_hand_side_values_face(n_face_q_points,
- Vector<double>(dim));
+ const EquationData::RightHandSide<dim> right_hand_side;
+ std::vector<Vector<double> > right_hand_side_values(n_q_points,
+ Vector<double>(dim));
+ std::vector<Vector<double> > right_hand_side_values_face(n_face_q_points,
+ Vector<double>(dim));
- FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
- Vector<double> cell_rhs(dofs_per_cell);
+ FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs(dofs_per_cell);
- std::vector<unsigned int> local_dof_indices(dofs_per_cell);
+ std::vector<unsigned int> local_dof_indices(dofs_per_cell);
- typename DoFHandler<dim>::active_cell_iterator cell =
- dof_handler.begin_active(), endc = dof_handler.end();
+ typename DoFHandler<dim>::active_cell_iterator cell =
+ dof_handler.begin_active(), endc = dof_handler.end();
- const FEValuesExtractors::Vector displacement(0);
+ const FEValuesExtractors::Vector displacement(0);
- const double kappa = 1.0;
- for (; cell != endc; ++cell)
- if (cell->is_locally_owned())
- {
- fe_values.reinit(cell);
- cell_matrix = 0;
- cell_rhs = 0;
+ const double kappa = 1.0;
+ for (; cell != endc; ++cell)
+ if (cell->is_locally_owned())
+ {
+ fe_values.reinit(cell);
+ cell_matrix = 0;
+ cell_rhs = 0;
- right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
- right_hand_side_values);
+ right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
+ right_hand_side_values);
- std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
- fe_values[displacement].get_function_symmetric_gradients(u,
- strain_tensor);
+ std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
+ fe_values[displacement].get_function_symmetric_gradients(u,
+ strain_tensor);
- for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
- {
- SymmetricTensor<4, dim> stress_strain_tensor_linearized;
- SymmetricTensor<4, dim> stress_strain_tensor;
- SymmetricTensor<2, dim> stress_tensor;
+ for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+ {
+ SymmetricTensor<4, dim> stress_strain_tensor_linearized;
+ SymmetricTensor<4, dim> stress_strain_tensor;
+ SymmetricTensor<2, dim> stress_tensor;
- plast_lin_hard->linearized_plast_linear_hardening(
- stress_strain_tensor_linearized, stress_strain_tensor,
- strain_tensor[q_point]);
+ plast_lin_hard->linearized_plast_linear_hardening(
+ stress_strain_tensor_linearized, stress_strain_tensor,
+ strain_tensor[q_point]);
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- {
- stress_tensor = stress_strain_tensor_linearized
- * plast_lin_hard->get_strain(fe_values, i, q_point);
-
- for (unsigned int j = 0; j < dofs_per_cell; ++j)
- {
- cell_matrix(i, j) += (stress_tensor
- * plast_lin_hard->get_strain(fe_values, j, q_point)
- * fe_values.JxW(q_point));
- }
-
- // the linearized part a(v^i;v^i,v) of the rhs
- cell_rhs(i) += (stress_tensor * strain_tensor[q_point]
- * fe_values.JxW(q_point));
-
- // the residual part a(v^i;v) of the rhs
- cell_rhs(i) -= (strain_tensor[q_point]
- * stress_strain_tensor
- * plast_lin_hard->get_strain(fe_values, i, q_point)
- * fe_values.JxW(q_point));
-
- // the residual part F(v) of the rhs
- Tensor<1, dim> rhs_values;
- rhs_values = 0;
- cell_rhs(i) += (fe_values[displacement].value(i, q_point)
- * rhs_values * fe_values.JxW(q_point));
- }
- }
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ stress_tensor = stress_strain_tensor_linearized
+ * plast_lin_hard->get_strain(fe_values, i, q_point);
+
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ {
+ cell_matrix(i, j) += (stress_tensor
+ * plast_lin_hard->get_strain(fe_values, j, q_point)
+ * fe_values.JxW(q_point));
+ }
+
+ // the linearized part a(v^i;v^i,v) of the rhs
+ cell_rhs(i) += (stress_tensor * strain_tensor[q_point]
+ * fe_values.JxW(q_point));
+
+ // the residual part a(v^i;v) of the rhs
+ cell_rhs(i) -= (strain_tensor[q_point]
+ * stress_strain_tensor
+ * plast_lin_hard->get_strain(fe_values, i, q_point)
+ * fe_values.JxW(q_point));
+
+ // the residual part F(v) of the rhs
+ Tensor<1, dim> rhs_values;
+ 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<dim>::faces_per_cell; ++face)
+ {
+ if (cell->face(face)->at_boundary()
+ && cell->face(face)->boundary_indicator() == 1)
+ {
+ fe_values_face.reinit(cell, 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)
+ {
+ Tensor<1, dim> rhs_values;
+ rhs_values[2] = right_hand_side_values[q_point][2];
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ cell_rhs(i) += (fe_values_face[displacement].value(i,
+ q_point) * rhs_values
+ * fe_values_face.JxW(q_point));
+ }
+ }
+ }
+
+ cell->get_dof_indices(local_dof_indices);
+ constraints.distribute_local_to_global(cell_matrix, cell_rhs,
+ local_dof_indices, system_matrix_newton, system_rhs_newton,
+ true);
+
+ };
+
+ system_matrix_newton.compress(VectorOperation::add);
+ system_rhs_newton.compress(VectorOperation::add);
+ }
+
+ template <int dim>
+ void
+ PlasticityContactProblem<dim>::residual_nl_system (
+ TrilinosWrappers::MPI::Vector &u)
+ {
+ QGauss<dim> quadrature_formula(fe.degree + 1);
+ QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
+
+ FEValues<dim> fe_values(fe, quadrature_formula,
+ UpdateFlags(
+ update_values | update_gradients | update_q_points
+ | update_JxW_values));
- for (unsigned int face = 0;
- face < GeometryInfo<dim>::faces_per_cell; ++face)
- {
- if (cell->face(face)->at_boundary()
- && cell->face(face)->boundary_indicator() == 1)
- {
- fe_values_face.reinit(cell, 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)
- {
- Tensor<1, dim> rhs_values;
- rhs_values[2] = right_hand_side_values[q_point][2];
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- cell_rhs(i) += (fe_values_face[displacement].value(i,
- q_point) * rhs_values
- * fe_values_face.JxW(q_point));
- }
- }
- }
+ FEFaceValues<dim> fe_values_face(fe, face_quadrature_formula,
+ update_values | update_quadrature_points | update_JxW_values);
- cell->get_dof_indices(local_dof_indices);
- constraints.distribute_local_to_global(cell_matrix, cell_rhs,
- local_dof_indices, system_matrix_newton, system_rhs_newton,
- true);
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
+ const unsigned int n_face_q_points = face_quadrature_formula.size();
- };
+ const EquationData::RightHandSide<dim> right_hand_side;
+ std::vector<Vector<double> > right_hand_side_values(n_q_points,
+ Vector<double>(dim));
+ std::vector<Vector<double> > right_hand_side_values_face(n_face_q_points,
+ Vector<double>(dim));
- system_matrix_newton.compress(VectorOperation::add);
- system_rhs_newton.compress(VectorOperation::add);
- }
+ Vector<double> cell_rhs(dofs_per_cell);
- template <int dim>
- void
- PlasticityContactProblem<dim>::residual_nl_system (
- TrilinosWrappers::MPI::Vector &u)
- {
- QGauss<dim> quadrature_formula(fe.degree + 1);
- QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
+ std::vector<unsigned int> local_dof_indices(dofs_per_cell);
- FEValues<dim> fe_values(fe, quadrature_formula,
- UpdateFlags(
- update_values | update_gradients | update_q_points
- | update_JxW_values));
+ const FEValuesExtractors::Vector displacement(0);
- FEFaceValues<dim> fe_values_face(fe, face_quadrature_formula,
- update_values | update_quadrature_points | update_JxW_values);
+ typename DoFHandler<dim>::active_cell_iterator cell =
+ dof_handler.begin_active(), endc = dof_handler.end();
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
- const unsigned int n_q_points = quadrature_formula.size();
- const unsigned int n_face_q_points = face_quadrature_formula.size();
+ unsigned int elast_points = 0;
+ unsigned int plast_points = 0;
+ double yield = 0;
+ unsigned int cell_number = 0;
+ cell_constitution = 0;
- const EquationData::RightHandSide<dim> right_hand_side;
- std::vector<Vector<double> > right_hand_side_values(n_q_points,
- Vector<double>(dim));
- std::vector<Vector<double> > right_hand_side_values_face(n_face_q_points,
- Vector<double>(dim));
+ for (; cell != endc; ++cell)
+ if (cell->is_locally_owned())
+ {
+ fe_values.reinit(cell);
+ cell_rhs = 0;
- Vector<double> cell_rhs(dofs_per_cell);
+ right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
+ right_hand_side_values);
- std::vector<unsigned int> local_dof_indices(dofs_per_cell);
+ std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
+ fe_values[displacement].get_function_symmetric_gradients(u,
+ strain_tensor);
- const FEValuesExtractors::Vector displacement(0);
+ for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+ {
+ SymmetricTensor<4, dim> stress_strain_tensor;
+ SymmetricTensor<2, dim> stress_tensor;
- typename DoFHandler<dim>::active_cell_iterator cell =
- dof_handler.begin_active(), endc = dof_handler.end();
+ plast_lin_hard->plast_linear_hardening(stress_strain_tensor,
+ strain_tensor[q_point], elast_points, plast_points, yield);
- unsigned int elast_points = 0;
- unsigned int plast_points = 0;
- double yield = 0;
- unsigned int cell_number = 0;
- cell_constitution = 0;
+ cell_constitution(cell_number) += yield;
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ cell_rhs(i) -= (strain_tensor[q_point]
+ * stress_strain_tensor
+ * //(stress_tensor) *
+ plast_lin_hard->get_strain(fe_values, i, q_point)
+ * fe_values.JxW(q_point));
+
+ Tensor<1, dim> rhs_values;
+ rhs_values = 0;
+ cell_rhs(i) += ((fe_values[displacement].value(i, q_point)
+ * rhs_values) * fe_values.JxW(q_point));
+ };
+ };
- for (; cell != endc; ++cell)
- if (cell->is_locally_owned())
- {
- fe_values.reinit(cell);
- cell_rhs = 0;
+ for (unsigned int face = 0;
+ face < GeometryInfo<dim>::faces_per_cell; ++face)
+ {
+ if (cell->face(face)->at_boundary()
+ && cell->face(face)->boundary_indicator() == 1)
+ {
+ fe_values_face.reinit(cell, face);
- right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
- right_hand_side_values);
+ right_hand_side.vector_value_list(
+ fe_values_face.get_quadrature_points(),
+ right_hand_side_values_face);
- std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
- fe_values[displacement].get_function_symmetric_gradients(u,
- strain_tensor);
+ for (unsigned int q_point = 0; q_point < n_face_q_points;
+ ++q_point)
+ {
+ Tensor<1, dim> rhs_values;
+ rhs_values[2] = right_hand_side_values[q_point][2];
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ cell_rhs(i) += (fe_values_face[displacement].value(i,
+ q_point) * rhs_values
+ * fe_values_face.JxW(q_point));
+ }
+ }
+ }
- for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
- {
- SymmetricTensor<4, dim> stress_strain_tensor;
- SymmetricTensor<2, dim> stress_tensor;
+ cell->get_dof_indices(local_dof_indices);
+ constraints_dirichlet_hanging_nodes.distribute_local_to_global(
+ cell_rhs, local_dof_indices, system_rhs_newton);
- plast_lin_hard->plast_linear_hardening(stress_strain_tensor,
- strain_tensor[q_point], elast_points, plast_points, yield);
+ for (unsigned int i = 0; i < dofs_per_cell; i++)
+ system_rhs_lambda(local_dof_indices[i]) += cell_rhs(i);
- cell_constitution(cell_number) += yield;
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- {
- cell_rhs(i) -= (strain_tensor[q_point]
- * stress_strain_tensor
- * //(stress_tensor) *
- plast_lin_hard->get_strain(fe_values, i, q_point)
- * fe_values.JxW(q_point));
-
- Tensor<1, dim> rhs_values;
- 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<dim>::faces_per_cell; ++face)
- {
- if (cell->face(face)->at_boundary()
- && cell->face(face)->boundary_indicator() == 1)
- {
- fe_values_face.reinit(cell, 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)
- {
- Tensor<1, dim> rhs_values;
- rhs_values[2] = right_hand_side_values[q_point][2];
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- cell_rhs(i) += (fe_values_face[displacement].value(i,
- q_point) * rhs_values
- * fe_values_face.JxW(q_point));
- }
- }
- }
+ cell_number += 1;
+ }
+ else
+ {
+ cell_constitution(cell_number) = 0;
+ cell_number += 1;
+ };
- cell->get_dof_indices(local_dof_indices);
- constraints_dirichlet_hanging_nodes.distribute_local_to_global(
- cell_rhs, local_dof_indices, system_rhs_newton);
+ cell_constitution /= n_q_points;
+ cell_constitution.compress(VectorOperation::add);
+ system_rhs_newton.compress(VectorOperation::add);
+ system_rhs_lambda.compress(VectorOperation::add);
- for (unsigned int i = 0; i < dofs_per_cell; i++)
- system_rhs_lambda(local_dof_indices[i]) += cell_rhs(i);
+// constraints_hanging_nodes.condense(system_rhs_lambda);
- cell_number += 1;
- }
- else
- {
- cell_constitution(cell_number) = 0;
- cell_number += 1;
- };
-
- cell_constitution /= n_q_points;
- cell_constitution.compress(VectorOperation::add);
- system_rhs_newton.compress(VectorOperation::add);
- system_rhs_lambda.compress(VectorOperation::add);
-
-// constraints_hanging_nodes.condense(system_rhs_lambda);
-
- unsigned int sum_elast_points = Utilities::MPI::sum(elast_points,
- mpi_communicator);
- unsigned int sum_plast_points = Utilities::MPI::sum(plast_points,
- mpi_communicator);
- pcout << " Number of elastic quadrature points: " << sum_elast_points
+ unsigned int sum_elast_points = Utilities::MPI::sum(elast_points,
+ mpi_communicator);
+ 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
<< std::endl;
- }
+ }
template <int dim>
- void
- PlasticityContactProblem<dim>::assemble_mass_matrix_diagonal (
- TrilinosWrappers::SparseMatrix &mass_matrix)
- {
- QGaussLobatto<dim - 1> face_quadrature_formula(fe.degree + 1);
+ void
+ PlasticityContactProblem<dim>::assemble_mass_matrix_diagonal (
+ TrilinosWrappers::SparseMatrix &mass_matrix)
+ {
+ QGaussLobatto<dim - 1> face_quadrature_formula(fe.degree + 1);
- FEFaceValues<dim> fe_values_face(fe, face_quadrature_formula,
- update_values | update_quadrature_points | update_JxW_values);
+ FEFaceValues<dim> fe_values_face(fe, face_quadrature_formula,
+ update_values | update_quadrature_points | update_JxW_values);
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
- const unsigned int n_face_q_points = face_quadrature_formula.size();
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_face_q_points = face_quadrature_formula.size();
- FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
- Tensor<1, dim, double> ones(dim);
- for (unsigned i = 0; i < dim; i++)
- ones[i] = 1.0;
+ FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+ Tensor<1, dim, double> ones(dim);
+ for (unsigned i = 0; i < dim; i++)
+ ones[i] = 1.0;
- std::vector<unsigned int> local_dof_indices(dofs_per_cell);
+ std::vector<unsigned int> local_dof_indices(dofs_per_cell);
- const FEValuesExtractors::Vector displacement(0);
+ const FEValuesExtractors::Vector displacement(0);
- typename DoFHandler<dim>::active_cell_iterator cell =
- dof_handler.begin_active(), endc = dof_handler.end();
+ typename DoFHandler<dim>::active_cell_iterator cell =
+ dof_handler.begin_active(), endc = dof_handler.end();
- for (; cell != endc; ++cell)
- if (cell->is_locally_owned())
- for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
- ++face)
- if (cell->face(face)->at_boundary()
- && cell->face(face)->boundary_indicator() == 1)
- {
- fe_values_face.reinit(cell, face);
- cell_matrix = 0;
+ for (; cell != endc; ++cell)
+ if (cell->is_locally_owned())
+ for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+ ++face)
+ if (cell->face(face)->at_boundary()
+ && cell->face(face)->boundary_indicator() == 1)
+ {
+ fe_values_face.reinit(cell, face);
+ cell_matrix = 0;
- for (unsigned int q_point = 0; q_point < n_face_q_points;
- ++q_point)
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- cell_matrix(i, i) += (fe_values_face[displacement].value(i,
- q_point) * ones * fe_values_face.JxW(q_point));
+ for (unsigned int q_point = 0; q_point < n_face_q_points;
+ ++q_point)
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ cell_matrix(i, i) += (fe_values_face[displacement].value(i,
+ q_point) * ones * fe_values_face.JxW(q_point));
- cell->get_dof_indices(local_dof_indices);
+ cell->get_dof_indices(local_dof_indices);
-// constraints_dirichlet_hanging_nodes.distribute_local_to_global(
-// cell_matrix, local_dof_indices, mass_matrix);
+// constraints_dirichlet_hanging_nodes.distribute_local_to_global(
+// cell_matrix, local_dof_indices, mass_matrix);
- for (unsigned int i = 0; i < dofs_per_cell; i++)
- mass_matrix.add(local_dof_indices[i], local_dof_indices[i],
- cell_matrix(i, i));
- }
- mass_matrix.compress(VectorOperation::add);
- }
+ for (unsigned int i = 0; i < dofs_per_cell; i++)
+ mass_matrix.add(local_dof_indices[i], local_dof_indices[i],
+ cell_matrix(i, i));
+ }
+ mass_matrix.compress(VectorOperation::add);
+ }
// @sect4{PlasticityContactProblem::update_solution_and_constraints}
// Projection and updating of the active set
// for the dofs which penetrates the obstacle.
template <int dim>
- void
- PlasticityContactProblem<dim>::update_solution_and_constraints ()
- {
- const EquationData::Obstacle<dim> obstacle(input_obstacle,
- (obstacle_filename != ""), (base_mesh == "box" ? 1.0 : 0.5));
- std::vector<bool> vertex_touched(dof_handler.n_dofs(), false);
-
- typename DoFHandler<dim>::active_cell_iterator cell =
- dof_handler.begin_active(), endc = dof_handler.end();
-
- TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton);
- distributed_solution = solution;
- TrilinosWrappers::MPI::Vector lambda(solution);
- lambda = resid_vector;
- TrilinosWrappers::MPI::Vector diag_mass_matrix_vector_relevant(solution);
- diag_mass_matrix_vector_relevant = diag_mass_matrix_vector;
-
- constraints.reinit(locally_relevant_dofs);
- active_set.clear();
- IndexSet active_set_locally_owned;
- active_set_locally_owned.set_size(locally_owned_dofs.size());
- const double c = 100.0 * e_modul;
-
- Quadrature<dim - 1> face_quadrature(fe.get_unit_face_support_points());
- FEFaceValues<dim> fe_values_face(fe, face_quadrature,
- update_quadrature_points);
-
- const unsigned int dofs_per_face = fe.dofs_per_face;
- const unsigned int n_face_q_points = face_quadrature.size();
+ void
+ PlasticityContactProblem<dim>::update_solution_and_constraints ()
+ {
+ const EquationData::Obstacle<dim> obstacle(input_obstacle,
+ (obstacle_filename != ""), (base_mesh == "box" ? 1.0 : 0.5));
+ std::vector<bool> vertex_touched(dof_handler.n_dofs(), false);
+
+ typename DoFHandler<dim>::active_cell_iterator cell =
+ dof_handler.begin_active(), endc = dof_handler.end();
+
+ TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton);
+ distributed_solution = solution;
+ TrilinosWrappers::MPI::Vector lambda(solution);
+ lambda = resid_vector;
+ TrilinosWrappers::MPI::Vector diag_mass_matrix_vector_relevant(solution);
+ diag_mass_matrix_vector_relevant = diag_mass_matrix_vector;
+
+ constraints.reinit(locally_relevant_dofs);
+ active_set.clear();
+ IndexSet active_set_locally_owned;
+ active_set_locally_owned.set_size(locally_owned_dofs.size());
+ const double c = 100.0 * e_modul;
+
+ Quadrature<dim - 1> face_quadrature(fe.get_unit_face_support_points());
+ FEFaceValues<dim> fe_values_face(fe, face_quadrature,
+ update_quadrature_points);
+
+ const unsigned int dofs_per_face = fe.dofs_per_face;
+ const unsigned int n_face_q_points = face_quadrature.size();
+
+ // pcout<< "dofs_per_face = " << dofs_per_face
+ // << "n_face_q_points = " << n_face_q_points
+ // <<std::endl;
+ unsigned int counter_hanging_nodes = 0;
+ for (; cell != endc; ++cell)
+ if (!cell->is_artificial())
+ for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+ ++face)
+ if (cell->face(face)->at_boundary()
+ && cell->face(face)->boundary_indicator() == 1)
+ {
+ fe_values_face.reinit(cell, face);
+ std::vector<unsigned int> dof_indices(dofs_per_face);
+ cell->face(face)->get_dof_indices(dof_indices);
- // pcout<< "dofs_per_face = " << dofs_per_face
- // << "n_face_q_points = " << n_face_q_points
- // <<std::endl;
- unsigned int counter_hanging_nodes = 0;
- for (; cell != endc; ++cell)
- if (!cell->is_artificial())
- for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
- ++face)
- if (cell->face(face)->at_boundary()
- && cell->face(face)->boundary_indicator() == 1)
- {
- fe_values_face.reinit(cell, face);
- std::vector<unsigned int> dof_indices(dofs_per_face);
- cell->face(face)->get_dof_indices(dof_indices);
+ for (unsigned int q_point = 0; q_point < n_face_q_points;
+ ++q_point)
+ {
+ unsigned int component = fe.face_system_to_component_index(
+ q_point).first;
- for (unsigned int q_point = 0; q_point < n_face_q_points;
- ++q_point)
- {
- unsigned int component = fe.face_system_to_component_index(
- q_point).first;
-
- if (component == 2)
- {
- unsigned int index_z = dof_indices[q_point];
-
- if (vertex_touched[index_z] == false)
- vertex_touched[index_z] = true;
- else
- continue;
-
- // the local row where
- Point<dim> point(
- fe_values_face.quadrature_point(q_point));
-
- double obstacle_value = obstacle.value(point, 2);
- double solution_index_z = solution(index_z);
- double gap = obstacle_value - point(2);
-
- if (lambda(index_z)
- / diag_mass_matrix_vector_relevant(index_z)
- + c * (solution_index_z - gap) > 0
- && !(constraints_hanging_nodes.is_constrained(
- index_z)))
- {
- constraints.add_line(index_z);
- constraints.set_inhomogeneity(index_z, gap);
- distributed_solution(index_z) = gap;
-
- if (locally_owned_dofs.is_element(index_z))
- {
- active_set_locally_owned.add_index(index_z);
- if (locally_relevant_dofs.is_element(index_z))
- active_set.add_index(index_z);
- }
-
- }
- else if (lambda(index_z)
- / diag_mass_matrix_vector_relevant(index_z)
- + c * (solution_index_z - gap) > 0
- && constraints_hanging_nodes.is_constrained(
- index_z))
- {
- if (locally_owned_dofs.is_element(index_z))
- {
- counter_hanging_nodes += 1;
-
-// std::cout << "index_z = " << index_z
-// << ", lambda = " << lambda (index_z)
-// << ", solution_index_z - gap = " << solution_index_z - gap
-// << ", diag_mass_matrix_vector_relevant = " << diag_mass_matrix_vector_relevant (index_z)
-// << ", x = " << point(0)
-// << ", y = " << point(1)
-// << std::endl;
- }
- }
- }
- }
- }
- distributed_solution.compress(VectorOperation::insert);
+ if (component == 2)
+ {
+ unsigned int index_z = dof_indices[q_point];
+
+ if (vertex_touched[index_z] == false)
+ vertex_touched[index_z] = true;
+ else
+ continue;
+
+ // the local row where
+ Point<dim> point(
+ fe_values_face.quadrature_point(q_point));
+
+ double obstacle_value = obstacle.value(point, 2);
+ double solution_index_z = solution(index_z);
+ double gap = obstacle_value - point(2);
+
+ if (lambda(index_z)
+ / diag_mass_matrix_vector_relevant(index_z)
+ + c * (solution_index_z - gap) > 0
+ && !(constraints_hanging_nodes.is_constrained(
+ index_z)))
+ {
+ constraints.add_line(index_z);
+ constraints.set_inhomogeneity(index_z, gap);
+ distributed_solution(index_z) = gap;
+
+ if (locally_owned_dofs.is_element(index_z))
+ {
+ active_set_locally_owned.add_index(index_z);
+ if (locally_relevant_dofs.is_element(index_z))
+ active_set.add_index(index_z);
+ }
+
+ }
+ else if (lambda(index_z)
+ / diag_mass_matrix_vector_relevant(index_z)
+ + c * (solution_index_z - gap) > 0
+ && constraints_hanging_nodes.is_constrained(
+ index_z))
+ {
+ if (locally_owned_dofs.is_element(index_z))
+ {
+ counter_hanging_nodes += 1;
+
+// std::cout << "index_z = " << index_z
+// << ", lambda = " << lambda (index_z)
+// << ", solution_index_z - gap = " << solution_index_z - gap
+// << ", diag_mass_matrix_vector_relevant = " << diag_mass_matrix_vector_relevant (index_z)
+// << ", x = " << point(0)
+// << ", y = " << point(1)
+// << std::endl;
+ }
+ }
+ }
+ }
+ }
+ distributed_solution.compress(VectorOperation::insert);
- unsigned int sum_contact_constraints = Utilities::MPI::sum(
- active_set_locally_owned.n_elements(), mpi_communicator);
- pcout << " Size of active set: " << sum_contact_constraints
+ unsigned int sum_contact_constraints = Utilities::MPI::sum(
+ active_set_locally_owned.n_elements(), mpi_communicator);
+ pcout << " Size of active set: " << sum_contact_constraints
<< std::endl;
- unsigned int sum_contact_hanging_nodes = Utilities::MPI::sum(
- counter_hanging_nodes, mpi_communicator);
- pcout << " Number of hanging nodes in contact: "
+ unsigned int sum_contact_hanging_nodes = Utilities::MPI::sum(
+ counter_hanging_nodes, mpi_communicator);
+ pcout << " Number of hanging nodes in contact: "
<< sum_contact_hanging_nodes << std::endl;
- solution = distributed_solution;
+ solution = distributed_solution;
- constraints.close();
+ constraints.close();
- // constraints_dirichlet_hanging_nodes.print (std::cout);
+ // constraints_dirichlet_hanging_nodes.print (std::cout);
- constraints.merge(constraints_dirichlet_hanging_nodes);
+ constraints.merge(constraints_dirichlet_hanging_nodes);
- //constraints.print (std::cout);
- }
+ //constraints.print (std::cout);
+ }
// @sect4{PlasticityContactProblem::dirichlet_constraints}
// the Dirichlet boundary values as well as the
// hanging nodes constraints.
template <int dim>
- void
- PlasticityContactProblem<dim>::dirichlet_constraints ()
- {
- /* boundary_indicators:
- _______
- / 1 /|
- /______ / |
- 8| | 8|
- | 8 | /
- |_______|/
- 6
- */
-
- constraints_dirichlet_hanging_nodes.reinit(locally_relevant_dofs);
- constraints_dirichlet_hanging_nodes.merge(constraints_hanging_nodes);
-
- // interpolate all components of the solution
- VectorTools::interpolate_boundary_values(dof_handler,
- base_mesh == "box" ? 6 : 0, EquationData::BoundaryValues<dim>(),
- constraints_dirichlet_hanging_nodes, ComponentMask());
-
- // interpolate x- and y-components of the
- // solution (this is a bit mask, so apply
- // operator| )
- FEValuesExtractors::Scalar x_displacement(0);
- FEValuesExtractors::Scalar y_displacement(1);
- VectorTools::interpolate_boundary_values(dof_handler, 8,
- EquationData::BoundaryValues<dim>(),
- constraints_dirichlet_hanging_nodes,
- (fe.component_mask(x_displacement) | fe.component_mask(y_displacement)));
- constraints_dirichlet_hanging_nodes.close();
- }
+ void
+ PlasticityContactProblem<dim>::dirichlet_constraints ()
+ {
+ /* boundary_indicators:
+ _______
+ / 1 /|
+ /______ / |
+ 8| | 8|
+ | 8 | /
+ |_______|/
+ 6
+ */
+
+ constraints_dirichlet_hanging_nodes.reinit(locally_relevant_dofs);
+ constraints_dirichlet_hanging_nodes.merge(constraints_hanging_nodes);
+
+ // interpolate all components of the solution
+ VectorTools::interpolate_boundary_values(dof_handler,
+ base_mesh == "box" ? 6 : 0, EquationData::BoundaryValues<dim>(),
+ constraints_dirichlet_hanging_nodes, ComponentMask());
+
+ // interpolate x- and y-components of the
+ // solution (this is a bit mask, so apply
+ // operator| )
+ FEValuesExtractors::Scalar x_displacement(0);
+ FEValuesExtractors::Scalar y_displacement(1);
+ VectorTools::interpolate_boundary_values(dof_handler, 8,
+ EquationData::BoundaryValues<dim>(),
+ constraints_dirichlet_hanging_nodes,
+ (fe.component_mask(x_displacement) | fe.component_mask(y_displacement)));
+ constraints_dirichlet_hanging_nodes.close();
+ }
// @sect4{PlasticityContactProblem::solve}
// value gamma the linear system becomes
// almost semi definite but still symmetric.
template <int dim>
- void
- PlasticityContactProblem<dim>::solve ()
- {
- TimerOutput::Scope t(computing_timer, "Solve");
+ void
+ PlasticityContactProblem<dim>::solve ()
+ {
+ TimerOutput::Scope t(computing_timer, "Solve");
- TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton);
- distributed_solution = solution;
+ TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton);
+ distributed_solution = solution;
- constraints_hanging_nodes.set_zero(distributed_solution);
- constraints_hanging_nodes.set_zero(system_rhs_newton);
- distributed_solution.compress(VectorOperation::insert);
- system_rhs_newton.compress(VectorOperation::insert);
+ constraints_hanging_nodes.set_zero(distributed_solution);
+ constraints_hanging_nodes.set_zero(system_rhs_newton);
+ distributed_solution.compress(VectorOperation::insert);
+ system_rhs_newton.compress(VectorOperation::insert);
- {
- TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner");
- preconditioner_u.initialize(system_matrix_newton, additional_data);
- }
+ {
+ TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner");
+ preconditioner_u.initialize(system_matrix_newton, additional_data);
+ }
- {
- TimerOutput::Scope t(computing_timer, "Solve: iterate");
-
- PrimitiveVectorMemory<TrilinosWrappers::MPI::Vector> mem;
- TrilinosWrappers::MPI::Vector tmp(system_rhs_newton);
- // 1e-4 seems to be the fasted option altogether, but to get more
- // reproducible parallel benchmark results, we use a small residual:
- double relative_accuracy = 1e-8;
- if (output_dir.compare("its/") == 0)
- relative_accuracy = 1e-4;
-
- const double solver_tolerance = relative_accuracy
- * system_matrix_newton.residual(tmp, distributed_solution,
- system_rhs_newton);
-
- SolverControl solver_control(system_matrix_newton.m(),
- solver_tolerance);
- SolverBicgstab<TrilinosWrappers::MPI::Vector> solver(solver_control,
- mem/*,
+ {
+ TimerOutput::Scope t(computing_timer, "Solve: iterate");
+
+ PrimitiveVectorMemory<TrilinosWrappers::MPI::Vector> mem;
+ TrilinosWrappers::MPI::Vector tmp(system_rhs_newton);
+ // 1e-4 seems to be the fasted option altogether, but to get more
+ // reproducible parallel benchmark results, we use a small residual:
+ double relative_accuracy = 1e-8;
+ if (output_dir.compare("its/") == 0)
+ relative_accuracy = 1e-4;
+
+ const double solver_tolerance = relative_accuracy
+ * system_matrix_newton.residual(tmp, distributed_solution,
+ system_rhs_newton);
+
+ SolverControl solver_control(system_matrix_newton.m(),
+ solver_tolerance);
+ SolverBicgstab<TrilinosWrappers::MPI::Vector> solver(solver_control,
+ mem/*,
SolverFGMRES<TrilinosWrappers::MPI::Vector>::
AdditionalData(30, true)*/);
- solver.solve(system_matrix_newton, distributed_solution,
- system_rhs_newton, preconditioner_u);
+ solver.solve(system_matrix_newton, distributed_solution,
+ system_rhs_newton, preconditioner_u);
- pcout << " Error: " << solver_control.initial_value()
- << " -> " << solver_control.last_value() << " in "
- << solver_control.last_step() << " Bicgstab iterations."
- << std::endl;
+ pcout << " Error: " << solver_control.initial_value()
+ << " -> " << solver_control.last_value() << " in "
+ << solver_control.last_step() << " Bicgstab iterations."
+ << std::endl;
- number_iterations += solver_control.last_step();
- }
+ number_iterations += solver_control.last_step();
+ }
- constraints.distribute(distributed_solution);
+ constraints.distribute(distributed_solution);
- solution = distributed_solution;
- }
+ solution = distributed_solution;
+ }
// @sect4{PlasticityContactProblem::solve_newton}
// will be used only if necessary. To obtain a good and reasonable
// starting value we solve an elastic problem in very first step (j=1).
template <int dim>
- void
- PlasticityContactProblem<dim>::solve_newton ()
- {
- TimerOutput::Scope t(computing_timer, "solve newton setup");
+ void
+ PlasticityContactProblem<dim>::solve_newton ()
+ {
+ TimerOutput::Scope t(computing_timer, "solve newton setup");
- double resid = 0;
- double resid_old = 100000;
- TrilinosWrappers::MPI::Vector old_solution(system_rhs_newton);
- TrilinosWrappers::MPI::Vector res(system_rhs_newton);
- TrilinosWrappers::MPI::Vector tmp_vector(system_rhs_newton);
+ double resid = 0;
+ double resid_old = 100000;
+ TrilinosWrappers::MPI::Vector old_solution(system_rhs_newton);
+ TrilinosWrappers::MPI::Vector res(system_rhs_newton);
+ TrilinosWrappers::MPI::Vector tmp_vector(system_rhs_newton);
- std::vector < std::vector<bool> > constant_modes;
- DoFTools::extract_constant_modes(dof_handler, ComponentMask(),
- constant_modes);
+ std::vector < std::vector<bool> > constant_modes;
+ DoFTools::extract_constant_modes(dof_handler, ComponentMask(),
+ constant_modes);
- double sigma_hlp = sigma_0;
+ double sigma_hlp = sigma_0;
- additional_data.constant_modes = constant_modes;
- additional_data.elliptic = true;
- additional_data.n_cycles = 1;
- additional_data.w_cycle = false;
- additional_data.output_details = false;
- additional_data.smoother_sweeps = 2;
- additional_data.aggregation_threshold = 1e-2;
+ additional_data.constant_modes = constant_modes;
+ additional_data.elliptic = true;
+ additional_data.n_cycles = 1;
+ additional_data.w_cycle = false;
+ additional_data.output_details = false;
+ additional_data.smoother_sweeps = 2;
+ additional_data.aggregation_threshold = 1e-2;
- IndexSet active_set_old(active_set);
+ IndexSet active_set_old(active_set);
- t.stop(); // stop newton setup timer
+ t.stop(); // stop newton setup timer
- unsigned int j = 1;
- unsigned int number_assemble_system = 0;
- for (; j <= 100; j++)
- {
- if (transfer_solution)
- {
- if (transfer_solution && j == 1 && cycle == 0)
- plast_lin_hard->set_sigma_0(1e+10);
- else if (transfer_solution && (j == 2 || cycle > 0))
- plast_lin_hard->set_sigma_0(sigma_hlp);
- }
- else
- {
- if (j == 1)
- plast_lin_hard->set_sigma_0(1e+10);
- else
- plast_lin_hard->set_sigma_0(sigma_hlp);
- }
+ unsigned int j = 1;
+ unsigned int number_assemble_system = 0;
+ for (; j <= 100; j++)
+ {
+ if (transfer_solution)
+ {
+ if (transfer_solution && j == 1 && cycle == 0)
+ plast_lin_hard->set_sigma_0(1e+10);
+ else if (transfer_solution && (j == 2 || cycle > 0))
+ plast_lin_hard->set_sigma_0(sigma_hlp);
+ }
+ else
+ {
+ if (j == 1)
+ plast_lin_hard->set_sigma_0(1e+10);
+ else
+ plast_lin_hard->set_sigma_0(sigma_hlp);
+ }
- pcout << " " << std::endl;
- pcout << " Newton iteration " << j << std::endl;
- pcout << " Updating active set..." << std::endl;
+ pcout << " " << std::endl;
+ pcout << " Newton iteration " << j << std::endl;
+ pcout << " Updating active set..." << std::endl;
- {
- TimerOutput::Scope t(computing_timer, "update active set");
- update_solution_and_constraints();
- }
+ {
+ TimerOutput::Scope t(computing_timer, "update active set");
+ update_solution_and_constraints();
+ }
- pcout << " Assembling system... " << std::endl;
- system_matrix_newton = 0;
- system_rhs_newton = 0;
- assemble_nl_system(solution); //compute Newton-Matrix
-
- number_assemble_system += 1;
-
- pcout << " Solving system... " << std::endl;
- solve();
-
- TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton);
- distributed_solution = solution;
-
- // We handle a highly nonlinear problem so we have to damp
- // the Newtons method. We refer that we iterate the new solution
- // in each Newton step and not only the solution update.
- // Since the solution set is a convex set and not a space we
- // compute for the damping a linear combination of the
- // previous and the current solution to guarantee that the
- // damped solution is in our solution set again.
- // At most we apply 10 damping steps.
- bool damped = false;
- tmp_vector = old_solution;
- double a = 0;
- for (unsigned int i = 0; (i < 5) && (!damped); i++)
- {
- a = std::pow(0.5, static_cast<double>(i));
- old_solution = tmp_vector;
- old_solution.sadd(1 - a, a, distributed_solution);
- old_solution.compress(VectorOperation::add);
+ pcout << " Assembling system... " << std::endl;
+ system_matrix_newton = 0;
+ system_rhs_newton = 0;
+ assemble_nl_system(solution); //compute Newton-Matrix
+
+ number_assemble_system += 1;
+
+ pcout << " Solving system... " << std::endl;
+ solve();
+
+ TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton);
+ distributed_solution = solution;
+
+ // We handle a highly nonlinear problem so we have to damp
+ // the Newtons method. We refer that we iterate the new solution
+ // in each Newton step and not only the solution update.
+ // Since the solution set is a convex set and not a space we
+ // compute for the damping a linear combination of the
+ // previous and the current solution to guarantee that the
+ // damped solution is in our solution set again.
+ // At most we apply 10 damping steps.
+ bool damped = false;
+ tmp_vector = old_solution;
+ double a = 0;
+ for (unsigned int i = 0; (i < 5) && (!damped); i++)
+ {
+ a = std::pow(0.5, static_cast<double>(i));
+ old_solution = tmp_vector;
+ old_solution.sadd(1 - a, a, distributed_solution);
+ old_solution.compress(VectorOperation::add);
- TimerOutput::Scope t(computing_timer, "Residual and lambda");
+ TimerOutput::Scope t(computing_timer, "Residual and lambda");
- system_rhs_newton = 0;
- system_rhs_lambda = 0;
+ system_rhs_newton = 0;
+ system_rhs_lambda = 0;
- solution = old_solution;
- residual_nl_system(solution);
- res = system_rhs_newton;
+ solution = old_solution;
+ residual_nl_system(solution);
+ res = system_rhs_newton;
- const unsigned int start_res = (res.local_range().first),
- end_res = (res.local_range().second);
- for (unsigned int n = start_res; n < end_res; ++n)
- if (constraints.is_inhomogeneously_constrained(n))
- res(n) = 0;
+ const unsigned int start_res = (res.local_range().first),
+ end_res = (res.local_range().second);
+ for (unsigned int n = start_res; n < end_res; ++n)
+ if (constraints.is_inhomogeneously_constrained(n))
+ res(n) = 0;
- res.compress(VectorOperation::insert);
+ res.compress(VectorOperation::insert);
- resid = res.l2_norm();
+ resid = res.l2_norm();
- if (resid < resid_old)
- damped = true;
+ if (resid < resid_old)
+ damped = true;
- pcout << " Residual of the non-contact part of the system: "
+ pcout << " Residual of the non-contact part of the system: "
<< resid << std::endl
<< " with a damping parameter alpha = " << a
<< std::endl;
- // The previous iteration of step 0 is the solution of an elastic problem.
- // So a linear combination of a plastic and an elastic solution makes no sense
- // since the elastic solution is not in the convex set of the plastic solution.
- if (!transfer_solution && j == 2)
- break;
- if (transfer_solution && j == 2 && cycle == 0)
- break;
- }
+ // The previous iteration of step 0 is the solution of an elastic problem.
+ // So a linear combination of a plastic and an elastic solution makes no sense
+ // since the elastic solution is not in the convex set of the plastic solution.
+ if (!transfer_solution && j == 2)
+ break;
+ if (transfer_solution && j == 2 && cycle == 0)
+ break;
+ }
- resid_old = resid;
+ resid_old = resid;
- resid_vector = system_rhs_lambda;
- resid_vector.compress(VectorOperation::insert);
+ resid_vector = system_rhs_lambda;
+ resid_vector.compress(VectorOperation::insert);
- int is_my_set_changed = (active_set == active_set_old) ? 0 : 1;
- int num_changed = Utilities::MPI::sum(is_my_set_changed,
- MPI_COMM_WORLD);
- if (num_changed == 0)
- {
- pcout << " Active set did not change!" << std::endl;
- if (output_dir.compare("its/") != 0 && resid < 1e-7)
- break;
- else if (output_dir.compare("its/") == 0 && resid < 1e-10)
- break;
- }
- active_set_old = active_set;
- }
+ int is_my_set_changed = (active_set == active_set_old) ? 0 : 1;
+ int num_changed = Utilities::MPI::sum(is_my_set_changed,
+ MPI_COMM_WORLD);
+ if (num_changed == 0)
+ {
+ pcout << " Active set did not change!" << std::endl;
+ if (output_dir.compare("its/") != 0 && resid < 1e-7)
+ break;
+ else if (output_dir.compare("its/") == 0 && resid < 1e-10)
+ break;
+ }
+ active_set_old = active_set;
+ }
- pcout << "" << std::endl << " Number of assembled systems = "
+ pcout << "" << std::endl << " Number of assembled systems = "
<< number_assemble_system << std::endl
<< " Number of Solver-Iterations = " << number_iterations
<< std::endl;
- }
+ }
// @sect3{The <code>refine_grid</code> function}
template <int dim>
- void
- PlasticityContactProblem<dim>::refine_grid ()
- {
- if (refinement_strategy == RefinementStrategy::refine_global)
- {
- triangulation.refine_global(1);
- }
- else
- {
- Vector<float> estimated_error_per_cell(
- triangulation.n_active_cells());
- KellyErrorEstimator<dim>::estimate(dof_handler,
- QGauss<dim - 1>(fe.degree + 2), typename FunctionMap<dim>::type(),
- solution, estimated_error_per_cell);
+ void
+ PlasticityContactProblem<dim>::refine_grid ()
+ {
+ if (refinement_strategy == RefinementStrategy::refine_global)
+ {
+ triangulation.refine_global(1);
+ }
+ else
+ {
+ Vector<float> estimated_error_per_cell(
+ triangulation.n_active_cells());
+ KellyErrorEstimator<dim>::estimate(dof_handler,
+ QGauss<dim - 1>(fe.degree + 2), typename FunctionMap<dim>::type(),
+ solution, estimated_error_per_cell);
- parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(
- triangulation, estimated_error_per_cell, 0.3, 0.03);
+ parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(
+ triangulation, estimated_error_per_cell, 0.3, 0.03);
- triangulation.prepare_coarsening_and_refinement();
- if (transfer_solution)
- soltrans->prepare_for_coarsening_and_refinement(solution);
+ triangulation.prepare_coarsening_and_refinement();
+ if (transfer_solution)
+ soltrans->prepare_for_coarsening_and_refinement(solution);
- triangulation.execute_coarsening_and_refinement();
- }
- }
+ triangulation.execute_coarsening_and_refinement();
+ }
+ }
// @sect3{The <code>move_mesh</code> function}
template <int dim>
- void
- PlasticityContactProblem<dim>::move_mesh (
- const TrilinosWrappers::MPI::Vector &_complete_displacement) const
- {
- std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
+ void
+ PlasticityContactProblem<dim>::move_mesh (
+ const TrilinosWrappers::MPI::Vector &_complete_displacement) const
+ {
+ std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
- for (typename DoFHandler<dim>::active_cell_iterator cell =
- dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
- if (cell->is_locally_owned())
- for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
- ++v)
- {
- if (vertex_touched[cell->vertex_index(v)] == false)
- {
- vertex_touched[cell->vertex_index(v)] = true;
+ for (typename DoFHandler<dim>::active_cell_iterator cell =
+ dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
+ if (cell->is_locally_owned())
+ for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
+ ++v)
+ {
+ if (vertex_touched[cell->vertex_index(v)] == false)
+ {
+ vertex_touched[cell->vertex_index(v)] = true;
- Point<dim> vertex_displacement;
- for (unsigned int d = 0; d < dim; ++d)
- {
- if (_complete_displacement(cell->vertex_dof_index(v, d))
- != 0)
- vertex_displacement[d] = _complete_displacement(
- cell->vertex_dof_index(v, d));
- }
+ Point<dim> vertex_displacement;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ if (_complete_displacement(cell->vertex_dof_index(v, d))
+ != 0)
+ vertex_displacement[d] = _complete_displacement(
+ cell->vertex_dof_index(v, d));
+ }
- cell->vertex(v) += vertex_displacement;
- }
- }
- }
+ cell->vertex(v) += vertex_displacement;
+ }
+ }
+ }
// @sect4{PlasticityContactProblem::output_results}
template <int dim>
- void
- PlasticityContactProblem<dim>::output_results (
- const std::string &title)
- {
- move_mesh(solution);
-
- // Calculation of the contact forces
- TrilinosWrappers::MPI::Vector lambda(solution);
- TrilinosWrappers::MPI::Vector distributed_lambda(system_rhs_newton);
- const unsigned int start_res = (resid_vector.local_range().first),
- end_res = (resid_vector.local_range().second);
- for (unsigned int n = start_res; n < end_res; ++n)
- if (constraints.is_inhomogeneously_constrained(n))
- distributed_lambda(n) = resid_vector(n) / diag_mass_matrix_vector(n);
- distributed_lambda.compress(VectorOperation::insert);
- constraints_hanging_nodes.distribute(distributed_lambda);
- lambda = distributed_lambda;
- TrilinosWrappers::MPI::Vector resid_vector_relevant(solution);
- TrilinosWrappers::MPI::Vector distributed_resid_vector(resid_vector);
- constraints_hanging_nodes.distribute(distributed_resid_vector);
- resid_vector_relevant = distributed_resid_vector;
-
- DataOut<dim> data_out;
-
- data_out.attach_dof_handler(dof_handler);
-
- const std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation(
- dim, DataComponentInterpretation::component_is_part_of_vector);
- data_out.add_data_vector(solution,
- std::vector < std::string > (dim, "Displacement"),
- DataOut<dim>::type_dof_data, data_component_interpretation);
- data_out.add_data_vector(lambda,
- std::vector < std::string > (dim, "ContactForce"),
- DataOut<dim>::type_dof_data, data_component_interpretation);
- data_out.add_data_vector(active_set,
- std::vector < std::string > (dim, "ActiveSet"),
- DataOut<dim>::type_dof_data, data_component_interpretation);
- data_out.add_data_vector(resid_vector_relevant,
- std::vector < std::string > (dim, "Residual"),
- DataOut<dim>::type_dof_data, data_component_interpretation);
-
- Vector<float> subdomain(triangulation.n_active_cells());
- for (unsigned int i = 0; i < subdomain.size(); ++i)
- subdomain(i) = triangulation.locally_owned_subdomain();
- data_out.add_data_vector(subdomain, "subdomain");
-
- data_out.add_data_vector(cell_constitution, "CellConstitution");
-
- data_out.build_patches();
-
- const std::string filename =
- (output_dir + title + "-"
- + Utilities::int_to_string(
- triangulation.locally_owned_subdomain(), 4));
-
- std::ofstream output_vtu((filename + ".vtu").c_str());
- data_out.write_vtu(output_vtu);
- pcout << output_dir + title << ".pvtu" << std::endl;
-
- if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
- {
- std::vector < std::string > filenames;
- for (unsigned int i = 0;
- i < Utilities::MPI::n_mpi_processes(mpi_communicator); ++i)
- filenames.push_back(
- title + "-" + Utilities::int_to_string(i, 4) + ".vtu");
-
- std::ofstream master_output((output_dir + title + ".pvtu").c_str());
- data_out.write_pvtu_record(master_output, filenames);
- }
+ void
+ PlasticityContactProblem<dim>::output_results (
+ const std::string &title)
+ {
+ move_mesh(solution);
+
+ // Calculation of the contact forces
+ TrilinosWrappers::MPI::Vector lambda(solution);
+ TrilinosWrappers::MPI::Vector distributed_lambda(system_rhs_newton);
+ const unsigned int start_res = (resid_vector.local_range().first),
+ end_res = (resid_vector.local_range().second);
+ for (unsigned int n = start_res; n < end_res; ++n)
+ if (constraints.is_inhomogeneously_constrained(n))
+ distributed_lambda(n) = resid_vector(n) / diag_mass_matrix_vector(n);
+ distributed_lambda.compress(VectorOperation::insert);
+ constraints_hanging_nodes.distribute(distributed_lambda);
+ lambda = distributed_lambda;
+ TrilinosWrappers::MPI::Vector resid_vector_relevant(solution);
+ TrilinosWrappers::MPI::Vector distributed_resid_vector(resid_vector);
+ constraints_hanging_nodes.distribute(distributed_resid_vector);
+ resid_vector_relevant = distributed_resid_vector;
+
+ DataOut<dim> data_out;
+
+ data_out.attach_dof_handler(dof_handler);
+
+ const std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation(
+ dim, DataComponentInterpretation::component_is_part_of_vector);
+ data_out.add_data_vector(solution,
+ std::vector < std::string > (dim, "Displacement"),
+ DataOut<dim>::type_dof_data, data_component_interpretation);
+ data_out.add_data_vector(lambda,
+ std::vector < std::string > (dim, "ContactForce"),
+ DataOut<dim>::type_dof_data, data_component_interpretation);
+ data_out.add_data_vector(active_set,
+ std::vector < std::string > (dim, "ActiveSet"),
+ DataOut<dim>::type_dof_data, data_component_interpretation);
+ data_out.add_data_vector(resid_vector_relevant,
+ std::vector < std::string > (dim, "Residual"),
+ DataOut<dim>::type_dof_data, data_component_interpretation);
+
+ Vector<float> subdomain(triangulation.n_active_cells());
+ for (unsigned int i = 0; i < subdomain.size(); ++i)
+ subdomain(i) = triangulation.locally_owned_subdomain();
+ data_out.add_data_vector(subdomain, "subdomain");
+
+ data_out.add_data_vector(cell_constitution, "CellConstitution");
+
+ data_out.build_patches();
+
+ const std::string filename =
+ (output_dir + title + "-"
+ + Utilities::int_to_string(
+ triangulation.locally_owned_subdomain(), 4));
+
+ std::ofstream output_vtu((filename + ".vtu").c_str());
+ data_out.write_vtu(output_vtu);
+ pcout << output_dir + title << ".pvtu" << std::endl;
+
+ if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+ {
+ std::vector < std::string > filenames;
+ for (unsigned int i = 0;
+ i < Utilities::MPI::n_mpi_processes(mpi_communicator); ++i)
+ filenames.push_back(
+ title + "-" + Utilities::int_to_string(i, 4) + ".vtu");
+
+ std::ofstream master_output((output_dir + title + ".pvtu").c_str());
+ data_out.write_pvtu_record(master_output, filenames);
+ }
- TrilinosWrappers::MPI::Vector tmp(solution);
- tmp *= -1;
- move_mesh(tmp);
- }
+ TrilinosWrappers::MPI::Vector tmp(solution);
+ tmp *= -1;
+ move_mesh(tmp);
+ }
// @sect4{PlasticityContactProblem::output_contact_force}
// a processor does not own the cell with the point we have to
// catch these cases.
template <int dim>
- void
- PlasticityContactProblem<dim>::output_contact_force (
- const unsigned int cycle)
+ void
+ PlasticityContactProblem<dim>::output_contact_force (
+ const unsigned int cycle)
+ {
+ Functions::FEFieldFunction<dim, DoFHandler<dim>,
+ TrilinosWrappers::MPI::Vector> solution_function(dof_handler,
+ solution);
+ std::cout.precision(10);
+
+ Vector<double> solution_p1(dim);
+ std::vector<Tensor<1, dim> > solution_gradient_p1(dim);
+
+ // Here we calculate the contact pressure as a vector lambda.
+ // If a dof is element of the active set lambda contains the
+ // nonlinear residual this dof divided by the according entry
+ // of the mass matrix. In all other dofs lambda will be set to
+ // zero.
+ TrilinosWrappers::MPI::Vector lambda(solution);
+ TrilinosWrappers::MPI::Vector distributed_lambda(system_rhs_newton);
+ const unsigned int start_res = (resid_vector.local_range().first),
+ end_res = (resid_vector.local_range().second);
+ for (unsigned int n = start_res; n < end_res; ++n)
+ if (constraints.is_inhomogeneously_constrained(n))
+ distributed_lambda(n) = resid_vector(n) / diag_mass_matrix_vector(n);
+ else
+ distributed_lambda(n) = 0;
+ distributed_lambda.compress(VectorOperation::insert);
+ constraints_hanging_nodes.distribute(distributed_lambda);
+ lambda = distributed_lambda;
+ Functions::FEFieldFunction<dim, DoFHandler<dim>,
+ TrilinosWrappers::MPI::Vector> lambda_function(dof_handler, lambda);
+
+ // Here we try to find the MPI-process which owns the cell
+ // with the point_of_interest. If it is the wrong MPI-process
+ // we catch this case and set point_found to false.
+ const Point<dim> point_of_interest(0.49, 0.5001, 1.0);
+ Vector<double> contact_pressure_in_point(dim);
+ bool point_found = true;
+
+ MPI_Barrier(MPI_COMM_WORLD);
+ try
+ {
+ lambda_function.vector_value(point_of_interest,
+ contact_pressure_in_point);
+ }
+ catch (const typename Functions::FEFieldFunction<dim, DoFHandler<dim>,
+ TrilinosWrappers::MPI::Vector>::ExcPointNotAvailableHere &)
+ {
+ point_found = false;
+ }
+
+ if (point_found == true)
+ {
+ std::cout << "PoI contact pressure: " << contact_pressure_in_point(2)
+ << std::endl;
+ }
+
+ // To obtain the contact force we have to compute an integral of the contact pressure
+ // in z-direction over the whole contact area. To be accurate enough we use the
+ // Gaussian quadrature rule with fe.degree + 1.
+ double contact_force = 0.0;
{
- Functions::FEFieldFunction<dim, DoFHandler<dim>,
- TrilinosWrappers::MPI::Vector> solution_function(dof_handler,
- solution);
- std::cout.precision(10);
-
- Vector<double> solution_p1(dim);
- std::vector<Tensor<1, dim> > solution_gradient_p1(dim);
-
- // Here we calculate the contact pressure as a vector lambda.
- // If a dof is element of the active set lambda contains the
- // nonlinear residual this dof divided by the according entry
- // of the mass matrix. In all other dofs lambda will be set to
- // zero.
- TrilinosWrappers::MPI::Vector lambda(solution);
- TrilinosWrappers::MPI::Vector distributed_lambda(system_rhs_newton);
- const unsigned int start_res = (resid_vector.local_range().first),
- end_res = (resid_vector.local_range().second);
- for (unsigned int n = start_res; n < end_res; ++n)
- if (constraints.is_inhomogeneously_constrained(n))
- distributed_lambda(n) = resid_vector(n) / diag_mass_matrix_vector(n);
- else
- distributed_lambda(n) = 0;
- distributed_lambda.compress(VectorOperation::insert);
- constraints_hanging_nodes.distribute(distributed_lambda);
- lambda = distributed_lambda;
- Functions::FEFieldFunction<dim, DoFHandler<dim>,
- TrilinosWrappers::MPI::Vector> lambda_function(dof_handler, lambda);
-
- // Here we try to find the MPI-process which owns the cell
- // with the point_of_interest. If it is the wrong MPI-process
- // we catch this case and set point_found to false.
- const Point<dim> point_of_interest(0.49, 0.5001, 1.0);
- Vector<double> contact_pressure_in_point(dim);
- bool point_found = true;
-
- MPI_Barrier(MPI_COMM_WORLD);
- try
- {
- lambda_function.vector_value(point_of_interest,
- contact_pressure_in_point);
- }
- catch (const typename Functions::FEFieldFunction<dim, DoFHandler<dim>,
- TrilinosWrappers::MPI::Vector>::ExcPointNotAvailableHere &)
- {
- point_found = false;
- }
+ QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
- if (point_found == true)
- {
- std::cout << "PoI contact pressure: " << contact_pressure_in_point(2)
- << std::endl;
- }
+ FEFaceValues<dim> fe_values_face(fe, face_quadrature_formula,
+ update_values | update_quadrature_points | update_JxW_values);
- // To obtain the contact force we have to compute an integral of the contact pressure
- // in z-direction over the whole contact area. To be accurate enough we use the
- // Gaussian quadrature rule with fe.degree + 1.
- double contact_force = 0.0;
- {
- QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
+ const unsigned int n_face_q_points = face_quadrature_formula.size();
- FEFaceValues<dim> fe_values_face(fe, face_quadrature_formula,
- update_values | update_quadrature_points | update_JxW_values);
+ const FEValuesExtractors::Vector displacement(0);
- const unsigned int n_face_q_points = face_quadrature_formula.size();
+ typename DoFHandler<dim>::active_cell_iterator cell =
+ dof_handler.begin_active(), endc = dof_handler.end();
+ for (; cell != endc; ++cell)
+ if (cell->is_locally_owned())
+ for (unsigned int face = 0;
+ face < GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face(face)->at_boundary()
+ && cell->face(face)->boundary_indicator() == 1)
+ {
+ fe_values_face.reinit(cell, face);
- const FEValuesExtractors::Vector displacement(0);
+ std::vector<Tensor<1, dim> > lambda_values(n_face_q_points);
+ fe_values_face[displacement].get_function_values(lambda,
+ lambda_values);
- typename DoFHandler<dim>::active_cell_iterator cell =
- dof_handler.begin_active(), endc = dof_handler.end();
- for (; cell != endc; ++cell)
- if (cell->is_locally_owned())
- for (unsigned int face = 0;
- face < GeometryInfo<dim>::faces_per_cell; ++face)
- if (cell->face(face)->at_boundary()
- && cell->face(face)->boundary_indicator() == 1)
+ for (unsigned int q_point = 0; q_point < n_face_q_points;
+ ++q_point)
{
- fe_values_face.reinit(cell, face);
-
- std::vector<Tensor<1, dim> > lambda_values(n_face_q_points);
- fe_values_face[displacement].get_function_values(lambda,
- lambda_values);
-
- for (unsigned int q_point = 0; q_point < n_face_q_points;
- ++q_point)
- {
- contact_force += lambda_values[q_point][2]
- * fe_values_face.JxW(q_point);
- }
+ contact_force += lambda_values[q_point][2]
+ * fe_values_face.JxW(q_point);
}
- contact_force = Utilities::MPI::sum(contact_force, MPI_COMM_WORLD);
- pcout << "Contact force = " << contact_force << std::endl;
- }
- MPI_Barrier(MPI_COMM_WORLD);
+ }
+ contact_force = Utilities::MPI::sum(contact_force, MPI_COMM_WORLD);
+ pcout << "Contact force = " << contact_force << std::endl;
}
+ MPI_Barrier(MPI_COMM_WORLD);
+ }
// @sect4{PlasticityContactProblem::run}
template <int dim>
- void
- PlasticityContactProblem<dim>::run ()
- {
+ void
+ PlasticityContactProblem<dim>::run ()
+ {
- if (obstacle_filename != "")
- {
- pcout << "Read the obstacle from '" << obstacle_filename << "' ... "
+ if (obstacle_filename != "")
+ {
+ pcout << "Read the obstacle from '" << obstacle_filename << "' ... "
<< std::flush;
- input_obstacle.reset(new Input<dim>(obstacle_filename.c_str()));
- pcout << "done." << std::endl;
- }
+ input_obstacle.reset(new Input<dim>(obstacle_filename.c_str()));
+ pcout << "done." << std::endl;
+ }
- computing_timer.reset();
- for (cycle = 0; cycle < n_cycles; ++cycle)
+ computing_timer.reset();
+ for (cycle = 0; cycle < n_cycles; ++cycle)
+ {
{
- {
- TimerOutput::Scope t(computing_timer, "Setup");
+ TimerOutput::Scope t(computing_timer, "Setup");
- pcout << std::endl;
- pcout << "Cycle " << cycle << ':' << std::endl;
-
- if (cycle == 0)
- {
- make_grid();
- }
- else
- {
- TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
- if (transfer_solution)
- soltrans.reset(
- new parallel::distributed::SolutionTransfer<dim,
- TrilinosWrappers::MPI::Vector>(dof_handler));
- refine_grid();
- }
+ pcout << std::endl;
+ pcout << "Cycle " << cycle << ':' << std::endl;
- setup_system();
+ if (cycle == 0)
+ {
+ make_grid();
+ }
+ else
+ {
+ TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
+ if (transfer_solution)
+ soltrans.reset(
+ new parallel::distributed::SolutionTransfer<dim,
+ TrilinosWrappers::MPI::Vector>(dof_handler));
+ refine_grid();
+ }
- if (transfer_solution && cycle > 0)
- {
- TrilinosWrappers::MPI::Vector distributed_solution(
- system_rhs_newton);
- distributed_solution = solution;
- soltrans->interpolate(distributed_solution);
- solution = distributed_solution;
- residual_nl_system(solution);
- resid_vector = system_rhs_lambda;
- resid_vector.compress(VectorOperation::insert);
- }
+ setup_system();
+ if (transfer_solution && cycle > 0)
+ {
+ TrilinosWrappers::MPI::Vector distributed_solution(
+ system_rhs_newton);
+ distributed_solution = solution;
+ soltrans->interpolate(distributed_solution);
+ solution = distributed_solution;
+ residual_nl_system(solution);
+ resid_vector = system_rhs_lambda;
+ resid_vector.compress(VectorOperation::insert);
}
- solve_newton();
+ }
- if (true) //Utilities::MPI::n_mpi_processes(mpi_communicator) <= 64)
- {
- pcout << " Writing graphical output... " << std::flush;
+ solve_newton();
- TimerOutput::Scope t(computing_timer, "Graphical output");
+ if (true) //Utilities::MPI::n_mpi_processes(mpi_communicator) <= 64)
+ {
+ pcout << " Writing graphical output... " << std::flush;
- std::ostringstream filename_solution;
- filename_solution << "solution-";
- filename_solution << Utilities::int_to_string(cycle, 2);
- output_results(filename_solution.str());
- }
+ TimerOutput::Scope t(computing_timer, "Graphical output");
+
+ std::ostringstream filename_solution;
+ filename_solution << "solution-";
+ filename_solution << Utilities::int_to_string(cycle, 2);
+ output_results(filename_solution.str());
+ }
- computing_timer.print_summary();
- computing_timer.reset();
+ computing_timer.print_summary();
+ computing_timer.reset();
- Utilities::System::MemoryStats stats;
- Utilities::System::get_memory_stats(stats);
- pcout << "VMPEAK, Resident in kB: " << stats.VmSize << " "
+ Utilities::System::MemoryStats stats;
+ Utilities::System::get_memory_stats(stats);
+ pcout << "VMPEAK, Resident in kB: " << stats.VmSize << " "
<< stats.VmRSS << std::endl;
- if (base_mesh == "box")
- output_contact_force(cycle);
- }
- }
+ if (base_mesh == "box")
+ output_contact_force(cycle);
+ }
+ }
}
// @sect3{The <code>main</code> function}
int
main (
- int argc, char *argv[])
+ int argc, char *argv[])
{
using namespace dealii;
using namespace Step42;
prm.read_input(argv[1]);
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
- {
- PlasticityContactProblem<3> problem(prm);
- problem.run();
- }
+ {
+ PlasticityContactProblem<3> problem(prm);
+ problem.run();
+ }
return 0;
}
projects / dealii-svn.git / commitdiff