1250 *
for (
unsigned int i = 0; i < dim; ++i)
1252 *
for (
unsigned int k = 0; k < dim; ++k)
1254 * coupling[displacement<dim> + i]
1256 * coupling[displacement_multiplier<dim> + k]
1262 * coupling[density_lower_slack<dim>][density_lower_slack<dim>] =
1264 * coupling[density_lower_slack<dim>][density_upper_slack<dim>] =
1266 * coupling[density_upper_slack<dim>][density_lower_slack<dim>] =
1269 * coupling[density_lower_slack_multiplier<dim>]
1271 * coupling[density_lower_slack_multiplier<dim>]
1273 * coupling[density_upper_slack_multiplier<dim>]
1279 * Before we can create the sparsity pattern, we also have to
1280 *
set up constraints. Since
this program does not adaptively
1281 *
refine the mesh, the only constraint we have is
one that
1282 * couples all density variables to enforce the
volume
1283 * constraint. This will ultimately lead to a dense sub-block
1284 * of the
matrix, but there is little we can
do about that.
1294 * constraints.clear();
1295 * constraints.add_line(last_density_dof);
1296 *
for (
unsigned int i = 0; i < density_dofs.
n_elements() - 1; ++i)
1297 * constraints.add_entry(last_density_dof,
1300 * constraints.set_inhomogeneity(last_density_dof, 0);
1302 * constraints.close();
1306 * We can now
finally create the sparsity pattern
for the
1307 *
matrix, taking into account which variables couple with
1308 * which other variables, and the constraints we have on the
1316 * The only part of the
matrix we have not dealt with is the
1318 * (integral) operators
for which deal.II does not currently
1319 * have
functions. What we will ultimately need to
do is go
1320 * over all cells and couple the unfiltered density on
this
1321 * cell to all filtered densities of neighboring cells that
1322 * are less than a threshold distance away, and the other way
1323 * around;
for the moment, we are only concerned with building
1324 * the sparsity pattern that would correspond to
this kind of
1325 *
matrix, so we perform the equivalent
loop and where later
1326 * on we would write into an entry of the
matrix, we now
1327 * simply add an entry to the sparsity
matrix:
1330 * for (
const auto &cell : dof_handler.active_cell_iterators())
1332 *
const unsigned int i = cell->active_cell_index();
1333 *
for (
const auto &check_cell : find_relevant_neighbors(cell))
1335 *
const double distance =
1336 * cell->center().distance(check_cell->center());
1337 *
if (distance < filter_r)
1340 * .block(SolutionBlocks::unfiltered_density,
1341 * SolutionBlocks::unfiltered_density_multiplier)
1342 * .add(i, check_cell->active_cell_index());
1344 * .block(SolutionBlocks::unfiltered_density_multiplier,
1345 * SolutionBlocks::unfiltered_density)
1346 * .add(i, check_cell->active_cell_index());
1353 * Having so generated the
"dynamic" sparsity pattern, we can
1354 *
finally copy it to the structure that is used to associate
1355 * matrices with a sparsity pattern. Because the sparsity
1356 * pattern is large and complex, we also output it into a file
1357 * of its own
for visualization purposes -- in other words,
1358 *
for "visual debugging".
1361 * sparsity_pattern.copy_from(dsp);
1363 * std::ofstream out(
"sparsity.plt");
1364 * sparsity_pattern.print_gnuplot(out);
1366 * system_matrix.reinit(sparsity_pattern);
1371 * What is left is to correctly size the various vectors and
1372 * their blocks, as well as setting
initial guesses
for some
1373 * of the components of the (nonlinear) solution vector. We
1374 * here use the symbolic component names
for individual blocks
1375 * of the solution vector and,
for brevity, use the same trick
1376 * with `
using namespace` as above:
1379 * nonlinear_solution.reinit(block_sizes);
1380 * system_rhs.reinit(block_sizes);
1383 *
using namespace SolutionBlocks;
1384 * nonlinear_solution.block(density).add(density_ratio);
1385 * nonlinear_solution.block(unfiltered_density).add(density_ratio);
1386 * nonlinear_solution.block(unfiltered_density_multiplier)
1387 * .add(density_ratio);
1388 * nonlinear_solution.block(density_lower_slack).add(density_ratio);
1389 * nonlinear_solution.block(density_lower_slack_multiplier).add(50);
1390 * nonlinear_solution.block(density_upper_slack).add(1 - density_ratio);
1391 * nonlinear_solution.block(density_upper_slack_multiplier).add(50);
1399 * <a name=
"Creatingthefiltermatrix"></a>
1400 * <h3>Creating the filter
matrix</h3>
1404 * Next up, a function that is used once at the beginning of the
1405 * program: It creates a
matrix @f$H@f$ so that the filtered density
1406 * vector equals @f$H@f$ times the unfiltered density. The creation
1407 * of
this matrix is non-trivial, and it is used in every
1408 * iteration, and so rather than reforming it as we
do with the
1409 * Newton
matrix, it is made only once and stored separately.
1413 * The way
this matrix is computed follows the outline used above
1414 * already to form its sparsity pattern. We repeat
this process here
1415 *
for the sparsity pattern of
this separately formed
matrix, and
1416 * then actually build the
matrix itself. You may want to check the
1417 * definition of
this matrix in the introduction to
this program.
1420 *
template <
int dim>
1421 *
void SANDTopOpt<dim>::setup_filter_matrix()
1425 * The sparsity pattern of the filter has already been determined
1426 * and implemented in the setup_system() function. We
copy the
1427 * structure from the appropriate block and use it again here.
1433 * filter_sparsity_pattern.copy_from(
1434 * sparsity_pattern.block(SolutionBlocks::unfiltered_density,
1435 * SolutionBlocks::unfiltered_density_multiplier));
1436 * filter_matrix.
reinit(filter_sparsity_pattern);
1440 * Having so built the sparsity pattern, now we re-do all of
1441 * these loops to actually compute the necessary
values of the
1448 * for (const auto &cell : dof_handler.active_cell_iterators())
1450 *
const unsigned int i = cell->active_cell_index();
1451 *
for (
const auto &check_cell : find_relevant_neighbors(cell))
1453 *
const double distance =
1454 * cell->center().distance(check_cell->center());
1455 *
if (distance < filter_r)
1457 * filter_matrix.add(i,
1458 * check_cell->active_cell_index(),
1459 * filter_r - distance);
1467 * The
final step is to normalize the
matrix so that
for each
1468 * row, the
sum of entries equals
one.
1471 *
for (
unsigned int i = 0; i < filter_matrix.m(); ++i)
1473 *
double denominator = 0;
1475 * iter != filter_matrix.end(i);
1477 * denominator = denominator + iter->value();
1479 * iter != filter_matrix.end(i);
1481 * iter->value() = iter->value() / denominator;
1487 * This function is used
for building the filter
matrix. We create a
set of
1488 * all the cell iterators within a certain radius of the cell that is input.
1489 * These are the neighboring cells that will be relevant
for the filter.
1492 *
template <
int dim>
1493 * std::set<typename Triangulation<dim>::cell_iterator>
1494 * SANDTopOpt<dim>::find_relevant_neighbors(
1497 * std::set<unsigned int> neighbor_ids;
1498 * std::set<typename Triangulation<dim>::cell_iterator> cells_to_check;
1500 * neighbor_ids.insert(cell->active_cell_index());
1501 * cells_to_check.insert(cell);
1503 *
bool new_neighbors_found;
1506 * new_neighbors_found =
false;
1507 *
for (
const auto &check_cell :
1509 * cells_to_check.begin(), cells_to_check.end()))
1511 *
for (
const auto n : check_cell->face_indices())
1513 *
if (!(check_cell->face(n)->at_boundary()))
1515 *
const auto & neighbor = check_cell->neighbor(n);
1516 *
const double distance =
1517 * cell->center().distance(neighbor->center());
1518 *
if ((distance < filter_r) &&
1519 * !(neighbor_ids.count(neighbor->active_cell_index())))
1521 * cells_to_check.insert(neighbor);
1522 * neighbor_ids.insert(neighbor->active_cell_index());
1523 * new_neighbors_found =
true;
1529 *
while (new_neighbors_found);
1530 *
return cells_to_check;
1536 * <a name=
"AssemblingtheNewtonmatrix"></a>
1537 * <h3>Assembling the Newton
matrix</h3>
1541 * Whereas the setup_filter_matrix function built a
matrix that is the same as
1542 *
long as the mesh does not change (which we don
't do anyway in
1543 * this program), the next function builds the matrix to be solved
1544 * in each iteration. This is where the magic happens. The components
1545 * of the system of linear equations describing Newton's method
for
1546 * finding the solution of the KKT conditions are implemented here.
1550 * The top of the function is as in most of these
functions and just
1551 * sets up all sorts of variables necessary
for the actual assembly,
1552 * including a whole bunch of extractors. The entire
set up should
1553 * look familiar, though somewhat lengthier,
if you
've previously
1554 * looked at @ref step_22 "step-22".
1557 * template <int dim>
1558 * void SANDTopOpt<dim>::assemble_system()
1560 * TimerOutput::Scope t(timer, "assembly");
1562 * system_matrix = 0;
1566 * MappingQGeneric<dim> mapping(1);
1567 * QGauss<dim> quadrature_formula(fe.degree + 1);
1568 * QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
1569 * FEValues<dim> fe_values(mapping,
1571 * quadrature_formula,
1572 * update_values | update_gradients |
1573 * update_quadrature_points | update_JxW_values);
1574 * FEFaceValues<dim> fe_face_values(mapping,
1576 * face_quadrature_formula,
1577 * update_values | update_quadrature_points |
1578 * update_normal_vectors |
1579 * update_JxW_values);
1581 * const unsigned int dofs_per_cell = fe.dofs_per_cell;
1582 * const unsigned int n_q_points = quadrature_formula.size();
1584 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1585 * Vector<double> dummy_cell_rhs(dofs_per_cell);
1587 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1589 * std::vector<double> lambda_values(n_q_points);
1590 * std::vector<double> mu_values(n_q_points);
1591 * const Functions::ConstantFunction<dim> lambda(1.);
1592 * const Functions::ConstantFunction<dim> mu(1.);
1593 * std::vector<Tensor<1, dim>> rhs_values(n_q_points);
1597 * At this point, we apply the filter to the unfiltered
1598 * density, and apply the adjoint (transpose) operation to the
1599 * unfiltered density multiplier, both to the current best
1600 * guess for the nonlinear solution. We use this later to tell
1601 * us how far off our filtered density is from the filter
1602 * applied to the unfiltered density. That is because while at
1603 * the solution of the nonlinear problem, we have
1604 * @f$\rho=H\varrho@f$, but at intermediate iterations, we in
1605 * general have @f$\rho^k\neq H\varrho^k@f$ and the "residual"
1606 * @f$\rho^k-H\varrho^k@f$ will then appear as the right hand side
1607 * of one of the Newton update equations that we compute
1611 * BlockVector<double> filtered_unfiltered_density_solution =
1612 * nonlinear_solution;
1613 * BlockVector<double> filter_adjoint_unfiltered_density_multiplier_solution =
1614 * nonlinear_solution;
1616 * filter_matrix.vmult(filtered_unfiltered_density_solution.block(
1617 * SolutionBlocks::unfiltered_density),
1618 * nonlinear_solution.block(
1619 * SolutionBlocks::unfiltered_density));
1620 * filter_matrix.Tvmult(
1621 * filter_adjoint_unfiltered_density_multiplier_solution.block(
1622 * SolutionBlocks::unfiltered_density_multiplier),
1623 * nonlinear_solution.block(SolutionBlocks::unfiltered_density_multiplier));
1626 * std::vector<double> old_density_values(n_q_points);
1627 * std::vector<Tensor<1, dim>> old_displacement_values(n_q_points);
1628 * std::vector<double> old_displacement_divs(n_q_points);
1629 * std::vector<SymmetricTensor<2, dim>> old_displacement_symmgrads(n_q_points);
1630 * std::vector<Tensor<1, dim>> old_displacement_multiplier_values(n_q_points);
1631 * std::vector<double> old_displacement_multiplier_divs(n_q_points);
1632 * std::vector<SymmetricTensor<2, dim>> old_displacement_multiplier_symmgrads(
1634 * std::vector<double> old_lower_slack_multiplier_values(n_q_points);
1635 * std::vector<double> old_upper_slack_multiplier_values(n_q_points);
1636 * std::vector<double> old_lower_slack_values(n_q_points);
1637 * std::vector<double> old_upper_slack_values(n_q_points);
1638 * std::vector<double> old_unfiltered_density_values(n_q_points);
1639 * std::vector<double> old_unfiltered_density_multiplier_values(n_q_points);
1640 * std::vector<double> filtered_unfiltered_density_values(n_q_points);
1641 * std::vector<double> filter_adjoint_unfiltered_density_multiplier_values(
1644 * using namespace ValueExtractors;
1645 * for (const auto &cell : dof_handler.active_cell_iterators())
1649 * cell->get_dof_indices(local_dof_indices);
1651 * fe_values.reinit(cell);
1653 * lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
1654 * mu.value_list(fe_values.get_quadrature_points(), mu_values);
1658 * As part of the construction of our system matrix, we need to
1659 * retrieve values from our current guess at the solution.
1660 * The following lines of code retrieve the needed values.
1663 * fe_values[densities<dim>].get_function_values(nonlinear_solution,
1664 * old_density_values);
1665 * fe_values[displacements<dim>].get_function_values(
1666 * nonlinear_solution, old_displacement_values);
1667 * fe_values[displacements<dim>].get_function_divergences(
1668 * nonlinear_solution, old_displacement_divs);
1669 * fe_values[displacements<dim>].get_function_symmetric_gradients(
1670 * nonlinear_solution, old_displacement_symmgrads);
1671 * fe_values[displacement_multipliers<dim>].get_function_values(
1672 * nonlinear_solution, old_displacement_multiplier_values);
1673 * fe_values[displacement_multipliers<dim>].get_function_divergences(
1674 * nonlinear_solution, old_displacement_multiplier_divs);
1675 * fe_values[displacement_multipliers<dim>]
1676 * .get_function_symmetric_gradients(
1677 * nonlinear_solution, old_displacement_multiplier_symmgrads);
1678 * fe_values[density_lower_slacks<dim>].get_function_values(
1679 * nonlinear_solution, old_lower_slack_values);
1680 * fe_values[density_lower_slack_multipliers<dim>].get_function_values(
1681 * nonlinear_solution, old_lower_slack_multiplier_values);
1682 * fe_values[density_upper_slacks<dim>].get_function_values(
1683 * nonlinear_solution, old_upper_slack_values);
1684 * fe_values[density_upper_slack_multipliers<dim>].get_function_values(
1685 * nonlinear_solution, old_upper_slack_multiplier_values);
1686 * fe_values[unfiltered_densities<dim>].get_function_values(
1687 * nonlinear_solution, old_unfiltered_density_values);
1688 * fe_values[unfiltered_density_multipliers<dim>].get_function_values(
1689 * nonlinear_solution, old_unfiltered_density_multiplier_values);
1690 * fe_values[unfiltered_densities<dim>].get_function_values(
1691 * filtered_unfiltered_density_solution,
1692 * filtered_unfiltered_density_values);
1693 * fe_values[unfiltered_density_multipliers<dim>].get_function_values(
1694 * filter_adjoint_unfiltered_density_multiplier_solution,
1695 * filter_adjoint_unfiltered_density_multiplier_values);
1697 * for (const auto q_point : fe_values.quadrature_point_indices())
1701 * We need several more values corresponding to the test functions
1702 * coming from the first derivatives taken from the Lagrangian,
1703 * that is the @f$d_{\bullet}@f$ functions. These are calculated here:
1706 * for (const auto i : fe_values.dof_indices())
1708 * const SymmetricTensor<2, dim> displacement_phi_i_symmgrad =
1709 * fe_values[displacements<dim>].symmetric_gradient(i, q_point);
1710 * const double displacement_phi_i_div =
1711 * fe_values[displacements<dim>].divergence(i, q_point);
1713 * const SymmetricTensor<2, dim>
1714 * displacement_multiplier_phi_i_symmgrad =
1715 * fe_values[displacement_multipliers<dim>].symmetric_gradient(
1717 * const double displacement_multiplier_phi_i_div =
1718 * fe_values[displacement_multipliers<dim>].divergence(i,
1721 * const double density_phi_i =
1722 * fe_values[densities<dim>].value(i, q_point);
1723 * const double unfiltered_density_phi_i =
1724 * fe_values[unfiltered_densities<dim>].value(i, q_point);
1725 * const double unfiltered_density_multiplier_phi_i =
1726 * fe_values[unfiltered_density_multipliers<dim>].value(i,
1729 * const double lower_slack_multiplier_phi_i =
1730 * fe_values[density_lower_slack_multipliers<dim>].value(
1733 * const double lower_slack_phi_i =
1734 * fe_values[density_lower_slacks<dim>].value(i, q_point);
1736 * const double upper_slack_phi_i =
1737 * fe_values[density_upper_slacks<dim>].value(i, q_point);
1739 * const double upper_slack_multiplier_phi_i =
1740 * fe_values[density_upper_slack_multipliers<dim>].value(
1744 * for (const auto j : fe_values.dof_indices())
1748 * Finally, we need values that come from the second round
1749 * of derivatives taken from the Lagrangian,
1750 * the @f$c_{\bullet}@f$ functions. These are calculated here:
1753 * const SymmetricTensor<2, dim> displacement_phi_j_symmgrad =
1754 * fe_values[displacements<dim>].symmetric_gradient(j,
1756 * const double displacement_phi_j_div =
1757 * fe_values[displacements<dim>].divergence(j, q_point);
1759 * const SymmetricTensor<2, dim>
1760 * displacement_multiplier_phi_j_symmgrad =
1761 * fe_values[displacement_multipliers<dim>]
1762 * .symmetric_gradient(j, q_point);
1763 * const double displacement_multiplier_phi_j_div =
1764 * fe_values[displacement_multipliers<dim>].divergence(
1767 * const double density_phi_j =
1768 * fe_values[densities<dim>].value(j, q_point);
1770 * const double unfiltered_density_phi_j =
1771 * fe_values[unfiltered_densities<dim>].value(j, q_point);
1772 * const double unfiltered_density_multiplier_phi_j =
1773 * fe_values[unfiltered_density_multipliers<dim>].value(
1777 * const double lower_slack_phi_j =
1778 * fe_values[density_lower_slacks<dim>].value(j, q_point);
1780 * const double upper_slack_phi_j =
1781 * fe_values[density_upper_slacks<dim>].value(j, q_point);
1783 * const double lower_slack_multiplier_phi_j =
1784 * fe_values[density_lower_slack_multipliers<dim>].value(
1787 * const double upper_slack_multiplier_phi_j =
1788 * fe_values[density_upper_slack_multipliers<dim>].value(
1793 * This is where the actual work starts. In
1794 * the following, we will build all of the
1795 * terms of the matrix -- they are numerous
1796 * and not entirely self-explanatory, also
1797 * depending on the previous solutions and its
1798 * derivatives (which we have already
1799 * evaluated above and put into the variables
1800 * called `old_*`). To understand what each of
1801 * these terms corresponds to, you will want
1802 * to look at the explicit form of these terms
1803 * in the introduction above.
1807 * The right hand sides of the equations being
1808 * driven to 0 give all the KKT conditions
1809 * for finding a local minimum -- the descriptions of what
1810 * each individual equation are given with the computations
1811 * of the right hand side.
1818 * cell_matrix(i, j) +=
1819 * fe_values.JxW(q_point) *
1822 * -density_phi_i * unfiltered_density_multiplier_phi_j
1824 * + density_penalty_exponent *
1825 * (density_penalty_exponent - 1) *
1826 * std::pow(old_density_values[q_point],
1827 * density_penalty_exponent - 2) *
1828 * density_phi_i * density_phi_j *
1829 * (old_displacement_multiplier_divs[q_point] *
1830 * old_displacement_divs[q_point] *
1831 * lambda_values[q_point] +
1832 * 2 * mu_values[q_point] *
1833 * (old_displacement_symmgrads[q_point] *
1834 * old_displacement_multiplier_symmgrads[q_point]))
1836 * + density_penalty_exponent *
1837 * std::pow(old_density_values[q_point],
1838 * density_penalty_exponent - 1) *
1840 * (displacement_multiplier_phi_j_div *
1841 * old_displacement_divs[q_point] *
1842 * lambda_values[q_point] +
1843 * 2 * mu_values[q_point] *
1844 * (old_displacement_symmgrads[q_point] *
1845 * displacement_multiplier_phi_j_symmgrad))
1847 * + density_penalty_exponent *
1848 * std::pow(old_density_values[q_point],
1849 * density_penalty_exponent - 1) *
1851 * (displacement_phi_j_div *
1852 * old_displacement_multiplier_divs[q_point] *
1853 * lambda_values[q_point] +
1854 * 2 * mu_values[q_point] *
1855 * (old_displacement_multiplier_symmgrads[q_point] *
1856 * displacement_phi_j_symmgrad)));
1859 * cell_matrix(i, j) +=
1860 * fe_values.JxW(q_point) *
1861 * (density_penalty_exponent *
1862 * std::pow(old_density_values[q_point],
1863 * density_penalty_exponent - 1) *
1865 * (old_displacement_multiplier_divs[q_point] *
1866 * displacement_phi_i_div * lambda_values[q_point] +
1867 * 2 * mu_values[q_point] *
1868 * (old_displacement_multiplier_symmgrads[q_point] *
1869 * displacement_phi_i_symmgrad))
1871 * + std::pow(old_density_values[q_point],
1872 * density_penalty_exponent) *
1873 * (displacement_multiplier_phi_j_div *
1874 * displacement_phi_i_div * lambda_values[q_point] +
1875 * 2 * mu_values[q_point] *
1876 * (displacement_multiplier_phi_j_symmgrad *
1877 * displacement_phi_i_symmgrad))
1881 * /* Equation 3, which has to do with the filter and which is
1882 * * calculated elsewhere. */
1883 * cell_matrix(i, j) +=
1884 * fe_values.JxW(q_point) *
1885 * (-1 * unfiltered_density_phi_i *
1886 * lower_slack_multiplier_phi_j +
1887 * unfiltered_density_phi_i * upper_slack_multiplier_phi_j);
1890 * /* Equation 4: Primal feasibility */
1891 * cell_matrix(i, j) +=
1892 * fe_values.JxW(q_point) *
1895 * density_penalty_exponent *
1896 * std::pow(old_density_values[q_point],
1897 * density_penalty_exponent - 1) *
1899 * (old_displacement_divs[q_point] *
1900 * displacement_multiplier_phi_i_div *
1901 * lambda_values[q_point] +
1902 * 2 * mu_values[q_point] *
1903 * (old_displacement_symmgrads[q_point] *
1904 * displacement_multiplier_phi_i_symmgrad))
1906 * + std::pow(old_density_values[q_point],
1907 * density_penalty_exponent) *
1908 * (displacement_phi_j_div *
1909 * displacement_multiplier_phi_i_div *
1910 * lambda_values[q_point] +
1911 * 2 * mu_values[q_point] *
1912 * (displacement_phi_j_symmgrad *
1913 * displacement_multiplier_phi_i_symmgrad)));
1915 * /* Equation 5: Primal feasibility */
1916 * cell_matrix(i, j) +=
1917 * -1 * fe_values.JxW(q_point) *
1918 * lower_slack_multiplier_phi_i *
1919 * (unfiltered_density_phi_j - lower_slack_phi_j);
1921 * /* Equation 6: Primal feasibility */
1922 * cell_matrix(i, j) +=
1923 * -1 * fe_values.JxW(q_point) *
1924 * upper_slack_multiplier_phi_i *
1925 * (-1 * unfiltered_density_phi_j - upper_slack_phi_j);
1927 * /* Equation 7: Primal feasibility - the part with the filter
1928 * * is added later */
1929 * cell_matrix(i, j) += -1 * fe_values.JxW(q_point) *
1930 * unfiltered_density_multiplier_phi_i *
1933 * /* Equation 8: Complementary slackness */
1934 * cell_matrix(i, j) +=
1935 * fe_values.JxW(q_point) *
1936 * (lower_slack_phi_i * lower_slack_multiplier_phi_j
1938 * + lower_slack_phi_i * lower_slack_phi_j *
1939 * old_lower_slack_multiplier_values[q_point] /
1940 * old_lower_slack_values[q_point]);
1942 * /* Equation 9: Complementary slackness */
1943 * cell_matrix(i, j) +=
1944 * fe_values.JxW(q_point) *
1945 * (upper_slack_phi_i * upper_slack_multiplier_phi_j
1948 * + upper_slack_phi_i * upper_slack_phi_j *
1949 * old_upper_slack_multiplier_values[q_point] /
1950 * old_upper_slack_values[q_point]);
1957 * Now that we have everything assembled, all we have to
1958 * do is deal with the effect of (Dirichlet) boundary
1959 * conditions and other constraints. We incorporate the
1960 * former locally with just the contributions from the
1961 * current cell, and then let the AffineConstraint class
1962 * deal with the latter while copying contributions from
1963 * the current cell into the global linear system:
1966 * MatrixTools::local_apply_boundary_values(boundary_values,
1967 * local_dof_indices,
1972 * constraints.distribute_local_to_global(cell_matrix,
1973 * local_dof_indices,
1979 * Having accumulated all of the terms that belong
1980 * into the Newton matrix, we now also have to
1981 * compute the terms for the right hand side
1982 * (i.e., the negative residual). We already do this
1983 * in another function, and so we call that here:
1986 * system_rhs = calculate_test_rhs(nonlinear_solution);
1990 * Here we use the filter matrix we have already
1991 * constructed. We only need to integrate this filter applied
1992 * to test functions, which are piecewise constant, and so the
1993 * integration becomes a simple multiplication by the measure
1994 * of the cell. Iterating over the pre-made filter matrix
1995 * allows us to use the information about which cells are in
1996 * or out of the filter without repeatedly checking neighbor
2000 * for (const auto &cell : dof_handler.active_cell_iterators())
2002 * const unsigned int i = cell->active_cell_index();
2003 * for (typename SparseMatrix<double>::iterator iter =
2004 * filter_matrix.begin(i);
2005 * iter != filter_matrix.end(i);
2008 * const unsigned int j = iter->column();
2009 * const double value = iter->value() * cell->measure();
2012 * .block(SolutionBlocks::unfiltered_density_multiplier,
2013 * SolutionBlocks::unfiltered_density)
2014 * .add(i, j, value);
2016 * .block(SolutionBlocks::unfiltered_density,
2017 * SolutionBlocks::unfiltered_density_multiplier)
2018 * .add(j, i, value);
2027 * <a name="SolvingtheNewtonlinearsystem"></a>
2028 * <h3>Solving the Newton linear system</h3>
2035 * We will need to solve a linear system in each iteration. We use
2036 * a direct solver, for now -- this is clearly not an efficient
2037 * choice for a matrix that has so many non-zeroes, and it will
2038 * not scale to anything interesting. For "real" applications, we
2039 * will need an iterative solver but the complexity of the system
2040 * means that an iterative solver algorithm will take a good deal
2041 * of work. Because this is not the focus of the current program,
2042 * we simply stick with the direct solver we have here -- the
2043 * function follows the same structure as used in @ref step_29 "step-29".
2046 * template <int dim>
2047 * BlockVector<double> SANDTopOpt<dim>::solve()
2049 * TimerOutput::Scope t(timer, "solver");
2051 * BlockVector<double> linear_solution;
2052 * linear_solution.reinit(nonlinear_solution);
2054 * SparseDirectUMFPACK A_direct;
2055 * A_direct.initialize(system_matrix);
2056 * A_direct.vmult(linear_solution, system_rhs);
2058 * constraints.distribute(linear_solution);
2060 * return linear_solution;
2067 * <a name="Detailsoftheoptimizationalgorithm"></a>
2068 * <h3>Details of the optimization algorithm</h3>
2072 * The next several functions deal with specific parts of the
2073 * optimization algorithm, most notably with deciding whether the
2074 * direction computed by solving the linearized (Newton) system is
2075 * viable and, if so, how far we want to go in this direction.
2080 * <a name="Computingsteplengths"></a>
2081 * <h4>Computing step lengths</h4>
2085 * We start with a function that does a binary search to figure
2086 * out the maximum step that meets the dual feasibility -- that
2087 * is, how far can we go so that @f$s>0@f$ and @f$z>0@f$. The function
2088 * returns a pair of values, one each for the @f$s@f$ and @f$z@f$ slack
2092 * template <int dim>
2093 * std::pair<double, double> SANDTopOpt<dim>::calculate_max_step_size(
2094 * const BlockVector<double> &state,
2095 * const BlockVector<double> &step) const
2097 * double fraction_to_boundary;
2098 * const double min_fraction_to_boundary = .8;
2099 * const double max_fraction_to_boundary = 1. - 1e-5;
2101 * if (min_fraction_to_boundary < 1 - barrier_size)
2103 * if (1 - barrier_size < max_fraction_to_boundary)
2104 * fraction_to_boundary = 1 - barrier_size;
2106 * fraction_to_boundary = max_fraction_to_boundary;
2109 * fraction_to_boundary = min_fraction_to_boundary;
2111 * double step_size_s_low = 0;
2112 * double step_size_z_low = 0;
2113 * double step_size_s_high = 1;
2114 * double step_size_z_high = 1;
2115 * double step_size_s, step_size_z;
2117 * const int max_bisection_method_steps = 50;
2118 * for (unsigned int k = 0; k < max_bisection_method_steps; ++k)
2120 * step_size_s = (step_size_s_low + step_size_s_high) / 2;
2121 * step_size_z = (step_size_z_low + step_size_z_high) / 2;
2123 * const BlockVector<double> state_test_s =
2124 * (fraction_to_boundary * state) + (step_size_s * step);
2126 * const BlockVector<double> state_test_z =
2127 * (fraction_to_boundary * state) + (step_size_z * step);
2129 * const bool accept_s =
2130 * (state_test_s.block(SolutionBlocks::density_lower_slack)
2131 * .is_non_negative()) &&
2132 * (state_test_s.block(SolutionBlocks::density_upper_slack)
2133 * .is_non_negative());
2134 * const bool accept_z =
2135 * (state_test_z.block(SolutionBlocks::density_lower_slack_multiplier)
2136 * .is_non_negative()) &&
2137 * (state_test_z.block(SolutionBlocks::density_upper_slack_multiplier)
2138 * .is_non_negative());
2141 * step_size_s_low = step_size_s;
2143 * step_size_s_high = step_size_s;
2146 * step_size_z_low = step_size_z;
2148 * step_size_z_high = step_size_z;
2151 * return {step_size_s_low, step_size_z_low};
2158 * <a name="Computingresiduals"></a>
2159 * <h4>Computing residuals</h4>
2163 * The next function computes a right hand side vector linearized
2164 * around a "test solution vector" that we can use to look at the
2165 * magnitude of the KKT conditions. This is then used for testing
2166 * the convergence before shrinking the barrier size, as well as in the
2167 * calculation of the @f$l_1@f$ merit.
2171 * The function is lengthy and complicated, but it is really just a
2172 * copy of the right hand side part of what the `assemble_system()`
2173 * function above did.
2176 * template <int dim>
2177 * BlockVector<double> SANDTopOpt<dim>::calculate_test_rhs(
2178 * const BlockVector<double> &test_solution) const
2182 * We first create a zero vector with size and blocking of system_rhs
2185 * BlockVector<double> test_rhs;
2186 * test_rhs.reinit(system_rhs);
2188 * MappingQGeneric<dim> mapping(1);
2189 * const QGauss<dim> quadrature_formula(fe.degree + 1);
2190 * const QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
2191 * FEValues<dim> fe_values(mapping,
2193 * quadrature_formula,
2194 * update_values | update_gradients |
2195 * update_quadrature_points | update_JxW_values);
2196 * FEFaceValues<dim> fe_face_values(mapping,
2198 * face_quadrature_formula,
2199 * update_values | update_quadrature_points |
2200 * update_normal_vectors |
2201 * update_JxW_values);
2203 * const unsigned int dofs_per_cell = fe.dofs_per_cell;
2204 * const unsigned int n_q_points = quadrature_formula.size();
2206 * Vector<double> cell_rhs(dofs_per_cell);
2207 * FullMatrix<double> dummy_cell_matrix(dofs_per_cell, dofs_per_cell);
2209 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2211 * std::vector<double> lambda_values(n_q_points);
2212 * std::vector<double> mu_values(n_q_points);
2214 * const Functions::ConstantFunction<dim> lambda(1.), mu(1.);
2215 * std::vector<Tensor<1, dim>> rhs_values(n_q_points);
2218 * BlockVector<double> filtered_unfiltered_density_solution = test_solution;
2219 * BlockVector<double> filter_adjoint_unfiltered_density_multiplier_solution =
2221 * filtered_unfiltered_density_solution.block(
2222 * SolutionBlocks::unfiltered_density) = 0;
2223 * filter_adjoint_unfiltered_density_multiplier_solution.block(
2224 * SolutionBlocks::unfiltered_density_multiplier) = 0;
2226 * filter_matrix.vmult(filtered_unfiltered_density_solution.block(
2227 * SolutionBlocks::unfiltered_density),
2228 * test_solution.block(
2229 * SolutionBlocks::unfiltered_density));
2230 * filter_matrix.Tvmult(
2231 * filter_adjoint_unfiltered_density_multiplier_solution.block(
2232 * SolutionBlocks::unfiltered_density_multiplier),
2233 * test_solution.block(SolutionBlocks::unfiltered_density_multiplier));
2236 * std::vector<double> old_density_values(n_q_points);
2237 * std::vector<Tensor<1, dim>> old_displacement_values(n_q_points);
2238 * std::vector<double> old_displacement_divs(n_q_points);
2239 * std::vector<SymmetricTensor<2, dim>> old_displacement_symmgrads(n_q_points);
2240 * std::vector<Tensor<1, dim>> old_displacement_multiplier_values(n_q_points);
2241 * std::vector<double> old_displacement_multiplier_divs(n_q_points);
2242 * std::vector<SymmetricTensor<2, dim>> old_displacement_multiplier_symmgrads(
2244 * std::vector<double> old_lower_slack_multiplier_values(n_q_points);
2245 * std::vector<double> old_upper_slack_multiplier_values(n_q_points);
2246 * std::vector<double> old_lower_slack_values(n_q_points);
2247 * std::vector<double> old_upper_slack_values(n_q_points);
2248 * std::vector<double> old_unfiltered_density_values(n_q_points);
2249 * std::vector<double> old_unfiltered_density_multiplier_values(n_q_points);
2250 * std::vector<double> filtered_unfiltered_density_values(n_q_points);
2251 * std::vector<double> filter_adjoint_unfiltered_density_multiplier_values(
2254 * using namespace ValueExtractors;
2255 * for (const auto &cell : dof_handler.active_cell_iterators())
2259 * cell->get_dof_indices(local_dof_indices);
2261 * fe_values.reinit(cell);
2263 * lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
2264 * mu.value_list(fe_values.get_quadrature_points(), mu_values);
2266 * fe_values[densities<dim>].get_function_values(test_solution,
2267 * old_density_values);
2268 * fe_values[displacements<dim>].get_function_values(
2269 * test_solution, old_displacement_values);
2270 * fe_values[displacements<dim>].get_function_divergences(
2271 * test_solution, old_displacement_divs);
2272 * fe_values[displacements<dim>].get_function_symmetric_gradients(
2273 * test_solution, old_displacement_symmgrads);
2274 * fe_values[displacement_multipliers<dim>].get_function_values(
2275 * test_solution, old_displacement_multiplier_values);
2276 * fe_values[displacement_multipliers<dim>].get_function_divergences(
2277 * test_solution, old_displacement_multiplier_divs);
2278 * fe_values[displacement_multipliers<dim>]
2279 * .get_function_symmetric_gradients(
2280 * test_solution, old_displacement_multiplier_symmgrads);
2281 * fe_values[density_lower_slacks<dim>].get_function_values(
2282 * test_solution, old_lower_slack_values);
2283 * fe_values[density_lower_slack_multipliers<dim>].get_function_values(
2284 * test_solution, old_lower_slack_multiplier_values);
2285 * fe_values[density_upper_slacks<dim>].get_function_values(
2286 * test_solution, old_upper_slack_values);
2287 * fe_values[density_upper_slack_multipliers<dim>].get_function_values(
2288 * test_solution, old_upper_slack_multiplier_values);
2289 * fe_values[unfiltered_densities<dim>].get_function_values(
2290 * test_solution, old_unfiltered_density_values);
2291 * fe_values[unfiltered_density_multipliers<dim>].get_function_values(
2292 * test_solution, old_unfiltered_density_multiplier_values);
2293 * fe_values[unfiltered_densities<dim>].get_function_values(
2294 * filtered_unfiltered_density_solution,
2295 * filtered_unfiltered_density_values);
2296 * fe_values[unfiltered_density_multipliers<dim>].get_function_values(
2297 * filter_adjoint_unfiltered_density_multiplier_solution,
2298 * filter_adjoint_unfiltered_density_multiplier_values);
2300 * for (const auto q_point : fe_values.quadrature_point_indices())
2302 * for (const auto i : fe_values.dof_indices())
2304 * const SymmetricTensor<2, dim> displacement_phi_i_symmgrad =
2305 * fe_values[displacements<dim>].symmetric_gradient(i, q_point);
2306 * const double displacement_phi_i_div =
2307 * fe_values[displacements<dim>].divergence(i, q_point);
2309 * const SymmetricTensor<2, dim>
2310 * displacement_multiplier_phi_i_symmgrad =
2311 * fe_values[displacement_multipliers<dim>].symmetric_gradient(
2313 * const double displacement_multiplier_phi_i_div =
2314 * fe_values[displacement_multipliers<dim>].divergence(i,
2318 * const double density_phi_i =
2319 * fe_values[densities<dim>].value(i, q_point);
2320 * const double unfiltered_density_phi_i =
2321 * fe_values[unfiltered_densities<dim>].value(i, q_point);
2322 * const double unfiltered_density_multiplier_phi_i =
2323 * fe_values[unfiltered_density_multipliers<dim>].value(i,
2326 * const double lower_slack_multiplier_phi_i =
2327 * fe_values[density_lower_slack_multipliers<dim>].value(
2330 * const double lower_slack_phi_i =
2331 * fe_values[density_lower_slacks<dim>].value(i, q_point);
2333 * const double upper_slack_phi_i =
2334 * fe_values[density_upper_slacks<dim>].value(i, q_point);
2336 * const double upper_slack_multiplier_phi_i =
2337 * fe_values[density_upper_slack_multipliers<dim>].value(
2340 * /* Equation 1: This equation, along with equations
2341 * * 2 and 3, are the variational derivatives of the
2342 * * Lagrangian with respect to the decision
2343 * * variables - the density, displacement, and
2344 * * unfiltered density. */
2346 * -1 * fe_values.JxW(q_point) *
2347 * (density_penalty_exponent *
2348 * std::pow(old_density_values[q_point],
2349 * density_penalty_exponent - 1) *
2351 * (old_displacement_multiplier_divs[q_point] *
2352 * old_displacement_divs[q_point] *
2353 * lambda_values[q_point] +
2354 * 2 * mu_values[q_point] *
2355 * (old_displacement_symmgrads[q_point] *
2356 * old_displacement_multiplier_symmgrads[q_point])) -
2358 * old_unfiltered_density_multiplier_values[q_point]);
2360 * /* Equation 2; the boundary terms will be added further down
2363 * -1 * fe_values.JxW(q_point) *
2364 * (std::pow(old_density_values[q_point],
2365 * density_penalty_exponent) *
2366 * (old_displacement_multiplier_divs[q_point] *
2367 * displacement_phi_i_div * lambda_values[q_point] +
2368 * 2 * mu_values[q_point] *
2369 * (old_displacement_multiplier_symmgrads[q_point] *
2370 * displacement_phi_i_symmgrad)));
2374 * -1 * fe_values.JxW(q_point) *
2375 * (unfiltered_density_phi_i *
2376 * filter_adjoint_unfiltered_density_multiplier_values
2378 * unfiltered_density_phi_i *
2379 * old_upper_slack_multiplier_values[q_point] +
2380 * -1 * unfiltered_density_phi_i *
2381 * old_lower_slack_multiplier_values[q_point]);
2385 * /* Equation 4; boundary term will again be dealt
2386 * * with below. This equation being driven to 0
2387 * * ensures that the elasticity equation is met as
2388 * * a constraint. */
2389 * cell_rhs(i) += -1 * fe_values.JxW(q_point) *
2390 * (std::pow(old_density_values[q_point],
2391 * density_penalty_exponent) *
2392 * (old_displacement_divs[q_point] *
2393 * displacement_multiplier_phi_i_div *
2394 * lambda_values[q_point] +
2395 * 2 * mu_values[q_point] *
2396 * (displacement_multiplier_phi_i_symmgrad *
2397 * old_displacement_symmgrads[q_point])));
2399 * /* Equation 5: This equation sets the lower slack
2400 * * variable equal to the unfiltered density,
2401 * * giving a minimum density of 0. */
2402 * cell_rhs(i) += fe_values.JxW(q_point) *
2403 * (lower_slack_multiplier_phi_i *
2404 * (old_unfiltered_density_values[q_point] -
2405 * old_lower_slack_values[q_point]));
2407 * /* Equation 6: This equation sets the upper slack
2408 * * variable equal to one minus the unfiltered
2410 * cell_rhs(i) += fe_values.JxW(q_point) *
2411 * (upper_slack_multiplier_phi_i *
2412 * (1 - old_unfiltered_density_values[q_point] -
2413 * old_upper_slack_values[q_point]));
2415 * /* Equation 7: This is the difference between the
2416 * * density and the filter applied to the
2417 * * unfiltered density. This being driven to 0 by
2418 * * the Newton steps ensures that the filter is
2419 * * applied correctly. */
2420 * cell_rhs(i) += fe_values.JxW(q_point) *
2421 * (unfiltered_density_multiplier_phi_i *
2422 * (old_density_values[q_point] -
2423 * filtered_unfiltered_density_values[q_point]));
2425 * /* Equation 8: This along with equation 9 give the
2426 * * requirement that s*z = \alpha for the barrier
2427 * * size alpha, and gives complementary slackness
2428 * * from KKT conditions when \alpha goes to 0. */
2430 * -1 * fe_values.JxW(q_point) *
2431 * (lower_slack_phi_i *
2432 * (old_lower_slack_multiplier_values[q_point] -
2433 * barrier_size / old_lower_slack_values[q_point]));
2437 * -1 * fe_values.JxW(q_point) *
2438 * (upper_slack_phi_i *
2439 * (old_upper_slack_multiplier_values[q_point] -
2440 * barrier_size / old_upper_slack_values[q_point]));
2444 * for (const auto &face : cell->face_iterators())
2446 * if (face->at_boundary() &&
2447 * face->boundary_id() == BoundaryIds::down_force)
2449 * fe_face_values.reinit(cell, face);
2451 * for (const auto face_q_point :
2452 * fe_face_values.quadrature_point_indices())
2454 * for (const auto i : fe_face_values.dof_indices())
2456 * Tensor<1, dim> traction;
2457 * traction[1] = -1.;
2461 * (traction * fe_face_values[displacements<dim>].value(
2462 * i, face_q_point)) *
2463 * fe_face_values.JxW(face_q_point);
2467 * fe_face_values[displacement_multipliers<dim>].value(
2468 * i, face_q_point)) *
2469 * fe_face_values.JxW(face_q_point);
2475 * MatrixTools::local_apply_boundary_values(boundary_values,
2476 * local_dof_indices,
2477 * dummy_cell_matrix,
2481 * constraints.distribute_local_to_global(cell_rhs,
2482 * local_dof_indices,
2493 * <a name="Computingthemeritfunction"></a>
2494 * <h4>Computing the merit function</h4>
2498 * The algorithm we use herein uses a "watchdog" strategy to
2499 * determine where and how far to go from the current iterate. We
2500 * base the watchdog strategy on an exact @f$l_1@f$ merit function. This
2501 * function calculates the exact @f$l_1@f$ merit of a given, putative,
2506 * The merit function consists of the sum of the objective function
2507 * (which is simply an integral of external forces (on the boundary
2508 * of the domain) times the displacement values of a test solution
2509 * (typically, the current solution plus some multiple of the Newton
2510 * update), and the @f$l_1@f$ norms of the Lagrange multiplier
2511 * components of residual vectors. The following code computes these
2515 * template <int dim>
2516 * double SANDTopOpt<dim>::calculate_exact_merit(
2517 * const BlockVector<double> &test_solution)
2519 * TimerOutput::Scope t(timer, "merit function");
2523 * Start with computing the objective function:
2526 * double objective_function_merit = 0;
2528 * MappingQGeneric<dim> mapping(1);
2529 * const QGauss<dim> quadrature_formula(fe.degree + 1);
2530 * const QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
2531 * FEValues<dim> fe_values(mapping,
2533 * quadrature_formula,
2534 * update_values | update_gradients |
2535 * update_quadrature_points | update_JxW_values);
2536 * FEFaceValues<dim> fe_face_values(mapping,
2538 * face_quadrature_formula,
2540 * update_quadrature_points |
2541 * update_normal_vectors |
2542 * update_JxW_values);
2544 * const unsigned int n_face_q_points = face_quadrature_formula.size();
2546 * std::vector<Tensor<1, dim>> displacement_face_values(n_face_q_points);
2548 * for (const auto &cell : dof_handler.active_cell_iterators())
2550 * for (const auto &face : cell->face_iterators())
2552 * if (face->at_boundary() &&
2553 * face->boundary_id() == BoundaryIds::down_force)
2555 * fe_face_values.reinit(cell, face);
2556 * fe_face_values[ValueExtractors::displacements<dim>]
2557 * .get_function_values(test_solution,
2558 * displacement_face_values);
2559 * for (unsigned int face_q_point = 0;
2560 * face_q_point < n_face_q_points;
2563 * Tensor<1, dim> traction;
2564 * traction[1] = -1.;
2566 * objective_function_merit +=
2567 * (traction * displacement_face_values[face_q_point]) *
2568 * fe_face_values.JxW(face_q_point);
2575 * for (const auto &cell : triangulation.active_cell_iterators())
2577 * objective_function_merit =
2578 * objective_function_merit -
2579 * barrier_size * cell->measure() *
2580 * std::log(test_solution.block(
2581 * SolutionBlocks::density_lower_slack)[cell->active_cell_index()]);
2582 * objective_function_merit =
2583 * objective_function_merit -
2584 * barrier_size * cell->measure() *
2585 * std::log(test_solution.block(
2586 * SolutionBlocks::density_upper_slack)[cell->active_cell_index()]);
2591 * Then compute the residual and take the @f$l_1@f$ norms of the
2592 * components that correspond to Lagrange mulipliers. We add
2593 * those to the objective function computed above, and return
2594 * the sum at the bottom:
2597 * const BlockVector<double> test_rhs = calculate_test_rhs(test_solution);
2599 * const double elasticity_constraint_merit =
2600 * penalty_multiplier *
2601 * test_rhs.block(SolutionBlocks::displacement_multiplier).l1_norm();
2602 * const double filter_constraint_merit =
2603 * penalty_multiplier *
2604 * test_rhs.block(SolutionBlocks::unfiltered_density_multiplier).l1_norm();
2605 * const double lower_slack_merit =
2606 * penalty_multiplier *
2607 * test_rhs.block(SolutionBlocks::density_lower_slack_multiplier).l1_norm();
2608 * const double upper_slack_merit =
2609 * penalty_multiplier *
2610 * test_rhs.block(SolutionBlocks::density_upper_slack_multiplier).l1_norm();
2612 * const double total_merit =
2613 * objective_function_merit + elasticity_constraint_merit +
2614 * filter_constraint_merit + lower_slack_merit + upper_slack_merit;
2615 * return total_merit;
2623 * <a name="Findingasearchdirection"></a>
2624 * <h4>Finding a search direction</h4>
2628 * Next up is the function that actually computes a search direction
2629 * starting at the current state (passed as the first argument) and
2630 * returns the resulting vector. To this end, the function first
2631 * calls the functions that assemble the linear system that
2632 * corresponds to the Newton system, and that solve it.
2636 * This function also updates the penalty multiplier in the merit
2637 * function, and then returns the largest scaled feasible step.
2638 * It uses the `calculate_max_step_sizes()` function to find the
2639 * largest feasible step that satisfies @f$s>0@f$ and @f$z>0@f$.
2645 * template <int dim>
2646 * BlockVector<double> SANDTopOpt<dim>::find_max_step()
2648 * assemble_system();
2649 * BlockVector<double> step = solve();
2653 * Next we are going to update penalty_multiplier. In
2654 * essence, a larger penalty multiplier makes us consider the
2655 * constraints more. Looking at the Hessian and gradient with
2656 * respect to the step we want to take with our decision
2657 * variables, and comparing that to the norm of our constraint
2658 * error gives us a way to ensure that our merit function is
2659 * "exact" - that is, it has a minimum in the same location
2660 * that the objective function does. As our merit function is
2661 * exact for any penalty multiplier over some minimum value,
2662 * we only keep the computed value if it increases the penalty
2669 * const std::vector<unsigned int> decision_variables = {
2670 * SolutionBlocks::density,
2671 * SolutionBlocks::displacement,
2672 * SolutionBlocks::unfiltered_density,
2673 * SolutionBlocks::density_upper_slack,
2674 * SolutionBlocks::density_lower_slack};
2675 * double hess_part = 0;
2676 * double grad_part = 0;
2677 * for (const unsigned int decision_variable_i : decision_variables)
2679 * for (const unsigned int decision_variable_j : decision_variables)
2681 * Vector<double> temp_vector(step.block(decision_variable_i).size());
2682 * system_matrix.block(decision_variable_i, decision_variable_j)
2683 * .vmult(temp_vector, step.block(decision_variable_j));
2684 * hess_part += step.block(decision_variable_i) * temp_vector;
2686 * grad_part -= system_rhs.block(decision_variable_i) *
2687 * step.block(decision_variable_i);
2690 * const std::vector<unsigned int> equality_constraint_multipliers = {
2691 * SolutionBlocks::displacement_multiplier,
2692 * SolutionBlocks::unfiltered_density_multiplier,
2693 * SolutionBlocks::density_lower_slack_multiplier,
2694 * SolutionBlocks::density_upper_slack_multiplier};
2695 * double constraint_norm = 0;
2696 * for (unsigned int multiplier_i : equality_constraint_multipliers)
2697 * constraint_norm += system_rhs.block(multiplier_i).linfty_norm();
2700 * double test_penalty_multiplier;
2701 * if (hess_part > 0)
2702 * test_penalty_multiplier =
2703 * (grad_part + .5 * hess_part) / (.05 * constraint_norm);
2705 * test_penalty_multiplier = (grad_part) / (.05 * constraint_norm);
2707 * penalty_multiplier = std::max(penalty_multiplier, test_penalty_multiplier);
2711 * Based on all of this, we can now compute step sizes for the
2712 * primal and dual (Lagrange multiplier) variables. Once we
2713 * have these, we scale the components of the solution vector,
2714 * and that is what this function returns.
2717 * const std::pair<double, double> max_step_sizes =
2718 * calculate_max_step_size(nonlinear_solution, step);
2719 * const double step_size_s = max_step_sizes.first;
2720 * const double step_size_z = max_step_sizes.second;
2722 * step.block(SolutionBlocks::density) *= step_size_s;
2723 * step.block(SolutionBlocks::displacement) *= step_size_s;
2724 * step.block(SolutionBlocks::unfiltered_density) *= step_size_s;
2725 * step.block(SolutionBlocks::displacement_multiplier) *= step_size_z;
2726 * step.block(SolutionBlocks::unfiltered_density_multiplier) *= step_size_z;
2727 * step.block(SolutionBlocks::density_lower_slack) *= step_size_s;
2728 * step.block(SolutionBlocks::density_lower_slack_multiplier) *= step_size_z;
2729 * step.block(SolutionBlocks::density_upper_slack) *= step_size_s;
2730 * step.block(SolutionBlocks::density_upper_slack_multiplier) *= step_size_z;
2740 * <a name="Computingascaledstep"></a>
2741 * <h4>Computing a scaled step</h4>
2745 * The next function then implements a back-tracking algorithm for a
2746 * line search. It keeps shrinking step size until it finds a step
2747 * where the merit is decreased, and then returns the new location
2748 * based on the current state vector, and the direction to go into,
2749 * times the step length.
2752 * template <int dim>
2753 * BlockVector<double>
2754 * SANDTopOpt<dim>::compute_scaled_step(const BlockVector<double> &state,
2755 * const BlockVector<double> &max_step,
2756 * const double descent_requirement)
2758 * const double merit_derivative =
2759 * (calculate_exact_merit(state + 1e-4 * max_step) -
2760 * calculate_exact_merit(state)) /
2762 * double step_size = 1;
2763 * unsigned int max_linesearch_iterations = 10;
2764 * for (unsigned int k = 0; k < max_linesearch_iterations; ++k)
2766 * if (calculate_exact_merit(state + step_size * max_step) <
2767 * calculate_exact_merit(state) +
2768 * step_size * descent_requirement * merit_derivative)
2771 * step_size = step_size / 2;
2773 * return state + (step_size * max_step);
2780 * <a name="Checkingforconvergence"></a>
2781 * <h4>Checking for convergence</h4>
2785 * The final auxiliary function in this block is the one that checks
2786 * to see if the KKT conditions are sufficiently met so that the
2787 * overall algorithm can lower the barrier size. It does so by
2788 * computing the @f$l_1@f$ norm of the residual, which is what
2789 * `calculate_test_rhs()` computes.
2792 * template <int dim>
2793 * bool SANDTopOpt<dim>::check_convergence(const BlockVector<double> &state)
2795 * const BlockVector<double> test_rhs = calculate_test_rhs(state);
2796 * const double test_rhs_norm = test_rhs.l1_norm();
2798 * const double convergence_condition = 1e-2;
2799 * const double target_norm = convergence_condition * barrier_size;
2801 * std::cout << " Checking convergence. Current rhs norm is "
2802 * << test_rhs_norm << ", target is " << target_norm << std::endl;
2804 * return (test_rhs_norm < target_norm);
2811 * <a name="Postprocessingthesolution"></a>
2812 * <h3>Postprocessing the solution</h3>
2816 * The first of the postprocessing functions outputs information
2817 * in a VTU file for visualization. It looks long, but it's really
2818 * just the same as what was done in @ref step_22
"step-22",
for example, just
2819 * with (a lot) more solution variables:
2822 * template <int dim>
2823 *
void SANDTopOpt<dim>::output_results(
const unsigned int iteration)
const
2825 * std::vector<std::string> solution_names(1,
"density");
2826 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
2827 * data_component_interpretation(
2829 *
for (
unsigned int i = 0; i < dim; ++i)
2831 * solution_names.emplace_back(
"displacement");
2832 * data_component_interpretation.push_back(
2835 * solution_names.emplace_back(
"unfiltered_density");
2836 * data_component_interpretation.push_back(
2838 *
for (
unsigned int i = 0; i < dim; ++i)
2840 * solution_names.emplace_back(
"displacement_multiplier");
2841 * data_component_interpretation.push_back(
2844 * solution_names.emplace_back(
"unfiltered_density_multiplier");
2845 * data_component_interpretation.push_back(
2847 * solution_names.emplace_back(
"low_slack");
2848 * data_component_interpretation.push_back(
2850 * solution_names.emplace_back(
"low_slack_multiplier");
2851 * data_component_interpretation.push_back(
2853 * solution_names.emplace_back(
"high_slack");
2854 * data_component_interpretation.push_back(
2856 * solution_names.emplace_back(
"high_slack_multiplier");
2857 * data_component_interpretation.push_back(
2865 * data_component_interpretation);
2868 * std::ofstream output(
"solution" +
std::to_string(iteration) +
".vtu");
2877 * printing. [STL](https:
2878 * files are made up of triangles and normal vectors, and we will
2879 * use it to show all of those cells with a density
value larger
2881 * to @f$z=0.25@f$, and then generating two triangles
for each face of
2882 * the cells with a sufficiently large density
value. The triangle
2883 * nodes must go counter-clockwise when looking from the outside,
2884 * and the normal vectors must be unit vectors pointing outwards,
2885 * which
requires a few checks.
2888 *
template <
int dim>
2889 *
void SANDTopOpt<dim>::write_as_stl()
2891 *
static_assert(dim == 2,
2892 *
"This function is not implemented for anything "
2893 *
"other than the 2d case.");
2895 * std::ofstream stlfile;
2896 * stlfile.open(
"bridge.stl");
2898 * stlfile <<
"solid bridge\n" << std::scientific;
2899 *
double height = .25;
2901 *
for (
const auto &cell : dof_handler.active_cell_iterators())
2903 *
if (nonlinear_solution.block(
2904 * SolutionBlocks::density)[cell->active_cell_index()] > 0.5)
2908 * We have now found a cell with a density
value larger
2909 * than
zero. Let us start by writing out the bottom
2910 * and top faces. Owing to the ordering issue mentioned
2911 * above, we have to make sure that we understand
2912 * whether a cell has a right- or left-handed
2913 * coordinate system. We
do this by interrogating the
2914 * directions of the two edges starting at vertex 0 and
2915 * whether they form a right-handed coordinate system.
2923 * {{edge_directions[0][0], edge_directions[0][1]},
2924 * {edge_directions[1][0], edge_directions[1][1]}});
2925 *
const bool is_right_handed_cell = (
determinant(edge_tensor) > 0);
2927 *
if (is_right_handed_cell)
2930 * stlfile <<
" facet normal " << 0.000000e+00 <<
" "
2931 * << 0.000000e+00 <<
" " << -1.000000e+00 <<
"\n";
2932 * stlfile <<
" outer loop\n";
2933 * stlfile <<
" vertex " << cell->vertex(0)[0] <<
" "
2934 * << cell->vertex(0)[1] <<
" " << 0.000000e+00 <<
"\n";
2935 * stlfile <<
" vertex " << cell->vertex(2)[0] <<
" "
2936 * << cell->vertex(2)[1] <<
" " << 0.000000e+00 <<
"\n";
2937 * stlfile <<
" vertex " << cell->vertex(1)[0] <<
" "
2938 * << cell->vertex(1)[1] <<
" " << 0.000000e+00 <<
"\n";
2939 * stlfile <<
" endloop\n";
2940 * stlfile <<
" endfacet\n";
2941 * stlfile <<
" facet normal " << 0.000000e+00 <<
" "
2942 * << 0.000000e+00 <<
" " << -1.000000e+00 <<
"\n";
2943 * stlfile <<
" outer loop\n";
2944 * stlfile <<
" vertex " << cell->vertex(1)[0] <<
" "
2945 * << cell->vertex(1)[1] <<
" " << 0.000000e+00 <<
"\n";
2946 * stlfile <<
" vertex " << cell->vertex(2)[0] <<
" "
2947 * << cell->vertex(2)[1] <<
" " << 0.000000e+00 <<
"\n";
2948 * stlfile <<
" vertex " << cell->vertex(3)[0] <<
" "
2949 * << cell->vertex(3)[1] <<
" " << 0.000000e+00 <<
"\n";
2950 * stlfile <<
" endloop\n";
2951 * stlfile <<
" endfacet\n";
2954 * stlfile <<
" facet normal " << 0.000000e+00 <<
" "
2955 * << 0.000000e+00 <<
" " << 1.000000e+00 <<
"\n";
2956 * stlfile <<
" outer loop\n";
2957 * stlfile <<
" vertex " << cell->vertex(0)[0] <<
" "
2958 * << cell->vertex(0)[1] <<
" " << height <<
"\n";
2959 * stlfile <<
" vertex " << cell->vertex(1)[0] <<
" "
2960 * << cell->vertex(1)[1] <<
" " << height <<
"\n";
2961 * stlfile <<
" vertex " << cell->vertex(2)[0] <<
" "
2962 * << cell->vertex(2)[1] <<
" " << height <<
"\n";
2963 * stlfile <<
" endloop\n";
2964 * stlfile <<
" endfacet\n";
2965 * stlfile <<
" facet normal " << 0.000000e+00 <<
" "
2966 * << 0.000000e+00 <<
" " << 1.000000e+00 <<
"\n";
2967 * stlfile <<
" outer loop\n";
2968 * stlfile <<
" vertex " << cell->vertex(1)[0] <<
" "
2969 * << cell->vertex(1)[1] <<
" " << height <<
"\n";
2970 * stlfile <<
" vertex " << cell->vertex(3)[0] <<
" "
2971 * << cell->vertex(3)[1] <<
" " << height <<
"\n";
2972 * stlfile <<
" vertex " << cell->vertex(2)[0] <<
" "
2973 * << cell->vertex(2)[1] <<
" " << height <<
"\n";
2974 * stlfile <<
" endloop\n";
2975 * stlfile <<
" endfacet\n";
2980 * stlfile <<
" facet normal " << 0.000000e+00 <<
" "
2981 * << 0.000000e+00 <<
" " << -1.000000e+00 <<
"\n";
2982 * stlfile <<
" outer loop\n";
2983 * stlfile <<
" vertex " << cell->vertex(0)[0] <<
" "
2984 * << cell->vertex(0)[1] <<
" " << 0.000000e+00 <<
"\n";
2985 * stlfile <<
" vertex " << cell->vertex(1)[0] <<
" "
2986 * << cell->vertex(1)[1] <<
" " << 0.000000e+00 <<
"\n";
2987 * stlfile <<
" vertex " << cell->vertex(2)[0] <<
" "
2988 * << cell->vertex(2)[1] <<
" " << 0.000000e+00 <<
"\n";
2989 * stlfile <<
" endloop\n";
2990 * stlfile <<
" endfacet\n";
2991 * stlfile <<
" facet normal " << 0.000000e+00 <<
" "
2992 * << 0.000000e+00 <<
" " << -1.000000e+00 <<
"\n";
2993 * stlfile <<
" outer loop\n";
2994 * stlfile <<
" vertex " << cell->vertex(1)[0] <<
" "
2995 * << cell->vertex(1)[1] <<
" " << 0.000000e+00 <<
"\n";
2996 * stlfile <<
" vertex " << cell->vertex(3)[0] <<
" "
2997 * << cell->vertex(3)[1] <<
" " << 0.000000e+00 <<
"\n";
2998 * stlfile <<
" vertex " << cell->vertex(2)[0] <<
" "
2999 * << cell->vertex(2)[1] <<
" " << 0.000000e+00 <<
"\n";
3000 * stlfile <<
" endloop\n";
3001 * stlfile <<
" endfacet\n";
3004 * stlfile <<
" facet normal " << 0.000000e+00 <<
" "
3005 * << 0.000000e+00 <<
" " << 1.000000e+00 <<
"\n";
3006 * stlfile <<
" outer loop\n";
3007 * stlfile <<
" vertex " << cell->vertex(0)[0] <<
" "
3008 * << cell->vertex(0)[1] <<
" " << height <<
"\n";
3009 * stlfile <<
" vertex " << cell->vertex(2)[0] <<
" "
3010 * << cell->vertex(2)[1] <<
" " << height <<
"\n";
3011 * stlfile <<
" vertex " << cell->vertex(1)[0] <<
" "
3012 * << cell->vertex(1)[1] <<
" " << height <<
"\n";
3013 * stlfile <<
" endloop\n";
3014 * stlfile <<
" endfacet\n";
3015 * stlfile <<
" facet normal " << 0.000000e+00 <<
" "
3016 * << 0.000000e+00 <<
" " << 1.000000e+00 <<
"\n";
3017 * stlfile <<
" outer loop\n";
3018 * stlfile <<
" vertex " << cell->vertex(1)[0] <<
" "
3019 * << cell->vertex(1)[1] <<
" " << height <<
"\n";
3020 * stlfile <<
" vertex " << cell->vertex(2)[0] <<
" "
3021 * << cell->vertex(2)[1] <<
" " << height <<
"\n";
3022 * stlfile <<
" vertex " << cell->vertex(3)[0] <<
" "
3023 * << cell->vertex(3)[1] <<
" " << height <<
"\n";
3024 * stlfile <<
" endloop\n";
3025 * stlfile <<
" endfacet\n";
3030 * Next we need to deal with the four faces of the
3031 * cell, extended into the @f$z@f$ direction. However, we
3032 * only need to write these faces
if either the face
3033 * is on the domain boundary, or
if it is the
3034 *
interface between a cell with density greater than
3035 * 0.5, and a cell with a density less than 0.5.
3038 *
for (
unsigned int face_number = 0;
3039 * face_number < GeometryInfo<dim>::faces_per_cell;
3043 * cell->face(face_number);
3045 *
if ((face->at_boundary()) ||
3046 * (!face->at_boundary() &&
3047 * (nonlinear_solution.block(
3048 * 0)[cell->neighbor(face_number)->active_cell_index()] <
3052 * (face->center() - cell->center());
3053 *
const double normal_norm = normal_vector.
norm();
3054 *
if ((face->vertex(0)[0] - face->vertex(0)[0]) *
3055 * (face->vertex(1)[1] - face->vertex(0)[1]) *
3057 * (face->vertex(0)[1] - face->vertex(0)[1]) * (0 - 0) *
3058 * normal_vector[0] +
3060 * (face->vertex(1)[0] - face->vertex(0)[0]) *
3061 * normal_vector[1] -
3062 * (face->vertex(0)[0] - face->vertex(0)[0]) * (0 - 0) *
3063 * normal_vector[1] -
3064 * (face->vertex(0)[1] - face->vertex(0)[1]) *
3065 * (face->vertex(1)[0] - face->vertex(0)[0]) *
3066 * normal_vector[0] -
3068 * (face->vertex(1)[1] - face->vertex(0)[1]) * 0 >
3071 * stlfile <<
" facet normal "
3072 * << normal_vector[0] / normal_norm <<
" "
3073 * << normal_vector[1] / normal_norm <<
" "
3074 * << 0.000000e+00 <<
"\n";
3075 * stlfile <<
" outer loop\n";
3076 * stlfile <<
" vertex " << face->vertex(0)[0]
3077 * <<
" " << face->vertex(0)[1] <<
" "
3078 * << 0.000000e+00 <<
"\n";
3079 * stlfile <<
" vertex " << face->vertex(0)[0]
3080 * <<
" " << face->vertex(0)[1] <<
" " << height
3082 * stlfile <<
" vertex " << face->vertex(1)[0]
3083 * <<
" " << face->vertex(1)[1] <<
" "
3084 * << 0.000000e+00 <<
"\n";
3085 * stlfile <<
" endloop\n";
3086 * stlfile <<
" endfacet\n";
3087 * stlfile <<
" facet normal "
3088 * << normal_vector[0] / normal_norm <<
" "
3089 * << normal_vector[1] / normal_norm <<
" "
3090 * << 0.000000e+00 <<
"\n";
3091 * stlfile <<
" outer loop\n";
3092 * stlfile <<
" vertex " << face->vertex(0)[0]
3093 * <<
" " << face->vertex(0)[1] <<
" " << height
3095 * stlfile <<
" vertex " << face->vertex(1)[0]
3096 * <<
" " << face->vertex(1)[1] <<
" " << height
3098 * stlfile <<
" vertex " << face->vertex(1)[0]
3099 * <<
" " << face->vertex(1)[1] <<
" "
3100 * << 0.000000e+00 <<
"\n";
3101 * stlfile <<
" endloop\n";
3102 * stlfile <<
" endfacet\n";
3106 * stlfile <<
" facet normal "
3107 * << normal_vector[0] / normal_norm <<
" "
3108 * << normal_vector[1] / normal_norm <<
" "
3109 * << 0.000000e+00 <<
"\n";
3110 * stlfile <<
" outer loop\n";
3111 * stlfile <<
" vertex " << face->vertex(0)[0]
3112 * <<
" " << face->vertex(0)[1] <<
" "
3113 * << 0.000000e+00 <<
"\n";
3114 * stlfile <<
" vertex " << face->vertex(1)[0]
3115 * <<
" " << face->vertex(1)[1] <<
" "
3116 * << 0.000000e+00 <<
"\n";
3117 * stlfile <<
" vertex " << face->vertex(0)[0]
3118 * <<
" " << face->vertex(0)[1] <<
" " << height
3120 * stlfile <<
" endloop\n";
3121 * stlfile <<
" endfacet\n";
3122 * stlfile <<
" facet normal "
3123 * << normal_vector[0] / normal_norm <<
" "
3124 * << normal_vector[1] / normal_norm <<
" "
3125 * << 0.000000e+00 <<
"\n";
3126 * stlfile <<
" outer loop\n";
3127 * stlfile <<
" vertex " << face->vertex(0)[0]
3128 * <<
" " << face->vertex(0)[1] <<
" " << height
3130 * stlfile <<
" vertex " << face->vertex(1)[0]
3131 * <<
" " << face->vertex(1)[1] <<
" "
3132 * << 0.000000e+00 <<
"\n";
3133 * stlfile <<
" vertex " << face->vertex(1)[0]
3134 * <<
" " << face->vertex(1)[1] <<
" " << height
3136 * stlfile <<
" endloop\n";
3137 * stlfile <<
" endfacet\n";
3143 * stlfile <<
"endsolid bridge";
3151 * <a name=
"Therunfunctiondrivingtheoverallalgorithm"></a>
3152 * <h3>The
run() function driving the overall algorithm</h3>
3156 * This function finally provides the overall driver logic. It is,
3157 * in the grand scheme of things, a rather complicated function
3158 * primarily because the optimization algorithm is difficult: It
3159 * isn't just about finding a Newton direction like in @ref step_15 "step-15" and
3160 * then going a fixed distance in this direction any more, but
3161 * instead about (i) determining what the optimal
log-barrier
3162 * penalty parameter should be in the current step, (ii) a
3163 * complicated algorithm to determine how far we want to go, and
3164 * other ingredients. Let us see how we can break this down into
3165 * smaller chunks in the following documentation.
3169 * The function starts out simple enough with
first setting up the
3170 * mesh, the
DoFHandler, and then the various linear algebra objects
3171 * necessary for the following:
3174 * template <
int dim>
3175 *
void SANDTopOpt<dim>::
run()
3177 * std::cout <<
"filter r is: " << filter_r << std::endl;
3184 * dof_handler.distribute_dofs(fe);
3187 * setup_boundary_values();
3188 * setup_block_system();
3189 * setup_filter_matrix();
3194 * We then
set a number of parameters that affect the
3195 *
log-barrier and line search components of the optimization
3199 * barrier_size = 25;
3200 *
const double min_barrier_size = .0005;
3202 *
const unsigned int max_uphill_steps = 8;
3203 *
const double descent_requirement = .0001;
3208 * Now start the principal iteration. The overall algorithm
3209 * works by
using an outer
loop in which we
loop until either
3210 * (i) the
log-barrier parameter has become small enough, or (ii)
3211 * we have reached convergence. In any
case, we terminate
if
3212 *
end up with too large a number of iterations. This overall
3213 * structure is encoded as a `
do { ... }
while (...)`
loop
3214 * where the convergence condition is at the bottom.
3217 *
unsigned int iteration_number = 0;
3218 *
const unsigned int max_iterations = 10000;
3222 * std::cout <<
"Starting outer step in iteration " << iteration_number
3223 * <<
" with barrier parameter " << barrier_size << std::endl;
3227 * Within
this outer
loop, we have an inner
loop in which we
3228 *
try to find an update direction
using the watchdog
3229 * algorithm described in the introduction.
3233 * The
general idea of the watchdog algorithm itself is
3234 *
this: For a maximum of `max_uphill_steps` (i.e., a
loop
3235 * within the
"inner loop" mentioned above) attempts, we use
3236 * `find_max_step()` to compute a Newton update step, and
3237 * add these up in the `nonlinear_solution` vector. In each of
3238 * these attempts (starting from the place reached at the
3239 *
end of the previous attempt), we check whether we have
3240 * reached a target
value of the merit function described
3241 * above. The target
value is computed based on where
this
3242 * algorithm starts (the `nonlinear_solution` at the beginning of
3243 * the watchdog
loop, saves as `watchdog_state`) and the
3244 *
first proposed direction provided by `find_max_step()` in
3245 * the
first go-around of
this loop (the `k==0`
case).
3250 * std::cout <<
" Starting inner step in iteration "
3251 * << iteration_number
3252 * <<
" with merit function penalty multiplier "
3253 * << penalty_multiplier << std::endl;
3255 *
bool watchdog_step_found =
false;
3259 *
double target_merit = numbers::signaling_nan<double>();
3260 *
double merit_derivative = numbers::signaling_nan<double>();
3262 *
for (
unsigned int k = 0; k < max_uphill_steps; ++k)
3264 * ++iteration_number;
3269 * first_step = update_step;
3270 * merit_derivative =
3271 * ((calculate_exact_merit(watchdog_state +
3272 * .0001 * first_step) -
3273 * calculate_exact_merit(watchdog_state)) /
3275 * target_merit = calculate_exact_merit(watchdog_state) +
3276 * descent_requirement * merit_derivative;
3279 * nonlinear_solution += update_step;
3280 *
const double current_merit =
3281 * calculate_exact_merit(nonlinear_solution);
3283 * std::cout <<
" current watchdog state merit is: "
3284 * << current_merit <<
"; target merit is "
3285 * << target_merit << std::endl;
3287 *
if (current_merit < target_merit)
3289 * watchdog_step_found =
true;
3290 * std::cout <<
" found workable step after " << k + 1
3291 * <<
" iterations" << std::endl;
3299 * The next part of the algorithm then depends on
3300 * whether the watchdog
loop above succeeded. If it
3301 * did, then we are satisfied and no further action is
3302 * necessary: We just stay where we are. If, however,
3303 * we took the maximal number of unsuccessful steps in
3304 * the
loop above, then we need to
do something
else,
3305 * and
this is what the following code block does.
3309 * Specifically, from the
final (unsuccessful) state
3310 * of the
loop above, we seek
one more update
3311 * direction and take what is called a
"stretch
3312 * step". If that stretch state satisfies a condition
3313 * involving the merit function, then we go there. On
3314 * the other hand,
if the stretch state is also
3315 * unacceptable (as all of the watchdog steps above
3316 * were), then we discard all of the watchdog steps
3317 * taken above and start over again where we had
3318 * started the watchdog iterations -- that place was
3319 * stored in the `watchdog_state` variable above. More
3320 * specifically, the conditions below
first test
3321 * whether we take a step from `watchdog_state` in
3322 * direction `first_step`, or whether we can
do one
3323 * more update from the stretch state to find a
new
3324 * place. It is possible that neither of these is
3325 * actually better than the state we started from at
3326 * the beginning of the watchdog algorithm, but even
3327 *
if that is so, that place clearly was a difficult
3328 * place to be in, and getting away to start the next
3329 * iteration from another place might be a useful
3330 * strategy to eventually converge.
3334 * We keep repeating the watchdog steps above along
3335 * with the logic below until
this inner iteration is
3336 *
finally converged (or
if we
run up against the
3337 * maximal number of iterations -- where we count the
3338 * number of linear solves as iterations and increment
3339 * the counter every time we
call `find_max_step()`
3340 * since that is where the linear solve actually
3341 * happens). In any
case, at the
end of each of these
3342 * inner iterations we also output the solution in a
3343 * form suitable
for visualization.
3349 *
if (watchdog_step_found ==
false)
3351 * ++iteration_number;
3354 * compute_scaled_step(nonlinear_solution,
3356 * descent_requirement);
3360 * If we did not get a successful watchdog step,
3361 * we now need to decide between going back to
3362 * where we started, or
using the
final state. We
3363 * compare the merits of both of these locations,
3364 * and then take a scaled step from whichever
3365 * location is better. As the scaled step is
3366 * guaranteed to lower the merit, we will
end up
3367 * keeping
one of the two.
3370 *
if ((calculate_exact_merit(nonlinear_solution) <
3371 * calculate_exact_merit(watchdog_state)) ||
3372 * (calculate_exact_merit(stretch_state) < target_merit))
3374 * std::cout <<
" Taking scaled step from end of watchdog"
3376 * nonlinear_solution = stretch_state;
3381 * <<
" Taking scaled step from beginning of watchdog"
3383 *
if (calculate_exact_merit(stretch_state) >
3384 * calculate_exact_merit(watchdog_state))
3386 * nonlinear_solution =
3387 * compute_scaled_step(watchdog_state,
3389 * descent_requirement);
3393 * ++iteration_number;
3394 * nonlinear_solution = stretch_state;
3397 * nonlinear_solution =
3398 * compute_scaled_step(nonlinear_solution,
3400 * descent_requirement);
3405 * output_results(iteration_number);
3407 *
while ((iteration_number < max_iterations) &&
3408 * (check_convergence(nonlinear_solution) ==
false));
3413 * At the
end of the outer
loop, we have to update the
3414 * barrier parameter,
for which we use the following
3415 * formula. The rest of the function is then simply about
3416 * checking the outer
loop convergence condition, and
if
3417 * we decide to terminate computations, about writing the
3418 *
final "design" as an STL file
for use in 3
d printing,
3419 * and to output some timing information.
3422 *
const double barrier_size_multiplier = .8;
3423 *
const double barrier_size_exponent = 1.2;
3427 *
std::pow(barrier_size, barrier_size_exponent)),
3428 * min_barrier_size);
3430 * std::cout << std::endl;
3432 *
while (((barrier_size > min_barrier_size) ||
3433 * (check_convergence(nonlinear_solution) ==
false)) &&
3434 * (iteration_number < max_iterations));
3437 * timer.print_summary();
3444 * <a name=
"Themainfunction"></a>
3445 * <h3>The main function</h3>
3449 * The remainder of the code, the `main()` function, is as usual:
3456 * SAND::SANDTopOpt<2> elastic_problem_2d;
3457 * elastic_problem_2d.run();
3459 *
catch (std::exception &exc)
3461 * std::cerr << std::endl
3463 * <<
"----------------------------------------------------"
3465 * std::cerr <<
"Exception on processing: " << std::endl
3466 * << exc.what() << std::endl
3467 * <<
"Aborting!" << std::endl
3468 * <<
"----------------------------------------------------"
3475 * std::cerr << std::endl
3477 * <<
"----------------------------------------------------"
3479 * std::cerr <<
"Unknown exception!" << std::endl
3480 * <<
"Aborting!" << std::endl
3481 * <<
"----------------------------------------------------"
3488<a name=
"Results"></a><h1>Results</h1>
3490<a name=
"TestProblem"></a><h3>Test Problem</h3>
3492The algorithms used above are tested against a traditional topology optimization
3493 problem called the Messerschmitt-Bolkow-Blohm Beam (MBB Beam).
3495This problem considers the optimal 2-
d structure that can be built on a
3496rectangle 6 units wide, and 1 unit tall. The bottom corners are fixed in place
3497in the @f$y@f$ direction
using a
zero Dirichlet boundary condition, and a downward
3498force is applied in the
center of the top of the beam by enforcing a Neumann
3499boundary condition. The rest of the boundary is allowed to move, and has no
3500external force applied, which takes the form of a
zero Neumann boundary
3501condition. In essence, we are asking the following question: How should we
3502design a bridge in a way so that
if the bottom left and bottom right
point of
3503the bridge are on rollers that allow these points to move horizontally but not
3504vertically, and so that the displacement in response to the vertical force in
3507While the total
volume of the domain is 6 units, 3 units of material are allowed
for
3508the structure. Because of the symmetry of the problem, it could be posed on a
3509rectangle of width 3 and height 1 by cutting the original domain in half, and
3510using zero Dirichlet boundary conditions in the @f$x@f$ direction along the cut
3511edge. That said, symmetry of the solution is a good indicator that the program
3512is working as expected, so we solved the problem on the whole domain,
3513as shown below. @cite Bendse2004
3515<div style=
"text-align:center;">
3516 <img src=
"https://www.dealii.org/images/steps/developer/step-79.mbbgeometry.png"
3517 alt=
"The MBB problem domain and boundary conditions">
3521Using the program discussed above, we find the minimum
volume of the MBB Beam and the
3522individual components of the solution look as follows:
3524<div
class=
"onecolumn" style=
"width: 80%; text-align: center;">
3526 <img src=
"https://www.dealii.org/images/steps/developer/step-79.filtereddensity.png"
3527 alt=
"Filtered Density Solution">
3530 <img src=
"https://www.dealii.org/images/steps/developer/step-79.unfiltereddensity.png"
3531 alt=
"Unfiltered Density Solution">
3536These pictures show that what we find here is in accordance with what
one
3537typically sees in other publications on the topic @cite Bendse2004. Maybe more interestingly, the
3538result looks like a truss bridge (except that we apply the load at the top of
3539the trusses, rather than the bottom as in real truss bridges, akin to a
"deck
3540truss" bridge), suggesting that the designs that have been used in bridge-
3541building
for centuries are indeed based on ideas we can now show to be optimal
3545<a name=
"Possibilitiesforextensions"></a><h4>Possibilities
for extensions</h4>
3548The results shown above took around 75 iterations to find, which is quite
3549concerning given the expense in solving the large linear systems in each
3550iteration. Looking at the evolution, it does look as though the convergence has
3551moments of happening quickly and moments of happening slowly. We believe
this to
3552be due to both a lack of precision on when and how to decrease the boundary
3553values, as well as our choice of merit function being sub-optimal. In the future,
3554a LOQO barrier update replacing the monotone
reduction, as well as a Markov
3555Filter in place of a merit function will decrease the number of necessary
3556iterations significantly.
3558The barrier decrease is most sensitive in the middle of the convergence, which
3559is problematic, as it seems like we need it to decrease quickly, then slowly,
3562Secondly, the linear solver used here is just the sparse direct solver based on
3564but the formulation of the optimization problem detailed above has quite a large
3565number of variables and so the linear problem is not only large but also has a
3566lot of
nonzero entries in many rows, even on meshes that overall are still
3567relatively coarse. As a consequence, the solver time dominates the
3568computations, and more sophisticated approaches at solving the linear system
3572<a name=
"PlainProg"></a>
3573<h1> The plain program</h1>
3574@include
"step-79.cc"
std::vector< bool > component_mask
void attach_dof_handler(const DoFHandlerType &)
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation=std::vector< DataComponentInterpretation::DataComponentInterpretation >())
virtual void build_patches(const unsigned int n_subdivisions=0)
size_type n_elements() const
size_type nth_index_in_set(const size_type local_index) const
const_iterator begin() const
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
numbers::NumberTraits< Number >::real_type norm() const
VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &x)
__global__ void reduction(Number *result, const Number *v, const size_type N)
__global__ void set(Number *val, const Number s, const size_type N)
std::string to_string(const T &t)
void write_vtu(std::ostream &out) const
typename ActiveSelector::face_iterator face_iterator
void loop(ITERATOR begin, typename identity< ITERATOR >::type 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, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ component_is_part_of_vector
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void create_triangulation(Triangulation< dim, dim > &tria, const AdditionalData &additional_data=AdditionalData())
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
static const types::blas_int zero
@ matrix
Contents is actually a matrix.
@ general
No special properties.
static const types::blas_int one
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
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)
void run(const Iterator &begin, const typename identity< Iterator >::type &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 reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)