603 *
const unsigned int component = 0)
const override
615 *
class BoundaryValues :
public Function<dim>
619 *
const unsigned int component = 0) const override
632 * We describe the obstacle
function by a cascaded barrier (think: stair
637 * class Obstacle : public
Function<dim>
641 *
const unsigned int component = 0) const override
648 *
else if (p(0) >= -0.5 && p(0) < 0.0)
650 *
else if (p(0) >= 0.0 && p(0) < 0.5)
662 * <a name=
"ImplementationofthecodeObstacleProblemcodeclass"></a>
663 * <h3>Implementation of the <code>ObstacleProblem</code>
class</h3>
671 * <a name=
"ObstacleProblemObstacleProblem"></a>
672 * <h4>ObstacleProblem::ObstacleProblem</h4>
676 * To everyone who has taken a look at the
first few tutorial programs, the
677 * constructor is completely obvious:
681 * ObstacleProblem<dim>::ObstacleProblem()
690 * <a name=
"ObstacleProblemmake_grid"></a>
691 * <h4>ObstacleProblem::make_grid</h4>
695 * We solve our obstacle problem on the square @f$[-1,1]\times [-1,1]@f$ in
696 * 2D. This
function therefore just sets up
one of the simplest possible
701 *
void ObstacleProblem<dim>::make_grid()
706 * std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
716 * <a name=
"ObstacleProblemsetup_system"></a>
717 * <h4>ObstacleProblem::setup_system</h4>
721 * In
this first function of note, we
set up the degrees of freedom handler,
722 * resize vectors and matrices, and deal with the constraints. Initially,
723 * the constraints are, of course, only given by boundary values, so we
724 *
interpolate them towards the top of the
function.
728 *
void ObstacleProblem<dim>::setup_system()
730 * dof_handler.distribute_dofs(fe);
731 * active_set.set_size(dof_handler.n_dofs());
733 * std::cout <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
739 * BoundaryValues<dim>(),
741 * constraints.close();
746 * system_matrix.reinit(dsp);
747 * complete_system_matrix.reinit(dsp);
749 *
IndexSet solution_index_set = dof_handler.locally_owned_dofs();
750 * solution.reinit(solution_index_set, MPI_COMM_WORLD);
751 * system_rhs.reinit(solution_index_set, MPI_COMM_WORLD);
752 * complete_system_rhs.reinit(solution_index_set, MPI_COMM_WORLD);
753 * contact_force.reinit(solution_index_set, MPI_COMM_WORLD);
757 * The only other thing to
do here is to compute the factors in the @f$B@f$
758 *
matrix which is used to
scale the residual. As discussed in the
759 * introduction, we
'll use a little trick to make this mass matrix
760 * diagonal, and in the following then first compute all of this as a
761 * matrix and then extract the diagonal elements for later use:
764 * TrilinosWrappers::SparseMatrix mass_matrix;
765 * mass_matrix.reinit(dsp);
766 * assemble_mass_matrix_diagonal(mass_matrix);
767 * diagonal_of_mass_matrix.reinit(solution_index_set);
768 * for (unsigned int j = 0; j < solution.size(); j++)
769 * diagonal_of_mass_matrix(j) = mass_matrix.diag_element(j);
776 * <a name="ObstacleProblemassemble_system"></a>
777 * <h4>ObstacleProblem::assemble_system</h4>
781 * This function at once assembles the system matrix and right-hand-side and
782 * applied the constraints (both due to the active set as well as from
783 * boundary values) to our system. Otherwise, it is functionally equivalent
784 * to the corresponding function in, for example, @ref step_4 "step-4".
788 * void ObstacleProblem<dim>::assemble_system()
790 * std::cout << " Assembling system..." << std::endl;
795 * const QGauss<dim> quadrature_formula(fe.degree + 1);
796 * RightHandSide<dim> right_hand_side;
798 * FEValues<dim> fe_values(fe,
799 * quadrature_formula,
800 * update_values | update_gradients |
801 * update_quadrature_points | update_JxW_values);
803 * const unsigned int dofs_per_cell = fe.dofs_per_cell;
804 * const unsigned int n_q_points = quadrature_formula.size();
806 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
807 * Vector<double> cell_rhs(dofs_per_cell);
809 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
811 * for (const auto &cell : dof_handler.active_cell_iterators())
813 * fe_values.reinit(cell);
817 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
818 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
820 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
821 * cell_matrix(i, j) +=
822 * (fe_values.shape_grad(i, q_point) *
823 * fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point));
826 * (fe_values.shape_value(i, q_point) *
827 * right_hand_side.value(fe_values.quadrature_point(q_point)) *
828 * fe_values.JxW(q_point));
831 * cell->get_dof_indices(local_dof_indices);
833 * constraints.distribute_local_to_global(cell_matrix,
847 * <a name="ObstacleProblemassemble_mass_matrix_diagonal"></a>
848 * <h4>ObstacleProblem::assemble_mass_matrix_diagonal</h4>
852 * The next function is used in the computation of the diagonal mass matrix
853 * @f$B@f$ used to scale variables in the active set method. As discussed in the
854 * introduction, we get the mass matrix to be diagonal by choosing the
855 * trapezoidal rule for quadrature. Doing so we don't really need the triple
856 *
loop over quadrature points, indices @f$i@f$ and indices @f$j@f$ any more and
857 * can, instead, just use a
double loop. The rest of the
function is obvious
858 * given what we have discussed in many of the previous tutorial programs.
862 * Note that at the time
this function is called, the constraints
object
863 * only contains boundary
value constraints; we therefore
do not have to pay
864 * attention in the last
copy-local-to-global step to
preserve the values of
865 *
matrix entries that may later on be constrained by the active
set.
869 * Note also that the trick with the trapezoidal rule only works
if we have
870 * in fact @f$Q_1@f$ elements. For higher order elements,
one would need to use
871 * a quadrature formula that has quadrature points at all the support points
872 * of the finite element. Constructing such a quadrature formula isn
't
873 * really difficult, but not the point here, and so we simply assert at the
874 * top of the function that our implicit assumption about the finite element
875 * is in fact satisfied.
879 * void ObstacleProblem<dim>::assemble_mass_matrix_diagonal(
880 * TrilinosWrappers::SparseMatrix &mass_matrix)
882 * Assert(fe.degree == 1, ExcNotImplemented());
884 * const QTrapez<dim> quadrature_formula;
885 * FEValues<dim> fe_values(fe,
886 * quadrature_formula,
887 * update_values | update_JxW_values);
889 * const unsigned int dofs_per_cell = fe.dofs_per_cell;
890 * const unsigned int n_q_points = quadrature_formula.size();
892 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
893 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
895 * for (const auto &cell : dof_handler.active_cell_iterators())
897 * fe_values.reinit(cell);
900 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
901 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
902 * cell_matrix(i, i) +=
903 * (fe_values.shape_value(i, q_point) *
904 * fe_values.shape_value(i, q_point) * fe_values.JxW(q_point));
906 * cell->get_dof_indices(local_dof_indices);
908 * constraints.distribute_local_to_global(cell_matrix,
918 * <a name="ObstacleProblemupdate_solution_and_constraints"></a>
919 * <h4>ObstacleProblem::update_solution_and_constraints</h4>
923 * In a sense, this is the central function of this program. It updates the
924 * active set of constrained degrees of freedom as discussed in the
925 * introduction and computes an AffineConstraints object from it that can then
926 * be used to eliminate constrained degrees of freedom from the solution of
927 * the next iteration. At the same time we set the constrained degrees of
928 * freedom of the solution to the correct value, namely the height of the
933 * Fundamentally, the function is rather simple: We have to loop over all
934 * degrees of freedom and check the sign of the function @f$\Lambda^k_i +
935 * c([BU^k]_i - G_i) = \Lambda^k_i + cB_i(U^k_i - [g_h]_i)@f$ because in our
936 * case @f$G_i = B_i[g_h]_i@f$. To this end, we use the formula given in the
937 * introduction by which we can compute the Lagrange multiplier as the
938 * residual of the original linear system (given via the variables
939 * <code>complete_system_matrix</code> and <code>complete_system_rhs</code>.
940 * At the top of this function, we compute this residual using a function
941 * that is part of the matrix classes.
945 * void ObstacleProblem<dim>::update_solution_and_constraints()
947 * std::cout << " Updating active set..." << std::endl;
949 * const double penalty_parameter = 100.0;
951 * TrilinosWrappers::MPI::Vector lambda(
952 * complete_index_set(dof_handler.n_dofs()));
953 * complete_system_matrix.residual(lambda, solution, complete_system_rhs);
957 * compute contact_force[i] = - lambda[i] * diagonal_of_mass_matrix[i]
960 * contact_force = lambda;
961 * contact_force.scale(diagonal_of_mass_matrix);
962 * contact_force *= -1;
966 * The next step is to reset the active set and constraints objects and to
967 * start the loop over all degrees of freedom. This is made slightly more
968 * complicated by the fact that we can't just
loop over all elements of
969 * the solution vector since there is no way
for us then to find out what
970 * location a DoF is associated with; however, we need
this location to
971 * test whether the displacement of a DoF is larger or smaller than the
972 * height of the obstacle at
this location.
976 * We work around
this by looping over all cells and DoFs defined on each
977 * of these cells. We use here that the displacement is described
using a
978 * @f$Q_1@f$
function for which degrees of freedom are
always located on the
979 *
vertices of the cell; thus, we can get the index of each degree of
980 * freedom and its location by asking the vertex
for this information. On
981 * the other hand,
this clearly wouldn
't work for higher order elements,
982 * and so we add an assertion that makes sure that we only deal with
983 * elements for which all degrees of freedom are located in vertices to
984 * avoid tripping ourselves with non-functional code in case someone wants
985 * to play with increasing the polynomial degree of the solution.
989 * The price to pay for having to loop over cells rather than DoFs is that
990 * we may encounter some degrees of freedom more than once, namely each
991 * time we visit one of the cells adjacent to a given vertex. We will
992 * therefore have to keep track which vertices we have already touched and
993 * which we haven't so far. We
do so by
using an array of flags
994 * <code>dof_touched</code>:
997 * constraints.
clear();
998 * active_set.clear();
1000 *
const Obstacle<dim> obstacle;
1001 * std::vector<bool> dof_touched(dof_handler.n_dofs(),
false);
1003 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1004 *
for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
1006 *
Assert(dof_handler.get_fe().dofs_per_cell ==
1010 *
const unsigned int dof_index = cell->vertex_dof_index(v, 0);
1012 *
if (dof_touched[dof_index] ==
false)
1013 * dof_touched[dof_index] =
true;
1019 * Now that we know that we haven
't touched this DoF yet, let's get
1020 * the
value of the displacement
function there as well as the
value
1021 * of the obstacle
function and use
this to decide whether the
1022 * current DoF belongs to the active
set. For that we use the
1023 *
function given above and in the introduction.
1027 * If we decide that the DoF should be part of the active
set, we
1028 * add its index to the active
set, introduce an inhomogeneous
1030 * solution
value to the height of the obstacle. Finally, the
1031 * residual of the non-contact part of the system serves as an
1032 * additional control (the residual equals the remaining,
1033 * unaccounted forces, and should be
zero outside the contact zone),
1034 * so we
zero out the components of the residual vector (i.e., the
1035 * Lagrange multiplier lambda) that correspond to the area where the
1036 * body is in contact; at the
end of the
loop over all cells, the
1037 * residual will therefore only consist of the residual in the
1038 * non-contact zone. We output the
norm of
this residual along with
1039 * the size of the active
set after the
loop.
1042 *
const double obstacle_value = obstacle.value(cell->vertex(v));
1043 *
const double solution_value = solution(dof_index);
1045 *
if (
lambda(dof_index) + penalty_parameter *
1046 * diagonal_of_mass_matrix(dof_index) *
1047 * (solution_value - obstacle_value) <
1050 * active_set.add_index(dof_index);
1051 * constraints.add_line(dof_index);
1052 * constraints.set_inhomogeneity(dof_index, obstacle_value);
1054 * solution(dof_index) = obstacle_value;
1059 * std::cout <<
" Size of active set: " << active_set.n_elements()
1062 * std::cout <<
" Residual of the non-contact part of the system: "
1063 * <<
lambda.l2_norm() << std::endl;
1067 * In a
final step, we add to the
set of constraints on DoFs we have so
1068 * far from the active
set those that result from Dirichlet boundary
1069 * values, and close the constraints object:
1074 * BoundaryValues<dim>(),
1076 * constraints.close();
1082 * <a name=
"ObstacleProblemsolve"></a>
1083 * <h4>ObstacleProblem::solve</h4>
1087 * There is nothing to say really about the solve
function. In the context
1088 * of a Newton method, we are not typically interested in very high accuracy
1089 * (why ask
for a highly accurate solution of a linear problem that we know
1090 * only gives us an approximation of the solution of the nonlinear problem),
1092 * either an absolute tolerance is reached (
for which we choose @f$10^{-12}@f$)
1093 * or when the residual is reduced by a certain factor (here, @f$10^{-3}@f$).
1096 * template <int dim>
1097 *
void ObstacleProblem<dim>::solve()
1099 * std::cout <<
" Solving system..." << std::endl;
1106 * solver.solve(system_matrix, solution, system_rhs, precondition);
1107 * constraints.distribute(solution);
1109 * std::cout <<
" Error: " << reduction_control.initial_value() <<
" -> "
1110 * << reduction_control.last_value() <<
" in "
1111 * << reduction_control.last_step() <<
" CG iterations."
1119 * <a name=
"ObstacleProblemoutput_results"></a>
1120 * <h4>ObstacleProblem::output_results</h4>
1124 * We use the
vtk-format
for the output. The file contains the displacement
1125 * and a numerical representation of the active
set.
1128 *
template <
int dim>
1129 *
void ObstacleProblem<dim>::output_results(
const unsigned int iteration)
const
1131 * std::cout <<
" Writing graphical output..." << std::endl;
1134 * dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
1135 *
for (
const auto index : active_set)
1136 * active_set_vector[index] = 1.;
1147 * std::ofstream output_vtk(
"output_" +
1157 * <a name=
"ObstacleProblemrun"></a>
1162 * This is the
function which has the top-
level control over everything. It
1163 * is not very long, and in fact rather straightforward: in every iteration
1164 * of the active
set method, we
assemble the linear system, solve it, update
1165 * the active
set and
project the solution back to the feasible
set, and
1166 * then output the results. The iteration is terminated whenever the active
1167 *
set has not changed in the previous iteration.
1171 * The only trickier part is that we have to save the linear system (i.e.,
1172 * the
matrix and right hand side) after assembling it in the
first
1173 * iteration. The reason is that
this is the only step where we can access
1174 * the linear system as built without any of the contact constraints
1175 * active. We need
this to compute the residual of the solution at other
1176 * iterations, but in other iterations that linear system we form has the
1177 * rows and columns that correspond to constrained degrees of freedom
1178 * eliminated, and so we can no longer access the full residual of the
1179 * original equation.
1182 *
template <
int dim>
1188 *
IndexSet active_set_old(active_set);
1189 *
for (
unsigned int iteration = 0; iteration <= solution.size(); ++iteration)
1191 * std::cout <<
"Newton iteration " << iteration << std::endl;
1193 * assemble_system();
1195 *
if (iteration == 0)
1197 * complete_system_matrix.copy_from(system_matrix);
1198 * complete_system_rhs = system_rhs;
1202 * update_solution_and_constraints();
1203 * output_results(iteration);
1205 *
if (active_set == active_set_old)
1208 * active_set_old = active_set;
1210 * std::cout << std::endl;
1219 * <a name=
"Thecodemaincodefunction"></a>
1220 * <h3>The <code>main</code>
function</h3>
1224 * And
this is the main
function. It follows the pattern of all other main
1225 *
functions. The
call to initialize MPI exists because the Trilinos library
1226 * upon which we build our linear solvers in
this program requires it.
1229 *
int main(
int argc,
char *argv[])
1233 *
using namespace dealii;
1234 *
using namespace Step41;
1241 * This program can only be
run in serial. Otherwise,
throw an exception.
1246 *
"This program can only be run in serial, use ./step-41"));
1248 * ObstacleProblem<2> obstacle_problem;
1249 * obstacle_problem.run();
1251 *
catch (std::exception &exc)
1253 * std::cerr << std::endl
1255 * <<
"----------------------------------------------------"
1257 * std::cerr <<
"Exception on processing: " << std::endl
1258 * << exc.what() << std::endl
1259 * <<
"Aborting!" << std::endl
1260 * <<
"----------------------------------------------------"
1267 * std::cerr << std::endl
1269 * <<
"----------------------------------------------------"
1271 * std::cerr <<
"Unknown exception!" << std::endl
1272 * <<
"Aborting!" << std::endl
1273 * <<
"----------------------------------------------------"
1281 <a name=
"Results"></a><h1>Results</h1>
1284 Running the program produces output like
this:
1286 Number of active cells: 16384
1287 Total number of cells: 21845
1288 Number of degrees of freedom: 16641
1291 Assembling system...
1293 Error: 0.310059 -> 5.16619e-05 in 5 CG iterations.
1294 Updating active
set...
1295 Size of active
set: 13164
1296 Residual of the non-contact part of the system: 1.61863e-05
1297 Writing graphical output...
1300 Assembling system...
1302 Error: 1.11987 -> 0.00109377 in 6 CG iterations.
1303 Updating active
set...
1304 Size of active
set: 12363
1305 Residual of the non-contact part of the system: 3.9373
1306 Writing graphical output...
1311 Assembling system...
1313 Error: 0.00713308 -> 2.29249e-06 in 4 CG iterations.
1314 Updating active
set...
1315 Size of active
set: 5399
1316 Residual of the non-contact part of the system: 0.000957525
1317 Writing graphical output...
1320 Assembling system...
1322 Error: 0.000957525 -> 2.8033e-07 in 4 CG iterations.
1323 Updating active
set...
1324 Size of active
set: 5399
1325 Residual of the non-contact part of the system: 2.8033e-07
1326 Writing graphical output...
1329 The iterations
end once the active
set doesn
't change any more (it has
1330 5,399 constrained degrees of freedom at that point). The algebraic
1331 precondition is apparently working nicely since we only need 4-6 CG
1332 iterations to solve the linear system (although this also has a lot to
1333 do with the fact that we are not asking for very high accuracy of the
1336 More revealing is to look at a sequence of graphical output files
1337 (every third step is shown, with the number of the iteration in the
1340 <table align="center">
1346 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.00.png" alt="">
1349 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.00.png" alt="">
1352 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.00.png" alt="">
1360 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.03.png" alt="">
1363 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.03.png" alt="">
1366 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.03.png" alt="">
1374 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.06.png" alt="">
1377 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.06.png" alt="">
1380 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.06.png" alt="">
1388 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.09.png" alt="">
1391 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.09.png" alt="">
1394 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.09.png" alt="">
1402 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.12.png" alt="">
1405 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.12.png" alt="">
1408 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.12.png" alt="">
1416 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.15.png" alt="">
1419 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.15.png" alt="">
1422 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.15.png" alt="">
1430 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.18.png" alt="">
1433 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.18.png" alt="">
1436 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.18.png" alt="">
1441 The pictures show that in the first step, the solution (which has been
1442 computed without any of the constraints active) bends through so much
1443 that pretty much every interior point has to be bounced back to the
1444 stairstep function, producing a discontinuous solution. Over the
1445 course of the active set iterations, this unphysical membrane shape is
1446 smoothed out, the contact with the lower-most stair step disappears,
1447 and the solution stabilizes.
1449 In addition to this, the program also outputs the values of the
1450 Lagrange multipliers. Remember that these are the contact forces and
1451 so should only be positive on the contact set, and zero outside. If,
1452 on the other hand, a Lagrange multiplier is negative in the active
1453 set, then this degree of freedom must be removed from the active
1454 set. The following pictures show the multipliers in iterations 1, 9
1455 and 18, where we use red and browns to indicate positive values, and
1456 blue for negative values.
1458 <table align="center">
1461 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.01.png" alt="">
1464 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.09.png" alt="">
1467 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.18.png" alt="">
1483 It is easy to see that the positive values converge nicely to moderate
1484 values in the interior of the contact set and large upward forces at
1485 the edges of the steps, as one would expect (to support the large
1486 curvature of the membrane there); at the fringes of the active set,
1487 multipliers are initially negative, causing the set to shrink until,
1488 in iteration 18, there are no more negative multipliers and the
1489 algorithm has converged.
1493 <a name="extensions"></a>
1494 <a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1497 As with any of the programs of this tutorial, there are a number of
1498 obvious possibilities for extensions and experiments. The first one is
1499 clear: introduce adaptivity. Contact problems are prime candidates for
1500 adaptive meshes because the solution has lines along which it is less
1501 regular (the places where contact is established between membrane and
1502 obstacle) and other areas where the solution is very smooth (or, in
1503 the present context, constant wherever it is in contact with the
1504 obstacle). Adding this to the current program should not pose too many
1505 difficulties, but it is not trivial to find a good error estimator for
1508 A more challenging task would be an extension to 3d. The problem here
1509 is not so much to simply make everything run in 3d. Rather, it is that
1510 when a 3d body is deformed and gets into contact with an obstacle,
1511 then the obstacle does not act as a constraining body force within the
1512 domain as is the case here. Rather, the contact force only acts on the
1513 boundary of the object. The inequality then is not in the differential
1514 equation but in fact in the (Neumann-type) boundary conditions, though
1515 this leads to a similar kind of variational
1516 inequality. Mathematically, this means that the Lagrange multiplier
1517 only lives on the surface, though it can of course be extended by zero
1518 into the domain if that is convenient. As in the current program, one
1519 does not need to form and store this Lagrange multiplier explicitly.
1521 A further interesting problem for the 3d case is to consider contact problems
1522 with friction. In almost every mechanical process friction has a big influence.
1523 For the modelling we have to take into account tangential stresses at the contact
1524 surface. Also we have to observe that friction adds another nonlinearity to
1527 Another nontrivial modification is to implement a more complex constitutive
1528 law like nonlinear elasticity or elasto-plastic material behavior.
1529 The difficulty here is to handle the additional nonlinearity arising
1530 through the nonlinear constitutive law.
1533 <a name="PlainProg"></a>
1534 <h1> The plain program</h1>
1535 @include "step-41.cc"