599 *
const unsigned int component = 0)
const override
611 *
class BoundaryValues :
public Function<dim>
615 *
const unsigned int component = 0) const override
628 * We describe the obstacle function by a cascaded barrier (think: stair
633 * class Obstacle : public
Function<dim>
637 * const unsigned
int component = 0) const override
640 *
Assert(component == 0, ExcIndexRange(component, 0, 1));
644 *
else if (p[0] >= -0.5 && p[0] < 0.0)
646 *
else if (p[0] >= 0.0 && p[0] < 0.5)
658 * <a name=
"step_41-ImplementationofthecodeObstacleProblemcodeclass"></a>
659 * <h3>Implementation of the <code>ObstacleProblem</code>
class</h3>
667 * <a name=
"step_41-ObstacleProblemObstacleProblem"></a>
668 * <h4>ObstacleProblem::ObstacleProblem</h4>
672 * To everyone who has taken a look at the
first few tutorial programs, the
673 * constructor is completely obvious:
677 * ObstacleProblem<dim>::ObstacleProblem()
686 * <a name=
"step_41-ObstacleProblemmake_grid"></a>
687 * <h4>ObstacleProblem::make_grid</h4>
691 * We solve our obstacle problem on the square @f$[-1,1]\times [-1,1]@f$ in
692 * 2
d. This function therefore just sets up one of the simplest possible
697 *
void ObstacleProblem<dim>::make_grid()
702 * std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
712 * <a name=
"step_41-ObstacleProblemsetup_system"></a>
713 * <h4>ObstacleProblem::setup_system</h4>
717 * In
this first function of note, we set up the degrees of freedom handler,
718 * resize vectors and matrices, and deal with the constraints. Initially,
719 * the constraints are, of course, only given by boundary
values, so we
720 *
interpolate them towards the top of the function.
724 *
void ObstacleProblem<dim>::setup_system()
726 * dof_handler.distribute_dofs(fe);
727 * active_set.set_size(dof_handler.n_dofs());
729 * std::cout <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
735 * BoundaryValues<dim>(),
737 * constraints.close();
742 * system_matrix.reinit(dsp);
743 * complete_system_matrix.reinit(dsp);
745 *
IndexSet solution_index_set = dof_handler.locally_owned_dofs();
746 * solution.reinit(solution_index_set, MPI_COMM_WORLD);
747 * system_rhs.reinit(solution_index_set, MPI_COMM_WORLD);
748 * complete_system_rhs.reinit(solution_index_set, MPI_COMM_WORLD);
749 * contact_force.reinit(solution_index_set, MPI_COMM_WORLD);
753 * The only other thing to
do here is to compute the factors in the @f$B@f$
754 *
matrix which is used to
scale the residual. As discussed in the
755 * introduction, we
'll use a little trick to make this mass matrix
756 * diagonal, and in the following then first compute all of this as a
757 * matrix and then extract the diagonal elements for later use:
760 * TrilinosWrappers::SparseMatrix mass_matrix;
761 * mass_matrix.reinit(dsp);
762 * assemble_mass_matrix_diagonal(mass_matrix);
763 * diagonal_of_mass_matrix.reinit(solution_index_set);
764 * for (unsigned int j = 0; j < solution.size(); ++j)
765 * diagonal_of_mass_matrix(j) = mass_matrix.diag_element(j);
772 * <a name="step_41-ObstacleProblemassemble_system"></a>
773 * <h4>ObstacleProblem::assemble_system</h4>
777 * This function at once assembles the system matrix and right-hand-side and
778 * applied the constraints (both due to the active set as well as from
779 * boundary values) to our system. Otherwise, it is functionally equivalent
780 * to the corresponding function in, for example, @ref step_4 "step-4".
784 * void ObstacleProblem<dim>::assemble_system()
786 * std::cout << " Assembling system..." << std::endl;
791 * const QGauss<dim> quadrature_formula(fe.degree + 1);
792 * RightHandSide<dim> right_hand_side;
794 * FEValues<dim> fe_values(fe,
795 * quadrature_formula,
796 * update_values | update_gradients |
797 * update_quadrature_points | update_JxW_values);
799 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
800 * const unsigned int n_q_points = quadrature_formula.size();
802 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
803 * Vector<double> cell_rhs(dofs_per_cell);
805 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
807 * for (const auto &cell : dof_handler.active_cell_iterators())
809 * fe_values.reinit(cell);
813 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
814 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
816 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
817 * cell_matrix(i, j) +=
818 * (fe_values.shape_grad(i, q_point) *
819 * fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point));
822 * (fe_values.shape_value(i, q_point) *
823 * right_hand_side.value(fe_values.quadrature_point(q_point)) *
824 * fe_values.JxW(q_point));
827 * cell->get_dof_indices(local_dof_indices);
829 * constraints.distribute_local_to_global(cell_matrix,
843 * <a name="step_41-ObstacleProblemassemble_mass_matrix_diagonal"></a>
844 * <h4>ObstacleProblem::assemble_mass_matrix_diagonal</h4>
848 * The next function is used in the computation of the diagonal mass matrix
849 * @f$B@f$ used to scale variables in the active set method. As discussed in the
850 * introduction, we get the mass matrix to be diagonal by choosing the
851 * trapezoidal rule for quadrature. Doing so we don't really need the triple
852 *
loop over quadrature points, indices @f$i@f$ and indices @f$j@f$ any more and
853 * can, instead, just use a
double loop. The rest of the function is obvious
854 * given what we have discussed in many of the previous tutorial programs.
858 * Note that at the time
this function is called, the constraints
object
859 * only contains boundary
value constraints; we therefore
do not have to pay
861 *
matrix entries that may later on be constrained by the active set.
865 * Note also that the trick with the trapezoidal rule only works
if we have
866 * in fact @f$Q_1@f$ elements. For higher order elements, one would need to use
867 * a quadrature formula that has quadrature points at all the support points
868 * of the finite element. Constructing such a quadrature formula isn
't
869 * really difficult, but not the point here, and so we simply assert at the
870 * top of the function that our implicit assumption about the finite element
871 * is in fact satisfied.
875 * void ObstacleProblem<dim>::assemble_mass_matrix_diagonal(
876 * TrilinosWrappers::SparseMatrix &mass_matrix)
878 * Assert(fe.degree == 1, ExcNotImplemented());
880 * const QTrapezoid<dim> quadrature_formula;
881 * FEValues<dim> fe_values(fe,
882 * quadrature_formula,
883 * update_values | update_JxW_values);
885 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
886 * const unsigned int n_q_points = quadrature_formula.size();
888 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
889 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
891 * for (const auto &cell : dof_handler.active_cell_iterators())
893 * fe_values.reinit(cell);
896 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
897 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
898 * cell_matrix(i, i) +=
899 * (fe_values.shape_value(i, q_point) *
900 * fe_values.shape_value(i, q_point) * fe_values.JxW(q_point));
902 * cell->get_dof_indices(local_dof_indices);
904 * constraints.distribute_local_to_global(cell_matrix,
914 * <a name="step_41-ObstacleProblemupdate_solution_and_constraints"></a>
915 * <h4>ObstacleProblem::update_solution_and_constraints</h4>
919 * In a sense, this is the central function of this program. It updates the
920 * active set of constrained degrees of freedom as discussed in the
921 * introduction and computes an AffineConstraints object from it that can then
922 * be used to eliminate constrained degrees of freedom from the solution of
923 * the next iteration. At the same time we set the constrained degrees of
924 * freedom of the solution to the correct value, namely the height of the
929 * Fundamentally, the function is rather simple: We have to loop over all
930 * degrees of freedom and check the sign of the function @f$\Lambda^k_i +
931 * c([BU^k]_i - G_i) = \Lambda^k_i + cB_i(U^k_i - [g_h]_i)@f$ because in our
932 * case @f$G_i = B_i[g_h]_i@f$. To this end, we use the formula given in the
933 * introduction by which we can compute the Lagrange multiplier as the
934 * residual of the original linear system (given via the variables
935 * <code>complete_system_matrix</code> and <code>complete_system_rhs</code>.
936 * At the top of this function, we compute this residual using a function
937 * that is part of the matrix classes.
941 * void ObstacleProblem<dim>::update_solution_and_constraints()
943 * std::cout << " Updating active set..." << std::endl;
945 * const double penalty_parameter = 100.0;
947 * TrilinosWrappers::MPI::Vector lambda(
948 * complete_index_set(dof_handler.n_dofs()));
949 * complete_system_matrix.residual(lambda, solution, complete_system_rhs);
953 * compute contact_force[i] = - lambda[i] * diagonal_of_mass_matrix[i]
956 * contact_force = lambda;
957 * contact_force.scale(diagonal_of_mass_matrix);
958 * contact_force *= -1;
962 * The next step is to reset the active set and constraints objects and to
963 * start the loop over all degrees of freedom. This is made slightly more
964 * complicated by the fact that we can't just
loop over all elements of
965 * the solution vector since there is no way
for us then to find out what
966 * location a DoF is associated with; however, we need
this location to
967 * test whether the displacement of a DoF is larger or smaller than the
968 * height of the obstacle at
this location.
972 * We work around
this by looping over all cells and DoFs defined on each
973 * of these cells. We use here that the displacement is described
using a
974 * @f$Q_1@f$ function
for which degrees of freedom are
always located on the
975 *
vertices of the cell; thus, we can get the
index of each degree of
976 * freedom and its location by asking the vertex
for this information. On
977 * the other hand,
this clearly wouldn
't work for higher order elements,
978 * and so we add an assertion that makes sure that we only deal with
979 * elements for which all degrees of freedom are located in vertices to
980 * avoid tripping ourselves with non-functional code in case someone wants
981 * to play with increasing the polynomial degree of the solution.
985 * The price to pay for having to loop over cells rather than DoFs is that
986 * we may encounter some degrees of freedom more than once, namely each
987 * time we visit one of the cells adjacent to a given vertex. We will
988 * therefore have to keep track which vertices we have already touched and
989 * which we haven't so far. We
do so by
using an array of flags
990 * <code>dof_touched</code>:
993 * constraints.clear();
994 * active_set.clear();
996 *
const Obstacle<dim> obstacle;
997 * std::vector<bool> dof_touched(dof_handler.n_dofs(),
false);
999 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1002 *
Assert(dof_handler.get_fe().n_dofs_per_cell() == cell->n_vertices(),
1003 * ExcNotImplemented());
1005 *
const unsigned int dof_index = cell->vertex_dof_index(v, 0);
1007 *
if (dof_touched[dof_index] ==
false)
1008 * dof_touched[dof_index] =
true;
1014 * Now that we know that we haven
't touched this DoF yet, let's get
1015 * the
value of the displacement function there as well as the
value
1016 * of the obstacle function and use
this to decide whether the
1017 * current DoF belongs to the active set. For that we use the
1018 * function given above and in the introduction.
1022 * If we decide that the DoF should be part of the active set, we
1023 * add its
index to the active set, introduce an inhomogeneous
1025 * solution
value to the height of the obstacle. Finally, the
1026 * residual of the non-contact part of the system serves as an
1027 * additional control (the residual equals the remaining,
1028 * unaccounted forces, and should be zero outside the contact zone),
1029 * so we zero out the components of the residual vector (i.e., the
1030 * Lagrange multiplier lambda) that correspond to the area where the
1031 * body is in contact; at the
end of the
loop over all cells, the
1032 * residual will therefore only consist of the residual in the
1033 * non-contact zone. We output the
norm of
this residual along with
1034 * the size of the active set after the
loop.
1037 *
const double obstacle_value = obstacle.value(cell->vertex(v));
1038 *
const double solution_value = solution(dof_index);
1040 *
if (
lambda(dof_index) + penalty_parameter *
1041 * diagonal_of_mass_matrix(dof_index) *
1042 * (solution_value - obstacle_value) <
1045 * active_set.add_index(dof_index);
1046 * constraints.add_constraint(dof_index, {}, obstacle_value);
1048 * solution(dof_index) = obstacle_value;
1053 * std::cout <<
" Size of active set: " << active_set.n_elements()
1056 * std::cout <<
" Residual of the non-contact part of the system: "
1057 * <<
lambda.l2_norm() << std::endl;
1061 * In a
final step, we add to the set of constraints on DoFs we have so
1062 * far from the active set those that result from Dirichlet boundary
1063 *
values, and close the constraints object:
1068 * BoundaryValues<dim>(),
1070 * constraints.close();
1076 * <a name=
"step_41-ObstacleProblemsolve"></a>
1077 * <h4>ObstacleProblem::solve</h4>
1081 * There is
nothing to say really about the solve function. In the context
1082 * of a Newton method, we are not typically interested in very high accuracy
1083 * (why ask
for a highly accurate solution of a linear problem that we know
1084 * only gives us an approximation of the solution of the nonlinear problem),
1086 * either an absolute tolerance is reached (
for which we choose @f$10^{-12}@f$)
1087 * or when the residual is reduced by a certain factor (here, @f$10^{-3}@f$).
1090 * template <int dim>
1091 *
void ObstacleProblem<dim>::solve()
1093 * std::cout <<
" Solving system..." << std::endl;
1100 * solver.solve(system_matrix, solution, system_rhs, precondition);
1101 * constraints.distribute(solution);
1103 * std::cout <<
" Error: " << reduction_control.initial_value() <<
" -> "
1104 * << reduction_control.last_value() <<
" in "
1105 * << reduction_control.last_step() <<
" CG iterations."
1113 * <a name=
"step_41-ObstacleProblemoutput_results"></a>
1114 * <h4>ObstacleProblem::output_results</h4>
1118 * We use the
vtk-format
for the output. The file contains the displacement
1119 * and a numerical representation of the active set.
1122 *
template <
int dim>
1123 *
void ObstacleProblem<dim>::output_results(
const unsigned int iteration)
const
1125 * std::cout <<
" Writing graphical output..." << std::endl;
1128 * dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
1129 *
for (
const auto index : active_set)
1130 * active_set_vector[
index] = 1.;
1135 * data_out.add_data_vector(solution,
"displacement");
1136 * data_out.add_data_vector(active_set_vector,
"active_set");
1137 * data_out.add_data_vector(contact_force,
"lambda");
1139 * data_out.build_patches();
1141 * std::ofstream output_vtk(
"output_" +
1143 * data_out.write_vtk(output_vtk);
1151 * <a name=
"step_41-ObstacleProblemrun"></a>
1152 * <h4>ObstacleProblem::run</h4>
1156 * This is the function which has the top-
level control over everything. It
1157 * is not very long, and in fact rather straightforward: in every iteration
1158 * of the active set method, we
assemble the linear system, solve it, update
1159 * the active set and
project the solution back to the feasible set, and
1160 * then output the results. The iteration is terminated whenever the active
1161 * set has not changed in the previous iteration.
1165 * The only trickier part is that we have to save the linear system (i.e.,
1166 * the matrix and right hand side) after assembling it in the
first
1167 * iteration. The reason is that
this is the only step where we can access
1168 * the linear system as built without any of the contact constraints
1169 * active. We need
this to compute the residual of the solution at other
1170 * iterations, but in other iterations that linear system we form has the
1171 * rows and columns that correspond to constrained degrees of freedom
1172 * eliminated, and so we can no longer access the full residual of the
1173 * original equation.
1176 *
template <
int dim>
1177 *
void ObstacleProblem<dim>::run()
1182 *
IndexSet active_set_old(active_set);
1183 *
for (
unsigned int iteration = 0; iteration <= solution.size(); ++iteration)
1185 * std::cout <<
"Newton iteration " << iteration << std::endl;
1187 * assemble_system();
1189 *
if (iteration == 0)
1191 * complete_system_matrix.copy_from(system_matrix);
1192 * complete_system_rhs = system_rhs;
1196 * update_solution_and_constraints();
1197 * output_results(iteration);
1199 *
if (active_set == active_set_old)
1202 * active_set_old = active_set;
1204 * std::cout << std::endl;
1213 * <a name=
"step_41-Thecodemaincodefunction"></a>
1214 * <h3>The <code>main</code> function</h3>
1218 * And
this is the main function. It follows the pattern of all other main
1219 *
functions. The call to initialize
MPI exists because the Trilinos library
1220 * upon which we build our linear solvers in
this program
requires it.
1223 *
int main(
int argc,
char *argv[])
1227 *
using namespace dealii;
1228 *
using namespace Step41;
1235 * This program can only be
run in
serial. Otherwise,
throw an exception.
1240 *
"This program can only be run in serial, use ./step-41"));
1242 * ObstacleProblem<2> obstacle_problem;
1243 * obstacle_problem.run();
1245 *
catch (std::exception &exc)
1247 * std::cerr << std::endl
1249 * <<
"----------------------------------------------------"
1251 * std::cerr <<
"Exception on processing: " << std::endl
1252 * << exc.what() << std::endl
1253 * <<
"Aborting!" << std::endl
1254 * <<
"----------------------------------------------------"
1261 * std::cerr << std::endl
1263 * <<
"----------------------------------------------------"
1265 * std::cerr <<
"Unknown exception!" << std::endl
1266 * <<
"Aborting!" << std::endl
1267 * <<
"----------------------------------------------------"
1275<a name=
"step_41-Results"></a><h1>Results</h1>
1278Running the program produces output like
this:
1280Number of active cells: 16384
1281Total number of cells: 21845
1282Number of degrees of freedom: 16641
1285 Assembling system...
1287 Error: 0.310059 -> 5.16619e-05 in 5 CG iterations.
1288 Updating active set...
1289 Size of active set: 13164
1290 Residual of the non-contact part of the system: 1.61863e-05
1291 Writing graphical output...
1294 Assembling system...
1296 Error: 1.11987 -> 0.00109377 in 6 CG iterations.
1297 Updating active set...
1298 Size of active set: 12363
1299 Residual of the non-contact part of the system: 3.9373
1300 Writing graphical output...
1305 Assembling system...
1307 Error: 0.00713308 -> 2.29249e-06 in 4 CG iterations.
1308 Updating active set...
1309 Size of active set: 5399
1310 Residual of the non-contact part of the system: 0.000957525
1311 Writing graphical output...
1314 Assembling system...
1316 Error: 0.000957525 -> 2.8033e-07 in 4 CG iterations.
1317 Updating active set...
1318 Size of active set: 5399
1319 Residual of the non-contact part of the system: 2.8033e-07
1320 Writing graphical output...
1323The iterations
end once the active set doesn
't change any more (it has
13245,399 constrained degrees of freedom at that point). The algebraic
1325precondition is apparently working nicely since we only need 4-6 CG
1326iterations to solve the linear system (although this also has a lot to
1327do with the fact that we are not asking for very high accuracy of the
1330More revealing is to look at a sequence of graphical output files
1331(every third step is shown, with the number of the iteration in the
1334<table align="center">
1340 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.00.png" alt="">
1343 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.00.png" alt="">
1346 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.00.png" alt="">
1354 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.03.png" alt="">
1357 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.03.png" alt="">
1360 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.03.png" alt="">
1368 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.06.png" alt="">
1371 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.06.png" alt="">
1374 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.06.png" alt="">
1382 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.09.png" alt="">
1385 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.09.png" alt="">
1388 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.09.png" alt="">
1396 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.12.png" alt="">
1399 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.12.png" alt="">
1402 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.12.png" alt="">
1410 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.15.png" alt="">
1413 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.15.png" alt="">
1416 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.15.png" alt="">
1424 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.18.png" alt="">
1427 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.18.png" alt="">
1430 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.18.png" alt="">
1435The pictures show that in the first step, the solution (which has been
1436computed without any of the constraints active) bends through so much
1437that pretty much every interior point has to be bounced back to the
1438stairstep function, producing a discontinuous solution. Over the
1439course of the active set iterations, this unphysical membrane shape is
1440smoothed out, the contact with the lower-most stair step disappears,
1441and the solution stabilizes.
1443In addition to this, the program also outputs the values of the
1444Lagrange multipliers. Remember that these are the contact forces and
1445so should only be positive on the contact set, and zero outside. If,
1446on the other hand, a Lagrange multiplier is negative in the active
1447set, then this degree of freedom must be removed from the active
1448set. The following pictures show the multipliers in iterations 1, 9
1449and 18, where we use red and browns to indicate positive values, and
1450blue for negative values.
1452<table align="center">
1455 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.01.png" alt="">
1458 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.09.png" alt="">
1461 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.18.png" alt="">
1477It is easy to see that the positive values converge nicely to moderate
1478values in the interior of the contact set and large upward forces at
1479the edges of the steps, as one would expect (to support the large
1480curvature of the membrane there); at the fringes of the active set,
1481multipliers are initially negative, causing the set to shrink until,
1482in iteration 18, there are no more negative multipliers and the
1483algorithm has converged.
1487<a name="step-41-extensions"></a>
1488<a name="step_41-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1491As with any of the programs of this tutorial, there are a number of
1492obvious possibilities for extensions and experiments. The first one is
1493clear: introduce adaptivity. Contact problems are prime candidates for
1494adaptive meshes because the solution has lines along which it is less
1495regular (the places where contact is established between membrane and
1496obstacle) and other areas where the solution is very smooth (or, in
1497the present context, constant wherever it is in contact with the
1498obstacle). Adding this to the current program should not pose too many
1499difficulties, but it is not trivial to find a good error estimator for
1502A more challenging task would be an extension to 3d. The problem here
1503is not so much to simply make everything run in 3d. Rather, it is that
1504when a 3d body is deformed and gets into contact with an obstacle,
1505then the obstacle does not act as a constraining body force within the
1506domain as is the case here. Rather, the contact force only acts on the
1507boundary of the object. The inequality then is not in the differential
1508equation but in fact in the (Neumann-type) boundary conditions, though
1509this leads to a similar kind of variational
1510inequality. Mathematically, this means that the Lagrange multiplier
1511only lives on the surface, though it can of course be extended by zero
1512into the domain if that is convenient. As in the current program, one
1513does not need to form and store this Lagrange multiplier explicitly.
1515A further interesting problem for the 3d case is to consider contact problems
1516with friction. In almost every mechanical process friction has a big influence.
1517For the modelling we have to take into account tangential stresses at the contact
1518surface. Also we have to observe that friction adds another nonlinearity to
1521Another nontrivial modification is to implement a more complex constitutive
1522law like nonlinear elasticity or elasto-plastic material behavior.
1523The difficulty here is to handle the additional nonlinearity arising
1524through the nonlinear constitutive law.
1527<a name="step_41-PlainProg"></a>
1528<h1> The plain program</h1>
1529@include "step-41.cc"
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
#define AssertThrow(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
std::vector< value_type > preserve(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
@ matrix
Contents is actually a matrix.
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
std::vector< unsigned int > serial(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm comm)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
static const unsigned int invalid_unsigned_int
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation