600 *
const unsigned int component = 0)
const override
612 *
class BoundaryValues :
public Function<dim>
616 *
const unsigned int component = 0) const override
629 * We describe the obstacle function by a cascaded barrier (think: stair
634 * class Obstacle : public
Function<dim>
638 * const unsigned
int component = 0) const override
641 *
Assert(component == 0, ExcIndexRange(component, 0, 1));
645 *
else if (p(0) >= -0.5 && p(0) < 0.0)
647 *
else if (p(0) >= 0.0 && p(0) < 0.5)
659 * <a name=
"ImplementationofthecodeObstacleProblemcodeclass"></a>
660 * <h3>Implementation of the <code>ObstacleProblem</code>
class</h3>
668 * <a name=
"ObstacleProblemObstacleProblem"></a>
669 * <h4>ObstacleProblem::ObstacleProblem</h4>
673 * To everyone who has taken a look at the
first few tutorial programs, the
674 * constructor is completely obvious:
678 * ObstacleProblem<dim>::ObstacleProblem()
687 * <a name=
"ObstacleProblemmake_grid"></a>
688 * <h4>ObstacleProblem::make_grid</h4>
692 * We solve our obstacle problem on the square @f$[-1,1]\times [-1,1]@f$ in
693 * 2
d. This function therefore just sets up one of the simplest possible
698 *
void ObstacleProblem<dim>::make_grid()
703 * std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
713 * <a name=
"ObstacleProblemsetup_system"></a>
714 * <h4>ObstacleProblem::setup_system</h4>
718 * In
this first function of note, we
set up the degrees of freedom handler,
719 * resize vectors and matrices, and deal with the constraints. Initially,
720 * the constraints are, of course, only given by boundary
values, so we
721 *
interpolate them towards the top of the function.
725 *
void ObstacleProblem<dim>::setup_system()
727 * dof_handler.distribute_dofs(fe);
728 * active_set.set_size(dof_handler.n_dofs());
730 * std::cout <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
736 * BoundaryValues<dim>(),
738 * constraints.close();
743 * system_matrix.reinit(dsp);
744 * complete_system_matrix.reinit(dsp);
746 *
IndexSet solution_index_set = dof_handler.locally_owned_dofs();
747 * solution.reinit(solution_index_set, MPI_COMM_WORLD);
748 * system_rhs.reinit(solution_index_set, MPI_COMM_WORLD);
749 * complete_system_rhs.reinit(solution_index_set, MPI_COMM_WORLD);
750 * contact_force.reinit(solution_index_set, MPI_COMM_WORLD);
754 * The only other thing to
do here is to compute the factors in the @f$B@f$
755 *
matrix which is used to
scale the residual. As discussed in the
756 * introduction, we
'll use a little trick to make this mass matrix
757 * diagonal, and in the following then first compute all of this as a
758 * matrix and then extract the diagonal elements for later use:
761 * TrilinosWrappers::SparseMatrix mass_matrix;
762 * mass_matrix.reinit(dsp);
763 * assemble_mass_matrix_diagonal(mass_matrix);
764 * diagonal_of_mass_matrix.reinit(solution_index_set);
765 * for (unsigned int j = 0; j < solution.size(); ++j)
766 * diagonal_of_mass_matrix(j) = mass_matrix.diag_element(j);
773 * <a name="ObstacleProblemassemble_system"></a>
774 * <h4>ObstacleProblem::assemble_system</h4>
778 * This function at once assembles the system matrix and right-hand-side and
779 * applied the constraints (both due to the active set as well as from
780 * boundary values) to our system. Otherwise, it is functionally equivalent
781 * to the corresponding function in, for example, @ref step_4 "step-4".
785 * void ObstacleProblem<dim>::assemble_system()
787 * std::cout << " Assembling system..." << std::endl;
792 * const QGauss<dim> quadrature_formula(fe.degree + 1);
793 * RightHandSide<dim> right_hand_side;
795 * FEValues<dim> fe_values(fe,
796 * quadrature_formula,
797 * update_values | update_gradients |
798 * update_quadrature_points | update_JxW_values);
800 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
801 * const unsigned int n_q_points = quadrature_formula.size();
803 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
804 * Vector<double> cell_rhs(dofs_per_cell);
806 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
808 * for (const auto &cell : dof_handler.active_cell_iterators())
810 * fe_values.reinit(cell);
814 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
815 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
817 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
818 * cell_matrix(i, j) +=
819 * (fe_values.shape_grad(i, q_point) *
820 * fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point));
823 * (fe_values.shape_value(i, q_point) *
824 * right_hand_side.value(fe_values.quadrature_point(q_point)) *
825 * fe_values.JxW(q_point));
828 * cell->get_dof_indices(local_dof_indices);
830 * constraints.distribute_local_to_global(cell_matrix,
844 * <a name="ObstacleProblemassemble_mass_matrix_diagonal"></a>
845 * <h4>ObstacleProblem::assemble_mass_matrix_diagonal</h4>
849 * The next function is used in the computation of the diagonal mass matrix
850 * @f$B@f$ used to scale variables in the active set method. As discussed in the
851 * introduction, we get the mass matrix to be diagonal by choosing the
852 * trapezoidal rule for quadrature. Doing so we don't really need the triple
853 *
loop over quadrature points, indices @f$i@f$ and indices @f$j@f$ any more and
854 * can, instead, just use a
double loop. The rest of the function is obvious
855 * given what we have discussed in many of the previous tutorial programs.
859 * Note that at the time
this function is called, the constraints
object
860 * only contains boundary
value constraints; we therefore
do not have to pay
862 *
matrix entries that may later on be constrained by the active
set.
866 * Note also that the trick with the trapezoidal rule only works
if we have
867 * in fact @f$Q_1@f$ elements. For higher order elements, one would need to use
868 * a quadrature formula that has quadrature points at all the support points
869 * of the finite element. Constructing such a quadrature formula isn
't
870 * really difficult, but not the point here, and so we simply assert at the
871 * top of the function that our implicit assumption about the finite element
872 * is in fact satisfied.
876 * void ObstacleProblem<dim>::assemble_mass_matrix_diagonal(
877 * TrilinosWrappers::SparseMatrix &mass_matrix)
879 * Assert(fe.degree == 1, ExcNotImplemented());
881 * const QTrapezoid<dim> quadrature_formula;
882 * FEValues<dim> fe_values(fe,
883 * quadrature_formula,
884 * update_values | update_JxW_values);
886 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
887 * const unsigned int n_q_points = quadrature_formula.size();
889 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
890 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
892 * for (const auto &cell : dof_handler.active_cell_iterators())
894 * fe_values.reinit(cell);
897 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
898 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
899 * cell_matrix(i, i) +=
900 * (fe_values.shape_value(i, q_point) *
901 * fe_values.shape_value(i, q_point) * fe_values.JxW(q_point));
903 * cell->get_dof_indices(local_dof_indices);
905 * constraints.distribute_local_to_global(cell_matrix,
915 * <a name="ObstacleProblemupdate_solution_and_constraints"></a>
916 * <h4>ObstacleProblem::update_solution_and_constraints</h4>
920 * In a sense, this is the central function of this program. It updates the
921 * active set of constrained degrees of freedom as discussed in the
922 * introduction and computes an AffineConstraints object from it that can then
923 * be used to eliminate constrained degrees of freedom from the solution of
924 * the next iteration. At the same time we set the constrained degrees of
925 * freedom of the solution to the correct value, namely the height of the
930 * Fundamentally, the function is rather simple: We have to loop over all
931 * degrees of freedom and check the sign of the function @f$\Lambda^k_i +
932 * c([BU^k]_i - G_i) = \Lambda^k_i + cB_i(U^k_i - [g_h]_i)@f$ because in our
933 * case @f$G_i = B_i[g_h]_i@f$. To this end, we use the formula given in the
934 * introduction by which we can compute the Lagrange multiplier as the
935 * residual of the original linear system (given via the variables
936 * <code>complete_system_matrix</code> and <code>complete_system_rhs</code>.
937 * At the top of this function, we compute this residual using a function
938 * that is part of the matrix classes.
942 * void ObstacleProblem<dim>::update_solution_and_constraints()
944 * std::cout << " Updating active set..." << std::endl;
946 * const double penalty_parameter = 100.0;
948 * TrilinosWrappers::MPI::Vector lambda(
949 * complete_index_set(dof_handler.n_dofs()));
950 * complete_system_matrix.residual(lambda, solution, complete_system_rhs);
954 * compute contact_force[i] = - lambda[i] * diagonal_of_mass_matrix[i]
957 * contact_force = lambda;
958 * contact_force.scale(diagonal_of_mass_matrix);
959 * contact_force *= -1;
963 * The next step is to reset the active set and constraints objects and to
964 * start the loop over all degrees of freedom. This is made slightly more
965 * complicated by the fact that we can't just
loop over all elements of
966 * the solution vector since there is no way
for us then to find out what
967 * location a DoF is associated with; however, we need
this location to
968 * test whether the displacement of a DoF is larger or smaller than the
969 * height of the obstacle at
this location.
973 * We work around
this by looping over all cells and DoFs defined on each
974 * of these cells. We use here that the displacement is described
using a
975 * @f$Q_1@f$ function
for which degrees of freedom are
always located on the
976 *
vertices of the cell; thus, we can get the
index of each degree of
977 * freedom and its location by asking the vertex
for this information. On
978 * the other hand,
this clearly wouldn
't work for higher order elements,
979 * and so we add an assertion that makes sure that we only deal with
980 * elements for which all degrees of freedom are located in vertices to
981 * avoid tripping ourselves with non-functional code in case someone wants
982 * to play with increasing the polynomial degree of the solution.
986 * The price to pay for having to loop over cells rather than DoFs is that
987 * we may encounter some degrees of freedom more than once, namely each
988 * time we visit one of the cells adjacent to a given vertex. We will
989 * therefore have to keep track which vertices we have already touched and
990 * which we haven't so far. We
do so by
using an array of flags
991 * <code>dof_touched</code>:
994 * constraints.clear();
995 * active_set.clear();
997 *
const Obstacle<dim> obstacle;
998 * std::vector<bool> dof_touched(dof_handler.n_dofs(),
false);
1000 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1003 *
Assert(dof_handler.get_fe().n_dofs_per_cell() == cell->n_vertices(),
1004 * ExcNotImplemented());
1006 *
const unsigned int dof_index = cell->vertex_dof_index(v, 0);
1008 *
if (dof_touched[dof_index] ==
false)
1009 * dof_touched[dof_index] =
true;
1015 * Now that we know that we haven
't touched this DoF yet, let's get
1016 * the
value of the displacement function there as well as the
value
1017 * of the obstacle function and use
this to decide whether the
1018 * current DoF belongs to the active
set. For that we use the
1019 * function given above and in the introduction.
1023 * If we decide that the DoF should be part of the active
set, we
1024 * add its
index to the active
set, introduce an inhomogeneous
1026 * solution
value to the height of the obstacle. Finally, the
1027 * residual of the non-contact part of the system serves as an
1028 * additional control (the residual equals the remaining,
1029 * unaccounted forces, and should be zero outside the contact zone),
1030 * so we zero out the components of the residual vector (i.e., the
1031 * Lagrange multiplier lambda) that correspond to the area where the
1032 * body is in contact; at the
end of the
loop over all cells, the
1033 * residual will therefore only consist of the residual in the
1034 * non-contact zone. We output the
norm of
this residual along with
1035 * the size of the active
set after the
loop.
1038 *
const double obstacle_value = obstacle.value(cell->vertex(v));
1039 *
const double solution_value = solution(dof_index);
1041 *
if (
lambda(dof_index) + penalty_parameter *
1042 * diagonal_of_mass_matrix(dof_index) *
1043 * (solution_value - obstacle_value) <
1046 * active_set.add_index(dof_index);
1047 * constraints.add_line(dof_index);
1048 * constraints.set_inhomogeneity(dof_index, obstacle_value);
1050 * solution(dof_index) = obstacle_value;
1055 * std::cout <<
" Size of active set: " << active_set.n_elements()
1058 * std::cout <<
" Residual of the non-contact part of the system: "
1059 * <<
lambda.l2_norm() << std::endl;
1063 * In a
final step, we add to the
set of constraints on DoFs we have so
1064 * far from the active
set those that result from Dirichlet boundary
1065 *
values, and close the constraints object:
1070 * BoundaryValues<dim>(),
1072 * constraints.close();
1078 * <a name=
"ObstacleProblemsolve"></a>
1079 * <h4>ObstacleProblem::solve</h4>
1083 * There is
nothing to say really about the solve function. In the context
1084 * of a Newton method, we are not typically interested in very high accuracy
1085 * (why ask
for a highly accurate solution of a linear problem that we know
1086 * only gives us an approximation of the solution of the nonlinear problem),
1088 * either an absolute tolerance is reached (
for which we choose @f$10^{-12}@f$)
1089 * or when the residual is reduced by a certain factor (here, @f$10^{-3}@f$).
1092 * template <int dim>
1093 *
void ObstacleProblem<dim>::solve()
1095 * std::cout <<
" Solving system..." << std::endl;
1102 * solver.solve(system_matrix, solution, system_rhs, precondition);
1103 * constraints.distribute(solution);
1105 * std::cout <<
" Error: " << reduction_control.initial_value() <<
" -> "
1106 * << reduction_control.last_value() <<
" in "
1107 * << reduction_control.last_step() <<
" CG iterations."
1115 * <a name=
"ObstacleProblemoutput_results"></a>
1116 * <h4>ObstacleProblem::output_results</h4>
1120 * We use the
vtk-format
for the output. The file contains the displacement
1121 * and a numerical representation of the active
set.
1124 *
template <
int dim>
1125 *
void ObstacleProblem<dim>::output_results(
const unsigned int iteration)
const
1127 * std::cout <<
" Writing graphical output..." << std::endl;
1130 * dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
1131 *
for (
const auto index : active_set)
1132 * active_set_vector[
index] = 1.;
1137 * data_out.add_data_vector(solution,
"displacement");
1138 * data_out.add_data_vector(active_set_vector,
"active_set");
1139 * data_out.add_data_vector(contact_force,
"lambda");
1141 * data_out.build_patches();
1143 * std::ofstream output_vtk(
"output_" +
1145 * data_out.write_vtk(output_vtk);
1153 * <a name=
"ObstacleProblemrun"></a>
1154 * <h4>ObstacleProblem::run</h4>
1158 * This is the function which has the top-
level control over everything. It
1159 * is not very long, and in fact rather straightforward: in every iteration
1160 * of the active
set method, we
assemble the linear system, solve it, update
1161 * the active
set and
project the solution back to the feasible
set, and
1162 * then output the results. The iteration is terminated whenever the active
1163 *
set has not changed in the previous iteration.
1167 * The only trickier part is that we have to save the linear system (i.e.,
1168 * the matrix and right hand side) after assembling it in the
first
1169 * iteration. The reason is that
this is the only step where we can access
1170 * the linear system as built without any of the contact constraints
1171 * active. We need
this to compute the residual of the solution at other
1172 * iterations, but in other iterations that linear system we form has the
1173 * rows and columns that correspond to constrained degrees of freedom
1174 * eliminated, and so we can no longer access the full residual of the
1175 * original equation.
1178 *
template <
int dim>
1179 *
void ObstacleProblem<dim>::run()
1184 *
IndexSet active_set_old(active_set);
1185 *
for (
unsigned int iteration = 0; iteration <= solution.size(); ++iteration)
1187 * std::cout <<
"Newton iteration " << iteration << std::endl;
1189 * assemble_system();
1191 *
if (iteration == 0)
1193 * complete_system_matrix.copy_from(system_matrix);
1194 * complete_system_rhs = system_rhs;
1198 * update_solution_and_constraints();
1199 * output_results(iteration);
1201 *
if (active_set == active_set_old)
1204 * active_set_old = active_set;
1206 * std::cout << std::endl;
1215 * <a name=
"Thecodemaincodefunction"></a>
1216 * <h3>The <code>main</code> function</h3>
1220 * And
this is the main function. It follows the pattern of all other main
1221 *
functions. The
call to initialize MPI exists because the Trilinos library
1222 * upon which we build our linear solvers in
this program
requires it.
1225 *
int main(
int argc,
char *argv[])
1229 *
using namespace dealii;
1230 *
using namespace Step41;
1237 * This program can only be
run in
serial. Otherwise,
throw an exception.
1242 *
"This program can only be run in serial, use ./step-41"));
1244 * ObstacleProblem<2> obstacle_problem;
1245 * obstacle_problem.run();
1247 *
catch (std::exception &exc)
1249 * std::cerr << std::endl
1251 * <<
"----------------------------------------------------"
1253 * std::cerr <<
"Exception on processing: " << std::endl
1254 * << exc.what() << std::endl
1255 * <<
"Aborting!" << std::endl
1256 * <<
"----------------------------------------------------"
1263 * std::cerr << std::endl
1265 * <<
"----------------------------------------------------"
1267 * std::cerr <<
"Unknown exception!" << std::endl
1268 * <<
"Aborting!" << std::endl
1269 * <<
"----------------------------------------------------"
1277<a name=
"Results"></a><h1>Results</h1>
1280Running the program produces output like
this:
1282Number of active cells: 16384
1283Total number of cells: 21845
1284Number of degrees of freedom: 16641
1287 Assembling system...
1289 Error: 0.310059 -> 5.16619e-05 in 5 CG iterations.
1290 Updating active
set...
1291 Size of active
set: 13164
1292 Residual of the non-contact part of the system: 1.61863e-05
1293 Writing graphical output...
1296 Assembling system...
1298 Error: 1.11987 -> 0.00109377 in 6 CG iterations.
1299 Updating active
set...
1300 Size of active
set: 12363
1301 Residual of the non-contact part of the system: 3.9373
1302 Writing graphical output...
1307 Assembling system...
1309 Error: 0.00713308 -> 2.29249e-06 in 4 CG iterations.
1310 Updating active
set...
1311 Size of active
set: 5399
1312 Residual of the non-contact part of the system: 0.000957525
1313 Writing graphical output...
1316 Assembling system...
1318 Error: 0.000957525 -> 2.8033e-07 in 4 CG iterations.
1319 Updating active
set...
1320 Size of active
set: 5399
1321 Residual of the non-contact part of the system: 2.8033e-07
1322 Writing graphical output...
1325The iterations
end once the active
set doesn
't change any more (it has
13265,399 constrained degrees of freedom at that point). The algebraic
1327precondition is apparently working nicely since we only need 4-6 CG
1328iterations to solve the linear system (although this also has a lot to
1329do with the fact that we are not asking for very high accuracy of the
1332More revealing is to look at a sequence of graphical output files
1333(every third step is shown, with the number of the iteration in the
1336<table align="center">
1342 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.00.png" alt="">
1345 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.00.png" alt="">
1348 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.00.png" alt="">
1356 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.03.png" alt="">
1359 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.03.png" alt="">
1362 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.03.png" alt="">
1370 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.06.png" alt="">
1373 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.06.png" alt="">
1376 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.06.png" alt="">
1384 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.09.png" alt="">
1387 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.09.png" alt="">
1390 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.09.png" alt="">
1398 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.12.png" alt="">
1401 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.12.png" alt="">
1404 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.12.png" alt="">
1412 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.15.png" alt="">
1415 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.15.png" alt="">
1418 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.15.png" alt="">
1426 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.18.png" alt="">
1429 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.18.png" alt="">
1432 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.18.png" alt="">
1437The pictures show that in the first step, the solution (which has been
1438computed without any of the constraints active) bends through so much
1439that pretty much every interior point has to be bounced back to the
1440stairstep function, producing a discontinuous solution. Over the
1441course of the active set iterations, this unphysical membrane shape is
1442smoothed out, the contact with the lower-most stair step disappears,
1443and the solution stabilizes.
1445In addition to this, the program also outputs the values of the
1446Lagrange multipliers. Remember that these are the contact forces and
1447so should only be positive on the contact set, and zero outside. If,
1448on the other hand, a Lagrange multiplier is negative in the active
1449set, then this degree of freedom must be removed from the active
1450set. The following pictures show the multipliers in iterations 1, 9
1451and 18, where we use red and browns to indicate positive values, and
1452blue for negative values.
1454<table align="center">
1457 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.01.png" alt="">
1460 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.09.png" alt="">
1463 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.18.png" alt="">
1479It is easy to see that the positive values converge nicely to moderate
1480values in the interior of the contact set and large upward forces at
1481the edges of the steps, as one would expect (to support the large
1482curvature of the membrane there); at the fringes of the active set,
1483multipliers are initially negative, causing the set to shrink until,
1484in iteration 18, there are no more negative multipliers and the
1485algorithm has converged.
1489<a name="extensions"></a>
1490<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1493As with any of the programs of this tutorial, there are a number of
1494obvious possibilities for extensions and experiments. The first one is
1495clear: introduce adaptivity. Contact problems are prime candidates for
1496adaptive meshes because the solution has lines along which it is less
1497regular (the places where contact is established between membrane and
1498obstacle) and other areas where the solution is very smooth (or, in
1499the present context, constant wherever it is in contact with the
1500obstacle). Adding this to the current program should not pose too many
1501difficulties, but it is not trivial to find a good error estimator for
1504A more challenging task would be an extension to 3d. The problem here
1505is not so much to simply make everything run in 3d. Rather, it is that
1506when a 3d body is deformed and gets into contact with an obstacle,
1507then the obstacle does not act as a constraining body force within the
1508domain as is the case here. Rather, the contact force only acts on the
1509boundary of the object. The inequality then is not in the differential
1510equation but in fact in the (Neumann-type) boundary conditions, though
1511this leads to a similar kind of variational
1512inequality. Mathematically, this means that the Lagrange multiplier
1513only lives on the surface, though it can of course be extended by zero
1514into the domain if that is convenient. As in the current program, one
1515does not need to form and store this Lagrange multiplier explicitly.
1517A further interesting problem for the 3d case is to consider contact problems
1518with friction. In almost every mechanical process friction has a big influence.
1519For the modelling we have to take into account tangential stresses at the contact
1520surface. Also we have to observe that friction adds another nonlinearity to
1523Another nontrivial modification is to implement a more complex constitutive
1524law like nonlinear elasticity or elasto-plastic material behavior.
1525The difficulty here is to handle the additional nonlinearity arising
1526through the nonlinear constitutive law.
1529<a name="PlainProg"></a>
1530<h1> The plain program</h1>
1531@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())
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
#define AssertThrow(cond, exc)
void loop(ITERATOR begin, std_cxx20::type_identity_t< ITERATOR > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), 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)
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
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