diff --git a/docs/problems/thermoelastic2d.md b/docs/problems/thermoelastic2d.md index c8aa303c..df3d0b68 100644 --- a/docs/problems/thermoelastic2d.md +++ b/docs/problems/thermoelastic2d.md @@ -19,6 +19,7 @@ The design space is then defined by a 2D array representing density values (para The objective of this problem is to minimize total compliance C under a volume fraction constraint V by placing a thermally conductive material. Total compliance is defined as a linear combination of thermal compliance and structural compliance. The weight configuration parameter defines the relative importance of thermal compliance and structural compliance in the total compliance calculation, where a weight of 1 corresponds to a purely structural problem and a weight of 0 corresponds to a purely thermal problem. +The target volume fraction `volfrac` is enforced as a design constraint rather than reported as an objective, following the Beams2D convention. ## Conditions @@ -32,7 +33,7 @@ The optimization process itself operates by calculating the sensitivities of the The optimization loop terminates when either an upper bound of the number of iterations has been reached or if the magnitude of the gradient update is below some threshold. ## Dataset -The dataset linked to this problem is on huggingface [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/thermoelastic_2d_v1). +The dataset linked to this problem is on huggingface [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/thermoelastic_2d_v2). This dataset contains a set of 1000 optimized thermoelastic designs in a 64x64 domain, where each design is optimized for a unique set of conditions. Each datapoint's conditions are randomly generated by arbitrarily placing: a single loaded element along the bottom boundary, two fixed elements (fixed in both the x and y direction) along the left and top boundary, and heatsink elements along the right boundary. Furthermore, values for the volume fraction constraint are randomly selected in the range $[0.2, 0.5]$. @@ -45,11 +46,10 @@ Relevant datapoint fields include: - `force_elements_x`: Encodes a binary NxN matrix specifying elements that have a structural load in the x-direction. - `force_elements_y`: Encodes a binary NxN matrix specifying elements that have a structural load in the y-direction. - `heatsink_elements`: Encodes a binary NxN matrix specifying elements that have a heat sink. -- `volume_fraction_error`: The volume fraction error with respect to the target volume fraction - `structural_compliance`: The structural compliance of the optimized design - `thermal_compliance`: The thermal compliance of the optimized design - `nelx`: The number of elements in the x-direction - `nely`: The number of elements in the y-direction -- `volume_fraction_target`: The volume fraction target of the optimized design +- `volfrac`: The volume fraction target of the optimized design - `rmin`: The filter size used in the optimization routine - `weight`: The domain weighting used in the optimization routine diff --git a/docs/problems/thermoelastic3d.md b/docs/problems/thermoelastic3d.md index cfaf7b72..96d2f4b1 100644 --- a/docs/problems/thermoelastic3d.md +++ b/docs/problems/thermoelastic3d.md @@ -20,6 +20,7 @@ The design space is then defined by a 3D array representing density values (para The objective of this problem is to minimize total compliance C under a volume fraction constraint V by placing a thermally conductive material. Total compliance is defined as a linear combination of thermal compliance and structural compliance. The weight configuration parameter defines the relative importance of thermal compliance and structural compliance in the total compliance calculation, where a weight of 1 corresponds to a purely structural problem and a weight of 0 corresponds to a purely thermal problem. +The target volume fraction `volfrac` is enforced as a design constraint rather than reported as an objective, following the Beams2D convention. ## Conditions @@ -33,7 +34,7 @@ The optimization process itself operates by calculating the sensitivities of the The optimization loop terminates when either an upper bound of the number of iterations has been reached or if the magnitude of the gradient update is below some threshold. ## Dataset -The dataset linked to this problem is on huggingface [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/thermoelastic_3d_v0). +The dataset linked to this problem is on huggingface [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/thermoelastic_3d_v1). This dataset contains a set of 100 optimized thermoelastic designs in a 16x16x16 domain, where each design is optimized for a unique set of conditions. Each datapoint's conditions are randomly generated by arbitrarily placing: a single loaded element along the bottom boundary, two fixed elements (fixed in the x, y, and z direction) along the left and top boundary, and heatsink elements along the right boundary. Furthermore, values for the volume fraction constraint are randomly selected in the range $[0.2, 0.5]$. @@ -45,7 +46,6 @@ Relevant datapoint fields include: - `force_elements_y`: Encodes a binary NxNxN matrix specifying elements that have a structural load in the y-direction. - `force_elements_z`: Encodes a binary NxNxN matrix specifying elements that have a structural load in the z-direction. - `heatsink_elements`: Encodes a binary NxNxN matrix specifying elements that have a heat sink. -- `volume_fraction`: The volume fraction value of the optimized design - `structural_compliance`: The structural compliance of the optimized design - `thermal_compliance`: The thermal compliance of the optimized design - `nelx`: The number of elements in the x-direction diff --git a/engibench/problems/beams3d/model/fem_model.py b/engibench/problems/beams3d/model/fem_model.py index 8b77a751..1b2dfdde 100644 --- a/engibench/problems/beams3d/model/fem_model.py +++ b/engibench/problems/beams3d/model/fem_model.py @@ -235,10 +235,8 @@ def run( # noqa: PLR0915 f0val = f0valm if self.eval_only: - vf_error = abs(np.mean(x) - volfrac) return { "structural_compliance": float(f0valm), - "volume_fraction": vf_error, } obj_values = np.array([f0valm]) @@ -323,12 +321,10 @@ def run( # noqa: PLR0915 opti_steps[-1].obj_values_update = np.zeros_like(opti_steps[-1].obj_values) print("3D structural optimization finished.") - vf_error = abs(np.mean(x) - volfrac) return { "design": x, "bcs": bcs, "structural_compliance": float(f0valm), - "volume_fraction": vf_error, "opti_steps": opti_steps, } diff --git a/engibench/problems/beams3d/v0.py b/engibench/problems/beams3d/v0.py index 5850719e..68d8b311 100644 --- a/engibench/problems/beams3d/v0.py +++ b/engibench/problems/beams3d/v0.py @@ -14,6 +14,7 @@ from engibench.constraint import bounded from engibench.constraint import constraint +from engibench.constraint import Criticality from engibench.constraint import greater_than from engibench.constraint import IMPL from engibench.constraint import THEORY @@ -52,6 +53,16 @@ def _fixed_elements(nelx: int, nely: int, nelz: int) -> npt.NDArray[np.int64]: return out +@constraint(categories=THEORY, criticality=Criticality.Warning) +def volume_fraction_bound(design: npt.NDArray, volfrac: float) -> None: + """Constraint for volume fraction of the design.""" + actual_volfrac = design.mean() + tolerance = 0.01 + assert abs(actual_volfrac - volfrac) <= tolerance, ( + f"Volume fraction of the design {actual_volfrac:.4f} does not match target {volfrac:.4f} specified in the conditions. While the optimizer might fix it, this is likely to affect objective values as the initial design is not feasible given the constraints." + ) + + class Beams3D(Problem[npt.NDArray]): """3D structural topology optimization problem.""" @@ -75,6 +86,7 @@ class Conditions: """Fractional y-position of the vertical load on the top face.""" conditions = Conditions() + design_constraints = (volume_fraction_bound,) design_space = spaces.Box(low=0.0, high=1.0, shape=(NELY, NELX, NELZ), dtype=np.float32) dataset_id = f"IDEALLab/beams_3d_{NELX}_v0" container_id = None diff --git a/engibench/problems/thermoelastic2d/model/fea_model.py b/engibench/problems/thermoelastic2d/model/fea_model.py index 8af20b77..9b0e6c27 100644 --- a/engibench/problems/thermoelastic2d/model/fea_model.py +++ b/engibench/problems/thermoelastic2d/model/fea_model.py @@ -161,7 +161,7 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str Args: bcs (dict[str, any]): A dictionary containing boundary conditions and problem parameters. Expected keys include: - - 'volume_fraction_target' (float): Target volume fraction. + - 'volfrac' (float): Target volume fraction. - 'fixed_elements' (np.ndarray): NxN binary array encoding the location of fixed elements. - 'force_elements_x' (np.ndarray): NxN binary array encoding the location of loaded elements in the x direction. - 'force_elements_y' (np.ndarray): NxN binary array encoding the location of loaded elements in the y direction. @@ -174,10 +174,10 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str Dict[str, Any]: A dictionary containing the optimization results. The dictionary includes: - 'design' (np.ndarray): Final design layout. - 'bcs' (Dict[str, Any]): The input boundary conditions. - - 'sc' (float): Structural cost component. - - 'tc' (float): Thermal cost component. - - 'vf' (float): Volume fraction error. - If self.eval_only is True, returns a dictionary with keys 'sc', 'tc', and 'vf' only. + - 'structural_compliance' (float): Structural compliance. + - 'thermal_compliance' (float): Thermal compliance. + - 'opti_steps' (list[OptiStep]): The optimization history. + If self.eval_only is True, returns a dictionary with keys 'structural_compliance' and 'thermal_compliance' only. """ # WEIGHTING w1 = bcs.get("weight", 0.5) @@ -188,7 +188,7 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str nelx = fe_h - 1 nely = fe_w - 1 - volfrac = bcs["volume_fraction_target"] + volfrac = bcs["volfrac"] n = nely * nelx # Total number of elements # OptiSteps records @@ -326,16 +326,13 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str f0val = (f0valm * w1) + (f0valt * w2) if self.eval_only is True: - vf_error = np.abs(np.mean(x) - volfrac) return { "structural_compliance": f0valm, "thermal_compliance": f0valt, - "volume_fraction_error": vf_error, } # OptiStep Information - vf_error = np.abs(np.mean(x) - volfrac) - obj_values = np.array([f0valm, f0valt, vf_error]) + obj_values = np.array([f0valm, f0valt]) x_curr = x.copy() # Design variables before the gradient update (nely, nelx) df0dx = df0dx_mat.reshape(nely * nelx, 1) @@ -414,14 +411,12 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str extra_iter = True print("Optimization finished...") - vf_error = np.abs(np.mean(x) - volfrac) return { "design": x, "bcs": bcs, "structural_compliance": f0valm, "thermal_compliance": f0valt, - "volume_fraction_error": vf_error, "opti_steps": opti_steps, } @@ -440,7 +435,7 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str "fixed_elements": [lci[21], lci[32], lci[43]], "force_elements_y": [bri[31]], "heatsink_elements": [lci[31], lci[32], lci[33]], - "volume_fraction_target": 0.2, + "volfrac": 0.2, "rmin": 1.1, "weight": 1.0, # 1.0 for pure structural, 0.0 for pure thermal } diff --git a/engibench/problems/thermoelastic2d/v0.py b/engibench/problems/thermoelastic2d/v0.py index e79438b3..1f0e1aa7 100644 --- a/engibench/problems/thermoelastic2d/v0.py +++ b/engibench/problems/thermoelastic2d/v0.py @@ -14,6 +14,7 @@ from engibench.constraint import bounded from engibench.constraint import constraint +from engibench.constraint import Criticality from engibench.constraint import IMPL from engibench.constraint import THEORY from engibench.core import ObjectiveDirection @@ -33,6 +34,16 @@ HEATSINK_ELEMENTS = indices_to_binary_matrix([LCI[31], LCI[32], LCI[33]], NELX + 1, NELY + 1) +@constraint(categories=THEORY, criticality=Criticality.Warning) +def volume_fraction_bound(design: npt.NDArray, volfrac: float) -> None: + """Constraint for volume fraction of the design.""" + actual_volfrac = design.mean() + tolerance = 0.01 + assert abs(actual_volfrac - volfrac) <= tolerance, ( + f"Volume fraction of the design {actual_volfrac:.4f} does not match target {volfrac:.4f} specified in the conditions. While the optimizer might fix it, this is likely to affect objective values as the initial design is not feasible given the constraints." + ) + + class ThermoElastic2D(Problem[npt.NDArray]): r"""Truss 2D integer optimization problem. @@ -43,7 +54,6 @@ class ThermoElastic2D(Problem[npt.NDArray]): objectives: tuple[tuple[str, ObjectiveDirection], ...] = ( ("structural_compliance", ObjectiveDirection.MINIMIZE), ("thermal_compliance", ObjectiveDirection.MINIMIZE), - ("volume_fraction_error", ObjectiveDirection.MINIMIZE), ) @dataclass @@ -66,7 +76,7 @@ class Conditions: default_factory=lambda: HEATSINK_ELEMENTS ) """Binary NxN matrix specifying elements that have a heat sink""" - volume_fraction_target: Annotated[float, bounded(lower=0.0, upper=1.0).category(THEORY)] = 0.3 + volfrac: Annotated[float, bounded(lower=0.0, upper=1.0).category(THEORY)] = 0.3 """Target volume fraction for the volume fraction constraint""" rmin: Annotated[ float, bounded(lower=1.0).category(THEORY), bounded(lower=0.0, upper=3.0).warning().category(IMPL) @@ -76,8 +86,9 @@ class Conditions: """Control which objective is optimized for. 1.0 is pure structural optimization, while 0.0 is pure thermal optimization""" conditions = Conditions() + design_constraints = (volume_fraction_bound,) design_space = spaces.Box(low=0.0, high=1.0, shape=(NELX, NELY), dtype=np.float32) - dataset_id = "IDEALLab/thermoelastic_2d_v1" + dataset_id = "IDEALLab/thermoelastic_2d_v2" container_id = None @dataclass @@ -130,9 +141,7 @@ def simulate_verbose(self, design: npt.NDArray, config: dict[str, Any] | None = boundary_dict[key] = value results = FeaModel(plot=False, eval_only=True).run(boundary_dict, x_init=design) - return SimulationResult( - np.array([results["structural_compliance"], results["thermal_compliance"], results["volume_fraction_error"]]) - ) + return SimulationResult(np.array([results["structural_compliance"], results["thermal_compliance"]])) def optimize( self, starting_point: npt.NDArray, config: dict[str, Any] | None = None diff --git a/engibench/problems/thermoelastic3d/model/fem_model.py b/engibench/problems/thermoelastic3d/model/fem_model.py index ca63f2db..db0a563d 100644 --- a/engibench/problems/thermoelastic3d/model/fem_model.py +++ b/engibench/problems/thermoelastic3d/model/fem_model.py @@ -198,10 +198,10 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str Dict[str, Any]: A dictionary containing the optimization results. The dictionary includes: - 'design' (np.ndarray): Final design layout. - 'bcs' (Dict[str, Any]): The input boundary conditions. - - 'sc' (float): Structural cost component. - - 'tc' (float): Thermal cost component. - - 'vf' (float): Volume fraction error. - If self.eval_only is True, returns a dictionary with keys 'sc', 'tc', and 'vf' only. + - 'structural_compliance' (float): Structural compliance. + - 'thermal_compliance' (float): Thermal compliance. + - 'opti_steps' (list[OptiStep]): The optimization history. + If self.eval_only is True, returns a dictionary with keys 'structural_compliance' and 'thermal_compliance' only. """ # Weighting w1 = bcs.get("weight", 0.5) # structural @@ -368,14 +368,11 @@ def node_id(ix: int, iy: int, iz: int) -> int: f0val = (f0valm * w1) + (f0valt * w2) if self.eval_only: - vf_error = abs(np.mean(x) - volfrac) return { "structural_compliance": float(f0valm), "thermal_compliance": float(f0valt), - "volume_fraction": vf_error, } - vf_error = np.abs(np.mean(x) - volfrac) - obj_values = np.array([f0valm, f0valt, vf_error]) + obj_values = np.array([f0valm, f0valt]) x_curr = x.copy() xval = x.reshape(n, 1) @@ -473,13 +470,11 @@ def node_id(ix: int, iy: int, iz: int) -> int: extra_iter = True print("3D optimization finished.") - vf_error = abs(np.mean(x) - volfrac) return { "design": x, "bcs": bcs, "structural_compliance": float(f0valm), "thermal_compliance": float(f0valt), - "volume_fraction": vf_error, "opti_steps": opti_steps, } diff --git a/engibench/problems/thermoelastic3d/v0.py b/engibench/problems/thermoelastic3d/v0.py index bd1fa06b..f74c2822 100644 --- a/engibench/problems/thermoelastic3d/v0.py +++ b/engibench/problems/thermoelastic3d/v0.py @@ -12,6 +12,7 @@ from engibench.constraint import bounded from engibench.constraint import constraint +from engibench.constraint import Criticality from engibench.constraint import IMPL from engibench.constraint import THEORY from engibench.core import ObjectiveDirection @@ -48,6 +49,16 @@ def _default_heatsink_elements(nelx: int, nely: int, nelz: int) -> npt.NDArray[n return out +@constraint(categories=THEORY, criticality=Criticality.Warning) +def volume_fraction_bound(design: npt.NDArray, volfrac: float) -> None: + """Constraint for volume fraction of the design.""" + actual_volfrac = design.mean() + tolerance = 0.01 + assert abs(actual_volfrac - volfrac) <= tolerance, ( + f"Volume fraction of the design {actual_volfrac:.4f} does not match target {volfrac:.4f} specified in the conditions. While the optimizer might fix it, this is likely to affect objective values as the initial design is not feasible given the constraints." + ) + + class ThermoElastic3D(Problem[npt.NDArray]): """Truss 3D integer optimization problem. @@ -58,7 +69,6 @@ class ThermoElastic3D(Problem[npt.NDArray]): objectives: tuple[tuple[str, ObjectiveDirection], ...] = ( ("structural_compliance", ObjectiveDirection.MINIMIZE), ("thermal_compliance", ObjectiveDirection.MINIMIZE), - ("volume_fraction", ObjectiveDirection.MINIMIZE), ) @dataclass @@ -87,8 +97,9 @@ class Conditions: weight: Annotated[float, bounded(lower=0.0, upper=1.0).category(THEORY)] = 0.5 """Control which objective is optimized for. 1.0 is pure structural optimization, while 0.0 is pure thermal optimization""" + design_constraints = (volume_fraction_bound,) design_space = spaces.Box(low=0.0, high=1.0, shape=(NELX, NELY, NELZ), dtype=np.float32) - dataset_id = "IDEALLab/thermoelastic_3d_v0" + dataset_id = "IDEALLab/thermoelastic_3d_v1" container_id = None @dataclass @@ -206,9 +217,7 @@ def simulate_verbose(self, design: npt.NDArray, config: dict[str, Any] | None = boundary_dict[key] = value results = FeaModel3D(eval_only=True).run(boundary_dict, x_init=design) - return SimulationResult( - np.array([results["structural_compliance"], results["thermal_compliance"], results["volume_fraction"]]) - ) + return SimulationResult(np.array([results["structural_compliance"], results["thermal_compliance"]])) def optimize( self, starting_point: npt.NDArray, config: dict[str, Any] | None = None diff --git a/tests/reference/simulate/problems.thermoelastic2d.v0.ThermoElastic2D.json b/tests/reference/simulate/problems.thermoelastic2d.v0.ThermoElastic2D.json index 5fbff0e4..16b1421b 100644 --- a/tests/reference/simulate/problems.thermoelastic2d.v0.ThermoElastic2D.json +++ b/tests/reference/simulate/problems.thermoelastic2d.v0.ThermoElastic2D.json @@ -1,8 +1,7 @@ { "performance": [ 104399698.05000441, - 7603.173244221379, - 0.05000413414090871 + 7603.173244221379 ], "rtol": 5e-05 } diff --git a/tests/reference/simulate/problems.thermoelastic3d.v0.ThermoElastic3D.json b/tests/reference/simulate/problems.thermoelastic3d.v0.ThermoElastic3D.json index f01e1fe6..ab8f24ab 100644 --- a/tests/reference/simulate/problems.thermoelastic3d.v0.ThermoElastic3D.json +++ b/tests/reference/simulate/problems.thermoelastic3d.v0.ThermoElastic3D.json @@ -1,8 +1,7 @@ { "performance": [ 2330.7670793294246, - 4050.8264802757444, - 0.1499882920390519 + 4050.8264802757444 ], "rtol": 1e-07 }