diff --git a/examples/ansys/pymapdl_fem_tess/tesseract_api.py b/examples/ansys/pymapdl_fem_tess/tesseract_api.py
new file mode 100644
index 0000000..33d04ad
--- /dev/null
+++ b/examples/ansys/pymapdl_fem_tess/tesseract_api.py
@@ -0,0 +1,689 @@
+# Copyright 2025 Pasteur Labs. All Rights Reserved.
+# SPDX-License-Identifier: Apache-2.0
+
+# create a temp directory to work in
+import os
+
+# Check if HOME is writable, if not, set it to /tmp
+if "HOME" not in os.environ or not os.access(os.environ.get("HOME", ""), os.W_OK):
+    os.environ["HOME"] = "/tmp"
+
+
+import logging
+import time
+from collections.abc import Callable
+from functools import wraps
+from typing import 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
+
+# Set up module logger
+logger = logging.getLogger(__name__)
+
+# Type variables for decorator
+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."""
+
+    @wraps(func)
+    def wrapper(*args: P.args, **kwargs: P.kwargs) -> T:
+        method_name = func.__name__
+        start_time = time.time()
+        logger.info(f"Starting {method_name}...")
+        result = func(*args, **kwargs)
+        end_time = time.time()
+        elapsed = end_time - start_time
+        logger.info(f"Completed {method_name} in {elapsed:.4f} seconds")
+        return result
+
+    return wrapper
+
+
+class InputSchema(BaseModel):
+    """Inputs for the tess_simp_compliance Tesseract.
+
+    Relies on an instance of ANSYS MAPDL, as in:
+    `/f/Ansys/ANSYS Inc/v241/ansys/bin/winx64/ANSYS241.exe -grpc -port 50052`
+    """
+
+    host: str = Field(description="The IP of your MAPDL grpc server.")
+    port: str = Field(description="The port of your MAPDL grpc server.")
+    rho: Differentiable[
+        Array[
+            (
+                None,
+                None,
+            ),
+            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."
+    )
+    Ly: float = Field(
+        default=30.0,
+        description=("Length of the simulation box in the y direction."),
+    )
+    Lz: float = Field(
+        default=30.0, description="Length of the simulation box in the z direction."
+    )
+
+    Nx: int = Field(
+        default=60,
+        description=("Number of elements in the x direction."),
+    )
+    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."),
+    )
+
+    E0: float = Field(
+        default=1.0,
+        description="Base Young's modulus in Pa for SIMP material interpolation.",
+    )
+
+    rho_min: float = Field(
+        default=1e-6,
+        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).",
+    )
+
+    log_level: str = Field(
+        default="INFO",
+        description="Logging level for output messages (DEBUG, INFO, WARNING, ERROR).",
+    )
+
+    vtk_output: str | None = Field(
+        default=None,
+        description="The path to write the results in VTK format.",
+    )
+
+
+class OutputSchema(BaseModel):
+    """Outputs for the tess_simp_compliance Tesseract."""
+
+    compliance: Differentiable[
+        Array[
+            (),
+            Float32,
+        ]
+    ] = Field(description="Compliance of the structure, a measure of stiffness")
+    strain_energy: Array[(None,), Float32] = Field(
+        description="The element-wise strain energy of the solution."
+    )
+    sensitivity: Array[(None,), Float32] = Field(
+        description="The deriviative of compliance with respect to each densitiy variable."
+    )
+
+
+class SIMPElasticity:
+    """SIMP-based topology optimization solver using ANSYS MAPDL.
+
+    This class encapsulates the full workflow for running a SIMP
+    (Solid Isotropic Material with Penalization) topology optimization
+    analysis using ANSYS MAPDL for finite element analysis.
+    """
+
+    def __init__(self, inputs: InputSchema, mapdl: Mapdl) -> None:
+        """Initialize the SIMP elasticity solver.
+
+        Args:
+            inputs: Input parameters for the simulation
+            mapdl: Active MAPDL instance
+        """
+        # Store inputs
+        self.inputs = inputs
+        self.mapdl = mapdl
+
+        # 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.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
+
+        # Configure logger
+        logger.setLevel(getattr(logging, self.log_level.upper(), logging.INFO))
+        # Add console handler if not already present
+        if not logger.handlers:
+            console_handler = logging.StreamHandler()
+            console_handler.setLevel(logging.DEBUG)
+            formatter = logging.Formatter(
+                "%(asctime)s - %(name)s - %(levelname)s - %(message)s"
+            )
+            console_handler.setFormatter(formatter)
+            logger.addHandler(console_handler)
+
+        # Initialize result storage
+        self.element_numbers = None  # Store actual MAPDL element numbers
+        self.n_elements = None
+        self.node_numbers = None  # Store actual MAPDL node numbers
+        self.n_nodes = None
+        self.displacement_constraints = None
+        self.nodal_displacement = None
+        self.nodal_force = None
+        self.strain_energy = None
+        self.compliance = None
+        self.sensitivity = None
+        self.pvmesh = None
+
+    @log_timing
+    def solve(self) -> OutputSchema:
+        """Run the complete SIMP analysis workflow.
+
+        Returns:
+            OutputSchema containing compliance, strain energy, and sensitivity
+        """
+        logger.info("Starting SIMP elasticity analysis...")
+
+        self._create_geometry()
+        self._define_element()
+        self._create_mesh()
+        self._define_simp_materials()
+        self._assign_materials_to_elements()
+        self._apply_boundary_conditions()
+        self._run_analysis()
+        self._extract_displacement_constraints()
+        self._extract_nodal_displacement()
+        self._extract_nodal_force()
+        self._extract_strain_energy()
+        self._calculate_compliance()
+        self._calculate_sensitivity()
+        if self.vtk_output:
+            self._create_pvmesh()
+
+        logger.info("SIMP analysis complete!")
+        logger.debug(f"MAPDL status: {self.mapdl}")
+
+        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("", 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)
+
+        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.
+
+        Uses the SIMP (Solid Isotropic Material with Penalization) method to
+        scale material properties based on density values. Creates a unique
+        material for each element.
+        """
+        # Flatten rho array and ensure it matches the number of elements
+        rho_flat = np.array(self.rho).flatten()
+
+        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})"
+            )
+
+        # Apply minimum density constraint to avoid singular stiffness matrix
+        # Calculate all material properties at once
+        # SIMP formula: E = E0 * rho^p
+        E_values = self.E0 * (self.rho_min + (1 - self.rho_min) * (rho_flat**self.p))
+        dens_values = self.rho_min + (1 - self.rho_min) * rho_flat
+
+        # Build all commands as a list for batch execution
+        commands = []
+        for i in range(len(rho_flat)):
+            mat_id = i + 1
+            commands.append(f"MP,EX,{mat_id},{E_values[i]}")
+            commands.append(f"MP,DENS,{mat_id},{dens_values[i]}")
+            commands.append(f"MP,NUXY,{mat_id},0.3")
+
+        # Execute all commands at once using input_strings
+        self.mapdl.input_strings(commands)
+
+    @log_timing
+    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.
+        """
+        self.mapdl.et(1, "SOLID186")
+
+    @log_timing
+    def _assign_materials_to_elements(self) -> None:
+        """Assign materials to elements after meshing using batch commands.
+
+        For a structured mesh with Nx x Ny x Nz elements, this assigns
+        material ID to each element based on the stored element numbers.
+        Material i+1 is assigned to the element at position i in the element_numbers array.
+        """
+        # Get all element numbers from the mesh
+        self.mapdl.allsel()
+
+        # Build all EMODIF commands as a list for batch execution
+        # Use actual element numbers instead of assuming sequential 1, 2, 3, ...
+        commands = []
+        for i in range(self.n_elements):
+            elem_id = self.element_numbers[i]  # Use actual element number from mesh
+            mat_id = i + 1  # Material index corresponds to density array index
+            commands.append(f"EMODIF,{elem_id},MAT,{mat_id}")
+
+        # Execute all commands at once using input_strings
+        self.mapdl.input_strings(commands)
+
+    @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.
+        """
+        # 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=self.Nx, kforc=1)
+
+        self.mapdl.lsel("s", "loc", "y", self.Ly)
+        self.mapdl.lsel("r", "loc", "z", 0)
+        self.mapdl.lesize("all", ndiv=self.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("s", "loc", "y", self.Ly)
+        self.mapdl.lsel("r", "loc", "z", self.Lz)
+        self.mapdl.lesize("all", ndiv=self.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.lsel("s", "loc", "x", self.Lx)
+        self.mapdl.lsel("r", "loc", "z", 0)
+        self.mapdl.lesize("all", ndiv=self.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("s", "loc", "x", self.Lx)
+        self.mapdl.lsel("r", "loc", "z", self.Lz)
+        self.mapdl.lesize("all", ndiv=self.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.lsel("s", "loc", "x", self.Lx)
+        self.mapdl.lsel("r", "loc", "y", 0)
+        self.mapdl.lesize("all", ndiv=self.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("s", "loc", "x", self.Lx)
+        self.mapdl.lsel("r", "loc", "y", self.Ly)
+        self.mapdl.lesize("all", ndiv=self.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.Nx * self.Ny * self.Nz
+        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 _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", self.Lx)
+
+        # Refine selection between y=0 and y=0.2 Ly
+        y_max = 0.2 * self.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
+        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
+
+        # 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 {self.total_force} N distributed over {num_nodes} nodes"
+        )
+
+    @log_timing
+    def _run_analysis(self) -> None:
+        """Run static structural analysis."""
+        self.mapdl.slashsolu()
+        self.mapdl.allsel()  # making sure all nodes and elements are selected.
+        self.mapdl.antype("STATIC")
+        self.mapdl.ematwrite("YES")
+        output = self.mapdl.solve()
+        self.mapdl.save("my_analysis")
+        self.mapdl.finish()
+
+        logger.debug(f"Analysis output: {output}")
+
+    @log_timing
+    def _extract_strain_energy(self) -> None:
+        """Extract strain energy for all elements and save to file.
+
+        Uses APDL commands to extract strain energy from element table and
+        write to a text file, then reorders to match input density array indexing.
+
+        The *VGET command retrieves data in element number order (sorted by element ID),
+        so we must reorder the results to match the element_numbers array order.
+        """
+        logger.debug("Extracting strain energy data...")
+
+        # Enter POST1 mode if not already there
+        self.mapdl.post1()
+        self.mapdl.set(1, 1)  # Load first load step, first substep
+
+        # Create element table for strain energy
+        self.mapdl.etable("SENE", "SENE")
+
+        # Create array parameter and populate with strain energy values
+        # *VGET retrieves element table values into an array parameter
+        self.mapdl.run(f"*DIM,strain_e,ARRAY,{self.n_elements},1")
+        self.mapdl.run("*VGET,strain_e,ELEM,1,ETAB,SENE, , ,2")
+
+        # Write strain energy to text file using non_interactive mode for file operations
+        # this is much faster than using a get command
+        with self.mapdl.non_interactive:
+            self.mapdl.run("*CFOPEN,strain_energy,txt")
+            self.mapdl.run("*VWRITE,strain_e(1,1)")
+            self.mapdl.run("(E15.7)")
+            self.mapdl.run("*CFCLOS")
+
+        # Download the file from the MAPDL working directory
+        self.mapdl.download("strain_energy.txt", ".")
+
+        # Read the data - this is in sorted element number order
+        strain_energy_sorted = np.loadtxt("strain_energy.txt")
+
+        # Create mapping from element numbers to their position in sorted order
+        # The *VGET command returns data sorted by element ID
+        sorted_element_numbers = np.sort(self.element_numbers)
+
+        # Create inverse mapping: for each element in self.element_numbers,
+        # find its position in the sorted array
+        # This tells us where to find each element's data in strain_energy_sorted
+        reorder_indices = np.searchsorted(sorted_element_numbers, self.element_numbers)
+
+        # Reorder strain energy to match the order of self.element_numbers
+        # (which matches the input density array order)
+        self.strain_energy = strain_energy_sorted[reorder_indices]
+
+        logger.debug(
+            f"Strain energy reordered: first 5 values = {self.strain_energy[:5]}"
+        )
+
+    @log_timing
+    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
+        self.nodal_displacement = nodal_displacement
+
+    @log_timing
+    def _extract_nodal_force(self) -> None:
+        nodal_force = np.zeros((self.n_nodes, 3))
+        # populate force using vectorized NumPy advanced indexing
+        nnum, dof_idx, f = self.mapdl.result.nodal_input_force(0)
+        node_indices = nnum - 1
+        dof_indices = dof_idx - 1
+        nodal_force[node_indices, dof_indices] = f
+        self.nodal_force = nodal_force
+
+    @log_timing
+    def _extract_displacement_constraints(self) -> None:
+        displacement_constraints = np.zeros((self.n_nodes, 3))
+        # Get a list of all nodal constraints and process with optimized mapping
+        nodal_constraints = self.mapdl.get_nodal_constrains()
+
+        constraint_map = {"UX": 0, "UY": 1, "UZ": 2}
+        for idx, field, _, _ in nodal_constraints:
+            dof = constraint_map.get(field)
+            if dof is not None:
+                node_num = int(idx) - 1
+                displacement_constraints[node_num, dof] = 1.0
+        self.displacement_constraints = displacement_constraints
+
+    @log_timing
+    def _calculate_compliance(self) -> None:
+        """Calculate compliance from nodal forces and displacements."""
+        self.compliance = np.dot(
+            self.nodal_force.flatten(), self.nodal_displacement.flatten()
+        )
+
+    @log_timing
+    def _calculate_sensitivity(self) -> None:
+        """Calculate sensitivity of compliance with respect to density.
+
+        For this special case, the adjoint solution is equal to -U,
+        so the sensitivity is equal to:
+             - U_e^T * (dK_e / drho) * U_e
+              = - (p/rho) * U_e^T * K_e * U_e
+              = - 2 * (p/rho) * strain_energy
+        """
+        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)
+
+    @log_timing
+    def _create_pvmesh(self) -> None:
+        """Create PyVista grid with analysis results.
+
+        Creates a PyVista grid containing the mesh geometry with density and
+        Young's modulus values stored as cell data, and displacement constraints,
+        applied forces, and nodal displacements stored as point data (node data).
+        """
+        logger.debug("Creating PyVista results grid...")
+
+        # Enter POST1 to access solution results
+        self.mapdl.post1()
+        self.mapdl.set(1, 1)
+
+        # Try to get mesh directly from result object to avoid disk I/O
+        # The result.grid property provides direct access to PyVista mesh
+        try:
+            self.pvmesh = self.mapdl.result.grid
+            logger.debug("Successfully loaded mesh directly from result.grid")
+        except (AttributeError, Exception) as e:
+            # Fallback to file-based approach if direct access fails
+            logger.debug(f"Direct grid access failed ({e}), using file-based approach")
+            self.mapdl.download_result(".")
+
+            show_progress = logger.level <= logging.INFO
+            self.mapdl.result.save_as_vtk("file.vtk", progress_bar=show_progress)
+
+            self.pvmesh = pv.read("file.vtk")
+
+        # add attributes
+        rho_flat = np.array(self.rho).flatten()
+        self.pvmesh.cell_data["density"] = self.__convert_celldata_for_pv(rho_flat)
+        self.pvmesh.cell_data["strain_energy"] = self.__convert_celldata_for_pv(
+            self.strain_energy
+        )
+        self.pvmesh.cell_data["sensitivity"] = self.__convert_celldata_for_pv(
+            self.sensitivity
+        )
+        self.pvmesh.point_data["displacement_constraints"] = (
+            self.displacement_constraints
+        )
+        self.pvmesh.point_data["nodal_displacement"] = self.nodal_displacement
+        self.pvmesh.point_data["nodal_force"] = self.nodal_force
+
+        # Export to VTK
+        self.pvmesh.save(self.vtk_output)
+
+        logger.info("Exported results to mesh_density.vtk")
+
+    def __convert_celldata_for_pv(self, elem_array: np.array) -> np.array:
+        """Convert cell_data to follow pyvista cell ordering."""
+        # PyVista cells may be ordered differently than our element_numbers array
+        # Get PyVista cell to MAPDL element number mapping
+        pv_elem_nums = self.pvmesh.cell_data["ansys_elem_num"]
+
+        # Create mapping: for each PyVista cell, find the corresponding density
+        array_for_pv = np.zeros(len(pv_elem_nums))
+        for pv_idx, mapdl_elem_num in enumerate(pv_elem_nums):
+            # Find where this MAPDL element number appears in self.element_numbers
+            array_idx = np.where(self.element_numbers == mapdl_elem_num)[0][0]
+            array_for_pv[pv_idx] = elem_array[array_idx]
+        return array_for_pv
+
+    def _build_output_schema(self) -> OutputSchema:
+        """Build and return the output schema."""
+        return OutputSchema(
+            compliance=self.compliance,
+            strain_energy=self.strain_energy.flatten(),
+            sensitivity=self.sensitivity,
+        )
+
+
+def apply(inputs: InputSchema) -> OutputSchema:
+    """Run the tess_simp_compliance Tesseract.
+
+    Args:
+        inputs: Input parameters for the SIMP elasticity analysis
+
+    Returns:
+        OutputSchema containing compliance, strain energy, and sensitivity
+    """
+    # Initialize MAPDL
+    mapdl = Mapdl(inputs.host, port=inputs.port)
+    mapdl.clear()
+
+    # Create solver instance and run analysis
+    solver = SIMPElasticity(inputs, mapdl)
+    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")}
diff --git a/examples/ansys/pymapdl_fem_tess/tesseract_config.yaml b/examples/ansys/pymapdl_fem_tess/tesseract_config.yaml
new file mode 100644
index 0000000..f2932ec
--- /dev/null
+++ b/examples/ansys/pymapdl_fem_tess/tesseract_config.yaml
@@ -0,0 +1,27 @@
+# Tesseract configuration file
+
+name: "tess_simp_compliance"
+version: "0.0.0"
+description: "Evaluate compliance using the SIMP material model"
+
+build_config:
+  # Base image to use for the container, must be Ubuntu or Debian-based
+  # with Python 3.9 or later
+  # base_image: "python:3.12-slim-bookworm"
+
+  # Platform to build the container for. In general, images can only be executed
+  # on the platform they were built for.
+  # target_platform: "native"
+
+  # Additional packages to install in the container (via apt-get)
+  # extra_packages:
+  #  - package_name
+
+  # Data to copy into the container, relative to the project root
+  # package_data:
+  #  - [path/to/source, path/to/destination]
+
+  # Additional Dockerfile commands to run during the build process
+  # custom_build_steps:
+  #  - |
+  #   RUN echo "Hello, World!"
diff --git a/examples/ansys/pymapdl_fem_tess/tesseract_requirements.txt b/examples/ansys/pymapdl_fem_tess/tesseract_requirements.txt
new file mode 100644
index 0000000..f6eb739
--- /dev/null
+++ b/examples/ansys/pymapdl_fem_tess/tesseract_requirements.txt
@@ -0,0 +1 @@
+ansys-mapdl-core