diff --git a/examples/ansys/pymapdl_tess/tesseract_api.py b/examples/ansys/pymapdl_tess/tesseract_api.py
index 33d04ad..b405da2 100644
--- a/examples/ansys/pymapdl_tess/tesseract_api.py
+++ b/examples/ansys/pymapdl_tess/tesseract_api.py
@@ -13,18 +13,17 @@
 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
+from tesseract_core.runtime import Array, Differentiable, Float32, Int32, ShapeDType
 
 # Set up module logger
 logger = logging.getLogger(__name__)
 
-# Type variables for decorator
 P = ParamSpec("P")
 T = TypeVar("T")
 
@@ -46,6 +45,21 @@ 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 +78,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."
-    )
-
-    Nx: int = Field(
-        default=60,
-        description=("Number of elements in the x direction."),
+    dirichlet_mask: Array[(None,), Int32] = Field(
+        description="Mask for Dirichlet boundary conditions",
     )
-    Ny: int = Field(
-        default=30,
-        description=("Number of elements in the y direction."),
+    dirichlet_values: Array[(None,), Float32] = Field(
+        description="Values for Dirichlet boundary conditions",
     )
-    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,18 +104,13 @@ 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).",
     )
 
     log_level: str = Field(
-        default="INFO",
+        default="WARNING",
         description="Logging level for output messages (DEBUG, INFO, WARNING, ERROR).",
     )
 
@@ -158,15 +158,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,8 +204,6 @@ def solve(self) -> OutputSchema:
         """
         logger.info("Starting SIMP elasticity analysis...")
 
-        self._create_geometry()
-        self._define_element()
         self._create_mesh()
         self._define_simp_materials()
         self._assign_materials_to_elements()
@@ -226,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("", 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.
@@ -261,7 +237,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
@@ -285,10 +261,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:
@@ -314,151 +293,114 @@ def _assign_materials_to_elements(self) -> None:
 
     @log_timing
     def _create_mesh(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].
         """
-        # 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()
+        self.mapdl.prep7()
+        self._define_element()
+
+        logger.info(f"Creating {self.hex_mesh.n_points} 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...")
 
-        # Mesh the volume using the swept mesh command
-        self.mapdl.vsweep("all")
+        # 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))}")
 
-        # Store element numbers in the order they were created
-        # This is crucial for maintaining consistent indexing between
-        # density input and strain energy output
+        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)
 
-        # Validate and log element numbering
-        n_elements = self.Nx * self.Ny * self.Nz
-        if len(self.element_numbers) != n_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:
             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"Node count mismatch: created {self.n_nodes}, expected {self.hex_mesh.n_points}"
             )
-        self.n_elements = n_elements
 
-        # Repeat validation for nodes
-        self.node_numbers = self.mapdl.mesh.nnum
-        self.n_nodes = len(self.node_numbers)
+        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}"
+            )
 
-        logger.debug(
-            f"Element numbers: min={self.element_numbers.min()}, "
-            f"max={self.element_numbers.max()}, "
-            f"first 5={self.element_numbers[:5]}"
+        # 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)
         )
-
-        # Check if element numbering is sequential starting from 1
-        is_sequential = np.array_equal(
-            self.element_numbers, np.arange(1, n_elements + 1)
+        is_sequential_elements = np.array_equal(
+            self.element_numbers, np.arange(1, self.n_elements + 1)
         )
-        if is_sequential:
-            logger.info("Element numbering is sequential (1, 2, 3, ...)")
-        else:
+
+        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 _apply_boundary_conditions(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.
         """
-        # 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()
+        commands = []
 
-        logger.info(
-            f"Applied total force of {self.total_force} N distributed over {num_nodes} nodes"
-        )
+        # Build Dirichlet boundary condition commands
+        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, strict=True):
+            # 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, strict=True):
+            # 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 _run_analysis(self) -> None:
@@ -533,7 +475,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
@@ -580,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:
@@ -674,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 4187cba..bc8c55c 100644
--- a/examples/ansys/pymapdl_tess/test.py
+++ b/examples/ansys/pymapdl_tess/test.py
@@ -1,83 +1,179 @@
 import os
 
 import numpy as np
-from dotenv import load_dotenv
+import pyvista as pv
+from tesseract_api import HexMesh
 from tesseract_core import Tesseract
 
-load_dotenv()
-
+# 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("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.")
-
-
-tess_simp_compliance = Tesseract.from_tesseract_api(
-    "tesseract_pymapdl/tess_simp_compiance/tesseract_api.py"
-)
-Lx, Ly, Lz = 3, 2, 1
-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
-rho = (np.arange(0, n_elem, 1) / n_elem).reshape((n_elem, 1))
-rho = 0.5 * np.ones((n_elem, 1))
-
-
-inputs = {
-    "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
+    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")
+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=(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
+    mesh = (
+        grid.cast_to_structured_grid()
+        .cast_to_explicit_structured_grid()
+        .cast_to_unstructured_grid()
+    )
+
+    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):
+    # 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
+    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)
+    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
+            ),
+        ),
+    )
+    # 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)
+    )
+    return dirichlet_mask, dirichlet_values, von_neumann_mask, von_neumann_values
+
+
+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 = 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)
+    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)
-    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)}")
+
+    # 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()