// This class has the the only purpose
// to read in data from a picture file
- // that has to be stored as a pbm ascii
+ // that has to be stored in pbm ascii
// format. This data will be bilinear
// interpolated and provides in this way
// a function which describes an obstacle.
pcout (std::cout,
(Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
obstacle_data (0),
- lx (1.0), // length of the cube in x direction
- ly (1.0), // length of the cube in y direction
- hx (0), //
+ hx (0),
hy (0),
nx (0),
ny (0)
MPI_Comm mpi_communicator;
ConditionalOStream pcout;
std::vector<double> obstacle_data;
- double lx, ly;
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 (int i, int j)
{
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 (double x,double y)
{
return val;
}
+ // As mentioned above this function reads in the
+ // obstacle datas 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 char* name)
{
// This class provides an interface
// for a constitutive law. In this
- // example we are using an elastic
- // plastic material with linear,
+ // example we are using an elasto
+ // plastic material behavior with linear,
// isotropic hardening.
-
+ // For gamma = 0 we obtain perfect elasto
+ // plasticity behavior.
template <int dim>
class ConstitutiveLaw
{
ConditionalOStream pcout;
};
+ // The constructor of the ConstitutiveLaw class sets the
+ // required material parameter for our deformable body:
+ // E -> elastic modulus
+ // nu -> Passion's number
+ // sigma_0 -> yield stress
+ // gamma -> hardening parameter.
+ // Also it supplies the stress strain tensor of forth order
+ // 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),
stress_strain_tensor_mu = 2*mu*(identity_tensor<dim>() - outer_product(unit_symmetric_tensor<dim>(), unit_symmetric_tensor<dim>())/3.0);
}
+ // Calculates the strain for the shape functions.
template <int dim>
inline
SymmetricTensor<2,dim> ConstitutiveLaw<dim>::get_strain (const FEValues<dim> &fe_values,
return tmp;
}
+ // This is the implemented constitutive law. It projects the
+ // deviator part of the stresses in a quadrature point back to
+ // the yield stress plus the linear isotropic hardening.
+ // Also we sum up the elastic and the plastic quadrature
+ // points.
template <int dim>
void ConstitutiveLaw<dim>::plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor,
const SymmetricTensor<2,dim> &strain_tensor,
}
}
+ // This function returns the linearized stress strain tensor.
+ // It contains the derivative of the nonlinear constitutive law.
template <int dim>
void ConstitutiveLaw<dim>::linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized,
SymmetricTensor<4,dim> &stress_strain_tensor,
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>
{
if (component == 1)
return_value = 0.0;
if (component == 2)
- // if ((p(0)-0.5)*(p(0)-0.5)+(p(1)-0.5)*(p(1)-0.5) < 0.2)
- // return_value = -5000;
- // else
return_value = 0.0;
- // for (unsigned int i=0; i<dim; ++i)
- // return_value += 4*std::pow(p(i), 4);
return return_value;
}
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>
{
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).
template <int dim>
class Obstacle : public Function<dim>
{
computing_timer.exit_section("Solve");
}
+ // @sect4{PlasticityContactProblem::solve_newton}
-
+ // In this function the damped Newton method is implemented.
+ // That means two nested loops: the outer loop for the newton
+ // iteration and the inner loop for the damping steps which
+ // 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 ()
{
unsigned int number_assemble_system = 0;
for (; j<=100; j++)
{
- // Solve an elastic problem to obtain a better start solution
if (j == 1 && cycle == 0)
plast_lin_hard->set_sigma_0 (1e+10);
else if (j == 2 || cycle > 0)
TrilinosWrappers::MPI::Vector distributed_solution (system_rhs_newton);
distributed_solution = solution;
- int damped = 0;
+
+ // 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;
-// constraints.distribute (old_solution);
double a = 0;
for (unsigned int i=0; (i<10)&&(!damped); i++)
{
a=std::pow(0.5, static_cast<double>(i));
old_solution = tmp_vector;
old_solution.sadd(1-a,a, distributed_solution);
-// constraints.distribute (old_solution);
old_solution.compress (VectorOperation::add);
computing_timer.enter_section("Residual and lambda");
resid = res.l2_norm ();
if (resid<resid_old)
- damped=1;
+ damped=true;
computing_timer.exit_section("Residual and lambda");
if (num_changed==0 && resid < 1e-8)
break;
active_set_old = active_set;
- } // End of active-set-loop
+ }
pcout << "" << std::endl
<< " Number of assembled systems = " << number_assemble_system
projects / dealii-svn.git / commitdiff