diff --git a/src/configuration.py b/src/configuration.py
index 87e4612..7260139 100644
--- a/src/configuration.py
+++ b/src/configuration.py
@@ -21,7 +21,6 @@ def __init__(self, params=None):
     @classmethod
     def from_parameters(cls, params: Dict[str, Any]) -> 'Configuration':
         """Create and fully initialize Configuration from parameter dictionary."""
-        #print("Initializing vessel configuration...")
         cls._validate_parameters(params)
         
         config = cls(params)
@@ -31,13 +30,12 @@ def from_parameters(cls, params: Dict[str, Any]) -> 'Configuration':
             layer = Layer.from_parameters(layer_data)
             config.add_layer(layer)
         
-        # Compute homeostatic stress (direct method - no simulation params needed)
         # TODO: integrate this into Layer.from_parameters.
-        #print("\nComputing homeostatic stresses...")
         for layer in config.layers:
             layer.compute_homeostatic_stress_direct()
-    
-        #print(f"\nConfiguration initialized with {len(config.layers)} layers")
+
+        config.print_state()       
+
         return config
 
     @staticmethod
@@ -316,3 +314,28 @@ def write_timestep_data(self, io_handler: 'IOHandler', timestep: int, time: floa
                         const_data,
                         layer_name=layer.name
                     )
+
+# =========================================================================
+# OUTPUT: PRINTING
+# =========================================================================
+
+    def print_state(self, timestep: int = None):
+        """Print configuration state.
+        
+        Prints state of all layers and their constituents at specified timestep.
+            
+        Args:
+            timestep: Timestep to print data from. If None, prints homeostatic state.
+            
+        Examples:
+            config.print_state()            # Homeostatic state
+            config.print_state(timestep=10) # At timestep 10
+        """
+        if timestep is None:
+            self.logger.section("HOMEOSTATIC CONFIGURATION", log_level='INFO')
+        else:
+            self.logger.section(f"CONFIGURATION AT TIMESTEP {timestep}", log_level='INFO')
+        
+        # Print each layer
+        for layer in self.layers:
+            layer.print_state(timestep=timestep, logger=self.logger)
\ No newline at end of file
diff --git a/src/constituent.py b/src/constituent.py
index 0fbe177..cf116be 100644
--- a/src/constituent.py
+++ b/src/constituent.py
@@ -1,6 +1,7 @@
 import math
 import warnings
 import numpy as np
+from custom_logging import get_logger, CustomLogger
 from typing import List, Dict, Any, Optional
 from abc import ABC, abstractmethod
 from kinetics import Kinetics  
@@ -26,9 +27,9 @@ def __init__(self, name):
         self.layer: Optional['Layer'] = None  # Reference to parent layer (set by Layer.add_constituent)
         self.kinetics: Optional[Kinetics] = None
         self.constitutive_model: Optional[ConstitutiveModel] = None # TODO: consider making this mandatory for SingleConstituent. Never initialized for MultiFiberFamilyConstituent.
+        self.logger: CustomLogger = get_logger(__name__)
         
         self.homeostatic_referential_density: Optional[float] = None
-        self.deposition_stretch_scalar: Optional[float] = None
         self.deposition_stretch: Optional[np.ndarray] = None
 
         # Active stress properties (optional - only for contractile constituents)
@@ -41,13 +42,14 @@ def __init__(self, name):
         self.sigma_hat_act_history = []  # [σ^_act(0), σ^_act(1), ...] (tensors)
         self.active_stress_history = []  # [σ_act(0), σ_act(1), ...] (tensors)
 
-        # Referential mass density evolution (rhoR_alpha) - timestep-indexed history
+        # Mass density and kinetic variables histories
         self.rhoR_alpha_history = []  # List of mass densities over time
         self.survival_history = [] # List of survival fractions over time (for each cohort)
         self.k_alpha_history = []  # List of degradation rates of survival function over time
         self.mR_alpha_history = []  # List of production rates over time
         self.stimulus_function_history = []  # List of stimulus function values over time #TODO: only stored for plot currently, streamline with Plotting class
 
+        # Deformation and stress histories
         self.F_alpha_history: List[List[np.ndarray]] = []
         self.sigma_hat_history: List[List[np.ndarray]] = []  # List of constituent-based partial stress over time TODO: perhaps make temporary array?
         self.stress_history: List[np.ndarray] = []
@@ -180,9 +182,7 @@ def __init__(self, name):
          
     @classmethod
     def from_parameters(cls, name: str, properties: Dict[str, Any], layer: 'Layer' = None) -> 'SingleConstituent':
-        """Create and fully initialize single constituent from parameters."""
-        print(f"    Initializing single constituent: {name}")
-        
+        """Create and fully initialize single constituent from parameters."""        
         constituent = cls(name)
 
         # Set layer reference immediately if provided
@@ -192,12 +192,10 @@ def from_parameters(cls, name: str, properties: Dict[str, Any], layer: 'Layer' =
         # Initialize homeostatic mass density (t=0)
         constituent.homeostatic_referential_density = properties['mass_fraction'] * 1050.0
         constituent.rhoR_alpha_history.append(constituent.homeostatic_referential_density)
-        print(f"        Homeostatic mass density (rhoR_alpha at t=0): {constituent.homeostatic_referential_density:.2f} kg/m³")
 
         # Initialize constitutive model
         if 'constitutive_model' in properties:
             constituent.constitutive_model = ConstitutiveModel.from_parameters(properties['constitutive_model'])
-            print(f"        Constitutive model: {constituent.constitutive_model}")
         
         # Initialize active properties if present
         if 'active_properties' in properties:
@@ -211,8 +209,6 @@ def from_parameters(cls, name: str, properties: Dict[str, Any], layer: 'Layer' =
         if 'kinetics' in properties:
             constituent.kinetics = Kinetics.from_parameters(properties['kinetics'])
             constituent._initialize_homeostatic_kinetics()
-        else:
-            print(f"        No kinetics (elastin-like constituent)")
         
         # TODO: Set tau_min based on constituent type
 
@@ -246,7 +242,7 @@ def _initialize_deposition_stretch(self, deposition_stretch_input) -> None:
                 f"deposition_stretch must be a scalar, got {type(deposition_stretch_input)}"
             )
         
-        self.deposition_stretch_scalar = float(deposition_stretch_input)
+        deposition_stretch_scalar = float(deposition_stretch_input)
 
         has_fiber = (
             self.constitutive_model is not None and 
@@ -287,23 +283,16 @@ def _initialize_deposition_stretch(self, deposition_stretch_input) -> None:
             
             # Create tensor with ones except fiber direction
             self.deposition_stretch = np.eye(3)
-            self.deposition_stretch[fiber_component, fiber_component] = self.deposition_stretch_scalar
-            
-            print(f"        Deposition stretch: λ_dep = {self.deposition_stretch_scalar:.2f} (fiber-only)")
-            print(f"        Fiber angle: {fiber_angle:.1f}° → {coord_name}")
-            print(f"        G_α = diag({self.deposition_stretch[0,0]:.2f}, "
-                f"{self.deposition_stretch[1,1]:.2f}, {self.deposition_stretch[2,2]:.2f})")
+            self.deposition_stretch[fiber_component, fiber_component] = deposition_stretch_scalar
             
         else:
             # MODE 1: ISOTROPIC 3D DEPOSITION #TODO: rename
             # Incompressibility: λ_r · λ_θ · λ_z = 1
             # If λ_θ = λ_z = λ, then λ_r = 1/λ²
-            lambda_r = 1.0 / (self.deposition_stretch_scalar ** 2)
-            self.deposition_stretch = np.diag([lambda_r, self.deposition_stretch_scalar, self.deposition_stretch_scalar])
-            
-            print(f"        Deposition stretch: λ_dep = {self.deposition_stretch_scalar:.2f} (isotropic 3D)")
-            print(f"        G_α = diag({lambda_r:.4f}, {self.deposition_stretch_scalar:.2f}, {self.deposition_stretch_scalar:.2f})")
+            lambda_r = 1.0 / (deposition_stretch_scalar ** 2)
+            self.deposition_stretch = np.diag([lambda_r, deposition_stretch_scalar, deposition_stretch_scalar])
 
+            
     def _initialize_homeostatic_kinetics(self) -> None:
         """Initialize kinetics histories at t=0 (homeostatic state).
         
@@ -313,9 +302,7 @@ def _initialize_homeostatic_kinetics(self) -> None:
         - Production = Degradation (steady state)
         - mR_h = rhoR_h * k_alpha_h
         - q(0,0) = 1.0 (just deposited material survives 100%)
-        """
-        print(f"        Homeostatic kinetics:")
-               
+        """              
         # Get homeostatic degradation rate from degradation function
         k_alpha_h = self.kinetics.degradation_function.k_alpha_h
         
@@ -331,11 +318,6 @@ def _initialize_homeostatic_kinetics(self) -> None:
         # Initialize survival history
         # At t=0, we have one cohort deposited at τ=0 with q(0,0) = 1.0
         self.survival_history.append([1.0])
-        
-        # Print summary
-        print(f"          k_alpha(0) = {k_alpha_h:.6f} 1/day")
-        print(f"          mR_alpha(0) = {mR_alpha_h:.6f} kg/(m³·day)")
-        print(f"          q(0,0) = 1.0")
 
     def _initialize_active_properties(self, active_props: Dict[str, float]) -> None:
         """Initialize active properties for contractile constituents (e.g., SMCs).
@@ -349,8 +331,6 @@ def _initialize_active_properties(self, active_props: Dict[str, float]) -> None:
                 - CB: Basal calcium concentration (-)
                 - CS: Shear stress sensitivity (optional, default=0)
         """
-        print("        Active properties:")
-        
         self.is_active = True
         
         # Unit conversions
@@ -378,14 +358,6 @@ def _initialize_active_properties(self, active_props: Dict[str, float]) -> None:
         # Active stress at homeostasis (compute using homeostatic conditions)
         sigma_hat_act_h = self._compute_active_stress_at_homeostasis()
         self.sigma_hat_act_history.append(sigma_hat_act_h)
-        
-        # Print summary
-        print(f"          T_act = {active_props['T_act_h']:.1f} kPa")
-        print(f"          k_act = {self.active_properties['k_act']:.6f} 1/day")
-        print(f"          λ_0 = {self.active_properties['lambda_0']:.2f}, λ_m = {self.active_properties['lambda_m']:.2f}")
-        print(f"          CB = {self.active_properties['CB']:.4f}")
-        print(f"          a_act(0) = {a_h*1000:.4f} mm")
-        print(f"          σ_hat_act(0) = {sigma_hat_act_h/1000:.2f} kPa")
 
     # =========================================================================
     # COMPUTATION: MASS DENSITY
@@ -850,6 +822,155 @@ def to_dict(self, timestep: int, time: float = None) -> Dict[str, Any]:
 
         return data
 
+    
+    # =========================================================================
+    # OUTPUT: PRINTING
+    # =========================================================================
+
+    def print_state(self, timestep: int = None, logger: Optional[CustomLogger] = None):
+        """Print constituent state.
+        
+        Args:
+            timestep: Timestep to print (None = homeostatic)
+            logger: Logger to use for printing. 
+                    If None, uses own logger (when called standalone).
+                    If provided, uses that logger (when called within parent context).
+            
+        TODO: To refactor individual print functions after refactoring of classes.
+        """
+        # Use provided logger (for nesting) or own logger (for standalone)
+        if logger is None:
+            logger = self.logger
+            
+        if timestep is None:
+            box_title = f"Constituent: {self.name}"
+        else:
+            box_title = f"Constituent: {self.name} (t={timestep})"
+        
+        with logger.box(box_title, width=50):
+            self._print_mass_fraction(logger, timestep)
+            self._print_deposition_stretch(logger)
+            self._print_constitutive_model(logger)
+            
+            if self.is_active:
+                self._print_active_properties(logger)
+            
+            if self.kinetics is not None:
+                self._print_kinetics(logger)
+
+    def _print_mass_fraction(self, logger: CustomLogger, timestep: int):
+        """Print mass fraction and density."""
+        if timestep is None:
+            mass_fraction = self.homeostatic_referential_density / self.layer.homeostatic_density
+            rho_alpha = self.homeostatic_referential_density
+        else:
+            rho_alpha = self.get_rhoR_alpha(timestep)
+            mass_fraction = rho_alpha / self.layer.get_density(timestep)
+        
+        logger.box_item_aligned("Mass fraction:", f"{mass_fraction:.2f} ({rho_alpha:.2f} kg/m³)", label_width=20)
+
+    def _print_deposition_stretch(self, logger: CustomLogger):
+        # TODO: Refactor to generalize for anisotropic deposition stretch in the future.
+        """Print deposition stretch."""
+        G = self.deposition_stretch
+        lambda_r = G[0, 0]
+        lambda_theta = G[1, 1]
+        lambda_z = G[2, 2]
+        
+        logger.box_item_aligned("Deposition stretch:", f"[{lambda_r:.2f}, {lambda_theta:.2f}, {lambda_z:.2f}]", label_width=20)
+
+    def _print_constitutive_model(self, logger: CustomLogger):
+        # TODO: Refactor to remove hardcoded parts, by automatically getting this info from constitutve_laws.py.
+        """Print constitutive model information."""
+        logger.box_section("Constitutive Model")
+        
+        # Get model type name
+        model_type = self.constitutive_model.__class__.__name__
+        logger.box_item_aligned("Type:", f"{model_type}", label_width=20)
+        
+        # Print parameters based on model type
+        if model_type == "NeoHookeanModel":
+            c_kPa = self.constitutive_model.c / 1000
+            logger.box_item_aligned(f"Parameters:", f"c = {c_kPa:.1f} kPa", label_width=20)
+        
+        elif model_type == "FungExponentialModel":
+            c1_kPa = self.constitutive_model.c1 / 1000
+            c2 = self.constitutive_model.c2
+            logger.box_item_aligned(f"Parameters:", f"c1 = {c1_kPa:.1f} kPa, c2 = {c2:.1f}", label_width=20)
+            
+            # Print fiber orientation if present
+            if hasattr(self.constitutive_model, 'fiber_angle_degrees') and self.constitutive_model.fiber_angle_degrees is not None:
+                angle = self.constitutive_model.fiber_angle_degrees
+                if angle == 90.0:
+                    orientation = "circumferential"
+                elif angle == 0.0:
+                    orientation = "axial"
+                else:
+                    orientation = f"{angle}°"
+                logger.box_item_aligned("Fiber orientation:", f"{angle}° ({orientation})", label_width=20)
+
+    def _print_active_properties(self, logger: CustomLogger):
+        # TODO: Refactor when active stress class has been implemented.
+        """Print active stress properties (for contractile constituents)."""
+        logger.box_section("Active Properties")
+        
+        # Convert to display units
+        T_act_kPa = self.active_properties['T_act'] / 1000
+        k_act = self.active_properties['k_act']
+        lambda_0 = self.active_properties['lambda_0']
+        lambda_m = self.active_properties['lambda_m']
+        CB = self.active_properties['CB']
+        CS = self.active_properties['CS']
+        
+        logger.box_item_aligned(f"T_act_h:", f"{T_act_kPa:.2f} kPa", label_width=20)
+        logger.box_item_aligned(f"k_act:", f"{k_act:.3f} 1/day", label_width=20)
+        logger.box_item_aligned(f"λ_0:", f"{lambda_0:.3f}", label_width=20)
+        logger.box_item_aligned(f"λ_m:", f"{lambda_m:.3f}", label_width=20)
+        logger.box_item_aligned(f"CB:", f"{CB:.3f}", label_width=20)
+        logger.box_item_aligned(f"CS:", f"{CS:.3f}", label_width=20)
+
+    def _print_kinetics(self, logger: CustomLogger):
+        # TODO: Refactor to remove hardcoded parts, by automatically getting this info from kinetics.
+        """Print kinetics information."""        
+        # =========================================================================
+        # PRODUCTION
+        # =========================================================================
+        prod_func = self.kinetics.production_function
+        prod_type = prod_func.__class__.__name__
+
+        logger.box_section("Kinetics - Mass Production")     
+        if prod_type == "LinearProductionRate":
+            gains_str = ", ".join([f"K_{k} = {v}" for k, v in prod_func.gain_params.items()])
+            logger.box_item_aligned("Stimulus function:", f"Linear", label_width=30)
+            logger.box_item_aligned("Gain parameters:", f"{gains_str}", label_width=30)
+        else:
+            # Other production types (future)
+            pass
+    
+        # =========================================================================
+        # DEGRADATION
+        # =========================================================================
+        deg_func = self.kinetics.degradation_function
+        k_alpha_h = deg_func.k_alpha_h
+        deg_type = deg_func.__class__.__name__
+        
+        logger.box_section("Kinetics - Mass Degradation") 
+        if deg_type == "ConstantDegradationRate":
+            logger.box_item_aligned(f"Survival function:", "Exponential", label_width=30)
+            logger.box_item_aligned(f"Degradation rate function:", "Constant", label_width=30)
+            logger.box_item_aligned("Homeostatic deg. rate:", f"k_alpha_h = {k_alpha_h:.4f} 1/day", label_width=30)
+        elif deg_type == "QuadraticDegradationRate":
+            logger.box_item_aligned(f"Survival function:", "Exponential", label_width=30)
+            # Print rate modulation if gains present
+            if deg_func.gain_params:
+                gains_str = ", ".join([f"K_{k} = {v}" for k, v in deg_func.gain_params.items()])
+                logger.box_item_aligned(f"Degradation rate function:", f"Quadratic", label_width=30)
+                logger.box_item_aligned("Homeostatic deg. rate:", f"k_alpha_h = {k_alpha_h:.4f} 1/day", label_width=30)
+                logger.box_item_aligned("Gain parameters:", f"{gains_str}", label_width=30)     
+        else:
+            # Other degradation types (future)
+            pass
+                            
 # =========================================================================
 # MULTI-FIBER FAMILY CONSTITUENT IMPLEMENTATION
 # =========================================================================
@@ -870,9 +991,7 @@ def __init__(self, name):
         
     @classmethod
     def from_parameters(cls, name, properties):
-        """Create and fully initialize multi-fiber family constituent from parameters."""
-        print(f"    Initializing multi-fiber family constituent: {name}")
-        
+        """Create and fully initialize multi-fiber family constituent from parameters."""      
         constituent = cls(name)
         
         # Store shared properties
@@ -946,8 +1065,6 @@ def _initialize_total_kinetics_histories(self) -> None:
             print(f"        No kinetics for multi-fiber constituent")
             return
         
-        print(f"        Total homeostatic kinetics:")
-        
         # Sum kinetics values from all fiber families at t=0
         total_k_alpha_h = sum(
             family.k_alpha_history[0] for family in self.fiber_families
@@ -962,10 +1079,6 @@ def _initialize_total_kinetics_histories(self) -> None:
         
         # Survival history for multi-fiber is aggregate (same as fiber families)
         self.survival_history.append([1.0])
-        
-        print(f"          k_alpha(0) = {total_k_alpha_h:.6f} 1/day")
-        print(f"          mR_alpha(0) = {total_mR_alpha_h:.6f} kg/(m³·day)")
-        print(f"          q(0,0) = 1.0")
 
     # =========================================================================
     # COMPUTATION: MASS DENSITY
diff --git a/src/constitutive_laws.py b/src/constitutive_laws.py
index 2ae32dc..cd443f6 100644
--- a/src/constitutive_laws.py
+++ b/src/constitutive_laws.py
@@ -107,9 +107,6 @@ def set_fiber_orientation(self, angle_degrees):
         
         # Store as unit vector
         self.fiber_orientation = a0 / np.linalg.norm(a0)
-        
-        print(f"        Fiber orientation: {angle_degrees}° → "
-              f"a0 = [{a0[0]:.3f}, {a0[1]:.3f}, {a0[2]:.3f}]")
     
     def is_isotropic(self):
         """Check if material is isotropic."""
diff --git a/src/custom_logging.py b/src/custom_logging.py
index 03bd858..93b5094 100644
--- a/src/custom_logging.py
+++ b/src/custom_logging.py
@@ -12,24 +12,211 @@
 from typing import Optional
 
 
-class IndentedLogger:
-    """Wrapper around Python logging with indentation support.
+class Box:
+    """Handles box-drawing formatting for structured output.
+    
+    Box structure visualization:
+        ┌─ Title ─────────────────────  ← box_start()
+        │  Content line                  ← format_line(spacing=2)
+        │    Indented content            ← format_line(spacing=4)
+        │    Label:      Value           ← format_line(spacing=4, aligned)
+        └─────────────────────────────  ← box_end()
+    """
+    
+    # Box drawing characters
+    VERTICAL = "│"
+    HORIZONTAL = "─"
+    TOP_LEFT = "┌"
+    BOTTOM_LEFT = "└"
+    
+    def __init__(self, indent: str = "", outer_box_prefix: str = ""):
+        """Initialize box formatter.
+        
+        Args:
+            indent: Base indentation string (e.g., "" for layer, "  " for constituent)
+            outer_box_prefix: Prefix from outer boxes (e.g., "│  " for nested boxes)
+        """
+        self.indent = indent
+        self.outer_box_prefix = outer_box_prefix
+    
+    def box_start(self, title: str, width: int = 60) -> str:
+        """Format box opening with title.
+        
+        Args:
+            title: Box title text
+            width: Total width of the box border
+            
+        Returns:
+            Formatted box header string
+            
+        Visual output:
+            ┌─ Title ─────────────────────────────────────
+        """
+        # Calculate separator to match box_end width exactly
+        # Total: outer_prefix + indent + "┌─ " + title + " " + separator
+        # We want total line length = outer_prefix + indent + width
+        used_chars = len(f"┌─ {title} ")
+        separator_length = width - used_chars
+        separator = self.HORIZONTAL * separator_length
+        
+        return f"{self.outer_box_prefix}{self.indent}{self.TOP_LEFT}{self.HORIZONTAL} {title} {separator}"
+    
+    def box_end(self, width: int = 60) -> str:
+        """Format box closing.
+        
+        Args:
+            width: Total width of the box border
+            
+        Returns:
+            Formatted box footer string
+            
+        Visual output:
+            └─────────────────────────────────────────────
+        """
+        # Match box_start width exactly
+        separator = self.HORIZONTAL * (width - 1)  # -1 for the └ character
+        return f"{self.outer_box_prefix}{self.indent}{self.BOTTOM_LEFT}{separator}"
+    
+    def format_line(self, text: str = "", spacing: int = 2, label: str = None, label_width: int = 12) -> str:
+        """Format a line within the box.
+        
+        Args:
+            text: Content to display (or value if label is provided)
+            spacing: Number of spaces after vertical bar (default: 2)
+            label: Optional label for aligned key-value pairs
+            label_width: Width for label column when using aligned mode
+            
+        Returns:
+            Formatted line string
+            
+        Visual output examples:
+            format_line("Content", spacing=2)
+            → │  Content
+            
+            format_line("Item", spacing=4)
+            → │    Item
+            
+            format_line("Value", spacing=4, label="Key:", label_width=12)
+            → │    Key:         Value
+            
+            format_line("", spacing=0)
+            → │
+        """
+        if not text and label is None:
+            # Blank line
+            return f"{self.outer_box_prefix}{self.indent}{self.VERTICAL}"
+        
+        spaces = " " * spacing
+        
+        if label is not None:
+            # Aligned key-value pair
+            return f"{self.outer_box_prefix}{self.indent}{self.VERTICAL}{spaces}{label:<{label_width}} {text}"
+        else:
+            # Plain content
+            return f"{self.outer_box_prefix}{self.indent}{self.VERTICAL}{spaces}{text}"
+
+
+class BoxContext:
+    """Context manager for box-drawing mode.
+    
+    Handles box lifecycle: creates box, prints header, manages indentation,
+    prints footer, and restores previous box state.
+    
+    Usage:
+        with logger.box("Layer: Cerebral artery"):
+            logger.box_line("Kinematics: ThinWall")
+    """
+    
+    def __init__(self, logger: 'CustomLogger', title: str, width: int = 60):
+        """Initialize box context.
+        
+        Args:
+            logger: Parent CustomLogger instance
+            title: Title for box header
+            width: Width of box border
+        """
+        self.logger = logger
+        self.title = title
+        self.width = width
+        self.previous_box = None
+    
+    def __enter__(self):
+        """Enter box context - print header and create box.
+        
+        Creates new box at current indentation level, prints header,
+        saves previous box, and increases indentation for box contents.
+        """
+        # Create new box with current indentation
+        indent = '  ' * self.logger.indent_level
+        
+        # Build outer box prefix (continuation of parent boxes)
+        outer_prefix = ""
+        if self.logger.indent_level > 0:
+            # We're inside another box - need to continue its vertical line
+            indent = "" # Don't add extra indent for nested box - it will be handled by outer box prefix
+            parent_box = self.logger.current_box
+            outer_prefix = f"{parent_box.outer_box_prefix}{parent_box.indent}{Box.VERTICAL}  "
+        
+        # Create new box with current indentation and outer prefix
+        new_box = Box(indent=indent, outer_box_prefix=outer_prefix)
+        
+        # Print box opening
+        header = new_box.box_start(self.title, self.width)
+        self.logger.logger.info(header)
+        
+        # Save previous box and set new one
+        self.previous_box = self.logger.current_box
+        self.logger.current_box = new_box
+        
+        # Increase indent for content inside box
+        self.logger.indent()
+        
+        return self.logger
+    
+    def __exit__(self, exc_type, exc_val, exc_tb):
+        """Exit box context - print footer and restore previous box.
+        
+        Decreases indentation, prints box footer, and restores the
+        previous box state.
+        """
+        # Decrease indent
+        self.logger.dedent()
+        
+        # Print box closing
+        footer = self.logger.current_box.box_end(self.width)
+        self.logger.logger.info(footer)
+        
+        # Restore previous box
+        self.logger.current_box = self.previous_box
+        
+        return False
+
+
+class CustomLogger:
+    """Wrapper around Python logging with indentation and box-drawing support.
     
     Provides:
     - Indentation support (indent/dedent methods)
     - Context manager for automatic indent/dedent
+    - Box-drawing utilities for structured output
     - Emoji prefixes for errors/warnings
     - Section headers
     
     Usage:
         logger = get_logger(__name__)
         
+        # Standard logging with auto-indentation
         logger.info("Starting computation")
-        with logger:  # Auto-indent
+        with logger:
             logger.debug("Step 1")
             logger.debug("Step 2")
-        # Auto-dedent
         logger.info("Complete")
+        
+        # Box-drawing for structured output
+        with logger.box("Layer: Arterial Wall"):
+            logger.box_line("Kinematics: ThinWall")
+            logger.box_section("Geometry")
+            logger.box_item("Inner radius: 1.4 mm")
     """
     
     def __init__(self, name: str):
@@ -41,9 +228,97 @@ def __init__(self, name: str):
         #TODO: consider adding log_level here?
         self.logger = logging.getLogger(name)
         self.indent_level = 0
+        self.current_box = Box(indent="")  # Default box for when not in box context
     
     # =========================================================================
-    # LOGGING METHODS
+    # BOX CONTEXT MANAGER
+    # =========================================================================
+    
+    def box(self, title: str, width: int = 60):
+        """Create box context manager for structured output.
+        
+        Args:
+            title: Title for box header
+            width: Width of box border (default: 60)
+            
+        Returns:
+            BoxContext for use with 'with' statement
+            
+        Example:
+            with logger.box("Layer: Cerebral artery"):
+                logger.box_line("Kinematics: ThinWall")
+                logger.box_section("Homeostatic Geometry")
+                logger.box_item("Inner radius: 1.4 mm")
+        """
+        return BoxContext(self, title, width)
+    
+    # =========================================================================
+    # BOX DRAWING METHODS
+    # =========================================================================
+    
+    def box_line(self, text: str = ""):
+        """Print content line within current box.
+        
+        Args:
+            text: Content to display (empty for blank line)
+            
+        Visual output:
+            │  Text content here
+            │  (blank if text is empty)
+        """
+        line = self.current_box.format_line(text, spacing=2)
+        self.logger.info(line)
+    
+    def box_section(self, title: str):
+        """Print section header within current box.
+        
+        Adds blank line before section title for visual separation.
+        
+        Args:
+            title: Section title
+            
+        Visual output:
+            │
+            │  Section Title
+        """
+        self.logger.info(self.current_box.format_line("", spacing=0))  # Blank line
+        self.logger.info(self.current_box.format_line(title, spacing=2))  # Section title
+    
+    def box_item(self, text: str, extra_indent: int = 2):
+        """Print indented item within a section.
+        
+        Args:
+            text: Item content
+            extra_indent: Additional indentation spaces (default: 2)
+            
+        Visual output (extra_indent=2):
+            │    Item content here
+        """
+        item = self.current_box.format_line(text, spacing=2 + extra_indent)
+        self.logger.info(item)
+    
+    def box_item_aligned(self, label: str, value: str, label_width: int = 12, extra_indent: int = 2):
+        """Print aligned key-value pair within box.
+        
+        Args:
+            label: Left-aligned label
+            value: Value to display
+            label_width: Width for label column (default: 12)
+            extra_indent: Additional indentation spaces (default: 2)
+            
+        Visual output (label_width=12, extra_indent=2):
+            │    Label:       Value here
+        """
+        item = self.current_box.format_line(
+            text=value, 
+            spacing=2 + extra_indent, 
+            label=label, 
+            label_width=label_width
+        )
+        self.logger.info(item)
+    
+    # =========================================================================
+    # STANDARD LOGGING METHODS
     # =========================================================================
     
     def debug(self, msg: str):
@@ -91,10 +366,16 @@ def critical(self, msg: str):
     # =========================================================================
     
     def section(self, title: str, log_level: str = 'INFO'):
-        """Print section header with indentation.
+        """Print section header with separator lines.
         
         Args:
             title: Section title
+            log_level: Logging level (default: 'INFO')
+            
+        Visual output:
+            ============================================================
+            Section Title
+            ============================================================
         """
         indent = '  ' * self.indent_level
         separator = '=' * 60
@@ -105,17 +386,17 @@ def section(self, title: str, log_level: str = 'INFO'):
             print(f"\n{indent}{separator}")
             print(f"{indent}{title}")
             print(f"{indent}{separator}")
-        
+    
     # =========================================================================
     # INDENTATION CONTROL
     # =========================================================================
     
     def indent(self):
-        """Increase indentation level."""
+        """Increase indentation level by one step (2 spaces)."""
         self.indent_level += 1
     
     def dedent(self):
-        """Decrease indentation level."""
+        """Decrease indentation level by one step (2 spaces)."""
         self.indent_level = max(0, self.indent_level - 1)
 
     def _format_msg(self, msg: str) -> str:
@@ -129,9 +410,9 @@ def _format_msg(self, msg: str) -> str:
         """
         indent = '  ' * self.indent_level
         return f"{indent}{msg}"
-        
+    
     # =========================================================================
-    # CONTEXT MANAGER
+    # CONTEXT MANAGER (for simple indentation)
     # =========================================================================
     
     def __enter__(self):
@@ -219,18 +500,18 @@ def init_logging(log_level: str, config_file: str = 'logging_config.yaml') -> No
         )
     
     # Set the root logger level (affects all loggers)
-    numeric_level = getattr(logging, log_level_upper)
+    numeric_level = getattr(logging, log_level.upper())
     logging.getLogger().setLevel(numeric_level)
 
 
-def get_logger(name: str) -> IndentedLogger:
-    """Get enhanced logger with indentation support.
+def get_logger(name: str) -> CustomLogger:
+    """Get enhanced logger with indentation and box-drawing support.
     
     Args:
         name: Logger name (typically __name__ from calling module)
         
     Returns:
-        IndentedLogger instance
+        CustomLogger instance
         
     Example:
         >>> logger = get_logger(__name__)
@@ -238,5 +519,8 @@ def get_logger(name: str) -> IndentedLogger:
         >>> with logger:
         ...     logger.debug("Detailed step")
         >>> logger.info("Process complete")
+        
+        >>> with logger.box("Configuration"):
+        ...     logger.box_line("Parameter: value")
     """
-    return IndentedLogger(name)
\ No newline at end of file
+    return CustomLogger(name)
\ No newline at end of file
diff --git a/src/deformation_kinematics.py b/src/deformation_kinematics.py
index 5c59382..9d8d24d 100644
--- a/src/deformation_kinematics.py
+++ b/src/deformation_kinematics.py
@@ -606,8 +606,8 @@ def compute_wss(
             raise ValueError(f"Blood viscosity must be positive, got {blood_viscosity}")
         
         # Poiseuille WSS formula
-        # TODO: unlcear why blood viscosity not used in Latorre 2018. Rethink units.
-        # wss = (4.0 * blood_viscosity * flow_rate) / (math.pi * inner_radius**3)
+        # TODO: unclear why blood viscosity not used in Latorre 2018.
+        #wss = (4.0 * blood_viscosity * flow_rate) / (math.pi * inner_radius**3)
         wss = flow_rate / inner_radius**3
 
         return wss
diff --git a/src/fixed_point_solver.py b/src/fixed_point_solver.py
index 7f3d3c2..aa56eae 100644
--- a/src/fixed_point_solver.py
+++ b/src/fixed_point_solver.py
@@ -292,7 +292,7 @@ def solve(
             if all_converged:
                 converged = True
                 self.logger.info(
-                    f"✓ Fixed-point converged in {iteration} iterations "
+                    f"✓ Fixed-point solver converged in {iteration} iterations "
                     f"(max Δρ/ρ = {max_change:.2e})"
                 )
                 break
diff --git a/src/layer.py b/src/layer.py
index 98cc4e0..84b14f5 100644
--- a/src/layer.py
+++ b/src/layer.py
@@ -1,4 +1,6 @@
+from turtle import title
 import numpy as np
+from custom_logging import get_logger, CustomLogger
 from typing import List, Dict, Any, Optional
 from constituent import Constituent
 from deformation_kinematics import DeformationKinematics, KinematicsFactory
@@ -16,6 +18,7 @@ def __init__(self, name):
         self.kinematics: Optional[DeformationKinematics] = None
         self.mechanics = Mechanics()
         self.perturbation_manager = PerturbationManager(scope='layer')
+        self.logger: CustomLogger = get_logger(__name__)
         
         # Homeostatic geometry (= reference configuration in G&R)
         self.homeostatic_inner_radius: Optional[float] = None
@@ -50,6 +53,8 @@ def __init__(self, name):
         self.flow_rate_history: List[float] = []
         self.wss_history: List[float] = []
 
+
+
     # =========================================================================
     # INITIALIZATION
     # =========================================================================
@@ -71,10 +76,8 @@ def from_parameters(cls, layer_data: Dict[str, Any]) -> 'Layer':
             
         Returns:
             Fully initialized Layer instance
-        """
-        print(f"\nInitializing layer: {layer_data['layer_name']}")
-        
-        # Step 1: Create empty layer (just name)
+        """      
+        # Step 1: Create empty layer
         layer = cls(name=layer_data['layer_name'])
 
         # STEP 2: Set kinematics from geometry type
@@ -95,9 +98,7 @@ def from_parameters(cls, layer_data: Dict[str, Any]) -> 'Layer':
         )
 
         # Step 4: Initialize constituents for this layer
-        constituents = layer_data['constituents']
-        print(f"  Constituents ({len(constituents)} total):")
-        
+        constituents = layer_data['constituents']        
         for name, properties in constituents.items():
             constituent = Constituent.from_parameters(name, properties, layer)
             layer.add_constituent(constituent)
@@ -110,9 +111,9 @@ def from_parameters(cls, layer_data: Dict[str, Any]) -> 'Layer':
     
     def set_kinematics(self, geometry_type: str) -> None:
         """Set kinematics for this layer."""
+        # TODO: validation should happen at the inputfile leve, not when setting kinematics.
         try:
             self.kinematics = KinematicsFactory.create(geometry_type)
-            print(f"    Kinematics set: {type(self.kinematics).__name__}")
         except ValueError as e:
             available_types = KinematicsFactory.available_types()
             raise ValueError(
@@ -145,9 +146,7 @@ def _initialize_homeostatic_geometry(
             rhoR_h: Homeostatic mass density (kg/m³)
         """
         self._ensure_kinematics_set()
-        
-        print("  Homeostatic geometry:")
-        
+    
         # Unit conversions
         mm_to_m = 1.0e-3
         kPa_to_Pa = 1.0e3
@@ -197,19 +196,6 @@ def _initialize_homeostatic_geometry(
         self.flow_rate_history = [self.homeostatic_flow_rate]
         self.wss_history = [self.homeostatic_wss]
 
-
-        # Print summary
-        print(f"    Inner radius (a_h): {a_h} mm → {self.homeostatic_inner_radius:.6f} m")
-        print(f"    Thickness (h_h): {h_h} mm → {self.homeostatic_thickness:.6f} m")
-        print(f"    Axial stretch (λ_z_h): {self.homeostatic_axial_stretch}")
-        print(f"    Reference density (ρ_h): {self.homeostatic_density} kg/m³")
-        print(f"    Loading:")
-        print(f"      P_h = {pressure_h} kPa → {self.homeostatic_pressure:.2f} Pa")
-        print(f"      Q_h = {flow_rate_h} m³/day")
-        print(f"      μ = {self.blood_viscosity*1000:.2f} mPa·s")        
-        print(f"    Wall shear stress (computed):")
-        print(f"      WSS_h = {self.homeostatic_wss:.4f} Pa")
-
     def _initialize_perturbations(self, perturbations_config: Dict[str, Any]):
         """Initialize perturbations from configuration.
         
@@ -220,8 +206,6 @@ def _initialize_perturbations(self, perturbations_config: Dict[str, Any]):
         for variable_name, pert_config in perturbations_config.items():
             perturbation = PerturbationFactory.create(pert_config)
             self.perturbation_manager.add_perturbation(variable_name, perturbation)
-
-            print(f"    Registered perturbation for '{variable_name}': {perturbation}")
        
     def _ensure_kinematics_set(self) -> None:
         """Ensure kinematics has been set before use."""
@@ -320,11 +304,10 @@ def get_homeostatic_geometry(self) -> Dict[str, float]:
         
         return {
             'inner_radius': self.homeostatic_inner_radius,
-            'thickness': self.homeostatic_thickness,
             'mid_radius': self.homeostatic_mid_radius,
             'outer_radius': self.homeostatic_outer_radius,
+            'thickness': self.homeostatic_thickness,
             'axial_stretch': self.homeostatic_axial_stretch,
-            'density': self.homeostatic_density
         }
 
     # ========================================================================
@@ -383,12 +366,10 @@ def get_geometry(self, timestep: int) -> Dict[str, float]:
         """Return layer geometry at specified timestep."""
         return {
             'inner_radius': self.get_inner_radius(timestep),
-            'thickness': self.get_thickness(timestep),
             'mid_radius': self.get_mid_radius(timestep),
             'outer_radius': self.get_outer_radius(timestep),
-            'axial_stretch': self.get_axial_stretch(timestep),
-            'density': self.get_density(timestep),
-            'stress': self.get_stress(timestep)
+            'thickness': self.get_thickness(timestep),
+            'axial_stretch': self.get_axial_stretch(timestep)
         }
 
     # ========================================================================
@@ -633,53 +614,39 @@ def compute_homeostatic_stress_direct(self) -> None:
 
         This avoids integration issues at t=0.
         """
-        print(f"  Computing homeostatic stress for layer '{self.name}' (direct method)")
-
         sigma_total = np.zeros((3, 3))
 
-        for constituent in self.constituents:
-            print(f"    {constituent.name}:")
-        
+        for constituent in self.constituents:        
             # Get constituent properties
             rho_alpha = constituent.homeostatic_referential_density
             rho_h = self.homeostatic_density
             G_alpha = constituent.deposition_stretch
         
-            print(f"      ρ_α = {rho_alpha:.2f} kg/m³, ρ_h = {rho_h:.2f} kg/m³")
-            print(f"      Mass fraction = {rho_alpha/rho_h:.4f}")
-            print(f"      G_α diagonal = [{G_alpha[0,0]:.3f}, {G_alpha[1,1]:.3f}, {G_alpha[2,2]:.3f}]")
-        
             # Compute F_α(0,0) = G_α
             F_alpha = G_alpha
-            print(f"      F_α diagonal = [{F_alpha[0,0]:.3f}, {F_alpha[1,1]:.3f}, {F_alpha[2,2]:.3f}]")
         
             # Compute C_α
             C_alpha = F_alpha.T @ F_alpha
-            print(f"      C_α diagonal = [{C_alpha[0,0]:.3f}, {C_alpha[1,1]:.3f}, {C_alpha[2,2]:.3f}]")
         
             # Get constitutive model
             model = constituent.constitutive_model
         
             # ✅ NEW: Compute S_hat_α from constitutive law
             S_hat_alpha = model.compute_PK2_stress(F_alpha)
-            print(f"      S_hat_α (PK2 stress) diagonal = [{S_hat_alpha[0,0]:.2f}, {S_hat_alpha[1,1]:.2f}, {S_hat_alpha[2,2]:.2f}] Pa")
         
             # ✅ NEW: Compute σ̂_α (push-forward)
             # TODO: compute from the ratio of J from layer. Currently hardcoded to 1 for testing.
             # J_alpha = np.linalg.det(F_alpha)
             J_alpha = 1
-            print(f"      J_α = {J_alpha:.4f}")
             
             sigma_hat_alpha = (1.0 / J_alpha) * (F_alpha @ S_hat_alpha @ F_alpha.T)
-            print(f"      σ̂_α (Cauchy stress) diagonal = [{sigma_hat_alpha[0,0]:.2f}, {sigma_hat_alpha[1,1]:.2f}, {sigma_hat_alpha[2,2]:.2f}] Pa")
 
             # Store sigma_hat_alpha at homeostasis for use in heredity integrals
             constituent.sigma_hat_history.append([sigma_hat_alpha])
 
             # Weight by mass fraction (PASSIVE STRESS)
             sigma_alpha_passive = (rho_alpha / rho_h) * sigma_hat_alpha
-            print(f"      σ_α (passive, weighted) diagonal = [{sigma_alpha_passive[0,0]/1000:.2f}, {sigma_alpha_passive[1,1]/1000:.2f}, {sigma_alpha_passive[2,2]/1000:.2f}] kPa")
-            
+ 
             # Start with passive stress
             sigma_alpha = sigma_alpha_passive.copy()
         
@@ -690,9 +657,6 @@ def compute_homeostatic_stress_direct(self) -> None:
                 # Weight by mass fraction
                 sigma_act_weighted = (rho_alpha / rho_h) * sigma_act
                 
-                print(f"      σ_act(0) = {sigma_act/1000:.2f} kPa (unweighted)")
-                print(f"      σ_act(0, weighted) = {sigma_act_weighted/1000:.2f} kPa")
-                
                 # Add to circumferential component
                 sigma_alpha[1, 1] += sigma_act_weighted
 
@@ -701,30 +665,12 @@ def compute_homeostatic_stress_direct(self) -> None:
 
             # Accumulate
             sigma_total += sigma_alpha
-        
-            print(f"      σ_α(0) TOTAL diagonal = [{sigma_alpha[0,0]/1000:.2f}, "
-                f"{sigma_alpha[1,1]/1000:.2f}, {sigma_alpha[2,2]/1000:.2f}] kPa")
-
-        # ✅ NEW: Print total before constraint
-        print(f"\n    Total stress BEFORE constraint:")
-        print(f"      σ_total diagonal = [{sigma_total[0,0]/1000:.2f}, "
-            f"{sigma_total[1,1]/1000:.2f}, {sigma_total[2,2]/1000:.2f}] kPa")
 
         # Apply kinematic constraint
         lagrange = self.kinematics.compute_lagrange_multiplier(sigma_total)
         
-        # ✅ NEW: Print Lagrange multiplier
-        print(f"    Lagrange multiplier:")
-        print(f"      λ diagonal = [{lagrange[0,0]/1000:.2f}, "
-            f"{lagrange[1,1]/1000:.2f}, {lagrange[2,2]/1000:.2f}] kPa")
-        
         sigma_constrained = sigma_total - lagrange
 
-        # ✅ NEW: Print total after constraint
-        print(f"    Total stress AFTER constraint:")
-        print(f"      σ_constrained diagonal = [{sigma_constrained[0,0]/1000:.2f}, "
-            f"{sigma_constrained[1,1]/1000:.2f}, {sigma_constrained[2,2]/1000:.2f}] kPa")
-
         # Store
         self.homeostatic_stress = sigma_constrained
         self.stress_history[0] = sigma_constrained.copy()
@@ -1159,6 +1105,138 @@ def to_dict(self, timestep: int, time: float) -> Dict[str, Any]:
             'wss': self.get_wss(timestep)
         }
 
+    # =========================================================================
+    # OUTPUT: PRINT SUMMARY
+    # =========================================================================
+
+    def print_state(self, timestep: int = None, logger: Optional[CustomLogger] = None):
+        """Print layer state.
+        
+        Args:
+            timestep: Timestep to print. If None, prints homeostatic state.
+            logger: Logger to use for printing.
+                If None, uses own logger (standalone mode).
+                If provided, uses that logger (nested mode).
+            
+        Examples:
+            layer.print_state()                    # Homeostatic, standalone
+            layer.print_state(timestep=10)         # At timestep 10, standalone
+            layer.print_state(timestep=10, logger=config.logger)  # Nested in config
+        """
+        # Use provided logger (for nesting) or own logger (for standalone)
+        if logger is None:
+            logger = self.logger
+    
+        if timestep is None:
+            box_title = f"Layer: {self.name}"
+        else:
+            box_title = f"Layer: {self.name} (timestep={timestep})"
+    
+        with logger.box(box_title):
+            self._print_content(timestep)
+    
+    def _print_content(self, timestep: int):
+        """Print layer content (shared by standalone and nested modes)."""
+        self.logger.box_line(f"Kinematics: {self.kinematics.__class__.__name__}")
+        self.logger.box_line(f"Constituents: {len(self.constituents)}")
+
+        self._print_mass(timestep)       
+        self._print_geometry(timestep)
+        self._print_loading(timestep)
+        self._print_stress(timestep)
+        
+        for constituent in self.constituents:
+            self.logger.box_line()
+            constituent.print_state(timestep=timestep, logger=self.logger)
+            
+
+    def _print_geometry(self, timestep: int):
+        """Print geometry section using kinematics metadata.
+        
+        TODO: Refactor kinematics to use registry pattern for variable metadata.
+        For now, this is a simplified implementation that will be replaced when
+        deformation_kinematics.py is restructured to define component symbols,
+        variable names, units, etc. in a systematic registry.
+        """
+        # Determine title
+        if timestep is None:
+            box_title = "Homeostatic Geometry"
+            geometry = self.get_homeostatic_geometry()
+        else:
+            box_title = f"Geometry (t={timestep})"
+            geometry = self.get_geometry(timestep)
+        
+        self.logger.box_section(box_title)
+        
+        # TODO: Replace with kinematics.get_geometry_display_info() once registry is implemented
+        # For now, hardcode thin-wall variables
+        
+        # Inner radius
+        a_mm = geometry['inner_radius'] * 1e3
+        h_mm = geometry['thickness'] * 1e3
+        lambda_z = geometry['axial_stretch']
+
+        if timestep is None:
+            self.logger.box_item_aligned("Inner radius (a_h):", f"{a_mm:.4f} mm", label_width=25)
+            self.logger.box_item_aligned("Thickness (h_h):", f"{h_mm:.4f} mm", label_width=25)
+            self.logger.box_item_aligned("Axial stretch (λ_z_h):", f"{lambda_z:.4f}", label_width=25)
+        else:
+            self.logger.box_item_aligned("Inner radius (a):", f"{a_mm:.4f} mm", label_width=25)
+            self.logger.box_item_aligned("Thickness (h):", f"{h_mm:.4f} mm", label_width=25)
+            self.logger.box_item_aligned("Axial stretch (λ_z):", f"{lambda_z:.4f}", label_width=25)
+            
+    def _print_loading(self, timestep: int):
+        """Print loading variables."""
+        # TODO: refactor
+        if timestep is None:
+            P_kPa = self.homeostatic_pressure * 1e-3
+            Q = self.homeostatic_flow_rate
+            wss_kPa = self.homeostatic_wss * 1e-3
+            box_title = "Homeostatic Loading"
+            suffix = "_h"
+        else:
+            P_kPa = self.get_pressure(timestep) * 1e-3
+            Q = self.get_flow_rate(timestep)
+            wss_kPa = self.get_wss(timestep) * 1e-3
+            box_title = f"Loading (timestep={timestep})"
+            suffix = ""
+        
+        self.logger.box_section(box_title)
+        self.logger.box_item_aligned(f"Pressure (P{suffix}):", f"{P_kPa:.2f} kPa", label_width=25)
+        self.logger.box_item_aligned(f"Flow rate (Q{suffix}):", f"{Q:.2f} m³/day", label_width=25)
+        # TODO: check wss computation and units (currently computing shear rate for Lattore2018)
+        #self.logger.box_item_aligned(f"WSS (τ_w{suffix}):", f"{wss_kPa:.2f} kPa", label_width=25)
+
+    def _print_mass(self, timestep: int):
+        """Print mass density."""
+        if timestep is None:
+            rho = self.homeostatic_density
+            box_title = "Homeostatic Mass Density"
+            suffix = "_h"
+        else:
+            rho = self.get_density(timestep)
+            box_title = f"Mass Density (timestep={timestep})"
+            suffix = ""
+        
+        self.logger.box_section(box_title)
+        self.logger.box_item_aligned(f"Referential mass density (ρ{suffix}):", f"{rho:.2f} kg/m³", label_width=25)
+        
+    def _print_stress(self, timestep: int):
+        """Print stress tensor."""
+        #TODO: Refactor to use kinematics.get_component_symbol(i) once registry is implemented
+        if timestep is None:
+            sigma = self.homeostatic_stress
+            box_title = "Homeostatic Intramural Stress (diagonal)"
+        else:
+            sigma = self.get_stress(timestep)
+            box_title = f"Intramural Stress (t={timestep}, diagonal)"
+        
+        self.logger.box_section(box_title)
+        for i in range(3):
+            component_symbol = ["rr", "θθ", "zz"][i]
+            stress_kPa = sigma[i, i] / 1000
+            self.logger.box_item_aligned(f"σ_{component_symbol}:", f"{stress_kPa:.2f} kPa", label_width=25)
+
     # def compute_homeostatic_stress_direct(self) -> None:
     #     """Compute homeostatic stress directly from constituent properties.
     
diff --git a/src/logging_config.yaml b/src/logging_config.yaml
index ae534ab..894a5b7 100644
--- a/src/logging_config.yaml
+++ b/src/logging_config.yaml
@@ -11,8 +11,8 @@ logging:
       format: '%(message)s'
     
     # With module name: [module_name] message
-    with_module:
-      format: '[%(name)s] %(message)s'
+    with_log_level:
+      format: '[%(levelname)s] %(message)s'
     
     # Full details: [LEVEL] module_name: message
     detailed:
@@ -26,7 +26,7 @@ logging:
     console:
       class: logging.StreamHandler
       level: DEBUG                    # Minimum level for this handler
-      formatter: detailed             # Use detailed format
+      formatter: with_log_level           
       stream: ext://sys.stdout        # Print to stdout
     
     # File handler - writes to file
@@ -34,7 +34,7 @@ logging:
       class: logging.FileHandler
       level: DEBUG                    # Log everything to file
       filename: simulation_debug.log  # Output file
-      formatter: detailed             # Use detailed format for file
+      formatter: with_log_level       
       mode: w                         # 'w' = overwrite, 'a' = append
 
   # =========================================================================
diff --git a/src/simulation.py b/src/simulation.py
index 944e11c..7365b6d 100644
--- a/src/simulation.py
+++ b/src/simulation.py
@@ -81,7 +81,7 @@ def run(self):
 
         self._setup_outputs()
 
-        self.logger.info("Starting G&R simulation")
+        self.logger.section("Starting G&R simulation")
         self.logger.info(f"Steps: {self.n_steps}, dt: {self.dt}")
 
         # Write initial state (t=0, homeostatic)