# Source code for gemseo.problems.topo_opt.topopt_initialize

# 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
"""Generation of the design space and disciplines of the topology optimization
problems."""
from __future__ import annotations

from numpy import arange
from numpy import concatenate
from numpy import dot
from numpy import fix
from numpy import full
from numpy import kron
from numpy import logical_and
from numpy import ones
from numpy import where
from numpy import zeros

from gemseo.algos.design_space import DesignSpace
from gemseo.core.discipline import MDODiscipline
from gemseo.problems.topo_opt.density_filter_disc import DensityFilter
from gemseo.problems.topo_opt.fea_disc import FininiteElementAnalysis
from gemseo.problems.topo_opt.material_model_interpolation_disc import (
MaterialModelInterpolation,
)
from gemseo.problems.topo_opt.volume_fraction_disc import VolumeFraction

[docs]def initialize_design_space_and_discipline_to(
problem: str,
n_x: int,
n_y: int,
e0: float,
nu: float,
penalty: float,
min_member_size: float,
vf0: float,
) -> tuple[DesignSpace, list[MDODiscipline]]:
"""Initialize design space and disciplines for 2D topology optimization problems.

Args:
problem: The problem name, one of "MBB", "L-Shape", "Short_Cantilever".
n_x: The number of elements in the x-direction.
n_y: The number of elements in the y-direction.
e0: The full material Young's modulus.
nu: The material Poisson's ratio.
penalty: The SIMP penalty coefficient.
min_member_size: The minimum structural member size.
vf0: The minimum structural element dimension
imposed in the topology optimization solution.

Returns:

- The design space.
- The disciplines.
"""
# Define the nodal coordinates
x = 0
yy = zeros((n_x + 1) * (n_y + 1))
xx = zeros((n_x + 1) * (n_y + 1))
for i in range(1, n_x + 2):
for j in range(1, n_y + 2):
yy[x] = j
xx[x] = i
x += 1
yy = n_y + 1 - yy
xx -= 1
# Compute the centroid coordinates
xc = zeros(n_x * n_y)
yc = zeros(n_x * n_y)
for xi in range(n_x):
xc[xi * n_y : (xi + 1) * n_y] = xi + 0.5
for yi in range(n_y):
yc[arange(yi, n_x * n_y, n_y)] = yi + 0.5
yc = n_y - yc
if "MBB" == problem:
excitation_node = 0  # Node where the force is applied
excitation_direction = 1  # 0 for x and 1 for y
amplitude = -1  # Amplitude of the force
fixednodes = concatenate(
([[where(xx == min(xx))], [(n_x + 1) * (n_y + 1) - 1]]), axis=None
)  # Fixed nodes
fixed_dir = concatenate(([[ones(n_y + 1)], ]), axis=None) - 1
emptyelts = []  # Mandatory empty elements
fullelts = []  # Mandatory full elements
elif "Short_Cantilever" == problem:
excitation_node = where(
logical_and(
(xx == max(xx)), (yy == fix(dot(0.5, min(yy)) + dot(0.5, max(yy))))
)
)[
0
]  # Nodes where the force is applied
excitation_direction = 1  # 0 for x and 1 for y
amplitude = -1  # Amplitude of the force
fixednodes = kron([1, 1], where(xx == min(xx)))  # Fixed nodes
fixed_dir = (
concatenate([[ones(n_y + 1)], [dot(2, ones(n_y + 1))]]).flatten() - 1
)
emptyelts = []  # Mandatory empty elements
fullelts = []  # Mandatory full elements
elif "L-Shape" == problem:
excitation_node = where(
logical_and(
(xx == max(xx)), (yy == fix(dot(0.5, min(yy)) + dot(0.5, max(yy))))
)
)[
0
]  # Nodes where the force is applied
excitation_direction = 1  # 0 for x and 1 for y
amplitude = -1  # Amplitude of the force
fixednodes = kron([1, 1], where(yy == max(yy)))  # Fixed nodes
fixed_dir = concatenate([[ones(n_x + 1)], [dot(2, ones(n_x + 1))]]).flatten()
emptyelts = where(
logical_and(xc >= (max(xx) + min(xx)) / 2, yc >= ((max(yy) + min(yy)) / 2))
)[
0
]  # Mandatory empty elements
fullelts = []  # Mandatory full element
else:
msg = "The examples covered by this function are MBB, L-Shape and Short_Cantilever."
raise NotImplementedError(msg)
initial_point = full((n_x * n_y,), vf0)
initial_point[emptyelts] = 0
initial_point[fullelts] = 1
ds = DesignSpace()
"x",
size=n_x * n_y,
l_b=zeros((n_x * n_y,)),
u_b=ones((n_x * n_y,)),
value=initial_point,
)
df = DensityFilter(n_x=n_x, n_y=n_y, min_member_size=min_member_size)
mmi = MaterialModelInterpolation(
e0=e0,
penalty=penalty,
n_x=n_x,
n_y=n_y,
empty_elements=emptyelts,
full_elements=fullelts,
)
fea = FininiteElementAnalysis(
nu=nu,
n_x=n_x,
n_y=n_y,
f_node=excitation_node,
f_direction=excitation_direction,
f_amplitude=amplitude,
fixed_nodes=fixednodes,
fixed_dir=fixed_dir,
)
vf = VolumeFraction(
n_x=n_x, n_y=n_y, empty_elements=emptyelts, full_elements=fullelts
)
disciplines = [df, mmi, fea, vf]
return ds, disciplines