1245 *
for (
unsigned int i = 0; i < dim; ++i)
1247 *
for (
unsigned int k = 0;
k < dim; ++
k)
1289 *
constraints.
clear();
1291 *
std::vector<std::pair<types::global_dof_index, double>>
1292 *
constraint_entries;
1293 *
constraint_entries.reserve(
density_dofs.n_elements() - 1);
1300 *
constraints.close();
1304 *
We can now
finally create
the sparsity pattern
for the
1328 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1330 *
const unsigned
int i = cell->active_cell_index();
1333 *
const double distance =
1334 *
cell->center().distance(
check_cell->center());
1356 *
for "visual debugging".
1359 *
sparsity_pattern.copy_from(
dsp);
1361 *
std::ofstream out(
"sparsity.plt");
1362 *
sparsity_pattern.print_gnuplot(out);
1364 *
system_matrix.reinit(sparsity_pattern);
1397 * <a name=
"step_79-Creatingthefiltermatrix"></a>
1418 *
template <
int dim>
1446 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1448 *
const unsigned int i = cell->active_cell_index();
1451 *
const double distance =
1452 *
cell->center().distance(
check_cell->center());
1490 *
template <
int dim>
1491 *
std::set<typename Triangulation<dim>::cell_iterator>
1496 *
std::set<typename Triangulation<dim>::cell_iterator>
cells_to_check;
1514 *
const double distance =
1515 *
cell->center().distance(neighbor->center());
1517 *
!(
neighbor_ids.count(neighbor->active_cell_index())))
1534 * <a name=
"step_79-AssemblingtheNewtonmatrix"></a>
1541 * this program), the next function builds the matrix to be solved
1542 * in each iteration. This is where the magic happens. The components
1543 * of the system of linear equations describing Newton's method
for
1552 * looked at @ref step_22 "step-22".
1555 * template <int dim>
1556 * void SANDTopOpt<dim>::assemble_system()
1558 * TimerOutput::Scope t(timer, "assembly");
1560 * system_matrix = 0;
1564 * const MappingQ<dim> mapping(1);
1565 * const QGauss<dim> quadrature_formula(fe.degree + 1);
1566 * const QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
1567 * FEValues<dim> fe_values(mapping,
1569 * quadrature_formula,
1570 * update_values | update_gradients |
1571 * update_quadrature_points | update_JxW_values);
1572 * FEFaceValues<dim> fe_face_values(mapping,
1574 * face_quadrature_formula,
1575 * update_values | update_quadrature_points |
1576 * update_normal_vectors |
1577 * update_JxW_values);
1579 * const unsigned int dofs_per_cell = fe.dofs_per_cell;
1580 * const unsigned int n_q_points = quadrature_formula.size();
1582 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1583 * Vector<double> dummy_cell_rhs(dofs_per_cell);
1585 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1587 * std::vector<double> lambda_values(n_q_points);
1588 * std::vector<double> mu_values(n_q_points);
1589 * const Functions::ConstantFunction<dim> lambda(1.);
1590 * const Functions::ConstantFunction<dim> mu(1.);
1591 * std::vector<Tensor<1, dim>> rhs_values(n_q_points);
1595 * At this point, we apply the filter to the unfiltered
1596 * density, and apply the adjoint (transpose) operation to the
1597 * unfiltered density multiplier, both to the current best
1598 * guess for the nonlinear solution. We use this later to tell
1599 * us how far off our filtered density is from the filter
1600 * applied to the unfiltered density. That is because while at
1601 * the solution of the nonlinear problem, we have
1602 * @f$\rho=H\varrho@f$, but at intermediate iterations, we in
1603 * general have @f$\rho^k\neq H\varrho^k@f$ and the "residual"
1604 * @f$\rho^k-H\varrho^k@f$ will then appear as the right hand side
1605 * of one of the Newton update equations that we compute
1609 * BlockVector<double> filtered_unfiltered_density_solution =
1610 * nonlinear_solution;
1611 * BlockVector<double> filter_adjoint_unfiltered_density_multiplier_solution =
1612 * nonlinear_solution;
1614 * filter_matrix.vmult(filtered_unfiltered_density_solution.block(
1615 * SolutionBlocks::unfiltered_density),
1616 * nonlinear_solution.block(
1617 * SolutionBlocks::unfiltered_density));
1618 * filter_matrix.Tvmult(
1619 * filter_adjoint_unfiltered_density_multiplier_solution.block(
1620 * SolutionBlocks::unfiltered_density_multiplier),
1621 * nonlinear_solution.block(SolutionBlocks::unfiltered_density_multiplier));
1624 * std::vector<double> old_density_values(n_q_points);
1625 * std::vector<Tensor<1, dim>> old_displacement_values(n_q_points);
1626 * std::vector<double> old_displacement_divs(n_q_points);
1627 * std::vector<SymmetricTensor<2, dim>> old_displacement_symmgrads(n_q_points);
1628 * std::vector<Tensor<1, dim>> old_displacement_multiplier_values(n_q_points);
1629 * std::vector<double> old_displacement_multiplier_divs(n_q_points);
1630 * std::vector<SymmetricTensor<2, dim>> old_displacement_multiplier_symmgrads(
1632 * std::vector<double> old_lower_slack_multiplier_values(n_q_points);
1633 * std::vector<double> old_upper_slack_multiplier_values(n_q_points);
1634 * std::vector<double> old_lower_slack_values(n_q_points);
1635 * std::vector<double> old_upper_slack_values(n_q_points);
1636 * std::vector<double> old_unfiltered_density_values(n_q_points);
1637 * std::vector<double> old_unfiltered_density_multiplier_values(n_q_points);
1638 * std::vector<double> filtered_unfiltered_density_values(n_q_points);
1639 * std::vector<double> filter_adjoint_unfiltered_density_multiplier_values(
1642 * using namespace ValueExtractors;
1643 * for (const auto &cell : dof_handler.active_cell_iterators())
1647 * cell->get_dof_indices(local_dof_indices);
1649 * fe_values.reinit(cell);
1651 * lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
1652 * mu.value_list(fe_values.get_quadrature_points(), mu_values);
1656 * As part of the construction of our system matrix, we need to
1657 * retrieve values from our current guess at the solution.
1658 * The following lines of code retrieve the needed values.
1661 * fe_values[densities<dim>].get_function_values(nonlinear_solution,
1662 * old_density_values);
1663 * fe_values[displacements<dim>].get_function_values(
1664 * nonlinear_solution, old_displacement_values);
1665 * fe_values[displacements<dim>].get_function_divergences(
1666 * nonlinear_solution, old_displacement_divs);
1667 * fe_values[displacements<dim>].get_function_symmetric_gradients(
1668 * nonlinear_solution, old_displacement_symmgrads);
1669 * fe_values[displacement_multipliers<dim>].get_function_values(
1670 * nonlinear_solution, old_displacement_multiplier_values);
1671 * fe_values[displacement_multipliers<dim>].get_function_divergences(
1672 * nonlinear_solution, old_displacement_multiplier_divs);
1673 * fe_values[displacement_multipliers<dim>]
1674 * .get_function_symmetric_gradients(
1675 * nonlinear_solution, old_displacement_multiplier_symmgrads);
1676 * fe_values[density_lower_slacks<dim>].get_function_values(
1677 * nonlinear_solution, old_lower_slack_values);
1678 * fe_values[density_lower_slack_multipliers<dim>].get_function_values(
1679 * nonlinear_solution, old_lower_slack_multiplier_values);
1680 * fe_values[density_upper_slacks<dim>].get_function_values(
1681 * nonlinear_solution, old_upper_slack_values);
1682 * fe_values[density_upper_slack_multipliers<dim>].get_function_values(
1683 * nonlinear_solution, old_upper_slack_multiplier_values);
1684 * fe_values[unfiltered_densities<dim>].get_function_values(
1685 * nonlinear_solution, old_unfiltered_density_values);
1686 * fe_values[unfiltered_density_multipliers<dim>].get_function_values(
1687 * nonlinear_solution, old_unfiltered_density_multiplier_values);
1688 * fe_values[unfiltered_densities<dim>].get_function_values(
1689 * filtered_unfiltered_density_solution,
1690 * filtered_unfiltered_density_values);
1691 * fe_values[unfiltered_density_multipliers<dim>].get_function_values(
1692 * filter_adjoint_unfiltered_density_multiplier_solution,
1693 * filter_adjoint_unfiltered_density_multiplier_values);
1695 * for (const auto q_point : fe_values.quadrature_point_indices())
1699 * We need several more values corresponding to the test functions
1700 * coming from the first derivatives taken from the Lagrangian,
1701 * that is the @f$d_{\bullet}@f$ functions. These are calculated here:
1704 * for (const auto i : fe_values.dof_indices())
1706 * const SymmetricTensor<2, dim> displacement_phi_i_symmgrad =
1707 * fe_values[displacements<dim>].symmetric_gradient(i, q_point);
1708 * const double displacement_phi_i_div =
1709 * fe_values[displacements<dim>].divergence(i, q_point);
1711 * const SymmetricTensor<2, dim>
1712 * displacement_multiplier_phi_i_symmgrad =
1713 * fe_values[displacement_multipliers<dim>].symmetric_gradient(
1715 * const double displacement_multiplier_phi_i_div =
1716 * fe_values[displacement_multipliers<dim>].divergence(i,
1719 * const double density_phi_i =
1720 * fe_values[densities<dim>].value(i, q_point);
1721 * const double unfiltered_density_phi_i =
1722 * fe_values[unfiltered_densities<dim>].value(i, q_point);
1723 * const double unfiltered_density_multiplier_phi_i =
1724 * fe_values[unfiltered_density_multipliers<dim>].value(i,
1727 * const double lower_slack_multiplier_phi_i =
1728 * fe_values[density_lower_slack_multipliers<dim>].value(
1731 * const double lower_slack_phi_i =
1732 * fe_values[density_lower_slacks<dim>].value(i, q_point);
1734 * const double upper_slack_phi_i =
1735 * fe_values[density_upper_slacks<dim>].value(i, q_point);
1737 * const double upper_slack_multiplier_phi_i =
1738 * fe_values[density_upper_slack_multipliers<dim>].value(
1742 * for (const auto j : fe_values.dof_indices())
1746 * Finally, we need values that come from the second round
1747 * of derivatives taken from the Lagrangian,
1748 * the @f$c_{\bullet}@f$ functions. These are calculated here:
1751 * const SymmetricTensor<2, dim> displacement_phi_j_symmgrad =
1752 * fe_values[displacements<dim>].symmetric_gradient(j,
1754 * const double displacement_phi_j_div =
1755 * fe_values[displacements<dim>].divergence(j, q_point);
1757 * const SymmetricTensor<2, dim>
1758 * displacement_multiplier_phi_j_symmgrad =
1759 * fe_values[displacement_multipliers<dim>]
1760 * .symmetric_gradient(j, q_point);
1761 * const double displacement_multiplier_phi_j_div =
1762 * fe_values[displacement_multipliers<dim>].divergence(
1765 * const double density_phi_j =
1766 * fe_values[densities<dim>].value(j, q_point);
1768 * const double unfiltered_density_phi_j =
1769 * fe_values[unfiltered_densities<dim>].value(j, q_point);
1770 * const double unfiltered_density_multiplier_phi_j =
1771 * fe_values[unfiltered_density_multipliers<dim>].value(
1775 * const double lower_slack_phi_j =
1776 * fe_values[density_lower_slacks<dim>].value(j, q_point);
1778 * const double upper_slack_phi_j =
1779 * fe_values[density_upper_slacks<dim>].value(j, q_point);
1781 * const double lower_slack_multiplier_phi_j =
1782 * fe_values[density_lower_slack_multipliers<dim>].value(
1785 * const double upper_slack_multiplier_phi_j =
1786 * fe_values[density_upper_slack_multipliers<dim>].value(
1791 * This is where the actual work starts. In
1792 * the following, we will build all of the
1793 * terms of the matrix -- they are numerous
1794 * and not entirely self-explanatory, also
1795 * depending on the previous solutions and its
1796 * derivatives (which we have already
1797 * evaluated above and put into the variables
1798 * called `old_*`). To understand what each of
1799 * these terms corresponds to, you will want
1800 * to look at the explicit form of these terms
1801 * in the introduction above.
1805 * The right hand sides of the equations being
1806 * driven to 0 give all the KKT conditions
1807 * for finding a local minimum -- the descriptions of what
1808 * each individual equation are given with the computations
1809 * of the right hand side.
1816 * cell_matrix(i, j) +=
1817 * fe_values.JxW(q_point) *
1820 * -density_phi_i * unfiltered_density_multiplier_phi_j
1822 * + density_penalty_exponent *
1823 * (density_penalty_exponent - 1) *
1824 * std::pow(old_density_values[q_point],
1825 * density_penalty_exponent - 2) *
1826 * density_phi_i * density_phi_j *
1827 * (old_displacement_multiplier_divs[q_point] *
1828 * old_displacement_divs[q_point] *
1829 * lambda_values[q_point] +
1830 * 2 * mu_values[q_point] *
1831 * (old_displacement_symmgrads[q_point] *
1832 * old_displacement_multiplier_symmgrads[q_point]))
1834 * + density_penalty_exponent *
1835 * std::pow(old_density_values[q_point],
1836 * density_penalty_exponent - 1) *
1838 * (displacement_multiplier_phi_j_div *
1839 * old_displacement_divs[q_point] *
1840 * lambda_values[q_point] +
1841 * 2 * mu_values[q_point] *
1842 * (old_displacement_symmgrads[q_point] *
1843 * displacement_multiplier_phi_j_symmgrad))
1845 * + density_penalty_exponent *
1846 * std::pow(old_density_values[q_point],
1847 * density_penalty_exponent - 1) *
1849 * (displacement_phi_j_div *
1850 * old_displacement_multiplier_divs[q_point] *
1851 * lambda_values[q_point] +
1852 * 2 * mu_values[q_point] *
1853 * (old_displacement_multiplier_symmgrads[q_point] *
1854 * displacement_phi_j_symmgrad)));
1857 * cell_matrix(i, j) +=
1858 * fe_values.JxW(q_point) *
1859 * (density_penalty_exponent *
1860 * std::pow(old_density_values[q_point],
1861 * density_penalty_exponent - 1) *
1863 * (old_displacement_multiplier_divs[q_point] *
1864 * displacement_phi_i_div * lambda_values[q_point] +
1865 * 2 * mu_values[q_point] *
1866 * (old_displacement_multiplier_symmgrads[q_point] *
1867 * displacement_phi_i_symmgrad))
1869 * + std::pow(old_density_values[q_point],
1870 * density_penalty_exponent) *
1871 * (displacement_multiplier_phi_j_div *
1872 * displacement_phi_i_div * lambda_values[q_point] +
1873 * 2 * mu_values[q_point] *
1874 * (displacement_multiplier_phi_j_symmgrad *
1875 * displacement_phi_i_symmgrad))
1879 * /* Equation 3, which has to do with the filter and which is
1880 * * calculated elsewhere. */
1881 * cell_matrix(i, j) +=
1882 * fe_values.JxW(q_point) *
1883 * (-1 * unfiltered_density_phi_i *
1884 * lower_slack_multiplier_phi_j +
1885 * unfiltered_density_phi_i * upper_slack_multiplier_phi_j);
1888 * /* Equation 4: Primal feasibility */
1889 * cell_matrix(i, j) +=
1890 * fe_values.JxW(q_point) *
1893 * density_penalty_exponent *
1894 * std::pow(old_density_values[q_point],
1895 * density_penalty_exponent - 1) *
1897 * (old_displacement_divs[q_point] *
1898 * displacement_multiplier_phi_i_div *
1899 * lambda_values[q_point] +
1900 * 2 * mu_values[q_point] *
1901 * (old_displacement_symmgrads[q_point] *
1902 * displacement_multiplier_phi_i_symmgrad))
1904 * + std::pow(old_density_values[q_point],
1905 * density_penalty_exponent) *
1906 * (displacement_phi_j_div *
1907 * displacement_multiplier_phi_i_div *
1908 * lambda_values[q_point] +
1909 * 2 * mu_values[q_point] *
1910 * (displacement_phi_j_symmgrad *
1911 * displacement_multiplier_phi_i_symmgrad)));
1913 * /* Equation 5: Primal feasibility */
1914 * cell_matrix(i, j) +=
1915 * -1 * fe_values.JxW(q_point) *
1916 * lower_slack_multiplier_phi_i *
1917 * (unfiltered_density_phi_j - lower_slack_phi_j);
1919 * /* Equation 6: Primal feasibility */
1920 * cell_matrix(i, j) +=
1921 * -1 * fe_values.JxW(q_point) *
1922 * upper_slack_multiplier_phi_i *
1923 * (-1 * unfiltered_density_phi_j - upper_slack_phi_j);
1925 * /* Equation 7: Primal feasibility - the part with the filter
1926 * * is added later */
1927 * cell_matrix(i, j) += -1 * fe_values.JxW(q_point) *
1928 * unfiltered_density_multiplier_phi_i *
1931 * /* Equation 8: Complementary slackness */
1932 * cell_matrix(i, j) +=
1933 * fe_values.JxW(q_point) *
1934 * (lower_slack_phi_i * lower_slack_multiplier_phi_j
1936 * + lower_slack_phi_i * lower_slack_phi_j *
1937 * old_lower_slack_multiplier_values[q_point] /
1938 * old_lower_slack_values[q_point]);
1940 * /* Equation 9: Complementary slackness */
1941 * cell_matrix(i, j) +=
1942 * fe_values.JxW(q_point) *
1943 * (upper_slack_phi_i * upper_slack_multiplier_phi_j
1946 * + upper_slack_phi_i * upper_slack_phi_j *
1947 * old_upper_slack_multiplier_values[q_point] /
1948 * old_upper_slack_values[q_point]);
1955 * Now that we have everything assembled, all we have to
1956 * do is deal with the effect of (Dirichlet) boundary
1957 * conditions and other constraints. We incorporate the
1958 * former locally with just the contributions from the
1959 * current cell, and then let the AffineConstraint class
1960 * deal with the latter while copying contributions from
1961 * the current cell into the global linear system:
1964 * MatrixTools::local_apply_boundary_values(boundary_values,
1965 * local_dof_indices,
1970 * constraints.distribute_local_to_global(cell_matrix,
1971 * local_dof_indices,
1977 * Having accumulated all of the terms that belong
1978 * into the Newton matrix, we now also have to
1979 * compute the terms for the right hand side
1980 * (i.e., the negative residual). We already do this
1981 * in another function, and so we call that here:
1984 * system_rhs = calculate_test_rhs(nonlinear_solution);
1988 * Here we use the filter matrix we have already
1989 * constructed. We only need to integrate this filter applied
1990 * to test functions, which are piecewise constant, and so the
1991 * integration becomes a simple multiplication by the measure
1992 * of the cell. Iterating over the pre-made filter matrix
1993 * allows us to use the information about which cells are in
1994 * or out of the filter without repeatedly checking neighbor
1998 * for (const auto &cell : dof_handler.active_cell_iterators())
2000 * const unsigned int i = cell->active_cell_index();
2001 * for (typename SparseMatrix<double>::iterator iter =
2002 * filter_matrix.begin(i);
2003 * iter != filter_matrix.end(i);
2006 * const unsigned int j = iter->column();
2007 * const double value = iter->value() * cell->measure();
2010 * .block(SolutionBlocks::unfiltered_density_multiplier,
2011 * SolutionBlocks::unfiltered_density)
2012 * .add(i, j, value);
2014 * .block(SolutionBlocks::unfiltered_density,
2015 * SolutionBlocks::unfiltered_density_multiplier)
2016 * .add(j, i, value);
2025 * <a name="step_79-SolvingtheNewtonlinearsystem"></a>
2026 * <h3>Solving the Newton linear system</h3>
2033 * We will need to solve a linear system in each iteration. We use
2034 * a direct solver, for now -- this is clearly not an efficient
2035 * choice for a matrix that has so many non-zeroes, and it will
2036 * not scale to anything interesting. For "real" applications, we
2037 * will need an iterative solver but the complexity of the system
2038 * means that an iterative solver algorithm will take a good deal
2039 * of work. Because this is not the focus of the current program,
2040 * we simply stick with the direct solver we have here -- the
2041 * function follows the same structure as used in @ref step_29 "step-29".
2044 * template <int dim>
2045 * BlockVector<double> SANDTopOpt<dim>::solve()
2047 * TimerOutput::Scope t(timer, "solver");
2049 * BlockVector<double> linear_solution;
2050 * linear_solution.reinit(nonlinear_solution);
2052 * SparseDirectUMFPACK A_direct;
2053 * A_direct.initialize(system_matrix);
2054 * A_direct.vmult(linear_solution, system_rhs);
2056 * constraints.distribute(linear_solution);
2058 * return linear_solution;
2065 * <a name="step_79-Detailsoftheoptimizationalgorithm"></a>
2066 * <h3>Details of the optimization algorithm</h3>
2070 * The next several functions deal with specific parts of the
2071 * optimization algorithm, most notably with deciding whether the
2072 * direction computed by solving the linearized (Newton) system is
2073 * viable and, if so, how far we want to go in this direction.
2078 * <a name="step_79-Computingsteplengths"></a>
2079 * <h4>Computing step lengths</h4>
2083 * We start with a function that does a binary search to figure
2084 * out the maximum step that meets the dual feasibility -- that
2085 * is, how far can we go so that @f$s>0@f$ and @f$z>0@f$. The function
2086 * returns a pair of values, one each for the @f$s@f$ and @f$z@f$ slack
2090 * template <int dim>
2091 * std::pair<double, double> SANDTopOpt<dim>::calculate_max_step_size(
2092 * const BlockVector<double> &state,
2093 * const BlockVector<double> &step) const
2095 * double fraction_to_boundary;
2096 * const double min_fraction_to_boundary = .8;
2097 * const double max_fraction_to_boundary = 1. - 1e-5;
2099 * if (min_fraction_to_boundary < 1 - barrier_size)
2101 * if (1 - barrier_size < max_fraction_to_boundary)
2102 * fraction_to_boundary = 1 - barrier_size;
2104 * fraction_to_boundary = max_fraction_to_boundary;
2107 * fraction_to_boundary = min_fraction_to_boundary;
2109 * double step_size_s_low = 0;
2110 * double step_size_z_low = 0;
2111 * double step_size_s_high = 1;
2112 * double step_size_z_high = 1;
2113 * double step_size_s, step_size_z;
2115 * const int max_bisection_method_steps = 50;
2116 * for (unsigned int k = 0; k < max_bisection_method_steps; ++k)
2118 * step_size_s = (step_size_s_low + step_size_s_high) / 2;
2119 * step_size_z = (step_size_z_low + step_size_z_high) / 2;
2121 * const BlockVector<double> state_test_s =
2122 * (fraction_to_boundary * state) + (step_size_s * step);
2124 * const BlockVector<double> state_test_z =
2125 * (fraction_to_boundary * state) + (step_size_z * step);
2127 * const bool accept_s =
2128 * (state_test_s.block(SolutionBlocks::density_lower_slack)
2129 * .is_non_negative()) &&
2130 * (state_test_s.block(SolutionBlocks::density_upper_slack)
2131 * .is_non_negative());
2132 * const bool accept_z =
2133 * (state_test_z.block(SolutionBlocks::density_lower_slack_multiplier)
2134 * .is_non_negative()) &&
2135 * (state_test_z.block(SolutionBlocks::density_upper_slack_multiplier)
2136 * .is_non_negative());
2139 * step_size_s_low = step_size_s;
2141 * step_size_s_high = step_size_s;
2144 * step_size_z_low = step_size_z;
2146 * step_size_z_high = step_size_z;
2149 * return {step_size_s_low, step_size_z_low};
2156 * <a name="step_79-Computingresiduals"></a>
2157 * <h4>Computing residuals</h4>
2161 * The next function computes a right hand side vector linearized
2162 * around a "test solution vector" that we can use to look at the
2163 * magnitude of the KKT conditions. This is then used for testing
2164 * the convergence before shrinking the barrier size, as well as in the
2165 * calculation of the @f$l_1@f$ merit.
2169 * The function is lengthy and complicated, but it is really just a
2170 * copy of the right hand side part of what the `assemble_system()`
2171 * function above did.
2174 * template <int dim>
2175 * BlockVector<double> SANDTopOpt<dim>::calculate_test_rhs(
2176 * const BlockVector<double> &test_solution) const
2180 * We first create a zero vector with size and blocking of system_rhs
2183 * BlockVector<double> test_rhs;
2184 * test_rhs.reinit(system_rhs);
2186 * const MappingQ<dim> mapping(1);
2187 * const QGauss<dim> quadrature_formula(fe.degree + 1);
2188 * const QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
2189 * FEValues<dim> fe_values(mapping,
2191 * quadrature_formula,
2192 * update_values | update_gradients |
2193 * update_quadrature_points | update_JxW_values);
2194 * FEFaceValues<dim> fe_face_values(mapping,
2196 * face_quadrature_formula,
2197 * update_values | update_quadrature_points |
2198 * update_normal_vectors |
2199 * update_JxW_values);
2201 * const unsigned int dofs_per_cell = fe.dofs_per_cell;
2202 * const unsigned int n_q_points = quadrature_formula.size();
2204 * Vector<double> cell_rhs(dofs_per_cell);
2205 * FullMatrix<double> dummy_cell_matrix(dofs_per_cell, dofs_per_cell);
2207 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2209 * std::vector<double> lambda_values(n_q_points);
2210 * std::vector<double> mu_values(n_q_points);
2212 * const Functions::ConstantFunction<dim> lambda(1.), mu(1.);
2213 * std::vector<Tensor<1, dim>> rhs_values(n_q_points);
2216 * BlockVector<double> filtered_unfiltered_density_solution = test_solution;
2217 * BlockVector<double> filter_adjoint_unfiltered_density_multiplier_solution =
2219 * filtered_unfiltered_density_solution.block(
2220 * SolutionBlocks::unfiltered_density) = 0;
2221 * filter_adjoint_unfiltered_density_multiplier_solution.block(
2222 * SolutionBlocks::unfiltered_density_multiplier) = 0;
2224 * filter_matrix.vmult(filtered_unfiltered_density_solution.block(
2225 * SolutionBlocks::unfiltered_density),
2226 * test_solution.block(
2227 * SolutionBlocks::unfiltered_density));
2228 * filter_matrix.Tvmult(
2229 * filter_adjoint_unfiltered_density_multiplier_solution.block(
2230 * SolutionBlocks::unfiltered_density_multiplier),
2231 * test_solution.block(SolutionBlocks::unfiltered_density_multiplier));
2234 * std::vector<double> old_density_values(n_q_points);
2235 * std::vector<Tensor<1, dim>> old_displacement_values(n_q_points);
2236 * std::vector<double> old_displacement_divs(n_q_points);
2237 * std::vector<SymmetricTensor<2, dim>> old_displacement_symmgrads(n_q_points);
2238 * std::vector<Tensor<1, dim>> old_displacement_multiplier_values(n_q_points);
2239 * std::vector<double> old_displacement_multiplier_divs(n_q_points);
2240 * std::vector<SymmetricTensor<2, dim>> old_displacement_multiplier_symmgrads(
2242 * std::vector<double> old_lower_slack_multiplier_values(n_q_points);
2243 * std::vector<double> old_upper_slack_multiplier_values(n_q_points);
2244 * std::vector<double> old_lower_slack_values(n_q_points);
2245 * std::vector<double> old_upper_slack_values(n_q_points);
2246 * std::vector<double> old_unfiltered_density_values(n_q_points);
2247 * std::vector<double> old_unfiltered_density_multiplier_values(n_q_points);
2248 * std::vector<double> filtered_unfiltered_density_values(n_q_points);
2249 * std::vector<double> filter_adjoint_unfiltered_density_multiplier_values(
2252 * using namespace ValueExtractors;
2253 * for (const auto &cell : dof_handler.active_cell_iterators())
2257 * cell->get_dof_indices(local_dof_indices);
2259 * fe_values.reinit(cell);
2261 * lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
2262 * mu.value_list(fe_values.get_quadrature_points(), mu_values);
2264 * fe_values[densities<dim>].get_function_values(test_solution,
2265 * old_density_values);
2266 * fe_values[displacements<dim>].get_function_values(
2267 * test_solution, old_displacement_values);
2268 * fe_values[displacements<dim>].get_function_divergences(
2269 * test_solution, old_displacement_divs);
2270 * fe_values[displacements<dim>].get_function_symmetric_gradients(
2271 * test_solution, old_displacement_symmgrads);
2272 * fe_values[displacement_multipliers<dim>].get_function_values(
2273 * test_solution, old_displacement_multiplier_values);
2274 * fe_values[displacement_multipliers<dim>].get_function_divergences(
2275 * test_solution, old_displacement_multiplier_divs);
2276 * fe_values[displacement_multipliers<dim>]
2277 * .get_function_symmetric_gradients(
2278 * test_solution, old_displacement_multiplier_symmgrads);
2279 * fe_values[density_lower_slacks<dim>].get_function_values(
2280 * test_solution, old_lower_slack_values);
2281 * fe_values[density_lower_slack_multipliers<dim>].get_function_values(
2282 * test_solution, old_lower_slack_multiplier_values);
2283 * fe_values[density_upper_slacks<dim>].get_function_values(
2284 * test_solution, old_upper_slack_values);
2285 * fe_values[density_upper_slack_multipliers<dim>].get_function_values(
2286 * test_solution, old_upper_slack_multiplier_values);
2287 * fe_values[unfiltered_densities<dim>].get_function_values(
2288 * test_solution, old_unfiltered_density_values);
2289 * fe_values[unfiltered_density_multipliers<dim>].get_function_values(
2290 * test_solution, old_unfiltered_density_multiplier_values);
2291 * fe_values[unfiltered_densities<dim>].get_function_values(
2292 * filtered_unfiltered_density_solution,
2293 * filtered_unfiltered_density_values);
2294 * fe_values[unfiltered_density_multipliers<dim>].get_function_values(
2295 * filter_adjoint_unfiltered_density_multiplier_solution,
2296 * filter_adjoint_unfiltered_density_multiplier_values);
2298 * for (const auto q_point : fe_values.quadrature_point_indices())
2300 * for (const auto i : fe_values.dof_indices())
2302 * const SymmetricTensor<2, dim> displacement_phi_i_symmgrad =
2303 * fe_values[displacements<dim>].symmetric_gradient(i, q_point);
2304 * const double displacement_phi_i_div =
2305 * fe_values[displacements<dim>].divergence(i, q_point);
2307 * const SymmetricTensor<2, dim>
2308 * displacement_multiplier_phi_i_symmgrad =
2309 * fe_values[displacement_multipliers<dim>].symmetric_gradient(
2311 * const double displacement_multiplier_phi_i_div =
2312 * fe_values[displacement_multipliers<dim>].divergence(i,
2316 * const double density_phi_i =
2317 * fe_values[densities<dim>].value(i, q_point);
2318 * const double unfiltered_density_phi_i =
2319 * fe_values[unfiltered_densities<dim>].value(i, q_point);
2320 * const double unfiltered_density_multiplier_phi_i =
2321 * fe_values[unfiltered_density_multipliers<dim>].value(i,
2324 * const double lower_slack_multiplier_phi_i =
2325 * fe_values[density_lower_slack_multipliers<dim>].value(
2328 * const double lower_slack_phi_i =
2329 * fe_values[density_lower_slacks<dim>].value(i, q_point);
2331 * const double upper_slack_phi_i =
2332 * fe_values[density_upper_slacks<dim>].value(i, q_point);
2334 * const double upper_slack_multiplier_phi_i =
2335 * fe_values[density_upper_slack_multipliers<dim>].value(
2338 * /* Equation 1: This equation, along with equations
2339 * * 2 and 3, are the variational derivatives of the
2340 * * Lagrangian with respect to the decision
2341 * * variables - the density, displacement, and
2342 * * unfiltered density. */
2344 * -1 * fe_values.JxW(q_point) *
2345 * (density_penalty_exponent *
2346 * std::pow(old_density_values[q_point],
2347 * density_penalty_exponent - 1) *
2349 * (old_displacement_multiplier_divs[q_point] *
2350 * old_displacement_divs[q_point] *
2351 * lambda_values[q_point] +
2352 * 2 * mu_values[q_point] *
2353 * (old_displacement_symmgrads[q_point] *
2354 * old_displacement_multiplier_symmgrads[q_point])) -
2356 * old_unfiltered_density_multiplier_values[q_point]);
2358 * /* Equation 2; the boundary terms will be added further down
2361 * -1 * fe_values.JxW(q_point) *
2362 * (std::pow(old_density_values[q_point],
2363 * density_penalty_exponent) *
2364 * (old_displacement_multiplier_divs[q_point] *
2365 * displacement_phi_i_div * lambda_values[q_point] +
2366 * 2 * mu_values[q_point] *
2367 * (old_displacement_multiplier_symmgrads[q_point] *
2368 * displacement_phi_i_symmgrad)));
2372 * -1 * fe_values.JxW(q_point) *
2373 * (unfiltered_density_phi_i *
2374 * filter_adjoint_unfiltered_density_multiplier_values
2376 * unfiltered_density_phi_i *
2377 * old_upper_slack_multiplier_values[q_point] +
2378 * -1 * unfiltered_density_phi_i *
2379 * old_lower_slack_multiplier_values[q_point]);
2383 * /* Equation 4; boundary term will again be dealt
2384 * * with below. This equation being driven to 0
2385 * * ensures that the elasticity equation is met as
2386 * * a constraint. */
2387 * cell_rhs(i) += -1 * fe_values.JxW(q_point) *
2388 * (std::pow(old_density_values[q_point],
2389 * density_penalty_exponent) *
2390 * (old_displacement_divs[q_point] *
2391 * displacement_multiplier_phi_i_div *
2392 * lambda_values[q_point] +
2393 * 2 * mu_values[q_point] *
2394 * (displacement_multiplier_phi_i_symmgrad *
2395 * old_displacement_symmgrads[q_point])));
2397 * /* Equation 5: This equation sets the lower slack
2398 * * variable equal to the unfiltered density,
2399 * * giving a minimum density of 0. */
2400 * cell_rhs(i) += fe_values.JxW(q_point) *
2401 * (lower_slack_multiplier_phi_i *
2402 * (old_unfiltered_density_values[q_point] -
2403 * old_lower_slack_values[q_point]));
2405 * /* Equation 6: This equation sets the upper slack
2406 * * variable equal to one minus the unfiltered
2408 * cell_rhs(i) += fe_values.JxW(q_point) *
2409 * (upper_slack_multiplier_phi_i *
2410 * (1 - old_unfiltered_density_values[q_point] -
2411 * old_upper_slack_values[q_point]));
2413 * /* Equation 7: This is the difference between the
2414 * * density and the filter applied to the
2415 * * unfiltered density. This being driven to 0 by
2416 * * the Newton steps ensures that the filter is
2417 * * applied correctly. */
2418 * cell_rhs(i) += fe_values.JxW(q_point) *
2419 * (unfiltered_density_multiplier_phi_i *
2420 * (old_density_values[q_point] -
2421 * filtered_unfiltered_density_values[q_point]));
2423 * /* Equation 8: This along with equation 9 give the
2424 * * requirement that s*z = \alpha for the barrier
2425 * * size alpha, and gives complementary slackness
2426 * * from KKT conditions when \alpha goes to 0. */
2428 * -1 * fe_values.JxW(q_point) *
2429 * (lower_slack_phi_i *
2430 * (old_lower_slack_multiplier_values[q_point] -
2431 * barrier_size / old_lower_slack_values[q_point]));
2435 * -1 * fe_values.JxW(q_point) *
2436 * (upper_slack_phi_i *
2437 * (old_upper_slack_multiplier_values[q_point] -
2438 * barrier_size / old_upper_slack_values[q_point]));
2442 * for (const auto &face : cell->face_iterators())
2444 * if (face->at_boundary() &&
2445 * face->boundary_id() == BoundaryIds::down_force)
2447 * fe_face_values.reinit(cell, face);
2449 * for (const auto face_q_point :
2450 * fe_face_values.quadrature_point_indices())
2452 * for (const auto i : fe_face_values.dof_indices())
2454 * Tensor<1, dim> traction;
2455 * traction[1] = -1.;
2459 * (traction * fe_face_values[displacements<dim>].value(
2460 * i, face_q_point)) *
2461 * fe_face_values.JxW(face_q_point);
2465 * fe_face_values[displacement_multipliers<dim>].value(
2466 * i, face_q_point)) *
2467 * fe_face_values.JxW(face_q_point);
2473 * MatrixTools::local_apply_boundary_values(boundary_values,
2474 * local_dof_indices,
2475 * dummy_cell_matrix,
2479 * constraints.distribute_local_to_global(cell_rhs,
2480 * local_dof_indices,
2491 * <a name="step_79-Computingthemeritfunction"></a>
2492 * <h4>Computing the merit function</h4>
2496 * The algorithm we use herein uses a "watchdog" strategy to
2497 * determine where and how far to go from the current iterate. We
2498 * base the watchdog strategy on an exact @f$l_1@f$ merit function. This
2499 * function calculates the exact @f$l_1@f$ merit of a given, putative,
2504 * The merit function consists of the sum of the objective function
2505 * (which is simply an integral of external forces (on the boundary
2506 * of the domain) times the displacement values of a test solution
2507 * (typically, the current solution plus some multiple of the Newton
2508 * update), and the @f$l_1@f$ norms of the Lagrange multiplier
2509 * components of residual vectors. The following code computes these
2513 * template <int dim>
2514 * double SANDTopOpt<dim>::calculate_exact_merit(
2515 * const BlockVector<double> &test_solution)
2517 * TimerOutput::Scope t(timer, "merit function");
2521 * Start with computing the objective function:
2524 * double objective_function_merit = 0;
2526 * const MappingQ<dim> mapping(1);
2527 * const QGauss<dim> quadrature_formula(fe.degree + 1);
2528 * const QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
2529 * FEValues<dim> fe_values(mapping,
2531 * quadrature_formula,
2532 * update_values | update_gradients |
2533 * update_quadrature_points | update_JxW_values);
2534 * FEFaceValues<dim> fe_face_values(mapping,
2536 * face_quadrature_formula,
2538 * update_quadrature_points |
2539 * update_normal_vectors |
2540 * update_JxW_values);
2542 * const unsigned int n_face_q_points = face_quadrature_formula.size();
2544 * std::vector<Tensor<1, dim>> displacement_face_values(n_face_q_points);
2546 * for (const auto &cell : dof_handler.active_cell_iterators())
2548 * for (const auto &face : cell->face_iterators())
2550 * if (face->at_boundary() &&
2551 * face->boundary_id() == BoundaryIds::down_force)
2553 * fe_face_values.reinit(cell, face);
2554 * fe_face_values[ValueExtractors::displacements<dim>]
2555 * .get_function_values(test_solution,
2556 * displacement_face_values);
2557 * for (unsigned int face_q_point = 0;
2558 * face_q_point < n_face_q_points;
2561 * Tensor<1, dim> traction;
2562 * traction[1] = -1.;
2564 * objective_function_merit +=
2565 * (traction * displacement_face_values[face_q_point]) *
2566 * fe_face_values.JxW(face_q_point);
2573 * for (const auto &cell : triangulation.active_cell_iterators())
2575 * objective_function_merit =
2576 * objective_function_merit -
2577 * barrier_size * cell->measure() *
2578 * std::log(test_solution.block(
2579 * SolutionBlocks::density_lower_slack)[cell->active_cell_index()]);
2580 * objective_function_merit =
2581 * objective_function_merit -
2582 * barrier_size * cell->measure() *
2583 * std::log(test_solution.block(
2584 * SolutionBlocks::density_upper_slack)[cell->active_cell_index()]);
2589 * Then compute the residual and take the @f$l_1@f$ norms of the
2590 * components that correspond to Lagrange multipliers. We add
2591 * those to the objective function computed above, and return
2592 * the sum at the bottom:
2595 * const BlockVector<double> test_rhs = calculate_test_rhs(test_solution);
2597 * const double elasticity_constraint_merit =
2598 * penalty_multiplier *
2599 * test_rhs.block(SolutionBlocks::displacement_multiplier).l1_norm();
2600 * const double filter_constraint_merit =
2601 * penalty_multiplier *
2602 * test_rhs.block(SolutionBlocks::unfiltered_density_multiplier).l1_norm();
2603 * const double lower_slack_merit =
2604 * penalty_multiplier *
2605 * test_rhs.block(SolutionBlocks::density_lower_slack_multiplier).l1_norm();
2606 * const double upper_slack_merit =
2607 * penalty_multiplier *
2608 * test_rhs.block(SolutionBlocks::density_upper_slack_multiplier).l1_norm();
2610 * const double total_merit =
2611 * objective_function_merit + elasticity_constraint_merit +
2612 * filter_constraint_merit + lower_slack_merit + upper_slack_merit;
2613 * return total_merit;
2621 * <a name="step_79-Findingasearchdirection"></a>
2622 * <h4>Finding a search direction</h4>
2626 * Next up is the function that actually computes a search direction
2627 * starting at the current state (passed as the first argument) and
2628 * returns the resulting vector. To this end, the function first
2629 * calls the functions that assemble the linear system that
2630 * corresponds to the Newton system, and that solve it.
2634 * This function also updates the penalty multiplier in the merit
2635 * function, and then returns the largest scaled feasible step.
2636 * It uses the `calculate_max_step_sizes()` function to find the
2637 * largest feasible step that satisfies @f$s>0@f$ and @f$z>0@f$.
2643 * template <int dim>
2644 * BlockVector<double> SANDTopOpt<dim>::find_max_step()
2646 * assemble_system();
2647 * BlockVector<double> step = solve();
2651 * Next we are going to update penalty_multiplier. In
2652 * essence, a larger penalty multiplier makes us consider the
2653 * constraints more. Looking at the Hessian and gradient with
2654 * respect to the step we want to take with our decision
2655 * variables, and comparing that to the norm of our constraint
2656 * error gives us a way to ensure that our merit function is
2657 * "exact" - that is, it has a minimum in the same location
2658 * that the objective function does. As our merit function is
2659 * exact for any penalty multiplier over some minimum value,
2660 * we only keep the computed value if it increases the penalty
2667 * const std::vector<unsigned int> decision_variables = {
2668 * SolutionBlocks::density,
2669 * SolutionBlocks::displacement,
2670 * SolutionBlocks::unfiltered_density,
2671 * SolutionBlocks::density_upper_slack,
2672 * SolutionBlocks::density_lower_slack};
2673 * double hess_part = 0;
2674 * double grad_part = 0;
2675 * for (const unsigned int decision_variable_i : decision_variables)
2677 * for (const unsigned int decision_variable_j : decision_variables)
2679 * Vector<double> temp_vector(step.block(decision_variable_i).size());
2680 * system_matrix.block(decision_variable_i, decision_variable_j)
2681 * .vmult(temp_vector, step.block(decision_variable_j));
2682 * hess_part += step.block(decision_variable_i) * temp_vector;
2684 * grad_part -= system_rhs.block(decision_variable_i) *
2685 * step.block(decision_variable_i);
2688 * const std::vector<unsigned int> equality_constraint_multipliers = {
2689 * SolutionBlocks::displacement_multiplier,
2690 * SolutionBlocks::unfiltered_density_multiplier,
2691 * SolutionBlocks::density_lower_slack_multiplier,
2692 * SolutionBlocks::density_upper_slack_multiplier};
2693 * double constraint_norm = 0;
2694 * for (const unsigned int multiplier_i : equality_constraint_multipliers)
2695 * constraint_norm += system_rhs.block(multiplier_i).linfty_norm();
2698 * double test_penalty_multiplier;
2699 * if (hess_part > 0)
2700 * test_penalty_multiplier =
2701 * (grad_part + .5 * hess_part) / (.05 * constraint_norm);
2703 * test_penalty_multiplier = (grad_part) / (.05 * constraint_norm);
2705 * penalty_multiplier = std::max(penalty_multiplier, test_penalty_multiplier);
2709 * Based on all of this, we can now compute step sizes for the
2710 * primal and dual (Lagrange multiplier) variables. Once we
2711 * have these, we scale the components of the solution vector,
2712 * and that is what this function returns.
2715 * const std::pair<double, double> max_step_sizes =
2716 * calculate_max_step_size(nonlinear_solution, step);
2717 * const double step_size_s = max_step_sizes.first;
2718 * const double step_size_z = max_step_sizes.second;
2720 * step.block(SolutionBlocks::density) *= step_size_s;
2721 * step.block(SolutionBlocks::displacement) *= step_size_s;
2722 * step.block(SolutionBlocks::unfiltered_density) *= step_size_s;
2723 * step.block(SolutionBlocks::displacement_multiplier) *= step_size_z;
2724 * step.block(SolutionBlocks::unfiltered_density_multiplier) *= step_size_z;
2725 * step.block(SolutionBlocks::density_lower_slack) *= step_size_s;
2726 * step.block(SolutionBlocks::density_lower_slack_multiplier) *= step_size_z;
2727 * step.block(SolutionBlocks::density_upper_slack) *= step_size_s;
2728 * step.block(SolutionBlocks::density_upper_slack_multiplier) *= step_size_z;
2738 * <a name="step_79-Computingascaledstep"></a>
2739 * <h4>Computing a scaled step</h4>
2743 * The next function then implements a back-tracking algorithm for a
2744 * line search. It keeps shrinking step size until it finds a step
2745 * where the merit is decreased, and then returns the new location
2746 * based on the current state vector, and the direction to go into,
2747 * times the step length.
2750 * template <int dim>
2751 * BlockVector<double>
2752 * SANDTopOpt<dim>::compute_scaled_step(const BlockVector<double> &state,
2753 * const BlockVector<double> &max_step,
2754 * const double descent_requirement)
2756 * const double merit_derivative =
2757 * (calculate_exact_merit(state + 1e-4 * max_step) -
2758 * calculate_exact_merit(state)) /
2760 * double step_size = 1;
2761 * unsigned int max_linesearch_iterations = 10;
2762 * for (unsigned int k = 0; k < max_linesearch_iterations; ++k)
2764 * if (calculate_exact_merit(state + step_size * max_step) <
2765 * calculate_exact_merit(state) +
2766 * step_size * descent_requirement * merit_derivative)
2769 * step_size = step_size / 2;
2771 * return state + (step_size * max_step);
2778 * <a name="step_79-Checkingforconvergence"></a>
2779 * <h4>Checking for convergence</h4>
2783 * The final auxiliary function in this block is the one that checks
2784 * to see if the KKT conditions are sufficiently met so that the
2785 * overall algorithm can lower the barrier size. It does so by
2786 * computing the @f$l_1@f$ norm of the residual, which is what
2787 * `calculate_test_rhs()` computes.
2790 * template <int dim>
2791 * bool SANDTopOpt<dim>::check_convergence(const BlockVector<double> &state)
2793 * const BlockVector<double> test_rhs = calculate_test_rhs(state);
2794 * const double test_rhs_norm = test_rhs.l1_norm();
2796 * const double convergence_condition = 1e-2;
2797 * const double target_norm = convergence_condition * barrier_size;
2799 * std::cout << " Checking convergence. Current rhs norm is "
2800 * << test_rhs_norm << ", target is " << target_norm << std::endl;
2802 * return (test_rhs_norm < target_norm);
2809 * <a name="step_79-Postprocessingthesolution"></a>
2810 * <h3>Postprocessing the solution</h3>
2814 * The first of the postprocessing functions outputs information
2815 * in a VTU file for visualization. It looks long, but it's
really
2824 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
2825 *
data_component_interpretation(
2827 *
for (
unsigned int i = 0; i < dim; ++i)
2830 *
data_component_interpretation.push_back(
2834 *
data_component_interpretation.push_back(
2836 *
for (
unsigned int i = 0; i < dim; ++i)
2839 *
data_component_interpretation.push_back(
2843 *
data_component_interpretation.push_back(
2846 *
data_component_interpretation.push_back(
2849 *
data_component_interpretation.push_back(
2852 *
data_component_interpretation.push_back(
2855 *
data_component_interpretation.push_back(
2859 *
data_out.attach_dof_handler(dof_handler);
2863 *
data_component_interpretation);
2864 *
data_out.build_patches();
2866 *
std::ofstream output(
"solution" + std::to_string(
iteration) +
".vtu");
2867 *
data_out.write_vtu(output);
2886 *
template <
int dim>
2889 *
static_assert(dim == 2,
2890 *
"This function is not implemented for anything "
2891 *
"other than the 2d case.");
2896 *
stlfile <<
"solid bridge\n" << std::scientific;
2897 *
double height = .25;
2899 *
for (
const auto &cell : dof_handler.active_cell_iterators())
2928 *
stlfile <<
" facet normal " << 0.000000e+00 <<
' '
2929 *
<< 0.000000e+00 <<
' ' << -1.000000e+00 <<
'\n';
2931 *
stlfile <<
" vertex " << cell->vertex(0)[0] <<
' '
2932 *
<< cell->vertex(0)[1] <<
' ' << 0.000000e+00 <<
'\n';
2933 *
stlfile <<
" vertex " << cell->vertex(2)[0] <<
' '
2934 *
<< cell->vertex(2)[1] <<
' ' << 0.000000e+00 <<
'\n';
2935 *
stlfile <<
" vertex " << cell->vertex(1)[0] <<
' '
2936 *
<< cell->vertex(1)[1] <<
' ' << 0.000000e+00 <<
'\n';
2939 *
stlfile <<
" facet normal " << 0.000000e+00 <<
' '
2940 *
<< 0.000000e+00 <<
' ' << -1.000000e+00 <<
'\n';
2942 *
stlfile <<
" vertex " << cell->vertex(1)[0] <<
' '
2943 *
<< cell->vertex(1)[1] <<
' ' << 0.000000e+00 <<
'\n';
2944 *
stlfile <<
" vertex " << cell->vertex(2)[0] <<
' '
2945 *
<< cell->vertex(2)[1] <<
' ' << 0.000000e+00 <<
'\n';
2946 *
stlfile <<
" vertex " << cell->vertex(3)[0] <<
' '
2947 *
<< cell->vertex(3)[1] <<
' ' << 0.000000e+00 <<
'\n';
2952 *
stlfile <<
" facet normal " << 0.000000e+00 <<
' '
2953 *
<< 0.000000e+00 <<
' ' << 1.000000e+00 <<
'\n';
2955 *
stlfile <<
" vertex " << cell->vertex(0)[0] <<
' '
2956 *
<< cell->vertex(0)[1] <<
' ' << height <<
'\n';
2957 *
stlfile <<
" vertex " << cell->vertex(1)[0] <<
' '
2958 *
<< cell->vertex(1)[1] <<
' ' << height <<
'\n';
2959 *
stlfile <<
" vertex " << cell->vertex(2)[0] <<
' '
2960 *
<< cell->vertex(2)[1] <<
' ' << height <<
'\n';
2963 *
stlfile <<
" facet normal " << 0.000000e+00 <<
' '
2964 *
<< 0.000000e+00 <<
' ' << 1.000000e+00 <<
'\n';
2966 *
stlfile <<
" vertex " << cell->vertex(1)[0] <<
' '
2967 *
<< cell->vertex(1)[1] <<
' ' << height <<
'\n';
2968 *
stlfile <<
" vertex " << cell->vertex(3)[0] <<
' '
2969 *
<< cell->vertex(3)[1] <<
' ' << height <<
'\n';
2970 *
stlfile <<
" vertex " << cell->vertex(2)[0] <<
' '
2971 *
<< cell->vertex(2)[1] <<
' ' << height <<
'\n';
2978 *
stlfile <<
" facet normal " << 0.000000e+00 <<
' '
2979 *
<< 0.000000e+00 <<
' ' << -1.000000e+00 <<
'\n';
2981 *
stlfile <<
" vertex " << cell->vertex(0)[0] <<
' '
2982 *
<< cell->vertex(0)[1] <<
' ' << 0.000000e+00 <<
'\n';
2983 *
stlfile <<
" vertex " << cell->vertex(1)[0] <<
' '
2984 *
<< cell->vertex(1)[1] <<
' ' << 0.000000e+00 <<
'\n';
2985 *
stlfile <<
" vertex " << cell->vertex(2)[0] <<
' '
2986 *
<< cell->vertex(2)[1] <<
' ' << 0.000000e+00 <<
'\n';
2989 *
stlfile <<
" facet normal " << 0.000000e+00 <<
' '
2990 *
<< 0.000000e+00 <<
' ' << -1.000000e+00 <<
'\n';
2992 *
stlfile <<
" vertex " << cell->vertex(1)[0] <<
' '
2993 *
<< cell->vertex(1)[1] <<
' ' << 0.000000e+00 <<
'\n';
2994 *
stlfile <<
" vertex " << cell->vertex(3)[0] <<
' '
2995 *
<< cell->vertex(3)[1] <<
' ' << 0.000000e+00 <<
'\n';
2996 *
stlfile <<
" vertex " << cell->vertex(2)[0] <<
' '
2997 *
<< cell->vertex(2)[1] <<
' ' << 0.000000e+00 <<
'\n';
3002 *
stlfile <<
" facet normal " << 0.000000e+00 <<
' '
3003 *
<< 0.000000e+00 <<
' ' << 1.000000e+00 <<
'\n';
3005 *
stlfile <<
" vertex " << cell->vertex(0)[0] <<
' '
3006 *
<< cell->vertex(0)[1] <<
' ' << height <<
'\n';
3007 *
stlfile <<
" vertex " << cell->vertex(2)[0] <<
' '
3008 *
<< cell->vertex(2)[1] <<
' ' << height <<
'\n';
3009 *
stlfile <<
" vertex " << cell->vertex(1)[0] <<
' '
3010 *
<< cell->vertex(1)[1] <<
' ' << height <<
'\n';
3013 *
stlfile <<
" facet normal " << 0.000000e+00 <<
' '
3014 *
<< 0.000000e+00 <<
' ' << 1.000000e+00 <<
'\n';
3016 *
stlfile <<
" vertex " << cell->vertex(1)[0] <<
' '
3017 *
<< cell->vertex(1)[1] <<
' ' << height <<
'\n';
3018 *
stlfile <<
" vertex " << cell->vertex(2)[0] <<
' '
3019 *
<< cell->vertex(2)[1] <<
' ' << height <<
'\n';
3020 *
stlfile <<
" vertex " << cell->vertex(3)[0] <<
' '
3021 *
<< cell->vertex(3)[1] <<
' ' << height <<
'\n';
3036 *
for (
unsigned int face_number = 0;
3041 *
cell->face(face_number);
3043 *
if ((face->at_boundary()) ||
3044 *
(!face->at_boundary() &&
3046 *
0)[cell->neighbor(face_number)->active_cell_index()] <
3050 *
(face->center() - cell->center());
3052 *
if ((face->vertex(0)[0] - face->vertex(0)[0]) *
3053 *
(face->vertex(1)[1] - face->vertex(0)[1]) *
3055 *
(face->vertex(0)[1] - face->vertex(0)[1]) * (0 - 0) *
3056 *
normal_vector[0] +
3058 *
(face->vertex(1)[0] - face->vertex(0)[0]) *
3059 *
normal_vector[1] -
3060 *
(face->vertex(0)[0] - face->vertex(0)[0]) * (0 - 0) *
3061 *
normal_vector[1] -
3062 *
(face->vertex(0)[1] - face->vertex(0)[1]) *
3063 *
(face->vertex(1)[0] - face->vertex(0)[0]) *
3064 *
normal_vector[0] -
3066 *
(face->vertex(1)[1] - face->vertex(0)[1]) * 0 >
3072 *
<< 0.000000e+00 <<
'\n';
3074 *
stlfile <<
" vertex " << face->vertex(0)[0]
3075 *
<<
' ' << face->vertex(0)[1] <<
' '
3076 *
<< 0.000000e+00 <<
'\n';
3077 *
stlfile <<
" vertex " << face->vertex(0)[0]
3078 *
<<
' ' << face->vertex(0)[1] <<
' ' << height
3080 *
stlfile <<
" vertex " << face->vertex(1)[0]
3081 *
<<
' ' << face->vertex(1)[1] <<
' '
3082 *
<< 0.000000e+00 <<
'\n';
3088 *
<< 0.000000e+00 <<
'\n';
3090 *
stlfile <<
" vertex " << face->vertex(0)[0]
3091 *
<<
' ' << face->vertex(0)[1] <<
' ' << height
3093 *
stlfile <<
" vertex " << face->vertex(1)[0]
3094 *
<<
' ' << face->vertex(1)[1] <<
' ' << height
3096 *
stlfile <<
" vertex " << face->vertex(1)[0]
3097 *
<<
' ' << face->vertex(1)[1] <<
' '
3098 *
<< 0.000000e+00 <<
'\n';
3107 *
<< 0.000000e+00 <<
'\n';
3109 *
stlfile <<
" vertex " << face->vertex(0)[0]
3110 *
<<
' ' << face->vertex(0)[1] <<
' '
3111 *
<< 0.000000e+00 <<
'\n';
3112 *
stlfile <<
" vertex " << face->vertex(1)[0]
3113 *
<<
' ' << face->vertex(1)[1] <<
' '
3114 *
<< 0.000000e+00 <<
'\n';
3115 *
stlfile <<
" vertex " << face->vertex(0)[0]
3116 *
<<
' ' << face->vertex(0)[1] <<
' ' << height
3123 *
<< 0.000000e+00 <<
'\n';
3125 *
stlfile <<
" vertex " << face->vertex(0)[0]
3126 *
<<
' ' << face->vertex(0)[1] <<
' ' << height
3128 *
stlfile <<
" vertex " << face->vertex(1)[0]
3129 *
<<
' ' << face->vertex(1)[1] <<
' '
3130 *
<< 0.000000e+00 <<
'\n';
3131 *
stlfile <<
" vertex " << face->vertex(1)[0]
3132 *
<<
' ' << face->vertex(1)[1] <<
' ' << height
3141 *
stlfile <<
"endsolid bridge";
3149 * <a name=
"step_79-Therunfunctiondrivingtheoverallalgorithm"></a>
3175 *
std::cout <<
"filter r is: " <<
filter_r << std::endl;
3182 *
dof_handler.distribute_dofs(fe);
3216 *
const unsigned int max_iterations = 10000;
3221 *
<<
" with barrier parameter " <<
barrier_size << std::endl;
3248 *
std::cout <<
" Starting inner step in iteration "
3250 *
<<
" with merit function penalty multiplier "
3257 *
double target_merit = numbers::signaling_nan<double>();
3281 *
std::cout <<
" current watchdog state merit is: "
3288 *
std::cout <<
" found workable step after " <<
k + 1
3289 *
<<
" iterations" << std::endl;
3372 *
std::cout <<
" Taking scaled step from end of watchdog"
3379 *
<<
" Taking scaled step from beginning of watchdog"
3428 *
std::cout << std::endl;
3441 * <a name=
"step_79-Themainfunction"></a>
3456 *
catch (std::exception &exc)
3458 *
std::cerr << std::endl
3460 *
<<
"----------------------------------------------------"
3462 *
std::cerr <<
"Exception on processing: " << std::endl
3463 *
<< exc.what() << std::endl
3464 *
<<
"Aborting!" << std::endl
3465 *
<<
"----------------------------------------------------"
3472 *
std::cerr << std::endl
3474 *
<<
"----------------------------------------------------"
3476 *
std::cerr <<
"Unknown exception!" << std::endl
3477 *
<<
"Aborting!" << std::endl
3478 *
<<
"----------------------------------------------------"
3513 <
img src=
"https://www.dealii.org/images/steps/developer/step-79.mbbgeometry.png"
3514 alt=
"The MBB problem domain and boundary conditions">
3521<
div class=
"onecolumn" style=
"width: 80%; text-align: center;">
3523 <
img src=
"https://www.dealii.org/images/steps/developer/step-79.filtereddensity.png"
3524 alt=
"Filtered Density Solution">
3527 <
img src=
"https://www.dealii.org/images/steps/developer/step-79.unfiltereddensity.png"
3528 alt=
"Unfiltered Density Solution">
3541<a name=
"step_79-Possibilitiesforextensions"></a><
h4>Possibilities
for extensions</
h4>
3568<a name=
"step_79-PlainProg"></a>
numbers::NumberTraits< Number >::real_type norm() const
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
typename ActiveSelector::face_iterator face_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, 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)
@ matrix
Contents is actually a matrix.
@ general
No special properties.
constexpr types::blas_int zero
constexpr 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)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::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)
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)