# 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
#
# 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))}