![]() |
deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
|
This program was contributed by Seyed Shahram Ghorashi <s.sh.ghorashi@gmail.com>.
It comes without any warranty or support by its authors or the authors of deal.II.
This program is part of the deal.II code gallery and consists of the following files (click to inspect):
The code deals with solving an elastoplasticity problem with linear isotropic hardening. At each load/displacement step, the error based on a prescribed quantity of interest (Goal-oriented error estimation) is computed by using the dual-weighted residual method.
Based on a prescribed error bound and estimated elementwise errors, the mesh is then refined/coarsened. Afterwards, the solution is projected to the new mesh and the analysis process is repeated.
The applied methodology and the solved numerical examples can be found in the following paper:
Ghorashi SSh, Rabczuk T.: "Goal-Oriented Error Estimation and Mesh Adaptivity in 3d Elastoplasticity Problems". International Journal of Fracture. Accepted. 2016.
The code contains several examples consisting of the three examples presented in the aforementioned paper. To run each of them you can switch to them, e.g.
Then by compiling it using the following commands
you can run it by typing
or
or
The three named input files are already in the current directory, and have the following options:
The set of include files is not much of a surprise any more at this time:
And here the only two new things among the header files: an include file in which symmetric tensors of rank 2 and 4 are implemented, as introduced in the introduction:
And a header that implements filters for iterators looping over all cells. We will use this when selecting only those cells for output that are owned by the present process in a parallel program:
This final include file provides the mkdir
function that we will use to create a directory for output files, if necessary:
Assert (input.n_levels() == 1, ExcMessage ("The input triangulations must be coarse meshes."));
if (dim == 2) { von_Mises_stress = std::sqrt( stress[0][0]*stress[0][0]
"Perforated_strip_tension" plane stress const double von_Mises_stress = std::sqrt( stress[0][0]*stress[0][0]
otherwise plane strain / 3d case
PointHistory
classAs was mentioned in the introduction, we have to store the old stress in quadrature point so that we can compute the residual forces at this point during the next time step. This alone would not warrant a structure with only one member, but in more complicated applications, we would have to store more information in quadrature points as well, such as the history variables of plasticity, etc. In essence, we have to store everything that affects the present state of the material here, which in plasticity is determined by the deformation history variables.
We will not give this class any meaningful functionality beyond being able to store data, i.e. there are no constructors, destructors, or other member functions. In such cases of ‘dumb’ classes, we usually opt to declare them as struct
rather than class
, to indicate that they are closer to C-style structures than C++-style classes.
ConstitutiveLaw
class templateThis class provides an interface for a constitutive law, i.e., for the relationship between strain \(\varepsilon(\mathbf u)\) and stress \(\sigma\). In this example we are using an elastoplastic material behavior with linear, isotropic hardening. Such materials are characterized by Young's modulus \(E\), Poisson's ratio \(\nu\), the initial yield stress \(\sigma_0\) and the isotropic hardening parameter \(\gamma\). For \(\gamma = 0\) we obtain perfect elastoplastic behavior.
As explained in the paper that describes this program, the first Newton steps are solved with a completely elastic material model to avoid having to deal with both nonlinearities (plasticity and contact) at once. To this end, this class has a function set_sigma_0()
that we use later on to simply set \(\sigma_0\) to a very large value – essentially guaranteeing that the actual stress will not exceed it, and thereby producing an elastic material. When we are ready to use a plastic model, we set \(\sigma_0\) back to its proper value, using the same function. As a result of this approach, we need to leave sigma_0
as the only non-const member variable of this class.
The constructor of the ConstitutiveLaw class sets the required material parameter for our deformable body. Material parameters for elastic isotropic media can be defined in a variety of ways, such as the pair \(E, \nu\) (elastic modulus and Poisson's number), using the Lame parameters \(\lambda,mu\) or several other commonly used conventions. Here, the constructor takes a description of material parameters in the form of \(E,\nu\), but since this turns out to these are not the coefficients that appear in the equations of the plastic projector, we immediately convert them into the more suitable set \(\kappa,\mu\) of bulk and shear moduli. In addition, the constructor takes \(\sigma_0\) (the yield stress absent any plastic strain) and \(\gamma\) (the hardening parameter) as arguments. In this constructor, we also compute the two principal components of the stress-strain relation and its linearization.
Plane stress kappa (((E*(1+2*nu)) / (std::pow((1+nu),2))) / (3 * (1 - 2 * (nu / (1+nu))))),
3d and plane strain
This is the principal component of the constitutive law. It projects the deviatoric part of the stresses in a quadrature point back to the yield stress (i.e., the original yield stress \(\sigma_0\) plus the term that describes linear isotropic hardening). We need this function to calculate the nonlinear residual in PlasticityContactProblem::residual_nl_system. The computations follow the formulas laid out in the introduction.
The function returns whether the quadrature point is plastic to allow for some statistics downstream on how many of the quadrature points are plastic and how many are elastic.
const SymmetricTensor<2, dim> deviator_stress_tensor = deviator(stress_tensor); const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
"Perforated_strip_tension" plane stress
otherwise plane strain / 3d case
This function returns the linearized stress strain tensor, linearized around the solution \(u^{i-1}\) of the previous Newton step \(i-1\). The parameter strain_tensor
(commonly denoted \(\varepsilon(u^{i-1})\)) must be passed as an argument, and serves as the linearization point. The function returns the derivative of the nonlinear constitutive law in the variable stress_strain_tensor, as well as the stress-strain tensor of the linearized problem in stress_strain_tensor_linearized. See PlasticityContactProblem::assemble_nl_system where this function is used.
Finally, below we will need a function that computes the rotation matrix induced by a displacement at a given point. In fact, of course, the displacement at a single point only has a direction and a magnitude, it is the change in direction and magnitude that induces rotations. In effect, the rotation matrix can be computed from the gradients of a displacement, or, more specifically, from the curl.
The formulas by which the rotation matrices are determined are a little awkward, especially in 3d. For 2d, there is a simpler way, so we implement this function twice, once for 2d and once for 3d, so that we can compile and use the program in both space dimensions if so desired – after all, deal.II is all about dimension independent programming and reuse of algorithm thoroughly tested with cheap computations in 2d, for the more expensive computations in 3d. Here is one case, where we have to implement different algorithms for 2d and 3d, but then can write the rest of the program in a way that is independent of the space dimension.
So, without further ado to the 2d implementation:
First, compute the curl of the velocity field from the gradients. Note that we are in 2d, so the rotation is a scalar:
From this, compute the angle of rotation:
And from this, build the antisymmetric rotation matrix:
The 3d case is a little more contrived:
Again first compute the curl of the velocity field. This time, it is a real vector:
From this vector, using its magnitude, compute the tangent of the angle of rotation, and from it the actual angle:
Now, here's one problem: if the angle of rotation is too small, that means that there is no rotation going on (for example a translational motion). In that case, the rotation matrix is the identity matrix.
The reason why we stress that is that in this case we have that tan_angle==0
. Further down, we need to divide by that number in the computation of the axis of rotation, and we would get into trouble when dividing doing so. Therefore, let's shortcut this and simply return the identity matrix if the angle of rotation is really small:
Otherwise compute the real rotation matrix. The algorithm for this is not exactly obvious, but can be found in a number of books, particularly on computer games where rotation is a very frequent operation. Online, you can find a description at http://www.makegames.com/3drotation/ and (this particular form, with the signs as here) at http://www.gamedev.net/reference/articles/article1199.asp:
The following should be relatively standard. We need classes for the boundary forcing term (which we here choose to be zero) and incremental boundary values.
BodyForce
classBody forces are generally mediated by one of the four basic physical types of forces: gravity, strong and weak interaction, and electromagnetism. Unless one wants to consider subatomic objects (for which quasistatic deformation is irrelevant and an inappropriate description anyway), only gravity and electromagnetic forces need to be considered. Let us, for simplicity assume that our body has a certain mass density, but is either non-magnetic and not electrically conducting or that there are no significant electromagnetic fields around. In that case, the body forces are simply rho g
, where rho
is the material density and g
is a vector in negative z-direction with magnitude 9.81 m/s^2. Both the density and g
are defined in the function, and we take as the density 7700 kg/m^3, a value commonly assumed for steel.
To be a little more general and to be able to do computations in 2d as well, we realize that the body force is always a function returning a dim
dimensional vector. We assume that gravity acts along the negative direction of the last, i.e. dim-1
th coordinate. The rest of the implementation of this function should be mostly self-explanatory given similar definitions in previous example programs. Note that the body force is independent of the location; to avoid compiler warnings about unused function arguments, we therefore comment out the name of the first argument of the vector_value
function:
IncrementalBoundaryValue
classIn addition to body forces, movement can be induced by boundary forces and forced boundary displacement. The latter case is equivalent to forces being chosen in such a way that they induce certain displacement.
For quasistatic displacement, typical boundary forces would be pressure on a body, or tangential friction against another body. We chose a somewhat simpler case here: we prescribe a certain movement of (parts of) the boundary, or at least of certain components of the displacement vector. We describe this by another vector-valued function that, for a given point on the boundary, returns the prescribed displacement.
Since we have a time-dependent problem, the displacement increment of the boundary equals the displacement accumulated during the length of the timestep. The class therefore has to know both the present time and the length of the present time step, and can then approximate the incremental displacement as the present velocity times the present timestep.
For the purposes of this program, we choose a simple form of boundary displacement: we displace the top boundary with constant velocity downwards. The rest of the boundary is either going to be fixed (and is then described using an object of type Functions::ZeroFunction
) or free (Neumann-type, in which case nothing special has to be done). The implementation of the class describing the constant downward motion should then be obvious using the knowledge we gained through all the previous example programs:
--------------------------— TimoshenkoBeam ------------------------------------—
compute traction on the right face of Timoshenko beam problem, t_bar
compute the fraction of imposed force
compute exact displacement of Timoshenko beam problem, u_bar
compute the fraction of imposed force
----------------------— Thick_tube_internal_pressure -------------------------------—
pressure (1.94e8),
compute traction on the inner boundary, t_bar
compute the fraction of imposed force
----------------------— Perforated_strip_tension -------------------------------—
compute the fraction of imposed force
impose displacement only on the top edge
compute the fraction of imposed displacement
----------------------— Cantiliver_beam_3d -------------------------------—
pressure should be imposed on the top surface, y = height
compute the fraction of imposed force
bound_size : size of the boundary, in 2d is the length and in the 3d case, area
bound_size : size of the boundary, in 2d is the length and in the 3d case, area
Mean stress at the specified domain is of interest. The interest domains are located on the bottom and top of the flanges close to the clamped face, z = 0 top domain: height/2 - thickness_flange <= y <= height/2 0 <= z <= 2 * thickness_flange bottom domain: -height/2 <= y <= -height/2 + thickness_flange 0 <= z <= 2 * thickness_flange
domain_size : size of the interested domain, in 2d is the area and in the 3d case, volume
top domain: height/2 - thickness_flange <= y <= height/2 0 <= z <= 2 * thickness_flange bottom domain: -height/2 <= y <= -height/2 + thickness_flange 0 <= z <= 2 * thickness_flange
Assemble right hand side of the dual problem when the quantity of interest is a nonlinear functional. In this case, the QoI should be linearized which depends on the solution of the primal problem. The extractor of the linearized QoI functional is the gradient of the the original QoI functional with the primal solution values.
bound_size : size of the boundary, in 2d is the length and in the 3d case, area
DualSolver class
constraints_hanging_nodes_dual.condense (sparsity_pattern_dual);
the boundary x = 0
the boundary y = 0
the boundary x = 0
the boundary y = 0
the boundary y = imposed incremental displacement
the boundary x = y = z = 0
+++ direct solver +++++++++
After the decomposition, we can use A_direct like a matrix representing the inverse of our system matrix, so to compute the solution we just have to multiply with the right hand side vector:
++++ iterative solver ++ CG ++++ doesn't work SolverControl solver_control (5000, 1e-12); SolverCG<> cg (solver_control);
PreconditionSSOR<> preconditioner; preconditioner.initialize(system_matrix_dual, 1.2);
cg.solve (system_matrix_dual, solution_dual, system_rhs_dual, preconditioner);
++++ iterative solver ++ BiCGStab ++++++ doesn't work SolverControl solver_control (5000, 1e-12); SolverBicgstab<> bicgstab (solver_control);
PreconditionJacobi<> preconditioner; preconditioner.initialize(system_matrix_dual, 1.0);
bicgstab.solve (system_matrix_dual, solution_dual, system_rhs_dual, preconditioner);
+++++++++++++++++++++++++++++++++++++++++++++++++
solve the dual problem
compuate the dual weights
estimate the error
--------------— estimate_some ----------------------—
------------— integrate_over_cell ----------------—
compute face_integrals
----------— integrate_over_regular_face --------—
----------— integrate_over_irregular_face ------—
PlasticityContactProblem
class templateThis is the main class of this program and supplies all functions and variables needed to describe the nonlinear contact problem. It is close to step-41 but with some additional features like handling hanging nodes, a Newton method, using Trilinos and p4est for parallel distributed computing. To deal with hanging nodes makes life a bit more complicated since we need another AffineConstraints object now. We create a Newton method for the active set method for the contact situation and to handle the nonlinear operator for the constitutive law.
The general layout of this class is very much like for most other tutorial programs. To make our life a bit easier, this class reads a set of input parameters from an input file. These parameters, using the ParameterHandler class, are declared in the declare_parameters
function (which is static so that it can be called before we even create an object of the current type), and a ParameterHandler object that has been used to read an input file will then be passed to the constructor of this class.
The remaining member functions are by and large as we have seen in several of the other tutorial programs, though with additions for the current nonlinear system. We will comment on their purpose as we get to them further below.
Next are three functions that handle the history variables stored in each quadrature point. The first one is called before the first timestep to set up a pristine state for the history variables. It only works on those quadrature points on cells that belong to the present processor:
The second one updates the history variables at the end of each timestep:
As far as member variables are concerned, we start with ones that we use to indicate the MPI universe this program runs on, and then two numbers telling us how many participating processors there are, and where in this world we are., a stream we use to let exactly one processor produce output to the console (see step-17) and a variable that is used to time the various sections of the program:
The next group describes the mesh and the finite element space. In particular, for this parallel program, the finite element space has associated with it variables that indicate which degrees of freedom live on the current processor (the index sets, see also step-40 and the Parallel computing with multiple processors using distributed memory documentation module) as well as a variety of constraints: those imposed by hanging nodes, by Dirichlet boundary conditions, and by the active set of contact nodes. Of the three AffineConstraints objects defined here, the first only contains hanging node constraints, the second also those associated with Dirichlet boundary conditions, and the third these plus the contact constraints.
The variable active_set
consists of those degrees of freedom constrained by the contact, and we use fraction_of_plastic_q_points_per_cell
to keep track of the fraction of quadrature points on each cell where the stress equals the yield stress. The latter is only used to create graphical output showing the plastic zone, but not for any further computation; the variable is a member variable of this class since the information is computed as a by-product of computing the residual, but is used only much later. (Note that the vector is a vector of length equal to the number of active cells on the local mesh; it is never used to exchange information between processors and can therefore be a regular deal.II vector.)
One difference of this program is that we declare the quadrature formula in the class declaration. The reason is that in all the other programs, it didn't do much harm if we had used different quadrature formulas when computing the matrix and the right hand side, for example. However, in the present case it does: we store information in the quadrature points, so we have to make sure all parts of the program agree on where they are and how many there are on each cell. Thus, let us first declare the quadrature formula that will be used throughout...
... and then also have a vector of history objects, one per quadrature point on those cells for which we are responsible (i.e. we don't store history data for quadrature points on cells that are owned by other processors).
The way this object is accessed is through a user pointer
that each cell, face, or edge holds: it is a void*
pointer that can be used by application programs to associate arbitrary data to cells, faces, or edges. What the program actually does with this data is within its own responsibility, the library just allocates some space for these pointers, and application programs can set and read the pointers for each of these objects.
The next block of variables corresponds to the solution and the linear systems we need to form. In particular, this includes the Newton matrix and right hand side; the vector that corresponds to the residual (i.e., the Newton right hand side) but from which we have not eliminated the various constraints and that is used to determine which degrees of freedom need to be constrained in the next iteration; and a vector that corresponds to the diagonal of the \(B\) matrix briefly mentioned in the introduction and discussed in the accompanying paper.
The next block of variables is then related to the time dependent nature of the problem: they denote the length of the time interval which we want to simulate, the present time and number of time step, and length of present timestep:
The next block contains the variables that describe the material response:
And then there is an assortment of other variables that are used to identify the mesh we are asked to build as selected by the parameter file, the obstacle that is being pushed into the deformable body, the mesh refinement strategy, whether to transfer the solution from one mesh to the next, and how many mesh refinement cycles to perform. As possible, we mark these kinds of variables as const
to help the reader identify which ones may or may not be modified later on (the output directory being an exception – it is never modified outside the constructor but it is awkward to initialize in the member-initializer-list following the colon in the constructor since there we have only one shot at setting it; the same is true for the mesh refinement criterion):
PlasticityContactProblem
classLet us start with the declaration of run-time parameters that can be selected in the input file. These values will be read back in the constructor of this class to initialize the member variables of this class:
PlasticityContactProblem
constructorGiven the declarations of member variables as well as the declarations of run-time parameters that are read from the input file, there is nothing surprising in this constructor. In the body we initialize the mesh refinement strategy and the output directory, creating such a directory if necessary.
The next block deals with constructing the starting mesh. We will use the following helper function and the first block of the make_grid()
to construct a mesh that corresponds to a half sphere. deal.II has a function that creates such a mesh, but it is in the wrong location and facing the wrong direction, so we need to shift and rotate it a bit before using it.
For later reference, as described in the documentation of GridGenerator::half_hyper_ball(), the flat surface of the halfsphere has boundary indicator zero, while the remainder has boundary indicator one.
give the indicators to boundaries for specification,
________100______ | | 0 | | 5 |________________| 100 0 to essential boundary conditions (left edge) which are as default 100 to the null boundaries (upper and lower edges) where we do not need to take care of them 5 to the natural boundaries (right edge) for imposing the traction force
Example 1 from the paper: Zhong Z., .... A new numerical method for determining collapse load-carrying capacity of structure made of elasto-plastic material, J. Cent. South Univ. (2014) 21: 398-404
give the indicators to boundaries for specification,
0 - inner boundary - natural boundary condition - impose the traction force 1 - outer boundary - free boundary - we do not need to take care of them 2 - left boundary - essential boundary condition - constrained to move along the x direction 3 - bottom boundary - essential boundary condition - constrained to move along the y direction
Example 2 from the paper: Zhong Z., .... A new numerical method for determining collapse load-carrying capacity of structure made of elasto-plastic material, J. Cent. South Univ. (2014) 21: 398-404
thickness = 0.01;
Modify the triangulation_1
distance of the face center from the center
Make the triangulation_2, a rectangular above the triangulation_1
make the triangulation_2d and refine it
Merge the two triangulation_1 and triangulation_2
Assign boundary indicators to the boundary faces
Extrude the triangulation_2d and make it 3d GridGenerator::extrude_triangulation(triangulation_2d, 2, thickness, triangulation);
Assign boundary indicators to the boundary faces
A rectangular tube made of Aluminium http://www.google.de/imgres?imgurl=http%3A%2F%2Fwww.americanaluminum.com%2Fimages%2Fstockshape-rectangletube.gif&imgrefurl=http%3A%2F%2Fwww.americanaluminum.com%2Fstandard%2FrectangleTube&h=280&w=300&tbnid=VPDNh4-DJz4wyM%3A&zoom=1&docid=9DoGJCkOeFqiSM&ei=L1AuVfG5GMvtO7DggdAF&tbm=isch&client=ubuntu&iact=rc&uact=3&dur=419&page=1&start=0&ndsp=33&ved=0CGYQrQMwFQ approximation of beam 17250 units are in meter
Make the triangulation_b, a rectangular at the bottom of rectangular tube
Make the triangulation_t, a rectangular at the top of rectangular tube
Make the triangulation_l, a rectangular at the left of rectangular tube
Make the triangulation_r, a rectangular at the right of rectangular tube
make the triangulation_2d
merging every two triangles to make triangulation_2d
Extrude the triangulation_2d and make it 3d
Assign boundary indicators to the boundary faces
The next piece in the puzzle is to set up the DoFHandler, resize vectors and take care of various other status variables such as index sets and constraint matrices.
In the following, each group of operations is put into a brace-enclosed block that is being timed by the variable declared at the top of the block (the constructor of the TimerOutput::Scope variable starts the timed section, the destructor that is called at the end of the block stops it again).
Finally, we set up sparsity patterns and matrices. We temporarily (ab)use the system matrix to also build the (diagonal) matrix that we use in eliminating degrees of freedom that are in contact with the obstacle, but we then immediately set the Newton matrix back to zero.
This function, broken out of the preceding one, computes the constraints associated with Dirichlet-type boundary conditions and puts them into the constraints_dirichlet_and_hanging_nodes
variable by merging with the constraints that come from hanging nodes.
As laid out in the introduction, we need to distinguish between two cases:
make_grid()
function, the former corresponds to boundary indicator 6, the latter to 8.the boundary x = 0
the boundary y = 0
the boundary x = 0
the boundary y = 0
the boundary y = imposed incremental displacement
the boundary x = y = z = 0
Given the complexity of the problem, it may come as a bit of a surprise that assembling the linear system we have to solve in each Newton iteration is actually fairly straightforward. The following function builds the Newton right hand side and Newton matrix. It looks fairly innocent because the heavy lifting happens in the call to ConstitutiveLaw::get_linearized_stress_strain_tensors()
and in particular in AffineConstraints<double>::distribute_local_to_global(), using the constraints we have previously computed.
std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
For assembling the local right hand side contributions, we need to access the prior linearized stress value in this quadrature point. To get it, we use the user pointer of this cell that points into the global array to the quadrature point data corresponding to the first quadrature point of the present cell, and then add an offset corresponding to the index of the quadrature point we presently consider:
In addition, we need the values of the external body forces at the quadrature points on this cell:
Having computed the stress-strain tensor and its linearization, we can now put together the parts of the matrix and right hand side. In both, we need the linearized stress-strain tensor times the symmetric gradient of \(\varphi_i\), i.e. the term \(I_\Pi\varepsilon(\varphi_i)\), so we introduce an abbreviation of this term. Recall that the matrix corresponds to the bilinear form \(A_{ij}=(I_\Pi\varepsilon(\varphi_i),\varepsilon(\varphi_j))\) in the notation of the accompanying publication, whereas the right hand side is \(F_i=([I_\Pi-P_\Pi C]\varepsilon(\varphi_i),\varepsilon(\mathbf u))\) where \(u\) is the current linearization points (typically the last solution). This might suggest that the right hand side will be zero if the material is completely elastic (where \(I_\Pi=P_\Pi\)) but this ignores the fact that the right hand side will also contain contributions from non-homogeneous constraints due to the contact.
The code block that follows this adds contributions that are due to boundary forces, should there be any.
The following function computes the nonlinear residual of the equation given the current solution (or any other linearization point). This is needed in the linear search algorithm where we need to try various linear combinations of previous and current (trial) solution to compute the (real, globalized) solution of the current Newton step.
That said, in a slight abuse of the name of the function, it actually does significantly more. For example, it also computes the vector that corresponds to the Newton residual but without eliminating constrained degrees of freedom. We need this vector to compute contact forces and, ultimately, to compute the next active set. Likewise, by keeping track of how many quadrature points we encounter on each cell that show plastic yielding, we also compute the fraction_of_plastic_q_points_per_cell
vector that we can later output to visualize the plastic zone. In both of these cases, the results are not necessary as part of the line search, and so we may be wasting a small amount of time computing them. At the same time, this information appears as a natural by-product of what we need to do here anyway, and we want to collect it once at the end of each Newton step, so we may as well do it here.
The actual implementation of this function should be rather obvious:
The last piece before we can discuss the actual Newton iteration on a single mesh is the solver for the linear systems. There are a couple of complications that slightly obscure the code, but mostly it is just setup then solve. Among the complications are:
The rest of the function is similar to step-40 and step-41 except that we use a BiCGStab solver instead of CG. This is due to the fact that for very small hardening parameters \(\gamma\), the linear system becomes almost semidefinite though still symmetric. BiCGStab appears to have an easier time with such linear systems.
----— Solver Bicgstab — Preconditioner AMG ----------------— TrilinosWrappers::PreconditionAMG preconditioner; { TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner");
std::vector<std::vector<bool> > constant_modes; DoFTools::extract_constant_modes(dof_handler, ComponentMask(), constant_modes);
TrilinosWrappers::PreconditionAMG::AdditionalData additional_data; 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;
preconditioner.initialize(newton_matrix, additional_data); }
{ TimerOutput::Scope t(computing_timer, "Solve: iterate");
TrilinosWrappers::MPI::Vector tmp(locally_owned_dofs, mpi_communicator);
// const double relative_accuracy = 1e-8; const double relative_accuracy = 1e-2; const double solver_tolerance = relative_accuracy
SolverControl solver_control(newton_matrix.m(), solver_tolerance); SolverBicgstab<TrilinosWrappers::MPI::Vector> solver(solver_control); solver.solve(newton_matrix, distributed_solution, newton_rhs, preconditioner);
pcout << " Error: " << solver_control.initial_value() << " -> " << solver_control.last_value() << " in " << solver_control.last_step() << " Bicgstab iterations." << std::endl; }
----— Solver CG — Preconditioner SSOR ----------------—
const double relative_accuracy = 1e-8;
SolverControl solver_control(newton_matrix.m(), solver_tolerance);
........................................................
This is, finally, the function that implements the damped Newton method on the current mesh. There are two nested loops: the outer loop for the Newton iteration and the inner loop for the line search which will be used only if necessary. To obtain a good and reasonable starting value we solve an elastic problem in very first Newton step on each mesh (or only on the first mesh if we transfer solutions between meshes). We do so by setting the yield stress to an unreasonably large value in these iterations and then setting it back to the correct value in subsequent iterations.
Other than this, the top part of this function should be reasonably obvious:
It gets a bit more hairy after we have computed the trial solution \(\tilde{\mathbf u}\) of the current Newton step. We handle a highly nonlinear problem so we have to damp Newton's method using a line search. To understand how we do this, recall that in our formulation, we compute a trial solution in each Newton step and not the update between old and new solution. Since the solution set is a convex set, we will use a line search that tries linear combinations of the previous and the trial solution to guarantee that the damped solution is in our solution set again. At most we apply 5 damping steps.
There are exceptions to when we use a line search. First, if this is the first Newton step on any mesh, then we don't have any point to compare the residual to, so we always accept a full step. Likewise, if this is the second Newton step on the first mesh (or the second on any mesh if we don't transfer solutions from mesh to mesh), then we have computed the first of these steps using just an elastic model (see how we set the yield stress sigma to an unreasonably large value above). In this case, the first Newton solution was a purely elastic one, the second one a plastic one, and any linear combination would not necessarily be expected to lie in the feasible set – so we just accept the solution we just got.
In either of these two cases, we bypass the line search and just update residual and other vectors as necessary.
The final step is to check for convergence. If the residual is less than a threshold of \(10^{-10}\), then we terminate the iteration on the current mesh: if (residual_norm < 1e-10)
make a non-parallel copy of tmp_solution
the dual function definition (it should be defined previously, e.g. input file)
......................................... Mean stress_yy over the bottom boundary
.........................................
displacement at Point A (x=0, y=height/2, z=length)
Mean stress at the specified domain is of interest. The interest domains are located on the bottom and top of the flanges close to the clamped face, z = 0 top domain: height/2 - thickness_flange <= y <= height/2 0 <= z <= 2 * thickness_flange bottom domain: -height/2 <= y <= -height/2 + thickness_flange 0 <= z <= 2 * thickness_flange
--------------— estimate_some ----------------------—
------------— integrate_over_cell ----------------—
compute face_integrals
----------— integrate_over_regular_face --------—
----------— integrate_over_irregular_face ------—
If you've made it this far into the deal.II tutorial, the following function refining the mesh should not pose any challenges to you any more. It refines the mesh, either globally or using the Kelly error estimator, and if so asked also transfers the solution from the previous to the next mesh. In the latter case, we also need to compute the active set and other quantities again, for which we need the information computed by compute_nonlinear_residual()
.
Make a field variable for history variables to be able to transfer the data to the quadrature points of the new mesh
Refine the mesh
const double refine_fraction_cells = .1, coarsen_fraction_cells = .3;
distributed_solution = solution;
distributed_incremental_displacement = incremental_displacement;
compute_nonlinear_residual(incremental_displacement);
stress
strain
stress
strain
Transfer the history data to the quadrature points of the new mesh In a final step, we have to get the data back from the now interpolated global field to the quadrature points on the new mesh. The following code will do that:
At the beginning of our computations, we needed to set up initial values of the history variables, such as the existing stresses in the material, that we store in each quadrature point. As mentioned above, we use the user_pointer
for this that is available in each cell.
To put this into larger perspective, we note that if we had previously available stresses in our model (which we assume do not exist for the purpose of this program), then we would need to interpolate the field of preexisting stresses to the quadrature points. Likewise, if we were to simulate elasto-plastic materials with hardening/softening, then we would have to store additional history variables like the present yield stress of the accumulated plastic strains in each quadrature points. Pre-existing hardening or weakening would then be implemented by interpolating these variables in the present function as well.
What we need to do here is to first count how many quadrature points are within the responsibility of this processor. This, of course, equals the number of cells that belong to this processor times the number of quadrature points our quadrature formula has on each cell.
For good measure, we also set all user pointers of all cells, whether ours of not, to the null pointer. This way, if we ever access the user pointer of a cell which we should not have accessed, a segmentation fault will let us know that this should not have happened:
Next, allocate as many quadrature objects as we need. Since the resize
function does not actually shrink the amount of allocated memory if the requested new size is smaller than the old size, we resort to a trick to first free all memory, and then reallocate it: we declare an empty vector as a temporary variable and then swap the contents of the old vector and this temporary variable. This makes sure that the quadrature_point_history
is now really empty, and we can let the temporary variable that now holds the previous contents of the vector go out of scope and be destroyed. In the next step. we can then re-allocate as many elements as we need, with the vector default-initializing the PointHistory
objects, which includes setting the stress variables to zero.
Finally loop over all cells again and set the user pointers from the cells that belong to the present processor to point to the first quadrature point objects corresponding to this cell in the vector of such objects:
At the end, for good measure make sure that our count of elements was correct and that we have both used up all objects we allocated previously, and not point to any objects beyond the end of the vector. Such defensive programming strategies are always good checks to avoid accidental errors and to guard against future changes to this function that forget to update all uses of a variable at the same time. Recall that constructs using the Assert
macro are optimized away in optimized mode, so do not affect the run time of optimized runs:
At the end of each time step, we should have computed an incremental displacement update so that the material in its new configuration accommodates for the difference between the external body and boundary forces applied during this time step minus the forces exerted through preexisting internal stresses. In order to have the preexisting stresses available at the next time step, we therefore have to update the preexisting stresses with the stresses due to the incremental displacement computed during the present time step. Ideally, the resulting sum of internal stresses would exactly counter all external forces. Indeed, a simple experiment can make sure that this is so: if we choose boundary conditions and body forces to be time independent, then the forcing terms (the sum of external forces and internal stresses) should be exactly zero. If you make this experiment, you will realize from the output of the norm of the right hand side in each time step that this is almost the case: it is not exactly zero, since in the first time step the incremental displacement and stress updates were computed relative to the undeformed mesh, which was then deformed. In the second time step, we again compute displacement and stress updates, but this time in the deformed mesh – there, the resulting updates are very small but not quite zero. This can be iterated, and in each such iteration the residual, i.e. the norm of the right hand side vector, is reduced; if one makes this little experiment, one realizes that the norm of this residual decays exponentially with the number of iterations, and after an initial very rapid decline is reduced by roughly a factor of about 3.5 in each iteration (for one testcase I looked at, other testcases, and other numbers of unknowns change the factor, but not the exponential decay).
In a sense, this can then be considered as a quasi-timestepping scheme to resolve the nonlinear problem of solving large-deformation elasticity on a mesh that is moved along in a Lagrangian manner.
Another complication is that the existing (old) stresses are defined on the old mesh, which we will move around after updating the stresses. If this mesh update involves rotations of the cell, then we need to also rotate the updated stress, since it was computed relative to the coordinate system of the old cell.
Thus, what we need is the following: on each cell which the present processor owns, we need to extract the old stress from the data stored with each quadrature point, compute the stress update, add the two together, and then rotate the result together with the incremental rotation computed from the incremental displacement at the present quadrature point. We will detail these steps below:
First, set up an FEValues
object by which we will evaluate the displacements and the gradients thereof at the quadrature points, together with a vector that will hold this information:
Then loop over all cells and do the job in the cells that belong to our subdomain:
Next, get a pointer to the quadrature point history data local to the present cell, and, as a defensive measure, make sure that this pointer is within the bounds of the global array:
Then initialize the FEValues
object on the present cell, and extract the strains of the displacement at the quadrature points
Then loop over the quadrature points of this cell:
The result of these operations is then written back into the original place:
The remaining three functions before we get to run()
have to do with generating output. The following one is an attempt at showing the deformed body in its deformed configuration. To this end, this function takes a displacement vector field and moves every vertex of the (local part) of the mesh by the previously computed displacement. We will call this function with the current displacement field before we generate graphical output, and we will call it again after generating graphical output with the negative displacement field to undo the changes to the mesh so made.
The function itself is pretty straightforward. All we have to do is keep track which vertices we have already touched, as we encounter the same vertices multiple times as we loop over cells.
Next is the function we use to actually generate graphical output. The function is a bit tedious, but not actually particularly complicated. It moves the mesh at the top (and moves it back at the end), then computes the contact forces along the contact surface. We can do so (as shown in the accompanying paper) by taking the untreated residual vector and identifying which degrees of freedom correspond to those with contact by asking whether they have an inhomogeneous constraints associated with them. As always, we need to be mindful that we can only write into completely distributed vectors (i.e., vectors without ghost elements) but that when we want to generate output, we need vectors that do indeed have ghost entries for all locally relevant degrees of freedom.
In the remainder of the function, we generate one VTU file on every processor, indexed by the subdomain id of this processor. On the first processor, we then also create a .pvtu
file that indexes all of the VTU files so that the entire set of output files can be read at once. These .pvtu
are used by Paraview to describe an entire parallel computation's output files. We then do the same again for the competitor of Paraview, the Visit visualization program, by creating a matching .visit
file.
produce eps files for mesh illustration
Extrapolate the stresses from Gauss point to the nodes
Then loop over the quadrature points of this cell:
Save stresses on nodes by nodal averaging construct a DoFHandler object based on FE_Q with 1 degree of freedom in order to compute stresses on nodes (by applying nodal averaging) Therefore, each vertex has one degree of freedom
begin check Point<dim> point1, point2; point1 = cell_1->vertex(v); point2 = dg_cell->vertex(v); AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError()); end check
Save figures of stresses
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Timoshenko beam
const double pressure (1.94e8), inner_radius (.1);
Plane stress const double mu (((e_modulus*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu))))); 3d and plane strain
make a non-parallel copy of solution
Compute stresses in the POLAR coordinates, 1- save it on Gauss points, 2- extrapolate them to nodes and taking their avarages (nodal avaraging)
we define a rotation matrix to be able to transform the stress from the Cartesian coordinate to the polar coordinate
Then loop over the quadrature points of this cell:
transform the stress from the Cartesian coordinate to the polar coordinate
rotation matrix
stress in polar coordinate
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ construct a DoFHandler object based on FE_Q with 1 degree of freedom in order to compute stresses on nodes (by applying nodal averaging) Therefore, each vertex has one degree of freedom
begin check Point<dim> point1, point2; point1 = cell_1->vertex(v); point2 = dg_cell->vertex(v); AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError()); end check
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
table_results_2: presenting the stress_rr and stress_tt on the nodes of bottom edge
table_results_3: presenting the mean stress_rr of the nodes on the inner radius
Plane stress const double mu (((e_modulus*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu))))); 3d and plane strain
table_results: Demonstrates the result of displacement at the top left corner versus imposed tension
make a non-parallel copy of solution
compute average sigma_yy on the bottom edge
table_results_2: Demonstrate the stress_yy on the nodes of bottom edge
if ( std::abs( (stress_yy_av/sigma_0) - .91 ) < .2 )
if ( true )
table_results_3: Demonstrate the Stress_mean (average tensile stress) on the bottom edge versus epsilon_yy on the bottom left corner
compute strain_yy_A Since the point A is the node on the bottom left corner, we need to work just with one element
Then loop over the quadrature points of this cell:
table_results: Demonstrates the result of displacement at the top front point, Point A
Quantity of interest: displacement at Point A (x=0, y=height/2, z=length)
make a non-parallel copy of solution
demonstrate the location and maximum von-Mises stress in the specified domain close to the clamped face, z = 0 top domain: height/2 - thickness_flange <= y <= height/2 0 <= z <= 2 * thickness_flange bottom domain: -height/2 <= y <= -height/2 + thickness_flange 0 <= z <= 2 * thickness_flange
Then loop over the quadrature points of this cell:
As in all other tutorial programs, the run()
function contains the overall logic. There is not very much to it here: in essence, it performs the loops over all mesh refinement cycles, and within each, hands things over to the Newton solver in solve_newton()
on the current mesh and calls the function that creates graphical output for the so-computed solution. It then outputs some statistics concerning both run times and memory consumption that has been collected over the course of computations on this mesh.
base_mesh == "Thick_tube_internal_pressure"
base_mesh == "Perforated_strip_tension"
---------------------— Refinement based on the relative error ----------------------------—
---------------------— Refinement based on the number of refinement -----------------------—
main
functionThere really isn't much to the main()
function. It looks like they always do: