diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml index 5305558..3098162 100644 --- a/.github/workflows/run_tests.yml +++ b/.github/workflows/run_tests.yml @@ -2,33 +2,42 @@ name: Tests on: push: - branches: [main, develop] + branches: ['*'] pull_request: - branches: [main, develop] jobs: build: - runs-on: ubuntu-latest + runs-on: ${{ matrix.platform }} strategy: - max-parallel: 6 fail-fast: false matrix: python-version: ['3.10', '3.11'] platform: [ubuntu-latest, macos-latest, windows-latest] steps: - uses: actions/checkout@v4 - with: - fetch-depth: 0 + - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v4 with: python-version: ${{ matrix.python-version }} - cache: pip + + - name: Get pip cache directory + id: pip-cache + run: echo "dir=$(pip cache dir)" >> $GITHUB_OUTPUT + shell: bash + + - name: Cache pip dependencies + uses: actions/cache@v3 + with: + path: ${{ steps.pip-cache.outputs.dir }} + key: ${{ runner.os }}-pip-${{ hashFiles('**/pyproject.toml') }} + restore-keys: ${{ runner.os }}-pip- - name: Install dependencies - run: | - python -m pip install --upgrade pip - python -m pip install .[dev] + run: pip install .[dev] + + - name: Type check with mypy + run: mypy tofea - name: Test with pytest run: pytest --cov tofea diff --git a/examples/compliance_2d.py b/examples/compliance_2d.py index 20ee070..ecd080e 100755 --- a/examples/compliance_2d.py +++ b/examples/compliance_2d.py @@ -3,76 +3,90 @@ import autograd.numpy as np import matplotlib.pyplot as plt import nlopt +import scipy.ndimage from autograd import value_and_grad +from autograd.extend import defvjp, primitive from tofea.fea2d import FEA2D_K -from tofea.topopt_helpers import simp_parametrization -max_its = 100 -volfrac = 0.5 -sigma = 0.5 -shape = (160, 80) -nelx, nely = shape -emin, emax = 1e-6, 1 +gaussian_filter = primitive(scipy.ndimage.gaussian_filter) +defvjp( + gaussian_filter, + lambda ans, x, *args, **kwargs: lambda g: gaussian_filter(g, *args, **kwargs), # noqa: ARG005 +) -dofs = np.arange(2 * (nelx + 1) * (nely + 1)).reshape(nelx + 1, nely + 1, 2) -fixed = np.zeros_like(dofs, dtype=bool) -load = np.zeros_like(dofs) -fixed[0, :, :] = 1 -load[-1, -1, 1] = 1 +def simp_parametrization(shape, sigma, vmin, vmax, penalty=3.0): + def _parametrization(x): + x = np.reshape(x, shape) + x = gaussian_filter(x, sigma) + x = vmin + (vmax - vmin) * x**penalty + return x -fem = FEA2D_K(fixed) -parametrization = simp_parametrization(shape, sigma, emin, emax) -x0 = np.full(shape, volfrac) + return _parametrization -def objective(x): - x = parametrization(x) - xp = np.pad(x, (0, 1), mode="constant", constant_values=emin) - xp = (xp - emin) / (emax - emin) - body = 1e-3 * np.stack([np.zeros_like(xp), xp], axis=-1) - x = fem.compliance(x, load + body) - return np.sum(x) +def main(): + max_its = 50 + volfrac = 0.3 + sigma = 1.0 + shape = (120, 60) + nelx, nely = shape + emin, emax = 1e-6, 1 + dofs = np.arange(2 * (nelx + 1) * (nely + 1)).reshape(nelx + 1, nely + 1, 2) + fixed = np.zeros_like(dofs, dtype=bool) + load = np.zeros_like(dofs) -def volume(x): - x = parametrization(x) - x = (x - emin) / (emax - emin) - return np.mean(x) + fixed[0, :, :] = 1 + load[-1, -1, 1] = 1 + fem = FEA2D_K(fixed) + parametrization = simp_parametrization(shape, sigma, emin, emax) + x0 = np.full(shape, volfrac) -plt.ion() -fig, ax = plt.subplots(1, 1) -im = ax.imshow(parametrization(x0).T, cmap="gray_r", vmin=emin, vmax=emax) -fig.tight_layout() + plt.ion() + fig, ax = plt.subplots(1, 1, tight_layout=True) + im = ax.imshow(parametrization(x0).T, cmap="gray_r", vmin=emin, vmax=emax) + @value_and_grad + def objective(x): + x = parametrization(x) + d = fem.displacement(x, load) + c = fem.compliance(x, d) + return c -def volume_constraint(x, gd): - v, g = value_and_grad(volume)(x) - if gd.size > 0: - gd[:] = g - return v - volfrac + @value_and_grad + def volume(x): + return np.mean(x) + def volume_constraint(x, gd): + v, g = volume(x) + if gd.size > 0: + gd[:] = g.ravel() + return v - volfrac -def nlopt_obj(x, gd): - c, dc = value_and_grad(objective)(x) + def nlopt_obj(x, gd): + v, g = objective(x) - if gd.size > 0: - gd[:] = dc.ravel() + if gd.size > 0: + gd[:] = g.ravel() - im.set_data(parametrization(x).T) - plt.pause(0.01) + im.set_data(parametrization(x).T) + plt.pause(0.01) - return c + return v + opt = nlopt.opt(nlopt.LD_MMA, x0.size) + opt.add_inequality_constraint(volume_constraint) + opt.set_min_objective(nlopt_obj) + opt.set_lower_bounds(0) + opt.set_upper_bounds(1) + opt.set_maxeval(max_its) + opt.optimize(x0.ravel()) -opt = nlopt.opt(nlopt.LD_CCSAQ, x0.size) -opt.add_inequality_constraint(volume_constraint, 1e-3) -opt.set_min_objective(nlopt_obj) -opt.set_lower_bounds(0) -opt.set_upper_bounds(1) -opt.set_maxeval(max_its) -opt.optimize(x0.ravel()) + plt.show(block=True) -plt.show(block=True) + +if __name__ == "__main__": + main() diff --git a/examples/heat_2d.py b/examples/heat_2d.py index b4d779a..4a25f51 100755 --- a/examples/heat_2d.py +++ b/examples/heat_2d.py @@ -3,72 +3,89 @@ import autograd.numpy as np import matplotlib.pyplot as plt import nlopt +import scipy.ndimage from autograd import value_and_grad +from autograd.extend import defvjp, primitive from tofea.fea2d import FEA2D_T -from tofea.topopt_helpers import simp_parametrization -max_its = 100 -volfrac = 0.5 -sigma = 0.5 -shape = (100, 100) -nelx, nely = shape -cmin, cmax = 1.0, 1e6 +gaussian_filter = primitive(scipy.ndimage.gaussian_filter) +defvjp( + gaussian_filter, + lambda ans, x, *args, **kwargs: lambda g: gaussian_filter(g, *args, **kwargs), # noqa: ARG005 +) -fixed = np.zeros((nelx + 1, nely + 1), dtype="?") -load = np.zeros_like(fixed) -fixed[40:60, -1] = 1 -load[:, :-10] = 1 +def simp_parametrization(shape, sigma, vmin, vmax, penalty=3.0): + def _parametrization(x): + x = np.reshape(x, shape) + x = gaussian_filter(x, sigma) + x = vmin + (vmax - vmin) * x**penalty + return x -fem = FEA2D_T(fixed) -parametrization = simp_parametrization(shape, sigma, cmin, cmax) -x0 = np.full(shape, volfrac) + return _parametrization -def objective(x): - x = parametrization(x) - x = fem.temperature(x, load) - return np.sum(x) +def main(): + max_its = 100 + volfrac = 0.5 + sigma = 1.0 + shape = (100, 100) + nelx, nely = shape + cmin, cmax = 1e-4, 1 + fixed = np.zeros((nelx + 1, nely + 1), dtype="?") + load = np.zeros_like(fixed) -def volume(x): - x = parametrization(x) - x = (x - cmin) / (cmax - cmin) - return np.mean(x) + fixed[shape[0] // 2 - 5 : shape[0] // 2 + 5, -1] = 1 + load[:, 0] = 1 + load[(0, -1), :] = 1 + fem = FEA2D_T(fixed) + parametrization = simp_parametrization(shape, sigma, cmin, cmax) + x0 = np.full(shape, volfrac) -plt.ion() -fig, ax = plt.subplots(1, 1) -im = ax.imshow(parametrization(x0).T, cmap="gray_r", vmin=cmin, vmax=cmax) -fig.tight_layout() + plt.ion() + fig, ax = plt.subplots(1, 1, tight_layout=True) + im = ax.imshow(parametrization(x0).T, cmap="gray_r", vmin=cmin, vmax=cmax) + @value_and_grad + def objective(x): + x = parametrization(x) + t = fem.temperature(x, load) + return np.mean(t) -def volume_constraint(x, gd): - v, g = value_and_grad(volume)(x) - if gd.size > 0: - gd[:] = g - return v - volfrac + @value_and_grad + def volume(x): + return np.mean(x) + def volume_constraint(x, gd): + v, g = volume(x) + if gd.size > 0: + gd[:] = g.ravel() + return v - volfrac -def nlopt_obj(x, gd): - c, dc = value_and_grad(objective)(x) + def nlopt_obj(x, gd): + v, g = objective(x) - if gd.size > 0: - gd[:] = dc.ravel() + if gd.size > 0: + gd[:] = g.ravel() - im.set_data(parametrization(x).T) - plt.pause(0.01) + im.set_data(parametrization(x).T) + plt.pause(0.01) - return c + return v + opt = nlopt.opt(nlopt.LD_MMA, x0.size) + opt.add_inequality_constraint(volume_constraint) + opt.set_min_objective(nlopt_obj) + opt.set_lower_bounds(0) + opt.set_upper_bounds(1) + opt.set_maxeval(max_its) + opt.optimize(x0.ravel()) -opt = nlopt.opt(nlopt.LD_CCSAQ, x0.size) -opt.add_inequality_constraint(volume_constraint, 1e-3) -opt.set_min_objective(nlopt_obj) -opt.set_lower_bounds(0) -opt.set_upper_bounds(1) -opt.set_maxeval(max_its) -opt.optimize(x0.ravel()) + plt.show(block=True) -plt.show(block=True) + +if __name__ == "__main__": + main() diff --git a/tests/test_elements.py b/tests/test_elements.py new file mode 100644 index 0000000..b4a77a8 --- /dev/null +++ b/tests/test_elements.py @@ -0,0 +1,26 @@ +import numpy as np +import pytest + +from tofea.elements import Q4Element_K, Q4Element_T + + +class TestQ4ElementK: + @pytest.fixture() + def q4element_k_instance(self): + return Q4Element_K(e=1.0, nu=1 / 3, dx=0.5, dy=0.5) + + def test_element(self, q4element_k_instance): + element = q4element_k_instance.element + assert isinstance(element, np.ndarray) + assert element.shape == (8, 8) + + +class TestQ4ElementT: + @pytest.fixture() + def q4element_t_instance(self): + return Q4Element_T(k=1.0, dx=0.5, dy=0.5) + + def test_element(self, q4element_t_instance): + element = q4element_t_instance.element + assert isinstance(element, np.ndarray) + assert element.shape == (4, 4) diff --git a/tests/test_fea2d.py b/tests/test_fea2d.py new file mode 100644 index 0000000..328074e --- /dev/null +++ b/tests/test_fea2d.py @@ -0,0 +1,56 @@ +import numpy as np +import pytest + +from tofea.fea2d import FEA2D_K, FEA2D_T + + +class TestFEA2DK: + @pytest.fixture() + def fea2d_k_instance(self): + fixed = np.zeros((5, 5, 2), dtype=bool) + return FEA2D_K(fixed) + + def test_shape(self, fea2d_k_instance): + assert fea2d_k_instance.shape == (4, 4) + + def test_dofs(self, fea2d_k_instance): + dofs = fea2d_k_instance.dofs + assert isinstance(dofs, np.ndarray) + assert dofs.shape == (5, 5, 2) + assert np.all(dofs == np.arange(50).reshape(5, 5, 2)) + + def test_fixdofs(self, fea2d_k_instance): + fixdofs = fea2d_k_instance.fixdofs + assert isinstance(fixdofs, np.ndarray) + assert fixdofs.size == 0 + + def test_freedofs(self, fea2d_k_instance): + freedofs = fea2d_k_instance.freedofs + assert isinstance(freedofs, np.ndarray) + assert np.all(freedofs == np.arange(50)) + + +class TestFEA2DT: + @pytest.fixture() + def fea2d_t_instance(self): + fixed = np.zeros((5, 5), dtype=bool) + return FEA2D_T(fixed) + + def test_shape(self, fea2d_t_instance): + assert fea2d_t_instance.shape == (4, 4) + + def test_dofs(self, fea2d_t_instance): + dofs = fea2d_t_instance.dofs + assert isinstance(dofs, np.ndarray) + assert dofs.shape == (5, 5) + assert np.all(dofs == np.arange(25).reshape(5, 5)) + + def test_fixdofs(self, fea2d_t_instance): + fixdofs = fea2d_t_instance.fixdofs + assert isinstance(fixdofs, np.ndarray) + assert fixdofs.size == 0 + + def test_freedofs(self, fea2d_t_instance): + freedofs = fea2d_t_instance.freedofs + assert isinstance(freedofs, np.ndarray) + assert np.all(freedofs == np.arange(25)) diff --git a/tofea/elements.py b/tofea/elements.py index 9eb397b..e04d3cc 100644 --- a/tofea/elements.py +++ b/tofea/elements.py @@ -1,13 +1,12 @@ from collections.abc import Iterable from dataclasses import dataclass from functools import cached_property +from typing import Any import numpy as np import sympy from numpy.typing import NDArray -__all__ = ["Q4Element_K", "Q4Element_T"] - @dataclass(frozen=True, slots=True, kw_only=True) class Element: @@ -22,7 +21,7 @@ def _b_entries( rule: Iterable[int], shape_funcs: Iterable[sympy.Expr], clist: Iterable[sympy.Symbol], - ) -> tuple[sympy.Expr]: + ) -> tuple[Any, ...]: shape_list = np.concatenate([x * np.asarray(rule) for x in shape_funcs]) return tuple(map(sympy.diff, shape_list, clist)) diff --git a/tofea/fea2d.py b/tofea/fea2d.py index a296df6..841ecf9 100644 --- a/tofea/fea2d.py +++ b/tofea/fea2d.py @@ -1,3 +1,5 @@ +from abc import ABC, abstractmethod +from dataclasses import dataclass from functools import cached_property import autograd.numpy as anp @@ -7,17 +9,51 @@ from tofea import DEFAULT_SOLVER from tofea.elements import Q4Element_K, Q4Element_T from tofea.primitives import solve_coo -from tofea.solvers import get_solver +from tofea.solvers import Solver, get_solver -class FEA2D: - def __init__(self, fixed: NDArray[np.bool_], solver: str = DEFAULT_SOLVER) -> None: - nx, ny = fixed.shape[:2] - dofs = np.arange(fixed.size, dtype=np.uint32).reshape(fixed.shape) - self.out_shape = (nx - 1, ny - 1) - self.fixdofs = dofs[fixed].ravel() - self.freedofs = dofs[~fixed].ravel() - self.solver = get_solver(solver) +@dataclass +class FEA2D(ABC): + fixed: NDArray[np.bool_] + solver: str = DEFAULT_SOLVER + dx: float = 0.5 + dy: float = 0.5 + + @property + @abstractmethod + def dof_dim(self) -> int: + ... + + @property + @abstractmethod + def element(self) -> NDArray: + ... + + @property + @abstractmethod + def dofmap(self) -> NDArray[np.uint32]: + ... + + @property + def shape(self) -> tuple[int, int]: + nx, ny = self.fixed.shape[:2] + return nx - 1, ny - 1 + + @cached_property + def dofs(self) -> NDArray[np.uint32]: + return np.arange(self.fixed.size, dtype=np.uint32).reshape(self.fixed.shape) + + @cached_property + def fixdofs(self) -> NDArray[np.uint32]: + return self.dofs[self.fixed].ravel() + + @cached_property + def freedofs(self) -> NDArray[np.uint32]: + return self.dofs[~self.fixed].ravel() + + @cached_property + def _solver(self) -> Solver: + return get_solver(self.solver) @cached_property def index_map(self) -> NDArray[np.uint32]: @@ -28,7 +64,7 @@ def index_map(self) -> NDArray[np.uint32]: @cached_property def e2sdofmap(self) -> NDArray[np.uint32]: - nx, ny = self.out_shape + nx, ny = self.shape x, y = np.unravel_index(np.arange(nx * ny), (nx, ny)) idxs = self.dof_dim * (y + x * (ny + 1)) return np.add(self.dofmap[None], idxs[:, None].astype(np.uint32)) @@ -52,21 +88,24 @@ def global_mat(self, x: NDArray) -> tuple[NDArray, NDArray]: def solve(self, x: NDArray, b: NDArray) -> NDArray: data, indices = self.global_mat(x) - u_nz = solve_coo(data, indices, b.ravel()[self.freedofs], self.solver) + u_nz = solve_coo(data, indices, b.ravel()[self.freedofs], self._solver) u = anp.concatenate([u_nz, np.zeros(len(self.fixdofs))])[self.index_map] return u +@dataclass class FEA2D_K(FEA2D): dof_dim: int = 2 + e: float = 1.0 + nu: float = 1 / 3 @cached_property def element(self) -> NDArray: - return Q4Element_K().element + return Q4Element_K(e=self.e, nu=self.nu, dx=self.dx, dy=self.dy).element @cached_property def dofmap(self) -> NDArray[np.uint32]: - _, nely = self.out_shape + _, nely = self.shape b = np.arange(2 * (nely + 1), 2 * (nely + 1) + 2) a = b + 2 return np.r_[2, 3, a, b, 0, 1].astype(np.uint32) @@ -74,25 +113,26 @@ def dofmap(self) -> NDArray[np.uint32]: def displacement(self, x: NDArray, b: NDArray) -> NDArray: return self.solve(x, b) - def compliance(self, x: NDArray, b: NDArray) -> NDArray: - displacement = self.displacement(x, b) - dofmap = np.reshape(self.e2sdofmap.T, (-1, *self.out_shape)) + def compliance(self, x: NDArray, displacement: NDArray) -> NDArray: + dofmap = np.reshape(self.e2sdofmap.T, (-1, *self.shape)) c = anp.einsum( "ixy,ij,jxy->xy", displacement[dofmap], self.element, displacement[dofmap] ) - return c + return anp.sum(x * c) +@dataclass class FEA2D_T(FEA2D): dof_dim: int = 1 + k: float = 1.0 @cached_property def element(self) -> NDArray: - return Q4Element_T().element + return Q4Element_T(k=self.k, dx=self.dx, dy=self.dy).element @cached_property def dofmap(self) -> NDArray[np.uint32]: - _, nely = self.out_shape + _, nely = self.shape return np.r_[1, (nely + 2), (nely + 1), 0].astype(np.uint32) def temperature(self, x: NDArray, b: NDArray) -> NDArray: diff --git a/tofea/topopt_helpers.py b/tofea/topopt_helpers.py deleted file mode 100644 index ccf5306..0000000 --- a/tofea/topopt_helpers.py +++ /dev/null @@ -1,44 +0,0 @@ -import autograd.numpy as np -import scipy.ndimage -from autograd.extend import defvjp, primitive - -gaussian_filter = primitive(scipy.ndimage.gaussian_filter) -defvjp( - gaussian_filter, - lambda ans, x, *args, **kwargs: lambda g: gaussian_filter(g, *args, **kwargs), # noqa: ARG005 -) - - -def sigmoid_projection(x, a, b=0.5): - num = np.tanh(a * b) + np.tanh(a * (x - b)) - denom = np.tanh(a * b) + np.tanh(a * (1 - b)) - return num / denom - - -def sigmoid_parametrization(shape, sigma, vmin, vmax, alpha=20, beta=0.5, flat=False): - def _parametrization(x): - x = np.reshape(x, shape) - x = gaussian_filter(x, sigma) - x = sigmoid_projection(x, alpha, beta) - x = vmin + (vmax - vmin) * x - return x.ravel() if flat else x - - return _parametrization - - -def simp_projection(x, vmin, vmax, penalty=3.0): - return vmin + x**penalty * (vmax - vmin) - - -def simp_parametrization(shape, sigma, vmin, vmax, penalty=3.0, flat=False): - def _parametrization(x): - x = np.reshape(x, shape) - x = gaussian_filter(x, sigma) - x = simp_projection(x, vmin, vmax, penalty) - return x.ravel() if flat else x - - return _parametrization - - -def gray_indicator(x): - return np.mean(4 * x * (1 - x))