From: Martin Kronbichler Date: Wed, 28 Oct 2020 13:38:53 +0000 (+0100) Subject: Merge pull request #11099 from jppelteret/step_44_update_assembly_01 X-Git-Tag: v9.3.0-rc1~960 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e038590ef8b38d256b1f3389e35367f634652048;p=dealii.git Merge pull request #11099 from jppelteret/step_44_update_assembly_01 Step 44: Unify assembly loop and introduce an optimisation --- e038590ef8b38d256b1f3389e35367f634652048 diff --cc examples/step-44/step-44.cc index ec6d15beb8,dfeeb33205..ca8515584e --- a/examples/step-44/step-44.cc +++ b/examples/step-44/step-44.cc @@@ -2387,155 -2236,168 +2236,164 @@@ namespace Step4 } } } + + // 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 - void Solid::make_constraints(const int &it_nr) + void Solid::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(n_components), - constraints, - fe.component_mask(x_displacement)); - else - VectorTools::interpolate_boundary_values( - dof_handler, - boundary_id, - Functions::ZeroFunction(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(n_components), - constraints, - fe.component_mask(y_displacement)); - else - VectorTools::interpolate_boundary_values( - dof_handler, - boundary_id, - Functions::ZeroFunction(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(n_components), - constraints, - (fe.component_mask(x_displacement) | - fe.component_mask(z_displacement))); - else - VectorTools::interpolate_boundary_values( - dof_handler, - boundary_id, - Functions::ZeroFunction(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(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(n_components), - constraints, - fe.component_mask(z_displacement)); - else - VectorTools::interpolate_boundary_values( - dof_handler, - boundary_id, - Functions::ZeroFunction(n_components), - constraints, - fe.component_mask(z_displacement)); + const int boundary_id = 2; + + VectorTools::interpolate_boundary_values( + dof_handler, + boundary_id, + Functions::ZeroFunction(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(n_components), - constraints, - (fe.component_mask(x_displacement) | - fe.component_mask(z_displacement))); - else - VectorTools::interpolate_boundary_values( - dof_handler, - boundary_id, - Functions::ZeroFunction(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(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(n_components), + constraints, + fe.component_mask(z_displacement)); + } + + { + const int boundary_id = 6; + + VectorTools::interpolate_boundary_values( + dof_handler, + boundary_id, + Functions::ZeroFunction(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(n_components), + constraints, + (fe.component_mask(x_displacement))); + } + { + const int boundary_id = 6; + + VectorTools::interpolate_boundary_values( + dof_handler, + boundary_id, + Functions::ZeroFunction(n_components), + constraints, + (fe.component_mask(x_displacement))); + } + } } else {