# 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
"""Finite element analysis (FEA) for 2D topology optimization problems."""
from __future__ import annotations
from typing import Sequence
import scipy
from numpy import arange
from numpy import array
from numpy import atleast_2d
from numpy import kron
from numpy import newaxis
from numpy import ones
from numpy import setdiff1d
from numpy import tile
from numpy import zeros
from gemseo.core.discipline import MDODiscipline
[docs]class FininiteElementAnalysis(MDODiscipline):
"""Finite Element Analysis for 2D topology optimization problems.
Take in input the Young Modulus vector E and computes in output the compliance, i.e.
twice the work of external forces.
"""
def __init__(
self,
nu: float = 0.3,
n_x: int = 100,
n_y: int = 100,
f_node: int | Sequence[int] = 101 * 101 - 1,
f_direction: int | Sequence[int] = 1,
f_amplitude: int | Sequence[int] = -1,
fixed_nodes: int | Sequence[int] | None = None,
fixed_dir: int | Sequence[int] | None = None,
name: str | None = None,
) -> None:
"""
Args:
nu: The material Poisson's ratio.
n_x: The number of elements in the x-direction.
n_y: The number of elements in the y-direction.
f_node: The indices of the nodes where the forces are applied.
f_direction: The force direction for each ``f_node``, either 0 for x or 1 for y.
f_amplitude: The force amplitude for each pair ``(f_node, f_direction)``.
fixed_nodes: The indices of the nodes where the structure is clamped.
If None, a default value is used.
fixed_dir: The clamped direction for each node, encode 0 for x and 1 for y.
If None, a default value is used.
name: The name of the discipline.
If None, use the class name.
"""
super().__init__(name=name)
if fixed_nodes is None:
fixed_nodes = tile(arange(101), 2)
if fixed_dir is None:
fixed_dir = array([0] * 101 + [1] * 101)
self.N_elements = n_x * n_y
self.N_nodes = (n_x + 1) * (n_y + 1)
self.N_DOFs = 2 * self.N_nodes
self.n_x = n_x
self.n_y = n_y
self.nu = nu
self.E = None
self.KE = None
self.iK = None
self.jK = None
self.edofMat = None
self.freedofs = None
self.f_node = f_node
self.f_direction = f_direction
self.f_amplitude = f_amplitude
self.fixednodes = fixed_nodes
self.fixed_dir = fixed_dir
self.prepare_fea()
self.input_grammar.update(["E"])
self.output_grammar.update(["compliance"])
self.default_inputs = {"E": ones(self.N_elements)}
def _run(self) -> None:
em = self.get_inputs_by_name("E")
sk = ((self.KE.flatten()[newaxis]).T * em).flatten(order="F")
k_mat = scipy.sparse.coo_matrix(
(sk, (self.iK, self.jK)), shape=(self.N_DOFs, self.N_DOFs)
).tocsc()
k_mat = k_mat[self.freedofs, :][:, self.freedofs]
u_vec = zeros((self.N_DOFs, 1))
f = zeros((self.N_DOFs, 1))
f[2 * self.f_node + self.f_direction, 0] = self.f_amplitude
u_vec[self.freedofs, 0] = scipy.sparse.linalg.spsolve(
k_mat, f[self.freedofs, 0]
)
# Objective function and sensitivity
ce = ones(self.N_elements)
ce[:] = (
(u_vec[self.edofMat].reshape(self.N_elements, 8) @ self.KE)
* u_vec[self.edofMat].reshape(self.N_elements, 8)
).sum(1)
self.local_data["compliance"] = array([(em * ce).sum()])
self._is_linearized = True
self._init_jacobian(with_zeros=True)
self.jac["compliance"] = {}
self.jac["compliance"]["E"] = atleast_2d(-ce)
[docs] def prepare_fea(self) -> None:
"""Prepare the Finite Element Analysis."""
self.KE = self.compute_elementary_stiffeness_matrix()
# FE: Build the index vectors for the for coo matrix format.
edof_mat = zeros((self.N_elements, 8), dtype=int)
for elx in range(self.n_x):
for ely in range(self.n_y):
el = ely + elx * self.n_y
n1 = (self.n_y + 1) * elx + ely
n2 = (self.n_y + 1) * (elx + 1) + ely
edof_mat[el, :] = array(
[
2 * n1 + 2,
2 * n1 + 3,
2 * n2 + 2,
2 * n2 + 3,
2 * n2,
2 * n2 + 1,
2 * n1,
2 * n1 + 1,
]
)
self.edofMat = edof_mat
# Construct the index pointers for the coo format
self.iK = kron(edof_mat, ones((8, 1))).flatten()
self.jK = kron(edof_mat, ones((1, 8))).flatten()
# Free DOFs
alldofs = array(range(0, 2 * self.N_nodes))
fixeddofs = 2 * self.fixednodes + self.fixed_dir
self.freedofs = setdiff1d(alldofs, fixeddofs)
[docs] def compute_elementary_stiffeness_matrix(
self,
) -> None: # noqa: D205,D212,D415
"""Compute the elementary stiffness matrix of 1x1 quadrilateral elements."""
em = 1.0
k = array(
[
1 / 2 - self.nu / 6,
1 / 8 + self.nu / 8,
-1 / 4 - self.nu / 12,
-1 / 8 + 3 * self.nu / 8,
-1 / 4 + self.nu / 12,
-1 / 8 - self.nu / 8,
self.nu / 6,
1 / 8 - 3 * self.nu / 8,
]
)
return (
em
/ (1 - self.nu**2)
* array(
[
[k[0], k[1], k[2], k[3], k[4], k[5], k[6], k[7]],
[k[1], k[0], k[7], k[6], k[5], k[4], k[3], k[2]],
[k[2], k[7], k[0], k[5], k[6], k[3], k[4], k[1]],
[k[3], k[6], k[5], k[0], k[7], k[2], k[1], k[4]],
[k[4], k[5], k[6], k[7], k[0], k[1], k[2], k[3]],
[k[5], k[4], k[3], k[2], k[1], k[0], k[7], k[6]],
[k[6], k[3], k[4], k[1], k[2], k[7], k[0], k[5]],
[k[7], k[2], k[1], k[4], k[3], k[6], k[5], k[0]],
]
)
)
