From 5d9a215fdcf71f8fc9a6423499b3f7067a79883c Mon Sep 17 00:00:00 2001 From: Nolan Black Date: Wed, 12 Nov 2025 14:38:17 -0500 Subject: [PATCH 1/7] fix test --- examples/ansys/pymapdl_tess/test.py | 28 +++++++++++++++++++--------- 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/examples/ansys/pymapdl_tess/test.py b/examples/ansys/pymapdl_tess/test.py index 4187cba..0600c3c 100644 --- a/examples/ansys/pymapdl_tess/test.py +++ b/examples/ansys/pymapdl_tess/test.py @@ -1,25 +1,35 @@ +""" +/usr/ansys_inc/v241/ansys/bin/ansys241 -port 50005 -grpc +""" import os import numpy as np -from dotenv import load_dotenv from tesseract_core import Tesseract -load_dotenv() +# load the MAPDL server's address from your environment host = os.getenv("MAPDL_HOST") if host is None: - raise ValueError("Unable to read $MAPDL_HOST from the environment.") + raise ValueError("Unable to read $MAPDL_HOST from the environment. " + \ + "Use 'export MAPDL_HOST=X.X.X.X' for local IP address of your MAPDL Instance.") port = os.getenv("MAPDL_PORT") if port is None: - raise ValueError("Unable to read $MAPDL_PORT from the environment.") + raise ValueError("Unable to read $MAPDL_PORT from the environment. " + \ + "Use 'export MAPDL_PORT=X' for the port of your MAPDL Instance.") -tess_simp_compliance = Tesseract.from_tesseract_api( - "tesseract_pymapdl/tess_simp_compiance/tesseract_api.py" -) +# Initialize Tesseract from api (for testing) +try: + tess_simp_compliance = Tesseract.from_tesseract_api( + "./tesseract_api.py" + ) +except RuntimeError as e: + raise RuntimeError("Unable to load tesseract from api. " \ + "Ensure that you have installed the build requirements using 'pip install -r tesseract_requirements.txt'") + Lx, Ly, Lz = 3, 2, 1 -Nx, Ny, Nz = 60, 40, 20 -# Nx, Ny, Nz = 6, 4, 2 +# Nx, Ny, Nz = 60, 40, 20 +Nx, Ny, Nz = 6, 4, 2 n_elem = Nx * Ny * Nz # Create a test density field varying from 0 to 1 From 4733b12caf339e317eeeb0612d5d12ddfe3e6bad Mon Sep 17 00:00:00 2001 From: Nolan Black Date: Wed, 12 Nov 2025 17:29:24 -0500 Subject: [PATCH 2/7] trying to fix hex meshing --- examples/ansys/pymapdl_tess/tesseract_api.py | 270 +++++++++++++------ examples/ansys/pymapdl_tess/test.py | 141 ++++++---- 2 files changed, 286 insertions(+), 125 deletions(-) diff --git a/examples/ansys/pymapdl_tess/tesseract_api.py b/examples/ansys/pymapdl_tess/tesseract_api.py index 33d04ad..664639b 100644 --- a/examples/ansys/pymapdl_tess/tesseract_api.py +++ b/examples/ansys/pymapdl_tess/tesseract_api.py @@ -11,6 +11,7 @@ import logging import time +import tempfile from collections.abc import Callable from functools import wraps from typing import ParamSpec, TypeVar @@ -18,17 +19,22 @@ import numpy as np import pyvista as pv from ansys.mapdl.core import Mapdl +from ansys.mapdl.reader import save_as_archive from pydantic import BaseModel, Field -from tesseract_core.runtime import Array, Differentiable, Float32 +from tesseract_core.runtime import Array, Differentiable, Float32, Int32 # Set up module logger logger = logging.getLogger(__name__) -# Type variables for decorator +# TODO delete me once debugging is over +Lx, Ly, Lz = 3, 2, 1 +Nx, Ny, Nz = 6, 4, 2 +# Nx, Ny, Nz = 60, 40, 20 +force_val = 0.1 +debug = True + P = ParamSpec("P") T = TypeVar("T") - - def log_timing(func: Callable[P, T]) -> Callable[P, T]: """Decorator to log the wall time of method execution.""" @@ -45,6 +51,19 @@ def wrapper(*args: P.args, **kwargs: P.kwargs) -> T: return wrapper +class HexMesh(BaseModel): + """Hexahedral mesh representation.""" + + points: Array[(None, 3), Float32] = Field(description="Array of vertex positions.") + faces: Array[(None, 8), Int32] = Field( + description="Array of hexahedral faces defined by indices into the points array." + ) + n_points: Int32 = Field( + default=0, description="Number of valid points in the points array." + ) + n_faces: Int32 = Field( + default=0, description="Number of valid faces in the faces array." + ) class InputSchema(BaseModel): """Inputs for the tess_simp_compliance Tesseract. @@ -64,29 +83,20 @@ class InputSchema(BaseModel): Float32, ] ] = Field(description="2D density field for topology optimization") - - Lx: float = Field( - default=60.0, description="Length of the simulation box in the x direction." + van_neumann_mask: Array[(None,), Int32] = Field( + description="Mask for van Neumann boundary conditions", ) - Ly: float = Field( - default=30.0, - description=("Length of the simulation box in the y direction."), + van_neumann_values: Array[(None, None), Float32] = Field( + description="Values for van Neumann boundary conditions", ) - Lz: float = Field( - default=30.0, description="Length of the simulation box in the z direction." + dirichlet_mask: Array[(None,), Int32] = Field( + description="Mask for Dirichlet boundary conditions", ) - - Nx: int = Field( - default=60, - description=("Number of elements in the x direction."), + dirichlet_values: Array[(None,), Float32] = Field( + description="Values for Dirichlet boundary conditions", ) - Ny: int = Field( - default=30, - description=("Number of elements in the y direction."), - ) - Nz: int = Field( - default=30, - description=("Number of elements in the z direction."), + hex_mesh: HexMesh = Field( + description="Hexahedral mesh representation of the geometry", ) E0: float = Field( @@ -99,11 +109,6 @@ class InputSchema(BaseModel): description="Minimum density value to avoid singular stiffness matrix in SIMP.", ) - total_force: float = Field( - default=5.0, - description="Total force magnitude in Newtons applied to the structure.", - ) - p: float = Field( default=3.0, description="SIMP penalty parameter for material interpolation (default: 3.0).", @@ -158,15 +163,14 @@ def __init__(self, inputs: InputSchema, mapdl: Mapdl) -> None: # Extract input parameters self.rho = inputs.rho - self.Lx = inputs.Lx - self.Ly = inputs.Ly - self.Lz = inputs.Lz - self.Nx = inputs.Nx - self.Ny = inputs.Ny - self.Nz = inputs.Nz + self.hex_mesh = inputs.hex_mesh + self.dirichlet_mask = inputs.dirichlet_mask + self.dirichlet_values = inputs.dirichlet_values + self.van_neumann_mask = inputs.van_neumann_mask + self.van_neumann_values = inputs.van_neumann_values + self.hex_mesh = inputs.hex_mesh self.E0 = inputs.E0 self.rho_min = inputs.rho_min - self.total_force = inputs.total_force self.p = inputs.p self.log_level = inputs.log_level self.vtk_output = inputs.vtk_output @@ -205,12 +209,21 @@ def solve(self) -> OutputSchema: """ logger.info("Starting SIMP elasticity analysis...") - self._create_geometry() + if False: + self._create_geometry() + else: + self._create_geometry2() self._define_element() - self._create_mesh() + if False: + self._create_mesh() + else: + self._create_mesh2() self._define_simp_materials() self._assign_materials_to_elements() - self._apply_boundary_conditions() + if debug: + self._apply_boundary_conditions() + else: + self._apply_boundary_conditions_2() self._run_analysis() self._extract_displacement_constraints() self._extract_nodal_displacement() @@ -226,6 +239,26 @@ def solve(self) -> OutputSchema: return self._build_output_schema() + @log_timing + def _create_geometry2(self) -> None: + """Create a rectangular block geometry. + + Creates 8 keypoints at the corners of a rectangular box and defines + a volume from these keypoints. + """ + # start pre-processor + self.mapdl.prep7() + + # create a pyvista mesh from self.hex_mesh + cells = np.concatenate([np.array([[8]] * self.hex_mesh.n_faces), self.hex_mesh.faces], 1) + cell_type = np.array([pv.CellType.HEXAHEDRON] * self.hex_mesh.n_faces) + mesh = pv.UnstructuredGrid(cells, cell_type, self.hex_mesh.points) + + # save the pyvista mesh, then reload it in mapdl + archive_filename = os.path.join(tempfile.gettempdir(), "tmp.cdb") + save_as_archive(archive_filename, mesh) + response = self.mapdl.cdread("db", archive_filename) + @log_timing def _create_geometry(self) -> None: """Create a rectangular block geometry. @@ -237,13 +270,13 @@ def _create_geometry(self) -> None: self.mapdl.prep7() k0 = self.mapdl.k("", 0, 0, 0) - k1 = self.mapdl.k("", self.Lx, 0, 0) - k2 = self.mapdl.k("", self.Lx, self.Ly, 0) - k3 = self.mapdl.k("", 0, self.Ly, 0) - k4 = self.mapdl.k("", 0, 0, self.Lz) - k5 = self.mapdl.k("", self.Lx, 0, self.Lz) - k6 = self.mapdl.k("", self.Lx, self.Ly, self.Lz) - k7 = self.mapdl.k("", 0, self.Ly, self.Lz) + k1 = self.mapdl.k("", Lx, 0, 0) + k2 = self.mapdl.k("", Lx, Ly, 0) + k3 = self.mapdl.k("", 0, Ly, 0) + k4 = self.mapdl.k("", 0, 0, Lz) + k5 = self.mapdl.k("", Lx, 0, Lz) + k6 = self.mapdl.k("", Lx, Ly, Lz) + k7 = self.mapdl.k("", 0, Ly, Lz) self.mapdl.v(k0, k1, k2, k3, k4, k5, k6, k7) @@ -285,10 +318,13 @@ def _define_simp_materials(self) -> None: def _define_element(self) -> None: """Define element type for structural analysis. - Sets element type 1 to SOLID186, a 20-node 3D structural solid element - with quadratic displacement behavior. + Sets element type 1 to SOLID185, a 8-node 3D structural solid element + with linear displacement behavior. + + Note: SOLID186, a 20-node 3D structural solid element + with quadratic displacement behavior, would be an improvement. """ - self.mapdl.et(1, "SOLID186") + self.mapdl.et(1, "SOLID185") @log_timing def _assign_materials_to_elements(self) -> None: @@ -313,6 +349,48 @@ def _assign_materials_to_elements(self) -> None: self.mapdl.input_strings(commands) @log_timing + def _create_mesh2(self) -> None: + """Create structured hexahedral mesh with specified divisions. + + Sets line divisions for each direction and generates a swept mesh using + VSWEEP. The resulting mesh will have Nx * Ny * Nz elements with uniform + element sizes in each direction. + """ + # Store element numbers in the order they were created + # This is crucial for maintaining consistent indexing between + # density input and strain energy output + self.element_numbers = self.mapdl.mesh.enum + + # Validate and log element numbering + n_elements = self.hex_mesh.n_faces + if len(self.element_numbers) != n_elements: + raise ValueError( + f"Number of created elements {len(self.element_numbers)} does not match " + f"expected {n_elements} (Nx={self.Nx}, Ny={self.Ny}, Nz={self.Nz})" + ) + self.n_elements = n_elements + + # Repeat validation for nodes + self.node_numbers = self.mapdl.mesh.nnum + self.n_nodes = len(self.node_numbers) + + logger.debug( + f"Element numbers: min={self.element_numbers.min()}, " + f"max={self.element_numbers.max()}, " + f"first 5={self.element_numbers[:5]}" + ) + + # Check if element numbering is sequential starting from 1 + is_sequential = np.array_equal( + self.element_numbers, np.arange(1, n_elements + 1) + ) + if is_sequential: + logger.info("Element numbering is sequential (1, 2, 3, ...)") + else: + logger.warning( + "Element numbering is NOT sequential - reordering will be applied" + ) + @log_timing def _create_mesh(self) -> None: """Create structured hexahedral mesh with specified divisions. @@ -324,53 +402,53 @@ def _create_mesh(self) -> None: # Lines parallel to X-axis (select by 2 constant coordinates) self.mapdl.lsel("s", "loc", "y", 0) self.mapdl.lsel("r", "loc", "z", 0) - self.mapdl.lesize("all", ndiv=self.Nx, kforc=1) + self.mapdl.lesize("all", ndiv=Nx, kforc=1) - self.mapdl.lsel("s", "loc", "y", self.Ly) + self.mapdl.lsel("s", "loc", "y", Ly) self.mapdl.lsel("r", "loc", "z", 0) - self.mapdl.lesize("all", ndiv=self.Nx, kforc=1) + self.mapdl.lesize("all", ndiv=Nx, kforc=1) self.mapdl.lsel("s", "loc", "y", 0) - self.mapdl.lsel("r", "loc", "z", self.Lz) - self.mapdl.lesize("all", ndiv=self.Nx, kforc=1) + self.mapdl.lsel("r", "loc", "z", Lz) + self.mapdl.lesize("all", ndiv=Nx, kforc=1) - self.mapdl.lsel("s", "loc", "y", self.Ly) - self.mapdl.lsel("r", "loc", "z", self.Lz) - self.mapdl.lesize("all", ndiv=self.Nx, kforc=1) + self.mapdl.lsel("s", "loc", "y", Ly) + self.mapdl.lsel("r", "loc", "z", Lz) + self.mapdl.lesize("all", ndiv=Nx, kforc=1) # Lines parallel to Y-axis (select by 2 constant coordinates) self.mapdl.lsel("s", "loc", "x", 0) self.mapdl.lsel("r", "loc", "z", 0) - self.mapdl.lesize("all", ndiv=self.Ny, kforc=1) + self.mapdl.lesize("all", ndiv=Ny, kforc=1) - self.mapdl.lsel("s", "loc", "x", self.Lx) + self.mapdl.lsel("s", "loc", "x", Lx) self.mapdl.lsel("r", "loc", "z", 0) - self.mapdl.lesize("all", ndiv=self.Ny, kforc=1) + self.mapdl.lesize("all", ndiv=Ny, kforc=1) self.mapdl.lsel("s", "loc", "x", 0) - self.mapdl.lsel("r", "loc", "z", self.Lz) - self.mapdl.lesize("all", ndiv=self.Ny, kforc=1) + self.mapdl.lsel("r", "loc", "z", Lz) + self.mapdl.lesize("all", ndiv=Ny, kforc=1) - self.mapdl.lsel("s", "loc", "x", self.Lx) - self.mapdl.lsel("r", "loc", "z", self.Lz) - self.mapdl.lesize("all", ndiv=self.Ny, kforc=1) + self.mapdl.lsel("s", "loc", "x", Lx) + self.mapdl.lsel("r", "loc", "z", Lz) + self.mapdl.lesize("all", ndiv=Ny, kforc=1) # Lines parallel to Z-axis (select by 2 constant coordinates) self.mapdl.lsel("s", "loc", "x", 0) self.mapdl.lsel("r", "loc", "y", 0) - self.mapdl.lesize("all", ndiv=self.Nz, kforc=1) + self.mapdl.lesize("all", ndiv=Nz, kforc=1) - self.mapdl.lsel("s", "loc", "x", self.Lx) + self.mapdl.lsel("s", "loc", "x", Lx) self.mapdl.lsel("r", "loc", "y", 0) - self.mapdl.lesize("all", ndiv=self.Nz, kforc=1) + self.mapdl.lesize("all", ndiv=Nz, kforc=1) self.mapdl.lsel("s", "loc", "x", 0) - self.mapdl.lsel("r", "loc", "y", self.Ly) - self.mapdl.lesize("all", ndiv=self.Nz, kforc=1) + self.mapdl.lsel("r", "loc", "y", Ly) + self.mapdl.lesize("all", ndiv=Nz, kforc=1) - self.mapdl.lsel("s", "loc", "x", self.Lx) - self.mapdl.lsel("r", "loc", "y", self.Ly) - self.mapdl.lesize("all", ndiv=self.Nz, kforc=1) + self.mapdl.lsel("s", "loc", "x", Lx) + self.mapdl.lsel("r", "loc", "y", Ly) + self.mapdl.lesize("all", ndiv=Nz, kforc=1) # Restore full selection of entities self.mapdl.allsel() @@ -384,7 +462,7 @@ def _create_mesh(self) -> None: self.element_numbers = self.mapdl.mesh.enum # Validate and log element numbering - n_elements = self.Nx * self.Ny * self.Nz + n_elements = self.hex_mesh.n_faces if len(self.element_numbers) != n_elements: raise ValueError( f"Number of created elements {len(self.element_numbers)} does not match " @@ -413,6 +491,38 @@ def _create_mesh(self) -> None: "Element numbering is NOT sequential - reordering will be applied" ) + @log_timing + def _apply_boundary_conditions_2(self) -> None: + """Apply boundary conditions, forces, and visualize. + + Fixes all nodes on the x=0 plane by constraining all translational + degrees of freedom (UX, UY, UZ). Then applies a concentrated force + in the negative z-direction to nodes at x=Lx, y=Ly, and z between + 0.4*Lz and 0.6*Lz. + """ + for i, x in zip(self.dirichlet_mask, self.dirichlet_values): + # Select node + self.mapdl.nsel("s", "NODE", "", i + 1) + + # set dirichlet conditions + self.mapdl.d("all", "all", x) + + # Reselect all nodes + self.mapdl.allsel() + + + for i, f in zip(self.van_neumann_mask, self.van_neumann_values): + # Select node + self.mapdl.nsel("s", "NODE", "", i + 1) + + # Apply force in negative z-direction to all selected nodes + self.mapdl.f("all", "fx", f[0]) + self.mapdl.f("all", "fy", f[1]) + self.mapdl.f("all", "fz", f[2]) + + # Reselect all nodes + self.mapdl.allsel() + @log_timing def _apply_boundary_conditions(self) -> None: """Apply boundary conditions, forces, and visualize. @@ -433,22 +543,22 @@ def _apply_boundary_conditions(self) -> None: # Apply force to nodes at x=Lx, y=0, z between 0.4*Lz and 0.6*Lz # Select nodes at x=Lx - self.mapdl.nsel("s", "loc", "x", self.Lx) + self.mapdl.nsel("s", "loc", "x", Lx) # Refine selection between y=0 and y=0.2 Ly - y_max = 0.2 * self.Ly + y_max = 0.2 * Ly self.mapdl.nsel("r", "loc", "y", 0, y_max) # Refine selection to z between 0.4*Lz and 0.6*Lz - z_min = 0.4 * self.Lz - z_max = 0.6 * self.Lz + z_min = 0.4 * Lz + z_max = 0.6 * Lz self.mapdl.nsel("r", "loc", "z", z_min, z_max) # Get number of selected nodes num_nodes = self.mapdl.mesh.n_node # Calculate force per node (negative for downward direction) - force_per_node = -self.total_force / num_nodes + force_per_node = -force_val / num_nodes # Apply force in negative z-direction to all selected nodes self.mapdl.f("all", "fz", force_per_node) @@ -457,7 +567,7 @@ def _apply_boundary_conditions(self) -> None: self.mapdl.allsel() logger.info( - f"Applied total force of {self.total_force} N distributed over {num_nodes} nodes" + f"Applied total force of {force_val} N distributed over {num_nodes} nodes" ) @log_timing @@ -533,7 +643,7 @@ def _extract_nodal_displacement(self) -> None: nodal_displacement = np.zeros((self.n_nodes, 3)) nnum, disp = self.mapdl.result.nodal_displacement(0) # 0 for first result set # populate nodal_displacement using vectorized NumPy indexing - nodal_displacement[nnum - 1] = disp + nodal_displacement[nnum - 1] = disp[:, 0:3] self.nodal_displacement = nodal_displacement @log_timing diff --git a/examples/ansys/pymapdl_tess/test.py b/examples/ansys/pymapdl_tess/test.py index 0600c3c..944bd9f 100644 --- a/examples/ansys/pymapdl_tess/test.py +++ b/examples/ansys/pymapdl_tess/test.py @@ -27,67 +27,118 @@ raise RuntimeError("Unable to load tesseract from api. " \ "Ensure that you have installed the build requirements using 'pip install -r tesseract_requirements.txt'") +# TODO +# we should probably create some re-usable mesh constructors... +# create a simple Hex mesh using Pyvista +import pyvista as pv Lx, Ly, Lz = 3, 2, 1 # Nx, Ny, Nz = 60, 40, 20 Nx, Ny, Nz = 6, 4, 2 +grid = pv.ImageData( + dimensions=np.array((Nx, Ny, Nz)) + 1, + #origin=(-Lx / 2, -Ly / 2, -Lz / 2), # The bottom left corner of the data set + origin=(0, 0, 0), # TODO + spacing=(Lx / Nx, Ly / Ny, Lz / Nz), # These are the cell sizes along each axis +) +# repeated casts will eventually expose cell_connectivitiy +mesh = grid.cast_to_structured_grid().cast_to_explicit_structured_grid().cast_to_unstructured_grid() + +from tesseract_api import HexMesh +hex_mesh = HexMesh( + points=mesh.points, + faces=mesh.cell_connectivity.reshape((Nx * Ny * Nz, 8)), + n_points=mesh.points.shape[0], + n_faces=Nx * Ny * Nz, +) + +# dirichlet condition +on_lhs = mesh.points[:, 0] <= -Lx / 2 +# TODO should this be an n_node vector? +# dirichlet_indices = np.where(on_lhs)[0] +# dirichlet_mask = np.zeros(hex_mesh.n_points) +# dirichlet_mask[dirichlet_indices] = 1 +dirichlet_mask = np.where(on_lhs)[0] +dirichlet_values = np.zeros((dirichlet_mask.shape[0])) + +# von Neumann condition +x_lim = Lx/2 +y_min = 0 # -Ly/2 # TODO +y_max = y_min + 0.2 * Ly +z_min = 0 # 0.0 - 0.1 * Lz +z_max = z_min + 0.2 * Lz +von_neumann = np.logical_and( + mesh.points[:, 0] >= x_lim, + np.logical_and(mesh.points[:,1] >= y_min, mesh.points[:,1] <= y_max), + np.logical_and(mesh.points[:,2] >= z_min, mesh.points[:,2] <= z_max), +) +# TODO should this be an n_node array? +von_neumann_mask = np.where(von_neumann)[0] +von_neumann_values = (0, 0.0, 0.1 / len(von_neumann_mask)) + np.zeros((von_neumann_mask.shape[0], 3)) -n_elem = Nx * Ny * Nz # Create a test density field varying from 0 to 1 +n_elem = Nx * Ny * Nz rho = (np.arange(0, n_elem, 1) / n_elem).reshape((n_elem, 1)) rho = 0.5 * np.ones((n_elem, 1)) - - + + inputs = { + "dirichlet_mask": dirichlet_mask, + "dirichlet_values": dirichlet_values, + "van_neumann_mask": von_neumann_mask, + "van_neumann_values": von_neumann_values, + "hex_mesh": dict(hex_mesh), "host": host, "port": port, "rho": rho, - "Lx": Lx, - "Ly": Ly, - "Lz": Lz, - "Nx": Nx, - "Ny": Ny, - "Nz": Nz, "E0": 1.0, "rho_min": 1e-6, - "total_force": 0.1, "log_level": "DEBUG", "vtk_output": "mesh_density.vtk", } -outputs = tess_simp_compliance.apply(inputs) - -# Verify relationship between compliance and strain energy -# For static analysis: Total Strain Energy = 0.5 * Compliance -strain_energy = outputs["strain_energy"] -compliance = outputs["compliance"] -total_strain_energy = np.sum(strain_energy) -print(f"\nCompliance: {compliance:.6e}") -print(f"Total Strain Energy: {total_strain_energy:.6e}") -print(f"0.5 * Compliance: {0.5 * compliance:.6e}") -print(f"Ratio (should be ~1.0): {total_strain_energy / (0.5 * compliance):.6f}") - -# Finite difference check -num_tests = 0 # set to 0 if you don't want to run this check -FD_delta = 1.0e-3 -f0 = outputs["compliance"] -sensitivity = outputs["sensitivity"] -FD_sensitivity = 0 * sensitivity -for i in range(num_tests): - print(i) - inputs["rho"][i] += FD_delta - outputs = tess_simp_compliance.apply(inputs) - fupp = outputs["compliance"] - FD_sensitivity[i] = (fupp - f0) / FD_delta - inputs["rho"][i] -= FD_delta +# cells = np.array([8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 8, 9, 10, 11, 12, 13, 14, 15]) +cells = np.concatenate([np.array([[8]] * hex_mesh.n_faces), hex_mesh.faces], 1) +print(cells) +celltypes = np.full(hex_mesh.n_faces, pv.CellType.HEXAHEDRON, dtype=np.uint8) +mesh2 = pv.UnstructuredGrid(cells.ravel(), celltypes, hex_mesh.points) +print(mesh2 == mesh) +print(mesh.cell_connectivity) +print(mesh2.cell_connectivity) +#outputs = tess_simp_compliance.apply(inputs) +# +## Verify relationship between compliance and strain energy +## For static analysis: Total Strain Energy = 0.5 * Compliance +#strain_energy = outputs["strain_energy"] +#compliance = outputs["compliance"] +#total_strain_energy = np.sum(strain_energy) +#print(f"\nCompliance: {compliance:.6e}") +#print(f"Total Strain Energy: {total_strain_energy:.6e}") +#print(f"0.5 * Compliance: {0.5 * compliance:.6e}") +#print(f"Ratio (should be ~1.0): {total_strain_energy / (0.5 * compliance):.6f}") -if num_tests > 0: - sens = sensitivity[0:num_tests] - FD_sens = FD_sensitivity[0:num_tests] - print(sens) - print(FD_sens) - errors = sens - FD_sens - print(errors) - rel_abs_error = np.abs(errors / sens) - print(rel_abs_error) - print(f"Should be under 1e-5: {np.max(rel_abs_error)}") +# # Finite difference check +# num_tests = 0 # set to 0 if you don't want to run this check +# FD_delta = 1.0e-3 +# f0 = outputs["compliance"] +# sensitivity = outputs["sensitivity"] +# FD_sensitivity = 0 * sensitivity +# for i in range(num_tests): +# print(i) +# inputs["rho"][i] += FD_delta +# outputs = tess_simp_compliance.apply(inputs) +# fupp = outputs["compliance"] +# FD_sensitivity[i] = (fupp - f0) / FD_delta +# inputs["rho"][i] -= FD_delta +# +# +# if num_tests > 0: +# sens = sensitivity[0:num_tests] +# FD_sens = FD_sensitivity[0:num_tests] +# print(sens) +# print(FD_sens) +# errors = sens - FD_sens +# print(errors) +# rel_abs_error = np.abs(errors / sens) +# print(rel_abs_error) +# print(f"Should be under 1e-5: {np.max(rel_abs_error)}") From 4ec7e45fa660175c9a4c0c185fbb5a36beaae839 Mon Sep 17 00:00:00 2001 From: Nolan Black Date: Mon, 17 Nov 2025 12:25:07 -0500 Subject: [PATCH 3/7] fixed generic hexmesh in inputschema, with ansys debug option --- examples/ansys/pymapdl_tess/tesseract_api.py | 181 ++++++++++--------- examples/ansys/pymapdl_tess/test.py | 109 +++++------ 2 files changed, 157 insertions(+), 133 deletions(-) diff --git a/examples/ansys/pymapdl_tess/tesseract_api.py b/examples/ansys/pymapdl_tess/tesseract_api.py index 664639b..6f347fb 100644 --- a/examples/ansys/pymapdl_tess/tesseract_api.py +++ b/examples/ansys/pymapdl_tess/tesseract_api.py @@ -28,10 +28,10 @@ # TODO delete me once debugging is over Lx, Ly, Lz = 3, 2, 1 +#Nx, Ny, Nz = 60, 40, 20 Nx, Ny, Nz = 6, 4, 2 -# Nx, Ny, Nz = 60, 40, 20 force_val = 0.1 -debug = True +use_ansys = False P = ParamSpec("P") T = TypeVar("T") @@ -209,18 +209,13 @@ def solve(self) -> OutputSchema: """ logger.info("Starting SIMP elasticity analysis...") - if False: - self._create_geometry() - else: - self._create_geometry2() - self._define_element() - if False: + if use_ansys: self._create_mesh() else: self._create_mesh2() self._define_simp_materials() self._assign_materials_to_elements() - if debug: + if use_ansys: self._apply_boundary_conditions() else: self._apply_boundary_conditions_2() @@ -239,26 +234,6 @@ def solve(self) -> OutputSchema: return self._build_output_schema() - @log_timing - def _create_geometry2(self) -> None: - """Create a rectangular block geometry. - - Creates 8 keypoints at the corners of a rectangular box and defines - a volume from these keypoints. - """ - # start pre-processor - self.mapdl.prep7() - - # create a pyvista mesh from self.hex_mesh - cells = np.concatenate([np.array([[8]] * self.hex_mesh.n_faces), self.hex_mesh.faces], 1) - cell_type = np.array([pv.CellType.HEXAHEDRON] * self.hex_mesh.n_faces) - mesh = pv.UnstructuredGrid(cells, cell_type, self.hex_mesh.points) - - # save the pyvista mesh, then reload it in mapdl - archive_filename = os.path.join(tempfile.gettempdir(), "tmp.cdb") - save_as_archive(archive_filename, mesh) - response = self.mapdl.cdread("db", archive_filename) - @log_timing def _create_geometry(self) -> None: """Create a rectangular block geometry. @@ -294,7 +269,7 @@ def _define_simp_materials(self) -> None: if len(rho_flat) != self.n_elements: raise ValueError( f"Density field size {len(rho_flat)} does not match " - f"number of elements {self.n_elements} (Nx={self.Nx}, Ny={self.Ny}, Nz={self.Nz})" + f"number of elements {self.n_elements}" ) # Apply minimum density constraint to avoid singular stiffness matrix @@ -350,46 +325,75 @@ def _assign_materials_to_elements(self) -> None: @log_timing def _create_mesh2(self) -> None: - """Create structured hexahedral mesh with specified divisions. + """Create hexahedral mesh directly using MAPDL commands. - Sets line divisions for each direction and generates a swept mesh using - VSWEEP. The resulting mesh will have Nx * Ny * Nz elements with uniform - element sizes in each direction. + Creates nodes and elements from hex_mesh using direct MAPDL commands (N and E). + This ensures perfect 1:1 correspondence: node i+1 = points[i], element i+1 = faces[i]. """ - # Store element numbers in the order they were created - # This is crucial for maintaining consistent indexing between - # density input and strain energy output - self.element_numbers = self.mapdl.mesh.enum + self.mapdl.prep7() + self._define_element() - # Validate and log element numbering - n_elements = self.hex_mesh.n_faces - if len(self.element_numbers) != n_elements: - raise ValueError( - f"Number of created elements {len(self.element_numbers)} does not match " - f"expected {n_elements} (Nx={self.Nx}, Ny={self.Ny}, Nz={self.Nz})" - ) - self.n_elements = n_elements + logger.info(f"Creating {self.hex_mesh.n_points} nodes...") - # Repeat validation for nodes + # Create nodes using batch commands for efficiency + node_commands = [] + for i in range(self.hex_mesh.n_points): + x, y, z = self.hex_mesh.points[i] + # Node number i+1 corresponds to point index i + node_commands.append(f"N,{i+1},{x},{y},{z}") + + self.mapdl.input_strings(node_commands) + + logger.info(f"Creating {self.hex_mesh.n_faces} hexahedral elements...") + + # Create elements using batch commands + element_commands = [] + for i in range(self.hex_mesh.n_faces): + # Get node indices for this element (0-based) + node_indices = self.hex_mesh.faces[i] + # Convert to node numbers (1-based) + node_nums = [int(idx + 1) for idx in node_indices] + # Element number i+1 corresponds to face index i + element_commands.append(f"E,{','.join(map(str, node_nums))}") + + self.mapdl.input_strings(element_commands) + + # Verify the mesh was created correctly + self.element_numbers = self.mapdl.mesh.enum self.node_numbers = self.mapdl.mesh.nnum + self.n_elements = len(self.element_numbers) self.n_nodes = len(self.node_numbers) - logger.debug( - f"Element numbers: min={self.element_numbers.min()}, " - f"max={self.element_numbers.max()}, " - f"first 5={self.element_numbers[:5]}" + logger.info( + f"Mesh created: {self.n_nodes} nodes, {self.n_elements} elements" ) - # Check if element numbering is sequential starting from 1 - is_sequential = np.array_equal( - self.element_numbers, np.arange(1, n_elements + 1) + # Validate counts + if self.n_nodes != self.hex_mesh.n_points: + raise ValueError( + f"Node count mismatch: created {self.n_nodes}, expected {self.hex_mesh.n_points}" + ) + + if self.n_elements != self.hex_mesh.n_faces: + raise ValueError( + f"Element count mismatch: created {self.n_elements}, expected {self.hex_mesh.n_faces}" + ) + + # Verify sequential numbering (should always be true with direct creation) + is_sequential_nodes = np.array_equal( + self.node_numbers, np.arange(1, self.n_nodes + 1) ) - if is_sequential: - logger.info("Element numbering is sequential (1, 2, 3, ...)") - else: + is_sequential_elements = np.array_equal( + self.element_numbers, np.arange(1, self.n_elements + 1) + ) + + if not is_sequential_nodes or not is_sequential_elements: logger.warning( - "Element numbering is NOT sequential - reordering will be applied" + f"Non-sequential numbering detected: nodes={is_sequential_nodes}, elements={is_sequential_elements}" ) + else: + logger.info("Verified: Node and element numbering is sequential (1, 2, 3, ...)") + @log_timing def _create_mesh(self) -> None: """Create structured hexahedral mesh with specified divisions. @@ -398,6 +402,12 @@ def _create_mesh(self) -> None: VSWEEP. The resulting mesh will have Nx * Ny * Nz elements with uniform element sizes in each direction. """ + # create the cube geometry + self._create_geometry() + + # assign element + self._define_element() + # Set element divisions for lines in each direction using ndiv parameter # Lines parallel to X-axis (select by 2 constant coordinates) self.mapdl.lsel("s", "loc", "y", 0) @@ -466,7 +476,7 @@ def _create_mesh(self) -> None: if len(self.element_numbers) != n_elements: raise ValueError( f"Number of created elements {len(self.element_numbers)} does not match " - f"expected {n_elements} (Nx={self.Nx}, Ny={self.Ny}, Nz={self.Nz})" + f"expected {n_elements}" ) self.n_elements = n_elements @@ -491,37 +501,46 @@ def _create_mesh(self) -> None: "Element numbering is NOT sequential - reordering will be applied" ) + @log_timing def _apply_boundary_conditions_2(self) -> None: - """Apply boundary conditions, forces, and visualize. + """Apply boundary conditions and forces using batch commands for efficiency. - Fixes all nodes on the x=0 plane by constraining all translational - degrees of freedom (UX, UY, UZ). Then applies a concentrated force - in the negative z-direction to nodes at x=Lx, y=Ly, and z between - 0.4*Lz and 0.6*Lz. + Since nodes are created with direct correspondence (node i+1 = points[i]), + we can directly reference node numbers in D and F commands. All commands + are batched using input_strings for optimal performance. """ - for i, x in zip(self.dirichlet_mask, self.dirichlet_values): - # Select node - self.mapdl.nsel("s", "NODE", "", i + 1) + commands = [] - # set dirichlet conditions - self.mapdl.d("all", "all", x) + # Build Dirichlet boundary condition commands + n_dirichlet = len(self.dirichlet_mask) + logger.info(f"Applying Dirichlet BCs to {n_dirichlet} nodes") - # Reselect all nodes - self.mapdl.allsel() + for i, x in zip(self.dirichlet_mask, self.dirichlet_values): + # Point index i corresponds to node number i+1 + node = i + 1 + # D,NODE,ALL,VALUE constrains all DOFs to value x + commands.append(f"D,{node},ALL,{x}") + # Build von Neumann boundary condition (force) commands + n_von_neumann = len(self.van_neumann_mask) + logger.info(f"Applying von Neumann BCs to {n_von_neumann} nodes") for i, f in zip(self.van_neumann_mask, self.van_neumann_values): - # Select node - self.mapdl.nsel("s", "NODE", "", i + 1) - - # Apply force in negative z-direction to all selected nodes - self.mapdl.f("all", "fx", f[0]) - self.mapdl.f("all", "fy", f[1]) - self.mapdl.f("all", "fz", f[2]) - - # Reselect all nodes - self.mapdl.allsel() + # Point index i corresponds to node number i+1 + node = i + 1 + # F,NODE,Lab,VALUE applies force component + # Only add commands for non-zero force components + if f[0] != 0: + commands.append(f"F,{node},FX,{f[0]}") + if f[1] != 0: + commands.append(f"F,{node},FY,{f[1]}") + if f[2] != 0: + commands.append(f"F,{node},FZ,{f[2]}") + + # Execute all boundary condition commands at once + logger.debug(f"Executing {len(commands)} boundary condition commands") + self.mapdl.input_strings(commands) @log_timing def _apply_boundary_conditions(self) -> None: diff --git a/examples/ansys/pymapdl_tess/test.py b/examples/ansys/pymapdl_tess/test.py index 944bd9f..4e9d7ea 100644 --- a/examples/ansys/pymapdl_tess/test.py +++ b/examples/ansys/pymapdl_tess/test.py @@ -32,7 +32,7 @@ # create a simple Hex mesh using Pyvista import pyvista as pv Lx, Ly, Lz = 3, 2, 1 -# Nx, Ny, Nz = 60, 40, 20 +Nx, Ny, Nz = 60, 40, 20 Nx, Ny, Nz = 6, 4, 2 grid = pv.ImageData( dimensions=np.array((Nx, Ny, Nz)) + 1, @@ -51,8 +51,8 @@ n_faces=Nx * Ny * Nz, ) -# dirichlet condition -on_lhs = mesh.points[:, 0] <= -Lx / 2 +# dirichlet condition (select nodes at x=0) +on_lhs = mesh.points[:, 0] <= 0 # TODO should this be an n_node vector? # dirichlet_indices = np.where(on_lhs)[0] # dirichlet_mask = np.zeros(hex_mesh.n_points) @@ -60,16 +60,18 @@ dirichlet_mask = np.where(on_lhs)[0] dirichlet_values = np.zeros((dirichlet_mask.shape[0])) -# von Neumann condition -x_lim = Lx/2 -y_min = 0 # -Ly/2 # TODO -y_max = y_min + 0.2 * Ly -z_min = 0 # 0.0 - 0.1 * Lz -z_max = z_min + 0.2 * Lz +# von Neumann condition (select nodes at x=Lx with constraints on y and z) +x_lim = Lx +y_min = 0 +y_max = 0.2 * Ly +z_min = 0.4 * Lz +z_max = 0.6 * Lz von_neumann = np.logical_and( mesh.points[:, 0] >= x_lim, - np.logical_and(mesh.points[:,1] >= y_min, mesh.points[:,1] <= y_max), - np.logical_and(mesh.points[:,2] >= z_min, mesh.points[:,2] <= z_max), + np.logical_and( + np.logical_and(mesh.points[:,1] >= y_min, mesh.points[:,1] <= y_max), + np.logical_and(mesh.points[:,2] >= z_min, mesh.points[:,2] <= z_max) + ) ) # TODO should this be an n_node array? von_neumann_mask = np.where(von_neumann)[0] @@ -98,47 +100,50 @@ # cells = np.array([8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 8, 9, 10, 11, 12, 13, 14, 15]) cells = np.concatenate([np.array([[8]] * hex_mesh.n_faces), hex_mesh.faces], 1) -print(cells) celltypes = np.full(hex_mesh.n_faces, pv.CellType.HEXAHEDRON, dtype=np.uint8) mesh2 = pv.UnstructuredGrid(cells.ravel(), celltypes, hex_mesh.points) -print(mesh2 == mesh) print(mesh.cell_connectivity) -print(mesh2.cell_connectivity) - -#outputs = tess_simp_compliance.apply(inputs) -# -## Verify relationship between compliance and strain energy -## For static analysis: Total Strain Energy = 0.5 * Compliance -#strain_energy = outputs["strain_energy"] -#compliance = outputs["compliance"] -#total_strain_energy = np.sum(strain_energy) -#print(f"\nCompliance: {compliance:.6e}") -#print(f"Total Strain Energy: {total_strain_energy:.6e}") -#print(f"0.5 * Compliance: {0.5 * compliance:.6e}") -#print(f"Ratio (should be ~1.0): {total_strain_energy / (0.5 * compliance):.6f}") - -# # Finite difference check -# num_tests = 0 # set to 0 if you don't want to run this check -# FD_delta = 1.0e-3 -# f0 = outputs["compliance"] -# sensitivity = outputs["sensitivity"] -# FD_sensitivity = 0 * sensitivity -# for i in range(num_tests): -# print(i) -# inputs["rho"][i] += FD_delta -# outputs = tess_simp_compliance.apply(inputs) -# fupp = outputs["compliance"] -# FD_sensitivity[i] = (fupp - f0) / FD_delta -# inputs["rho"][i] -= FD_delta -# -# -# if num_tests > 0: -# sens = sensitivity[0:num_tests] -# FD_sens = FD_sensitivity[0:num_tests] -# print(sens) -# print(FD_sens) -# errors = sens - FD_sens -# print(errors) -# rel_abs_error = np.abs(errors / sens) -# print(rel_abs_error) -# print(f"Should be under 1e-5: {np.max(rel_abs_error)}") + +outputs = tess_simp_compliance.apply(inputs) + +# Verify relationship between compliance and strain energy +# For static analysis: Total Strain Energy = 0.5 * Compliance +strain_energy = outputs["strain_energy"] +compliance = outputs["compliance"] +total_strain_energy = np.sum(strain_energy) +print(f"\nCompliance: {compliance:.6e}") +print(f"Total Strain Energy: {total_strain_energy:.6e}") +print(f"0.5 * Compliance: {0.5 * compliance:.6e}") +print(f"Ratio (should be ~1.0): {total_strain_energy / (0.5 * compliance):.6f}") + +# Finite difference check +num_tests = 5 # set to 0 if you don't want to run this check +FD_delta = 1.0e-3 +f0 = outputs["compliance"] +sensitivity = outputs["sensitivity"] +FD_sensitivity = 0 * sensitivity +for i in range(num_tests): + print(i) + inputs["rho"][i] += FD_delta + outputs = tess_simp_compliance.apply(inputs) + fupp = outputs["compliance"] + # FD_sensitivity[i] = (fupp - f0) / FD_delta + # inputs["rho"][i] -= FD_delta + + inputs["rho"][i] -= 2.0*FD_delta + outputs = tess_simp_compliance.apply(inputs) + fdown = outputs["compliance"] + FD_sensitivity[i] = (fupp - fdown) / FD_delta / 2.0 + inputs["rho"][i] += FD_delta + + +if num_tests > 0: + sens = sensitivity[0:num_tests] + FD_sens = FD_sensitivity[0:num_tests] + print(sens) + print(FD_sens) + errors = sens - FD_sens + print(errors) + rel_abs_error = np.abs(errors / sens) + print(rel_abs_error) + print(f"Should be under 1e-5: {np.max(rel_abs_error)}") From da32bb79203e7e8aadc1045b2ab1f9374657c959 Mon Sep 17 00:00:00 2001 From: Nolan Black Date: Mon, 17 Nov 2025 12:28:42 -0500 Subject: [PATCH 4/7] remove deprecated options --- examples/ansys/pymapdl_tess/tesseract_api.py | 215 ++----------------- 1 file changed, 14 insertions(+), 201 deletions(-) diff --git a/examples/ansys/pymapdl_tess/tesseract_api.py b/examples/ansys/pymapdl_tess/tesseract_api.py index 6f347fb..cd4807c 100644 --- a/examples/ansys/pymapdl_tess/tesseract_api.py +++ b/examples/ansys/pymapdl_tess/tesseract_api.py @@ -11,7 +11,6 @@ import logging import time -import tempfile from collections.abc import Callable from functools import wraps from typing import ParamSpec, TypeVar @@ -19,22 +18,16 @@ import numpy as np import pyvista as pv from ansys.mapdl.core import Mapdl -from ansys.mapdl.reader import save_as_archive from pydantic import BaseModel, Field from tesseract_core.runtime import Array, Differentiable, Float32, Int32 # Set up module logger logger = logging.getLogger(__name__) -# TODO delete me once debugging is over -Lx, Ly, Lz = 3, 2, 1 -#Nx, Ny, Nz = 60, 40, 20 -Nx, Ny, Nz = 6, 4, 2 -force_val = 0.1 -use_ansys = False - P = ParamSpec("P") T = TypeVar("T") + + def log_timing(func: Callable[P, T]) -> Callable[P, T]: """Decorator to log the wall time of method execution.""" @@ -51,6 +44,7 @@ def wrapper(*args: P.args, **kwargs: P.kwargs) -> T: return wrapper + class HexMesh(BaseModel): """Hexahedral mesh representation.""" @@ -65,6 +59,7 @@ class HexMesh(BaseModel): default=0, description="Number of valid faces in the faces array." ) + class InputSchema(BaseModel): """Inputs for the tess_simp_compliance Tesseract. @@ -209,16 +204,10 @@ def solve(self) -> OutputSchema: """ logger.info("Starting SIMP elasticity analysis...") - if use_ansys: - self._create_mesh() - else: - self._create_mesh2() + self._create_mesh() self._define_simp_materials() self._assign_materials_to_elements() - if use_ansys: - self._apply_boundary_conditions() - else: - self._apply_boundary_conditions_2() + self._apply_boundary_conditions() self._run_analysis() self._extract_displacement_constraints() self._extract_nodal_displacement() @@ -234,27 +223,6 @@ def solve(self) -> OutputSchema: return self._build_output_schema() - @log_timing - def _create_geometry(self) -> None: - """Create a rectangular block geometry. - - Creates 8 keypoints at the corners of a rectangular box and defines - a volume from these keypoints. - """ - # start pre-processor - self.mapdl.prep7() - - k0 = self.mapdl.k("", 0, 0, 0) - k1 = self.mapdl.k("", Lx, 0, 0) - k2 = self.mapdl.k("", Lx, Ly, 0) - k3 = self.mapdl.k("", 0, Ly, 0) - k4 = self.mapdl.k("", 0, 0, Lz) - k5 = self.mapdl.k("", Lx, 0, Lz) - k6 = self.mapdl.k("", Lx, Ly, Lz) - k7 = self.mapdl.k("", 0, Ly, Lz) - - self.mapdl.v(k0, k1, k2, k3, k4, k5, k6, k7) - @log_timing def _define_simp_materials(self) -> None: """Define materials using SIMP approach with batch commands for efficiency. @@ -324,7 +292,7 @@ def _assign_materials_to_elements(self) -> None: self.mapdl.input_strings(commands) @log_timing - def _create_mesh2(self) -> None: + def _create_mesh(self) -> None: """Create hexahedral mesh directly using MAPDL commands. Creates nodes and elements from hex_mesh using direct MAPDL commands (N and E). @@ -340,7 +308,7 @@ def _create_mesh2(self) -> None: for i in range(self.hex_mesh.n_points): x, y, z = self.hex_mesh.points[i] # Node number i+1 corresponds to point index i - node_commands.append(f"N,{i+1},{x},{y},{z}") + node_commands.append(f"N,{i + 1},{x},{y},{z}") self.mapdl.input_strings(node_commands) @@ -364,9 +332,7 @@ def _create_mesh2(self) -> None: self.n_elements = len(self.element_numbers) self.n_nodes = len(self.node_numbers) - logger.info( - f"Mesh created: {self.n_nodes} nodes, {self.n_elements} elements" - ) + logger.info(f"Mesh created: {self.n_nodes} nodes, {self.n_elements} elements") # Validate counts if self.n_nodes != self.hex_mesh.n_points: @@ -392,118 +358,12 @@ def _create_mesh2(self) -> None: f"Non-sequential numbering detected: nodes={is_sequential_nodes}, elements={is_sequential_elements}" ) else: - logger.info("Verified: Node and element numbering is sequential (1, 2, 3, ...)") - - @log_timing - def _create_mesh(self) -> None: - """Create structured hexahedral mesh with specified divisions. - - Sets line divisions for each direction and generates a swept mesh using - VSWEEP. The resulting mesh will have Nx * Ny * Nz elements with uniform - element sizes in each direction. - """ - # create the cube geometry - self._create_geometry() - - # assign element - self._define_element() - - # Set element divisions for lines in each direction using ndiv parameter - # Lines parallel to X-axis (select by 2 constant coordinates) - self.mapdl.lsel("s", "loc", "y", 0) - self.mapdl.lsel("r", "loc", "z", 0) - self.mapdl.lesize("all", ndiv=Nx, kforc=1) - - self.mapdl.lsel("s", "loc", "y", Ly) - self.mapdl.lsel("r", "loc", "z", 0) - self.mapdl.lesize("all", ndiv=Nx, kforc=1) - - self.mapdl.lsel("s", "loc", "y", 0) - self.mapdl.lsel("r", "loc", "z", Lz) - self.mapdl.lesize("all", ndiv=Nx, kforc=1) - - self.mapdl.lsel("s", "loc", "y", Ly) - self.mapdl.lsel("r", "loc", "z", Lz) - self.mapdl.lesize("all", ndiv=Nx, kforc=1) - - # Lines parallel to Y-axis (select by 2 constant coordinates) - self.mapdl.lsel("s", "loc", "x", 0) - self.mapdl.lsel("r", "loc", "z", 0) - self.mapdl.lesize("all", ndiv=Ny, kforc=1) - - self.mapdl.lsel("s", "loc", "x", Lx) - self.mapdl.lsel("r", "loc", "z", 0) - self.mapdl.lesize("all", ndiv=Ny, kforc=1) - - self.mapdl.lsel("s", "loc", "x", 0) - self.mapdl.lsel("r", "loc", "z", Lz) - self.mapdl.lesize("all", ndiv=Ny, kforc=1) - - self.mapdl.lsel("s", "loc", "x", Lx) - self.mapdl.lsel("r", "loc", "z", Lz) - self.mapdl.lesize("all", ndiv=Ny, kforc=1) - - # Lines parallel to Z-axis (select by 2 constant coordinates) - self.mapdl.lsel("s", "loc", "x", 0) - self.mapdl.lsel("r", "loc", "y", 0) - self.mapdl.lesize("all", ndiv=Nz, kforc=1) - - self.mapdl.lsel("s", "loc", "x", Lx) - self.mapdl.lsel("r", "loc", "y", 0) - self.mapdl.lesize("all", ndiv=Nz, kforc=1) - - self.mapdl.lsel("s", "loc", "x", 0) - self.mapdl.lsel("r", "loc", "y", Ly) - self.mapdl.lesize("all", ndiv=Nz, kforc=1) - - self.mapdl.lsel("s", "loc", "x", Lx) - self.mapdl.lsel("r", "loc", "y", Ly) - self.mapdl.lesize("all", ndiv=Nz, kforc=1) - - # Restore full selection of entities - self.mapdl.allsel() - - # Mesh the volume using the swept mesh command - self.mapdl.vsweep("all") - - # Store element numbers in the order they were created - # This is crucial for maintaining consistent indexing between - # density input and strain energy output - self.element_numbers = self.mapdl.mesh.enum - - # Validate and log element numbering - n_elements = self.hex_mesh.n_faces - if len(self.element_numbers) != n_elements: - raise ValueError( - f"Number of created elements {len(self.element_numbers)} does not match " - f"expected {n_elements}" - ) - self.n_elements = n_elements - - # Repeat validation for nodes - self.node_numbers = self.mapdl.mesh.nnum - self.n_nodes = len(self.node_numbers) - - logger.debug( - f"Element numbers: min={self.element_numbers.min()}, " - f"max={self.element_numbers.max()}, " - f"first 5={self.element_numbers[:5]}" - ) - - # Check if element numbering is sequential starting from 1 - is_sequential = np.array_equal( - self.element_numbers, np.arange(1, n_elements + 1) - ) - if is_sequential: - logger.info("Element numbering is sequential (1, 2, 3, ...)") - else: - logger.warning( - "Element numbering is NOT sequential - reordering will be applied" + logger.info( + "Verified: Node and element numbering is sequential (1, 2, 3, ...)" ) - @log_timing - def _apply_boundary_conditions_2(self) -> None: + def _apply_boundary_conditions(self) -> None: """Apply boundary conditions and forces using batch commands for efficiency. Since nodes are created with direct correspondence (node i+1 = points[i]), @@ -516,7 +376,7 @@ def _apply_boundary_conditions_2(self) -> None: n_dirichlet = len(self.dirichlet_mask) logger.info(f"Applying Dirichlet BCs to {n_dirichlet} nodes") - for i, x in zip(self.dirichlet_mask, self.dirichlet_values): + for i, x in zip(self.dirichlet_mask, self.dirichlet_values, strict=True): # Point index i corresponds to node number i+1 node = i + 1 # D,NODE,ALL,VALUE constrains all DOFs to value x @@ -526,7 +386,7 @@ def _apply_boundary_conditions_2(self) -> None: n_von_neumann = len(self.van_neumann_mask) logger.info(f"Applying von Neumann BCs to {n_von_neumann} nodes") - for i, f in zip(self.van_neumann_mask, self.van_neumann_values): + for i, f in zip(self.van_neumann_mask, self.van_neumann_values, strict=True): # Point index i corresponds to node number i+1 node = i + 1 # F,NODE,Lab,VALUE applies force component @@ -542,53 +402,6 @@ def _apply_boundary_conditions_2(self) -> None: logger.debug(f"Executing {len(commands)} boundary condition commands") self.mapdl.input_strings(commands) - @log_timing - def _apply_boundary_conditions(self) -> None: - """Apply boundary conditions, forces, and visualize. - - Fixes all nodes on the x=0 plane by constraining all translational - degrees of freedom (UX, UY, UZ). Then applies a concentrated force - in the negative z-direction to nodes at x=Lx, y=Ly, and z between - 0.4*Lz and 0.6*Lz. - """ - # Select all nodes on the x=0 plane - self.mapdl.nsel("s", "loc", "x", 0) - - # Apply displacement constraints (fix all DOFs: UX, UY, UZ) - self.mapdl.d("all", "all", 0) - - # Reselect all nodes - self.mapdl.allsel() - - # Apply force to nodes at x=Lx, y=0, z between 0.4*Lz and 0.6*Lz - # Select nodes at x=Lx - self.mapdl.nsel("s", "loc", "x", Lx) - - # Refine selection between y=0 and y=0.2 Ly - y_max = 0.2 * Ly - self.mapdl.nsel("r", "loc", "y", 0, y_max) - - # Refine selection to z between 0.4*Lz and 0.6*Lz - z_min = 0.4 * Lz - z_max = 0.6 * Lz - self.mapdl.nsel("r", "loc", "z", z_min, z_max) - - # Get number of selected nodes - num_nodes = self.mapdl.mesh.n_node - - # Calculate force per node (negative for downward direction) - force_per_node = -force_val / num_nodes - - # Apply force in negative z-direction to all selected nodes - self.mapdl.f("all", "fz", force_per_node) - - # Reselect all nodes - self.mapdl.allsel() - - logger.info( - f"Applied total force of {force_val} N distributed over {num_nodes} nodes" - ) - @log_timing def _run_analysis(self) -> None: """Run static structural analysis.""" From 9b6e8729c280064405ae86fad0ab46659da93117 Mon Sep 17 00:00:00 2001 From: Nolan Black Date: Mon, 17 Nov 2025 12:30:13 -0500 Subject: [PATCH 5/7] removed some unnec code --- examples/ansys/pymapdl_tess/test.py | 70 +++++++++++++++-------------- 1 file changed, 37 insertions(+), 33 deletions(-) diff --git a/examples/ansys/pymapdl_tess/test.py b/examples/ansys/pymapdl_tess/test.py index 4e9d7ea..c8858d2 100644 --- a/examples/ansys/pymapdl_tess/test.py +++ b/examples/ansys/pymapdl_tess/test.py @@ -1,49 +1,57 @@ -""" -/usr/ansys_inc/v241/ansys/bin/ansys241 -port 50005 -grpc -""" +"""/usr/ansys_inc/v241/ansys/bin/ansys241 -port 50005 -grpc""" + import os import numpy as np from tesseract_core import Tesseract - # load the MAPDL server's address from your environment host = os.getenv("MAPDL_HOST") if host is None: - raise ValueError("Unable to read $MAPDL_HOST from the environment. " + \ - "Use 'export MAPDL_HOST=X.X.X.X' for local IP address of your MAPDL Instance.") + raise ValueError( + "Unable to read $MAPDL_HOST from the environment. " + + "Use 'export MAPDL_HOST=X.X.X.X' for local IP address of your MAPDL Instance." + ) port = os.getenv("MAPDL_PORT") if port is None: - raise ValueError("Unable to read $MAPDL_PORT from the environment. " + \ - "Use 'export MAPDL_PORT=X' for the port of your MAPDL Instance.") + raise ValueError( + "Unable to read $MAPDL_PORT from the environment. " + + "Use 'export MAPDL_PORT=X' for the port of your MAPDL Instance." + ) # Initialize Tesseract from api (for testing) try: - tess_simp_compliance = Tesseract.from_tesseract_api( - "./tesseract_api.py" + tess_simp_compliance = Tesseract.from_tesseract_api("./tesseract_api.py") +except RuntimeError: + raise RuntimeError( + "Unable to load tesseract from api. " + "Ensure that you have installed the build requirements using 'pip install -r tesseract_requirements.txt'" ) -except RuntimeError as e: - raise RuntimeError("Unable to load tesseract from api. " \ - "Ensure that you have installed the build requirements using 'pip install -r tesseract_requirements.txt'") - + # TODO # we should probably create some re-usable mesh constructors... # create a simple Hex mesh using Pyvista import pyvista as pv + Lx, Ly, Lz = 3, 2, 1 Nx, Ny, Nz = 60, 40, 20 Nx, Ny, Nz = 6, 4, 2 grid = pv.ImageData( dimensions=np.array((Nx, Ny, Nz)) + 1, - #origin=(-Lx / 2, -Ly / 2, -Lz / 2), # The bottom left corner of the data set - origin=(0, 0, 0), # TODO - spacing=(Lx / Nx, Ly / Ny, Lz / Nz), # These are the cell sizes along each axis + # origin=(-Lx / 2, -Ly / 2, -Lz / 2), # The bottom left corner of the data set + origin=(0, 0, 0), # TODO + spacing=(Lx / Nx, Ly / Ny, Lz / Nz), # These are the cell sizes along each axis ) # repeated casts will eventually expose cell_connectivitiy -mesh = grid.cast_to_structured_grid().cast_to_explicit_structured_grid().cast_to_unstructured_grid() +mesh = ( + grid.cast_to_structured_grid() + .cast_to_explicit_structured_grid() + .cast_to_unstructured_grid() +) from tesseract_api import HexMesh + hex_mesh = HexMesh( points=mesh.points, faces=mesh.cell_connectivity.reshape((Nx * Ny * Nz, 8)), @@ -58,7 +66,7 @@ # dirichlet_mask = np.zeros(hex_mesh.n_points) # dirichlet_mask[dirichlet_indices] = 1 dirichlet_mask = np.where(on_lhs)[0] -dirichlet_values = np.zeros((dirichlet_mask.shape[0])) +dirichlet_values = np.zeros(dirichlet_mask.shape[0]) # von Neumann condition (select nodes at x=Lx with constraints on y and z) x_lim = Lx @@ -69,20 +77,22 @@ von_neumann = np.logical_and( mesh.points[:, 0] >= x_lim, np.logical_and( - np.logical_and(mesh.points[:,1] >= y_min, mesh.points[:,1] <= y_max), - np.logical_and(mesh.points[:,2] >= z_min, mesh.points[:,2] <= z_max) - ) + np.logical_and(mesh.points[:, 1] >= y_min, mesh.points[:, 1] <= y_max), + np.logical_and(mesh.points[:, 2] >= z_min, mesh.points[:, 2] <= z_max), + ), ) # TODO should this be an n_node array? von_neumann_mask = np.where(von_neumann)[0] -von_neumann_values = (0, 0.0, 0.1 / len(von_neumann_mask)) + np.zeros((von_neumann_mask.shape[0], 3)) +von_neumann_values = (0, 0.0, 0.1 / len(von_neumann_mask)) + np.zeros( + (von_neumann_mask.shape[0], 3) +) # Create a test density field varying from 0 to 1 n_elem = Nx * Ny * Nz rho = (np.arange(0, n_elem, 1) / n_elem).reshape((n_elem, 1)) rho = 0.5 * np.ones((n_elem, 1)) - - + + inputs = { "dirichlet_mask": dirichlet_mask, "dirichlet_values": dirichlet_values, @@ -98,12 +108,6 @@ "vtk_output": "mesh_density.vtk", } -# cells = np.array([8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 8, 9, 10, 11, 12, 13, 14, 15]) -cells = np.concatenate([np.array([[8]] * hex_mesh.n_faces), hex_mesh.faces], 1) -celltypes = np.full(hex_mesh.n_faces, pv.CellType.HEXAHEDRON, dtype=np.uint8) -mesh2 = pv.UnstructuredGrid(cells.ravel(), celltypes, hex_mesh.points) -print(mesh.cell_connectivity) - outputs = tess_simp_compliance.apply(inputs) # Verify relationship between compliance and strain energy @@ -117,7 +121,7 @@ print(f"Ratio (should be ~1.0): {total_strain_energy / (0.5 * compliance):.6f}") # Finite difference check -num_tests = 5 # set to 0 if you don't want to run this check +num_tests = 0 # set to 0 if you don't want to run this check FD_delta = 1.0e-3 f0 = outputs["compliance"] sensitivity = outputs["sensitivity"] @@ -130,7 +134,7 @@ # FD_sensitivity[i] = (fupp - f0) / FD_delta # inputs["rho"][i] -= FD_delta - inputs["rho"][i] -= 2.0*FD_delta + inputs["rho"][i] -= 2.0 * FD_delta outputs = tess_simp_compliance.apply(inputs) fdown = outputs["compliance"] FD_sensitivity[i] = (fupp - fdown) / FD_delta / 2.0 From dfc4732aa3d77338d4b27fef6b67299caf805850 Mon Sep 17 00:00:00 2001 From: Nolan Black Date: Mon, 17 Nov 2025 14:41:39 -0500 Subject: [PATCH 6/7] improved test, made notes for future ref --- examples/ansys/pymapdl_tess/tesseract_api.py | 79 ++++-- examples/ansys/pymapdl_tess/test.py | 274 ++++++++++--------- 2 files changed, 212 insertions(+), 141 deletions(-) diff --git a/examples/ansys/pymapdl_tess/tesseract_api.py b/examples/ansys/pymapdl_tess/tesseract_api.py index cd4807c..b405da2 100644 --- a/examples/ansys/pymapdl_tess/tesseract_api.py +++ b/examples/ansys/pymapdl_tess/tesseract_api.py @@ -13,13 +13,13 @@ import time from collections.abc import Callable from functools import wraps -from typing import ParamSpec, TypeVar +from typing import Any, ParamSpec, TypeVar import numpy as np import pyvista as pv from ansys.mapdl.core import Mapdl from pydantic import BaseModel, Field -from tesseract_core.runtime import Array, Differentiable, Float32, Int32 +from tesseract_core.runtime import Array, Differentiable, Float32, Int32, ShapeDType # Set up module logger logger = logging.getLogger(__name__) @@ -110,7 +110,7 @@ class InputSchema(BaseModel): ) log_level: str = Field( - default="INFO", + default="WARNING", description="Logging level for output messages (DEBUG, INFO, WARNING, ERROR).", ) @@ -522,9 +522,11 @@ def _calculate_sensitivity(self) -> None: inverse_rho = np.nan_to_num(1 / self.rho.flatten(), nan=0.0) self.sensitivity = -2.0 * self.p * inverse_rho * self.strain_energy.flatten() - # TODO improve this cache? - # stash the sensitivity s.t. it may be loaded in the vjp - np.save("sensitivity.npy", self.sensitivity) + # Cache sensitivity in temporary directory for use in VJP + import tempfile + + sensitivity_path = os.path.join(tempfile.gettempdir(), "sensitivity.npy") + np.save(sensitivity_path, self.sensitivity) @log_timing def _create_pvmesh(self) -> None: @@ -616,16 +618,55 @@ def apply(inputs: InputSchema) -> OutputSchema: return solver.solve() -# TODO -# def vector_jacobian_product( -# inputs: InputSchema, -# vjp_inputs: set[str], -# vjp_outputs: set[str], -# cotangent_vector: dict[str, Any], -# ): -# pass -# -# -# def abstract_eval(abstract_inputs): -# """Calculate output shape of apply from the shape of its inputs.""" -# return {"compliance": ShapeDType(shape=(), dtype="float32")} +def vector_jacobian_product( + inputs: InputSchema, + vjp_inputs: set[str], + vjp_outputs: set[str], + cotangent_vector: dict[str, Any], +) -> dict[str, Any]: + """Compute vector-Jacobian product for backpropagation through ANSYS solve. + + The sensitivity (dcompliance/drho) is already computed and cached during + the forward pass in _calculate_sensitivity(), so we just need to load it + and multiply by the upstream gradient. + + Args: + inputs: Original inputs to the forward pass (InputSchema) + vjp_inputs: Set of input names we need gradients for (e.g., {"rho"}) + vjp_outputs: Set of output names being differentiated (e.g., {"compliance"}) + cotangent_vector: Upstream gradients from the loss function + (e.g., {"compliance": dloss/dcompliance}) + + Returns: + Dictionary mapping input names to their gradients + (e.g., {"rho": dloss/∂rho}) + + """ + gradients = {} + assert vjp_inputs == {"rho"} + assert vjp_outputs == {"compliance"} + + # Load the cached sensitivity (∂compliance/∂rho) from temporary directory + # This was computed and saved in _calculate_sensitivity() + import tempfile + + sensitivity_path = os.path.join(tempfile.gettempdir(), "sensitivity.npy") + sensitivity = np.load(sensitivity_path) + + # Clean up the temporary file after loading + try: + os.unlink(sensitivity_path) + except FileNotFoundError: + pass # Already deleted, no problem + + grad_rho_flat = cotangent_vector["compliance"] * sensitivity + + # Reshape to match input rho shape + gradients["rho"] = grad_rho_flat.reshape(inputs.rho.shape) + + return gradients + + +def abstract_eval(abstract_inputs): + """Calculate output shape of apply from the shape of its inputs.""" + return {"compliance": ShapeDType(shape=(), dtype="float32")} diff --git a/examples/ansys/pymapdl_tess/test.py b/examples/ansys/pymapdl_tess/test.py index c8858d2..42c389f 100644 --- a/examples/ansys/pymapdl_tess/test.py +++ b/examples/ansys/pymapdl_tess/test.py @@ -1,11 +1,12 @@ -"""/usr/ansys_inc/v241/ansys/bin/ansys241 -port 50005 -grpc""" - import os import numpy as np +import pyvista as pv +from tesseract_api import HexMesh # TODO can we standardize these mesh definitions? from tesseract_core import Tesseract # load the MAPDL server's address from your environment +# usr/ansys_inc/v241/ansys/bin/ansys241 -port 50050 -grpc host = os.getenv("MAPDL_HOST") if host is None: raise ValueError( @@ -23,131 +24,160 @@ # Initialize Tesseract from api (for testing) try: tess_simp_compliance = Tesseract.from_tesseract_api("./tesseract_api.py") -except RuntimeError: +except RuntimeError as e: raise RuntimeError( "Unable to load tesseract from api. " "Ensure that you have installed the build requirements using 'pip install -r tesseract_requirements.txt'" + ) from e + + +def mesh_from_pyvista(Lx, Ly, Lz, Nx, Ny, Nz): + grid = pv.ImageData( + dimensions=np.array((Nx, Ny, Nz)) + 1, + # origin=(-Lx / 2, -Ly / 2, -Lz / 2), # The bottom left corner of the data set + origin=(0, 0, 0), # TODO + spacing=(Lx / Nx, Ly / Ny, Lz / Nz), # These are the cell sizes along each axis + ) + # repeated casts will eventually expose cell_connectivitiy + mesh = ( + grid.cast_to_structured_grid() + .cast_to_explicit_structured_grid() + .cast_to_unstructured_grid() ) -# TODO -# we should probably create some re-usable mesh constructors... -# create a simple Hex mesh using Pyvista -import pyvista as pv + hex_mesh = HexMesh( + points=mesh.points, + faces=mesh.cell_connectivity.reshape((Nx * Ny * Nz, 8)), + n_points=mesh.points.shape[0], + n_faces=Nx * Ny * Nz, + ) + return hex_mesh + + +def cantilever_bc(Lx, Ly, Lz, Nx, Ny, Nz, hex_mesh): + # dirichlet condition (select nodes at x=0) + on_lhs = hex_mesh.points[:, 0] <= 0 + # TODO should this be an n_node vector? + # dirichlet_indices = np.where(on_lhs)[0] + # dirichlet_mask = np.zeros(hex_mesh.n_points) + # dirichlet_mask[dirichlet_indices] = 1 + dirichlet_mask = np.where(on_lhs)[0] + dirichlet_values = np.zeros(dirichlet_mask.shape[0]) + + # von Neumann condition (select nodes at x=Lx with constraints on y and z) + x_lim = Lx + y_min = 0 + y_max = 0.2 * Ly + z_min = 0.4 * Lz + z_max = 0.6 * Lz + von_neumann = np.logical_and( + hex_mesh.points[:, 0] >= x_lim, + np.logical_and( + np.logical_and( + hex_mesh.points[:, 1] >= y_min, hex_mesh.points[:, 1] <= y_max + ), + np.logical_and( + hex_mesh.points[:, 2] >= z_min, hex_mesh.points[:, 2] <= z_max + ), + ), + ) + # TODO should this be an n_node array? + von_neumann_mask = np.where(von_neumann)[0] + von_neumann_values = np.array([0, 0.0, 0.1 / len(von_neumann_mask)]) + np.zeros( + (von_neumann_mask.shape[0], 3) + ) + return dirichlet_mask, dirichlet_values, von_neumann_mask, von_neumann_values -Lx, Ly, Lz = 3, 2, 1 -Nx, Ny, Nz = 60, 40, 20 -Nx, Ny, Nz = 6, 4, 2 -grid = pv.ImageData( - dimensions=np.array((Nx, Ny, Nz)) + 1, - # origin=(-Lx / 2, -Ly / 2, -Lz / 2), # The bottom left corner of the data set - origin=(0, 0, 0), # TODO - spacing=(Lx / Nx, Ly / Ny, Lz / Nz), # These are the cell sizes along each axis -) -# repeated casts will eventually expose cell_connectivitiy -mesh = ( - grid.cast_to_structured_grid() - .cast_to_explicit_structured_grid() - .cast_to_unstructured_grid() -) - -from tesseract_api import HexMesh - -hex_mesh = HexMesh( - points=mesh.points, - faces=mesh.cell_connectivity.reshape((Nx * Ny * Nz, 8)), - n_points=mesh.points.shape[0], - n_faces=Nx * Ny * Nz, -) - -# dirichlet condition (select nodes at x=0) -on_lhs = mesh.points[:, 0] <= 0 -# TODO should this be an n_node vector? -# dirichlet_indices = np.where(on_lhs)[0] -# dirichlet_mask = np.zeros(hex_mesh.n_points) -# dirichlet_mask[dirichlet_indices] = 1 -dirichlet_mask = np.where(on_lhs)[0] -dirichlet_values = np.zeros(dirichlet_mask.shape[0]) - -# von Neumann condition (select nodes at x=Lx with constraints on y and z) -x_lim = Lx -y_min = 0 -y_max = 0.2 * Ly -z_min = 0.4 * Lz -z_max = 0.6 * Lz -von_neumann = np.logical_and( - mesh.points[:, 0] >= x_lim, - np.logical_and( - np.logical_and(mesh.points[:, 1] >= y_min, mesh.points[:, 1] <= y_max), - np.logical_and(mesh.points[:, 2] >= z_min, mesh.points[:, 2] <= z_max), - ), -) -# TODO should this be an n_node array? -von_neumann_mask = np.where(von_neumann)[0] -von_neumann_values = (0, 0.0, 0.1 / len(von_neumann_mask)) + np.zeros( - (von_neumann_mask.shape[0], 3) -) - -# Create a test density field varying from 0 to 1 -n_elem = Nx * Ny * Nz -rho = (np.arange(0, n_elem, 1) / n_elem).reshape((n_elem, 1)) -rho = 0.5 * np.ones((n_elem, 1)) - - -inputs = { - "dirichlet_mask": dirichlet_mask, - "dirichlet_values": dirichlet_values, - "van_neumann_mask": von_neumann_mask, - "van_neumann_values": von_neumann_values, - "hex_mesh": dict(hex_mesh), - "host": host, - "port": port, - "rho": rho, - "E0": 1.0, - "rho_min": 1e-6, - "log_level": "DEBUG", - "vtk_output": "mesh_density.vtk", -} - -outputs = tess_simp_compliance.apply(inputs) - -# Verify relationship between compliance and strain energy -# For static analysis: Total Strain Energy = 0.5 * Compliance -strain_energy = outputs["strain_energy"] -compliance = outputs["compliance"] -total_strain_energy = np.sum(strain_energy) -print(f"\nCompliance: {compliance:.6e}") -print(f"Total Strain Energy: {total_strain_energy:.6e}") -print(f"0.5 * Compliance: {0.5 * compliance:.6e}") -print(f"Ratio (should be ~1.0): {total_strain_energy / (0.5 * compliance):.6f}") - -# Finite difference check -num_tests = 0 # set to 0 if you don't want to run this check -FD_delta = 1.0e-3 -f0 = outputs["compliance"] -sensitivity = outputs["sensitivity"] -FD_sensitivity = 0 * sensitivity -for i in range(num_tests): - print(i) - inputs["rho"][i] += FD_delta - outputs = tess_simp_compliance.apply(inputs) - fupp = outputs["compliance"] - # FD_sensitivity[i] = (fupp - f0) / FD_delta - # inputs["rho"][i] -= FD_delta - inputs["rho"][i] -= 2.0 * FD_delta +def sample_rho(hex_mesh): + # Create a test density field varying from 0 to 1 + n_elem = hex_mesh.n_faces + rho = (np.arange(0, n_elem, 1) / n_elem).reshape((n_elem, 1)) + rho = 0.5 * np.ones((n_elem, 1)) + return rho + + +def main(): + Lx, Ly, Lz = 3, 2, 1 + Nx, Ny, Nz = 60, 40, 20 + Nx, Ny, Nz = 6, 4, 2 + num_tests = 5 # set to 0 if you don't want to run this check + run_central_difference = True + + hex_mesh = mesh_from_pyvista(Lx, Ly, Lz, Nx, Ny, Nz) + dirichlet_mask, dirichlet_values, von_neumann_mask, von_neumann_values = ( + cantilever_bc(Lx, Ly, Lz, Nx, Ny, Nz, hex_mesh) + ) + rho = sample_rho(hex_mesh) + + inputs = { + "dirichlet_mask": dirichlet_mask, + "dirichlet_values": dirichlet_values, + "van_neumann_mask": von_neumann_mask, + "van_neumann_values": von_neumann_values, + "hex_mesh": dict(hex_mesh), + "host": host, + "port": port, + "rho": rho, + "E0": 1.0, + "rho_min": 1e-6, + "log_level": "DEBUG", + "vtk_output": "mesh_density.vtk", + } + outputs = tess_simp_compliance.apply(inputs) - fdown = outputs["compliance"] - FD_sensitivity[i] = (fupp - fdown) / FD_delta / 2.0 - inputs["rho"][i] += FD_delta - - -if num_tests > 0: - sens = sensitivity[0:num_tests] - FD_sens = FD_sensitivity[0:num_tests] - print(sens) - print(FD_sens) - errors = sens - FD_sens - print(errors) - rel_abs_error = np.abs(errors / sens) - print(rel_abs_error) - print(f"Should be under 1e-5: {np.max(rel_abs_error)}") + + # Verify relationship between compliance and strain energy + # For static analysis: Total Strain Energy = 0.5 * Compliance + strain_energy = outputs["strain_energy"] + compliance = outputs["compliance"] + total_strain_energy = np.sum(strain_energy) + print(f"\nCompliance: {compliance:.6e}") + print(f"Total Strain Energy: {total_strain_energy:.6e}") + print(f"0.5 * Compliance: {0.5 * compliance:.6e}") + print(f"Ratio (should be ~1.0): {total_strain_energy / (0.5 * compliance):.6f}") + + # a sample backwards pass + f0 = outputs["compliance"] + vjp = tess_simp_compliance.vector_jacobian_product( + inputs, {"rho"}, {"compliance"}, {"compliance": 1.0} + ) + sensitivity = vjp["rho"] + + # Finite difference check + FD_delta = 1.0e-4 + FD_sensitivity = 0 * sensitivity + for i in range(num_tests): + print(i) + inputs["rho"][i] += FD_delta + outputs = tess_simp_compliance.apply(inputs) + fupp = outputs["compliance"] + + if run_central_difference: + inputs["rho"][i] -= 2.0 * FD_delta + outputs = tess_simp_compliance.apply(inputs) + fdown = outputs["compliance"] + FD_sensitivity[i] = (fupp - fdown) / FD_delta / 2.0 + inputs["rho"][i] += FD_delta + + else: + FD_sensitivity[i] = (fupp - f0) / FD_delta + inputs["rho"][i] -= FD_delta + + if num_tests > 0: + sens = sensitivity[0:num_tests] + FD_sens = FD_sensitivity[0:num_tests] + print(sens) + print(FD_sens) + errors = sens - FD_sens + print(errors) + rel_abs_error = np.abs(errors / sens) + if run_central_difference: + print(f"Relative error should be near O({FD_delta})") + else: + print(f"Relative error should be O({FD_delta * 10})") + print(rel_abs_error) + + +if __name__ == "__main__": + main() From 4eb61784659418accb983545c27f5b714d6a195c Mon Sep 17 00:00:00 2001 From: Nolan Black Date: Tue, 18 Nov 2025 15:23:09 -0500 Subject: [PATCH 7/7] improve test comments --- examples/ansys/pymapdl_tess/test.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/examples/ansys/pymapdl_tess/test.py b/examples/ansys/pymapdl_tess/test.py index 42c389f..bc8c55c 100644 --- a/examples/ansys/pymapdl_tess/test.py +++ b/examples/ansys/pymapdl_tess/test.py @@ -2,7 +2,7 @@ import numpy as np import pyvista as pv -from tesseract_api import HexMesh # TODO can we standardize these mesh definitions? +from tesseract_api import HexMesh from tesseract_core import Tesseract # load the MAPDL server's address from your environment @@ -34,8 +34,7 @@ def mesh_from_pyvista(Lx, Ly, Lz, Nx, Ny, Nz): grid = pv.ImageData( dimensions=np.array((Nx, Ny, Nz)) + 1, - # origin=(-Lx / 2, -Ly / 2, -Lz / 2), # The bottom left corner of the data set - origin=(0, 0, 0), # TODO + origin=(0, 0, 0), spacing=(Lx / Nx, Ly / Ny, Lz / Nz), # These are the cell sizes along each axis ) # repeated casts will eventually expose cell_connectivitiy @@ -55,13 +54,10 @@ def mesh_from_pyvista(Lx, Ly, Lz, Nx, Ny, Nz): def cantilever_bc(Lx, Ly, Lz, Nx, Ny, Nz, hex_mesh): + # Create a dirichlet_mask of nodes indices associated with diricelet condition # dirichlet condition (select nodes at x=0) on_lhs = hex_mesh.points[:, 0] <= 0 - # TODO should this be an n_node vector? - # dirichlet_indices = np.where(on_lhs)[0] - # dirichlet_mask = np.zeros(hex_mesh.n_points) - # dirichlet_mask[dirichlet_indices] = 1 - dirichlet_mask = np.where(on_lhs)[0] + dirichlet_mask = np.where(on_lhs)[0] # size (num_dirichlet_nodes,) dirichlet_values = np.zeros(dirichlet_mask.shape[0]) # von Neumann condition (select nodes at x=Lx with constraints on y and z) @@ -81,7 +77,7 @@ def cantilever_bc(Lx, Ly, Lz, Nx, Ny, Nz, hex_mesh): ), ), ) - # TODO should this be an n_node array? + # A (num_von_neumann_nodes, n_dof) array von_neumann_mask = np.where(von_neumann)[0] von_neumann_values = np.array([0, 0.0, 0.1 / len(von_neumann_mask)]) + np.zeros( (von_neumann_mask.shape[0], 3) @@ -101,7 +97,7 @@ def main(): Lx, Ly, Lz = 3, 2, 1 Nx, Ny, Nz = 60, 40, 20 Nx, Ny, Nz = 6, 4, 2 - num_tests = 5 # set to 0 if you don't want to run this check + num_tests = 0 # set to 0 if you don't want to run this check run_central_difference = True hex_mesh = mesh_from_pyvista(Lx, Ly, Lz, Nx, Ny, Nz)