}
}
}
+
+ // Finally, we need to copy the lower half of the local matrix into the
+ // upper half:
+ for (const unsigned int i : scratch.fe_values.dof_indices())
+ for (const unsigned int j :
+ scratch.fe_values.dof_indices_starting_at(i + 1))
+ data.cell_matrix(i, j) = data.cell_matrix(j, i);
}
+
+
// @sect4{Solid::make_constraints}
// The constraints for this problem are simple to describe.
- // However, since we are dealing with an iterative Newton method,
- // it should be noted that any displacement constraints should only
- // be specified at the zeroth iteration and subsequently no
- // additional contributions are to be made since the constraints
- // are already exactly satisfied.
+ // In this particular example, the boundary values will be calculated for
+ // the two first iterations of Newton's algorithm. In general, one would
+ // build non-homogeneous constraints in the zeroth iteration (that is, when
+ // `apply_dirichlet_bc == true` in the code block that follows) and build
+ // only the corresponding homogeneous constraints in the following step. While
+ // the current example has only homogeneous constraints, previous experiences
+ // have shown that a common error is forgetting to add the extra condition
+ // when refactoring the code to specific uses. This could lead to errors that
+ // are hard to debug. In this spirit, we choose to make the code more verbose
+ // in terms of what operations are performed at each Newton step.
template <int dim>
- void Solid<dim>::make_constraints(const int &it_nr)
+ void Solid<dim>::make_constraints(const int it_nr)
{
- std::cout << " CST " << std::flush;
-
- // Since the constraints are different at different Newton iterations, we
- // need to clear the constraints matrix and completely rebuild
- // it. However, after the first iteration, the constraints remain the same
- // and we can simply skip the rebuilding step if we do not clear it.
- if (it_nr > 1)
- return;
- constraints.clear();
+ // Since we (a) are dealing with an iterative Newton method, (b) are using
+ // an incremental formulation for the displacement, and (c) apply the
+ // constraints to the incremental displacement field, any non-homogeneous
+ // constraints on the displacement update should only be specified at the
+ // zeroth iteration. No subsequent contributions are to be made since the
+ // constraints will be exactly satisfied after that iteration.
const bool apply_dirichlet_bc = (it_nr == 0);
- // In this particular example, the boundary values will be calculated for
- // the two first iterations of Newton's algorithm. In general, one would
- // build non-homogeneous constraints in the zeroth iteration (that is, when
- // `apply_dirichlet_bc == true`) and build only the corresponding
- // homogeneous constraints in the following step. While the current
- // example has only homogeneous constraints, previous experiences have
- // shown that a common error is forgetting to add the extra condition when
- // refactoring the code to specific uses. This could lead errors that are
- // hard to debug. In this spirit, we choose to make the code more verbose
- // in terms of what operations are performed at each Newton step.
- //
- // The boundary conditions for the indentation problem are as follows: On
- // the -x, -y and -z faces (IDs 0,2,4) we set up a symmetry condition to
- // allow only planar movement while the +x and +z faces (IDs 1,5) are
- // traction free. In this contrived problem, part of the +y face (ID 3) is
- // set to have no motion in the x- and z-component. Finally, as described
- // earlier, the other part of the +y face has an the applied pressure but
- // is also constrained in the x- and z-directions.
- //
- // In the following, we will have to tell the function interpolation
- // boundary values which components of the solution vector should be
- // constrained (i.e., whether it's the x-, y-, z-displacements or
- // combinations thereof). This is done using ComponentMask objects (see
- // @ref GlossComponentMask) which we can get from the finite element if we
- // provide it with an extractor object for the component we wish to
- // select. To this end we first set up such extractor objects and later
- // use it when generating the relevant component masks:
- const FEValuesExtractors::Scalar x_displacement(0);
- const FEValuesExtractors::Scalar y_displacement(1);
+ // Furthermore, after the first Newton iteration within a timestep, the
+ // constraints remain the same and we do not need to modify or rebuild them
+ // so long as we do not clear the @p constraints object.
+ if (it_nr > 1)
+ {
+ std::cout << " --- " << std::flush;
+ return;
+ }
- {
- const int boundary_id = 0;
-
- if (apply_dirichlet_bc == true)
- VectorTools::interpolate_boundary_values(
- dof_handler,
- boundary_id,
- Functions::ZeroFunction<dim>(n_components),
- constraints,
- fe.component_mask(x_displacement));
- else
- VectorTools::interpolate_boundary_values(
- dof_handler,
- boundary_id,
- Functions::ZeroFunction<dim>(n_components),
- constraints,
- fe.component_mask(x_displacement));
- }
- {
- const int boundary_id = 2;
-
- if (apply_dirichlet_bc == true)
- VectorTools::interpolate_boundary_values(
- dof_handler,
- boundary_id,
- Functions::ZeroFunction<dim>(n_components),
- constraints,
- fe.component_mask(y_displacement));
- else
- VectorTools::interpolate_boundary_values(
- dof_handler,
- boundary_id,
- Functions::ZeroFunction<dim>(n_components),
- constraints,
- fe.component_mask(y_displacement));
- }
+ std::cout << " CST " << std::flush;
- if (dim == 3)
+ if (apply_dirichlet_bc)
{
- const FEValuesExtractors::Scalar z_displacement(2);
+ // At the zeroth Newton iteration we wish to apply the full set of
+ // non-homogeneous and homogeneous constraints that represent the
+ // boundary conditions on the displacement increment. Since in general
+ // the constraints may be different at each time step, we need to clear
+ // the constraints matrix and completely rebuild it. An example case
+ // would be if a surface is accelerating; in such a scenario the change
+ // in displacement is non-constant between each time step.
+ constraints.clear();
+
+ // The boundary conditions for the indentation problem in 3D are as
+ // follows: On the -x, -y and -z faces (IDs 0,2,4) we set up a symmetry
+ // condition to allow only planar movement while the +x and +z faces
+ // (IDs 1,5) are traction free. In this contrived problem, part of the
+ // +y face (ID 3) is set to have no motion in the x- and z-component.
+ // Finally, as described earlier, the other part of the +y face has an
+ // the applied pressure but is also constrained in the x- and
+ // z-directions.
+ //
+ // In the following, we will have to tell the function interpolation
+ // boundary values which components of the solution vector should be
+ // constrained (i.e., whether it's the x-, y-, z-displacements or
+ // combinations thereof). This is done using ComponentMask objects (see
+ // @ref GlossComponentMask) which we can get from the finite element if we
+ // provide it with an extractor object for the component we wish to
+ // select. To this end we first set up such extractor objects and later
+ // use it when generating the relevant component masks:
+ const FEValuesExtractors::Scalar x_displacement(0);
+ const FEValuesExtractors::Scalar y_displacement(1);
{
- const int boundary_id = 3;
-
- if (apply_dirichlet_bc == true)
- VectorTools::interpolate_boundary_values(
- dof_handler,
- boundary_id,
- Functions::ZeroFunction<dim>(n_components),
- constraints,
- (fe.component_mask(x_displacement) |
- fe.component_mask(z_displacement)));
- else
- VectorTools::interpolate_boundary_values(
- dof_handler,
- boundary_id,
- Functions::ZeroFunction<dim>(n_components),
- constraints,
- (fe.component_mask(x_displacement) |
- fe.component_mask(z_displacement)));
+ const int boundary_id = 0;
+
+ VectorTools::interpolate_boundary_values(
+ dof_handler,
+ boundary_id,
+ Functions::ZeroFunction<dim>(n_components),
+ constraints,
+ fe.component_mask(x_displacement));
}
{
- const int boundary_id = 4;
-
- if (apply_dirichlet_bc == true)
- VectorTools::interpolate_boundary_values(
- dof_handler,
- boundary_id,
- Functions::ZeroFunction<dim>(n_components),
- constraints,
- fe.component_mask(z_displacement));
- else
- VectorTools::interpolate_boundary_values(
- dof_handler,
- boundary_id,
- Functions::ZeroFunction<dim>(n_components),
- constraints,
- fe.component_mask(z_displacement));
+ const int boundary_id = 2;
+
+ VectorTools::interpolate_boundary_values(
+ dof_handler,
+ boundary_id,
+ Functions::ZeroFunction<dim>(n_components),
+ constraints,
+ fe.component_mask(y_displacement));
}
- {
- const int boundary_id = 6;
-
- if (apply_dirichlet_bc == true)
- VectorTools::interpolate_boundary_values(
- dof_handler,
- boundary_id,
- Functions::ZeroFunction<dim>(n_components),
- constraints,
- (fe.component_mask(x_displacement) |
- fe.component_mask(z_displacement)));
- else
- VectorTools::interpolate_boundary_values(
- dof_handler,
- boundary_id,
- Functions::ZeroFunction<dim>(n_components),
- constraints,
- (fe.component_mask(x_displacement) |
- fe.component_mask(z_displacement)));
- }
+ if (dim == 3)
+ {
+ const FEValuesExtractors::Scalar z_displacement(2);
+
+ {
+ const int boundary_id = 3;
+
+ VectorTools::interpolate_boundary_values(
+ dof_handler,
+ boundary_id,
+ Functions::ZeroFunction<dim>(n_components),
+ constraints,
+ (fe.component_mask(x_displacement) |
+ fe.component_mask(z_displacement)));
+ }
+ {
+ const int boundary_id = 4;
+
+ VectorTools::interpolate_boundary_values(
+ dof_handler,
+ boundary_id,
+ Functions::ZeroFunction<dim>(n_components),
+ constraints,
+ fe.component_mask(z_displacement));
+ }
+
+ {
+ const int boundary_id = 6;
+
+ VectorTools::interpolate_boundary_values(
+ dof_handler,
+ boundary_id,
+ Functions::ZeroFunction<dim>(n_components),
+ constraints,
+ (fe.component_mask(x_displacement) |
+ fe.component_mask(z_displacement)));
+ }
+ }
+ else
+ {
+ {
+ const int boundary_id = 3;
+
+ VectorTools::interpolate_boundary_values(
+ dof_handler,
+ boundary_id,
+ Functions::ZeroFunction<dim>(n_components),
+ constraints,
+ (fe.component_mask(x_displacement)));
+ }
+ {
+ const int boundary_id = 6;
+
+ VectorTools::interpolate_boundary_values(
+ dof_handler,
+ boundary_id,
+ Functions::ZeroFunction<dim>(n_components),
+ constraints,
+ (fe.component_mask(x_displacement)));
+ }
+ }
}
else
{
projects / dealii.git / commitdiff