Source code for gemseo.problems.topo_opt.material_model_interpolation_disc

# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License version 3 as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
#    INITIAL AUTHORS - API and implementation and/or documentation
#        :author: Simone Coniglio
#    OTHER AUTHORS   - MACROSCOPIC CHANGES
"""A discipline for topology optimization material model interpolation."""
from __future__ import annotations

from typing import Sequence

from numpy import atleast_2d
from numpy import diag
from numpy import ones
from numpy import ones_like

from gemseo.core.discipline import MDODiscipline


[docs]class MaterialModelInterpolation(MDODiscipline): """Material Model Interpolation class for topology optimization problems. Compute the Young's modulus (E) and the material local density (rho) from filtered design variables xPhys with the SIMP (Solid Isotropic Material with Penalization) exponential method. """ def __init__( self, e0: float, penalty: float, n_x: int, n_y: int, empty_elements: Sequence[int], full_elements: Sequence[int], contrast: float = 1e9, ) -> None: # noqa: D205,D212,D415 """ Args: e0: The full material Young modulus. penalty: The SIMP penalty coefficient. n_x: The number of elements in the x-direction. n_y: The number of elements in the y-direction. empty_elements: The index of an empty element ids that are not part of the design space. full_elements: The index of full element ids that are not part of the design space. contrast: The ratio between the full material Young's modulus and void material Young's modulus. """ super().__init__() self.E0 = e0 self.penalty = penalty self.Emin = e0 / contrast self.empty_elements = empty_elements self.full_elements = full_elements self.N_elements = n_x * n_y self.input_grammar.update(["xPhys"]) self.output_grammar.update(["rho", "E"]) self.default_inputs = {"xPhys": ones(n_x * n_y)} def _run(self) -> None: xphys = self.get_inputs_by_name("xPhys") xphys[self.empty_elements] = 0 xphys[self.full_elements] = 1 xphys = xphys.flatten() rho = xphys[:] young_modulus = self.Emin + (self.E0 - self.Emin) * xphys**self.penalty self.local_data["E"] = young_modulus self.local_data["rho"] = rho self._is_linearized = True self._init_jacobian(with_zeros=True) dyoung_modulus_dxphys = ( self.penalty * xphys.ravel() ** (self.penalty - 1) * (self.E0 - self.Emin) ) dyoung_modulus_dxphys[self.empty_elements] = 0 dyoung_modulus_dxphys[self.full_elements] = 0 self.jac["E"] = {"xPhys": atleast_2d(diag(dyoung_modulus_dxphys))} drho_dxphys = ones_like(xphys) drho_dxphys[self.empty_elements] = 0 drho_dxphys[self.full_elements] = 0 self.jac["rho"] = {"xPhys": atleast_2d(diag(drho_dxphys))}