diff --git a/ExampleScript.py b/ExampleScript.py
index c0ea6d13..7b0bdd8c 100644
--- a/ExampleScript.py
+++ b/ExampleScript.py
@@ -1,36 +1,99 @@
-import parametric_stellarator
+import parastell
 import logging
 
 
 # Define plasma equilibrium VMEC file
 plas_eq = 'plas_eq.nc'
-# Define number of periods in stellarator plasma
-num_periods = 4
 # Define radial build
-radial_build = {
-    'sol': {'thickness': 10, 'h5m_tag': 'Vacuum'},
-    'first_wall': {'thickness': 5},
-    'blanket': {'thickness': 5},
-    'back_wall': {'thickness': 5},
-    'shield': {'thickness': 20},
-    'coolant_manifolds': {'thickness': 5},
-    'gap': {'thickness': 5, 'h5m_tag': 'Vacuum'},
-    # Note that some neutron transport codes (such as OpenMC) will interpret
-    # materials with "vacuum" in the name as void material
-    'vacuum_vessel': {'thickness': 20, 'h5m_tag': 'vv'}
+build = {
+    'phi_list': [0.0, 22.5, 45.0, 67.5, 90.0],
+    'theta_list': [0.0, 5.0, 90.0, 175.0, 180.0, 185.0, 270.0, 355.0, 360.0],
+    'wall_s': 1.2,
+    'radial_build': {
+        'first_wall': {
+            'thickness_matrix': [
+                [5, 5, 5, 5, 5, 5, 5, 5, 5],
+                [5, 5, 5, 5, 5, 5, 5, 5, 5],
+                [5, 5, 5, 5, 5, 5, 5, 5, 5],
+                [5, 5, 5, 5, 5, 5, 5, 5, 5],
+                [5, 5, 5, 5, 5, 5, 5, 5, 5]
+            ]
+        },
+        'breeder': {
+            'thickness_matrix': [
+                [100, 100, 30, 10, 10, 10, 30, 100, 100],
+                [30,  30,  10, 5,  5,  5,  20, 30,  30],
+                [25,  25,  25, 5,  5,  5,  25, 25,  25],
+                [30,  30,  20, 5,  5,  5,  10, 30,  30],
+                [100, 100, 30, 10, 10, 10, 30, 100, 100]
+            ]
+        },
+        'back_wall': {
+            'thickness_matrix': [
+                [5, 5, 5, 5, 5, 5, 5, 5, 5],
+                [5, 5, 5, 5, 5, 5, 5, 5, 5],
+                [5, 5, 5, 5, 5, 5, 5, 5, 5],
+                [5, 5, 5, 5, 5, 5, 5, 5, 5],
+                [5, 5, 5, 5, 5, 5, 5, 5, 5]
+            ]
+        },
+        'shield': {
+            'thickness_matrix': [
+                [25, 25, 25, 25, 25, 25, 25, 25, 25],
+                [25, 25, 25, 25, 25, 25, 25, 25, 25],
+                [25, 25, 25, 25, 25, 25, 25, 25, 25],
+                [25, 25, 25, 25, 25, 25, 25, 25, 25],
+                [25, 25, 25, 25, 25, 25, 25, 25, 25]
+            ]
+        },
+        'manifolds': {
+            'thickness_matrix': [
+                [50, 50, 15, 5,  5,  5,  15, 50, 50],
+                [20, 20, 5,  5,  5,  5,  15, 20, 20],
+                [15, 15, 15, 5,  5,  5,  15, 15, 15],
+                [20, 20, 15, 5,  5,  5,  5,  20, 20],
+                [50, 50, 15, 5,  5,  5,  15, 50, 50]
+            ]
+        },
+        'gap': {
+            'thickness_matrix': [
+                [5, 5, 5, 5, 5, 5, 5, 5, 5],
+                [5, 5, 5, 5, 5, 5, 5, 5, 5],
+                [5, 5, 5, 5, 5, 5, 5, 5, 5],
+                [5, 5, 5, 5, 5, 5, 5, 5, 5],
+                [5, 5, 5, 5, 5, 5, 5, 5, 5]
+            ],
+            'h5m_tag': 'Vacuum'
+        },
+        # Note that some neutron transport codes (such as OpenMC) will interpret
+        # materials with "vacuum" in the name as void material
+        'vacuum_vessel': {
+            'thickness_matrix': [
+                [15, 15, 15, 15, 15, 15, 15, 15, 15],
+                [15, 15, 15, 15, 15, 15, 15, 15, 15],
+                [15, 15, 15, 15, 15, 15, 15, 15, 15],
+                [15, 15, 15, 15, 15, 15, 15, 15, 15],
+                [15, 15, 15, 15, 15, 15, 15, 15, 15]
+            ],
+            'h5m_tag': 'vac_vessel'
+        }
+    }
 }
+# Define number of periods in stellarator plasma
+num_periods = 4
 # Define number of periods to generate
 gen_periods = 1
 # Define number of toroidal cross-sections to make
 num_phi = 60
 # Define number of poloidal points to include in each toroidal cross-section
-num_theta = 100
+num_theta = 60
 # Define magnet coil parameters
 magnets = {
     'file': 'coils.txt',
     'cross_section': ['circle', 20],
     'start': 3,
     'stop': None,
+    'sample': 6,
     'name': 'magnet_coils',
     'h5m_tag': 'magnets',
     'meshing': True
@@ -47,7 +110,8 @@
     'graveyard': False,
     'step_export': True,
     'h5m_export': 'Cubit',
-    'plas_h5m_tag': None,
+    'plas_h5m_tag': 'Vacuum',
+    'sol_h5m_tag': 'Vacuum',
     # Note the following export parameters are used only for Cubit H5M exports
     'facet_tol': 1,
     'len_tol': 5,
@@ -82,8 +146,8 @@
 logger.addHandler(f_handler)
 
 # Create stellarator
-parametric_stellarator.parametric_stellarator(
-    plas_eq, num_periods, radial_build, gen_periods, num_phi, num_theta,
+parastell.parastell(
+    plas_eq, num_periods, build, gen_periods, num_phi, num_theta,
     magnets = magnets, source = source,
     export = export, logger = logger
-    )
+)
diff --git a/README.md b/README.md
index 06f57c10..3174016b 100644
--- a/README.md
+++ b/README.md
@@ -1,10 +1,10 @@
-## parametric_stellarator
-Parametric 3D CAD modeling tool for stellarator fusion devices. This toolset takes VMEC plasma equilibrium data to build a blanket with components of uniform thickness from a user-specified radial build, as well as corresponding coil filament point-locus data to build magnet coils of user-specified cross-section.
+## ParaStell
+Parametric 3D CAD modeling tool for stellarator fusion devices. This toolset takes VMEC plasma equilibrium data to build intra-magnet components of varying thickness from a user-specified radial build, as well as corresponding coil filament point-locus data to build magnet coils of user-specified cross-section.
 
 # Dependencies
 This tool depends on:
 
-- [The VMEC tools](https://github.com/aaroncbader/pystell_uw/blob/master/read_vmec.py) developed by @aaroncbader 
+- [The VMEC tools](https://github.com/aaroncbader/pystell_uw/) developed by @aaroncbader 
 - [CadQuery](https://cadquery.readthedocs.io/en/latest/installation.html)
 - [MOAB](https://bitbucket.org/fathomteam/moab/src/master/)
 - [Coreform Cubit](https://coreform.com/products/downloads/) or [Cubit](https://cubit.sandia.gov/downloads/) by Sandia National Laboratories
diff --git a/magnet_coils.py b/magnet_coils.py
index 85bab4d8..9a412420 100644
--- a/magnet_coils.py
+++ b/magnet_coils.py
@@ -1,7 +1,6 @@
 import log
 import cubit
 import numpy as np
-import sys
 
 
 def unit_vector(vec):
@@ -277,7 +276,7 @@ def create_magnets(filaments, cross_section, meshing, logger):
     return vol_ids
 
 
-def extract_filaments(file, start, stop):
+def extract_filaments(file, start, stop, sample):
     """Extracts filament data from magnet coil data file.
     
     Arguments:
@@ -285,6 +284,7 @@ def extract_filaments(file, start, stop):
         start (int): index for line in data file where coil data begins.
         stop (int): index for line in data file where coil data ends (defaults
             to None).
+        sample (int): sampling modifier for filament points.
 
     Returns:
         filaments (list of list of list of float): list filament coordinates.
@@ -301,7 +301,7 @@ def extract_filaments(file, start, stop):
     coords = []
 
     # Extract magnet coil data
-    for line in data:
+    for i, line in enumerate(data):
         # Parse line in magnet coil data
         columns = line.strip().split()
 
@@ -317,10 +317,13 @@ def extract_filaments(file, start, stop):
         s = float(columns[3])
 
         # Coil current of zero signals end of filament
-        # If current is not zero, store coordinate in filament list
+        # If current is not zero, conditionally store coordinate in filament
+        # list
         if s != 0:
-            # Append coordinates to list
-            coords.append([x, y, z])
+            # Only store every five points
+            if i % sample == 0:
+                # Append coordinates to list
+                coords.append([x, y, z])
         # Otherwise, store filament coordinates but do not append final
         # filament point. In Cubit, continuous curves are created by setting
         # the initial and final vertex indices equal. This is handled in the
@@ -345,6 +348,9 @@ def magnet_coils(magnets, logger = None):
                 'cross_section': coil cross-section definition (list),
                 'start': starting line index for data in file (int),
                 'stop': stopping line index for data in file (int),
+                'sample': sampling modifier for filament points (int). For a
+                    user-supplied value of n, sample every n points in list of
+                    points in each filament,
                 'name': name to use for STEP export (str),
                 'h5m_tag': material tag to use in H5M neutronics model (str)
             }
@@ -369,7 +375,7 @@ def magnet_coils(magnets, logger = None):
     
     # Extract filament data
     filaments = extract_filaments(
-        magnets['file'], magnets['start'], magnets['stop']
+        magnets['file'], magnets['start'], magnets['stop'], magnets['sample']
     )
 
     # Generate magnet coil solids
diff --git a/parametric_stellarator.py b/parastell.py
similarity index 78%
rename from parametric_stellarator.py
rename to parastell.py
index b11869a8..cbce7e54 100644
--- a/parametric_stellarator.py
+++ b/parastell.py
@@ -1,5 +1,5 @@
-import magnet_coils
 import source_mesh
+import magnet_coils
 import log
 import read_vmec
 import cadquery as cq
@@ -7,6 +7,7 @@
 import cad_to_dagmc
 import numpy as np
 from scipy.optimize import root_scalar
+from scipy.interpolate import RegularGridInterpolator
 import os
 import sys
 from pymoab import core, types
@@ -291,11 +292,12 @@ def graveyard(vmec, offset, components, logger):
     return components
 
 
-def offset_point(vmec, zeta, theta, offset):
+def offset_point(vmec, s, zeta, theta, offset):
     """Stellarator offset surface root-finding problem.
 
     Arguments:
         vmec (object): plasma equilibrium object.
+        s (float): normalized magnetic closed flux surface value.
         zeta (float): toroidal angle being solved for (rad).
         theta (float): poloidal angle of interest (rad).
         offset (float): total offset of layer from plamsa (m).
@@ -306,14 +308,14 @@ def offset_point(vmec, zeta, theta, offset):
     """
     # Compute (x, y, z) point at edge of plasma for toroidal and poloidal
     # angles of interest
-    r = np.array(vmec.vmec2xyz(1.0, theta, zeta))
+    r = np.array(vmec.vmec2xyz(s, theta, zeta))
     
     # Define small number
     eta = 0.000001
 
     # Vary poloidal and toroidal angles by small amount
-    r_phi = np.array(vmec.vmec2xyz(1.0, theta, zeta + eta))
-    r_theta = np.array(vmec.vmec2xyz(1.0, theta + eta, zeta))
+    r_phi = np.array(vmec.vmec2xyz(s, theta, zeta + eta))
+    r_theta = np.array(vmec.vmec2xyz(s, theta + eta, zeta))
     
     # Take difference from plasma edge point
     delta_phi = r_phi - r
@@ -330,7 +332,7 @@ def offset_point(vmec, zeta, theta, offset):
     return pt
 
 
-def root_problem(zeta, vmec, theta, phi, offset):
+def root_problem(zeta, vmec, s, theta, phi, offset):
     """Stellarator offset surface root-finding problem. The algorithm finds the
     point on the plasma surface whose unit normal, multiplied by a factor of
     offset, reaches the desired point on the toroidal plane defined by phi.
@@ -338,6 +340,7 @@ def root_problem(zeta, vmec, theta, phi, offset):
     Arguments:
         zeta (float): toroidal angle being solved for (rad).
         vmec (object): plasma equilibrium object.
+        s (float): normalized magnetic closed flux surface value.
         theta (float): poloidal angle of interest (rad).
         phi (float): toroidal angle of interest (rad).
         offset (float): total offset of layer from plasma (m).
@@ -347,7 +350,7 @@ def root_problem(zeta, vmec, theta, phi, offset):
             angle of interest (root-finding problem definition).
     """
     # Define offset surface
-    x, y, z = offset_point(vmec, zeta, theta, offset)
+    x, y, z = offset_point(vmec, s, zeta, theta, offset)
 
     # Compute solved toroidal angle
     offset_phi = np.arctan2(y, x)
@@ -362,11 +365,12 @@ def root_problem(zeta, vmec, theta, phi, offset):
     return diff
 
 
-def offset_surface(vmec, theta, phi, offset, period_ext):
+def offset_surface(vmec, s, theta, phi, offset, period_ext):
     """Computes offset surface point.
 
     Arguments:
         vmec (object): plasma equilibrium object.
+        s (float): normalized magnetic closed flux surface value.
         theta (float): poloidal angle of interest (rad).
         phi (float): toroidal angle of interest (rad).
         offset (float): total offset of layer from plamsa (m).
@@ -376,58 +380,60 @@ def offset_surface(vmec, theta, phi, offset, period_ext):
         r (array): offset suface point (m).
     """
     # Conditionally offset poloidal profile
+    # Use VMEC plasma edge value for offset of 0
     if offset == 0:
         # Compute plasma edge point
-        r = vmec.vmec2xyz(1.0, theta, phi)
+        r = vmec.vmec2xyz(s, theta, phi)
+    
+    # Compute offset greater than zero
     elif offset > 0:
         # Root-solve for the toroidal angle at which the plasma
         # surface normal points to the toroidal angle of interest
         zeta = root_scalar(
-            root_problem, args = (vmec, theta, phi, offset), x0 = phi,
+            root_problem, args = (vmec, s, theta, phi, offset), x0 = phi,
             method = 'bisect',
             bracket = [phi - period_ext/2, phi + period_ext/2]
         )
         zeta = zeta.root
 
         # Compute offset surface point
-        r = offset_point(vmec, zeta, theta, offset)
+        r = offset_point(vmec, s, zeta, theta, offset)
+    
+    # Raise error for negative offset values
+    elif offset < 0:
+        raise ValueError(
+            'Offset must be greater than or equal to 0. Check thickness inputs '
+            'for negative values'
+        )
 
     return r
 
 
 def stellarator_torus(
-    vmec, num_periods, offset, cutter, gen_periods, num_phi, num_theta):
+    vmec, num_periods, s, cutter, gen_periods, phi_list_exp, theta_list_exp,
+    interpolator):
     """Creates a stellarator helical torus as a CadQuery object.
     
     Arguments:
         vmec (object): plasma equilibrium object.
         num_periods (int): number of periods in stellarator geometry.
-        offset (float): total offset of layer from plamsa (cm).
+        s (float): normalized magnetic closed flux surface value.
         cutter (object): CadQuery solid object used to cut period of
             stellarator torus.
         gen_periods (int): number of stellarator periods to build in model.
-        num_phi (int): number of phi geometric cross-sections to make.
-        num_theta (int): number of points defining the geometric cross-section.
-        logger (object): logger object.
+        phi_list_exp (list of float): interpolated list of toroidal angles
+            (deg).
+        theta_list_exp (list of float): interpolated list of poloidal angles
+            (deg).
+        interpolator (object): scipy.interpolate.RegularGridInterpolator object.
     
     Returns:
         torus (object): stellarator torus CadQuery solid object.
         cutter (object): updated cutting volume CadQuery solid object.
     """
-    # Check if offset distance is negative
-    if offset < 0:
-        raise ValueError('Offset must be greater than or equal to 0')
-    
-    # Determine toroidal extent of each period in radians
-    period_ext = 2*np.pi/num_periods
+    # Determine toroidal extent of each period in degrees
+    period_ext = 360.0/num_periods
     
-    # Define toroidal (phi) and poloidal (theta) arrays
-    phi_list = np.linspace(0, period_ext, num = num_phi + 1)
-    theta_list = np.linspace(0, 2*np.pi, num = num_theta + 1)[:-1]
-
-    # Convert toroidal extent of period to degrees
-    period_ext = np.rad2deg(period_ext)
-
     # Define initial angles defining each period needed
     initial_angles = np.linspace(
         period_ext, period_ext*(gen_periods - 1), num = gen_periods - 1
@@ -440,15 +446,23 @@ def stellarator_torus(
     period = cq.Workplane("XY")
 
     # Generate poloidal profiles
-    for phi in phi_list:
+    for phi in phi_list_exp:
         # Initialize points in poloidal profile
         pts = []
-        
+
+        # Convert toroidal (phi) angle from degrees to radians
+        phi = np.deg2rad(phi)
+
         # Compute array of points along poloidal profile
-        for theta in theta_list:
+        for theta in theta_list_exp[:-1]:
+            # Convert poloidal (theta) angle from degrees to radians
+            theta = np.deg2rad(theta)
+
+            # Interpolate offset according to toroidal and poloidal angles
+            offset = interpolator([np.rad2deg(phi), np.rad2deg(theta)])[0]
+
             # Compute offset surface point
-            x, y, z = offset_surface(vmec, theta, phi, offset, period_ext)
-            
+            x, y, z = offset_surface(vmec, s, theta, phi, offset, period_ext)
             # Convert from m to cm
             pt = (x*100, y*100, z*100)
             # Append point to poloidal profile
@@ -483,6 +497,44 @@ def stellarator_torus(
     return torus, cutter
 
 
+def expand_ang(ang_list, num_ang):
+    """Expands list of angles by linearly interpolating according to specified
+    number to include in stellarator build.
+
+    Arguments:
+        ang_list (list of float): user-supplied list of toroidal or poloidal
+            angles (deg).
+        num_ang (int): number of angles to include in stellarator build.
+    
+    Returns:
+        ang_list_exp (list of float): interpolated list of angles (deg).
+    """
+    # Initialize interpolated list of angles
+    ang_list_exp = []
+
+    # Compute total angular extent of supplied list
+    ang_ext = ang_list[-1] - ang_list[0]
+
+    # Compute average distance between angles to include in stellarator build
+    ang_diff_avg = ang_ext/num_ang
+    
+    # Loop over supplied angles
+    for ang, next_ang in zip(ang_list[:-1], ang_list[1:]):
+        # Compute number of angles to interpolate
+        n_ang = int(np.ceil((next_ang - ang)/ang_diff_avg))
+
+        # Interpolate angles and append to storage list
+        ang_list_exp = np.append(
+            ang_list_exp,
+            np.linspace(ang, next_ang, num = n_ang + 1)[:-1]
+        )
+
+    # Append final specified angle to storage list
+    ang_list_exp = np.append(ang_list_exp, ang_list[-1])
+
+    return ang_list_exp
+
+
 # Define default export dictionary
 export_def = {
     'exclude': [],
@@ -490,6 +542,7 @@ def stellarator_torus(
     'step_export': True,
     'h5m_export': None,
     'plas_h5m_tag': None,
+    'sol_h5m_tag': None,
     'facet_tol': None,
     'len_tol': None,
     'norm_tol': None,
@@ -501,8 +554,8 @@ def stellarator_torus(
 }
 
 
-def parametric_stellarator(
-    plas_eq, num_periods, radial_build, gen_periods, num_phi = 60,
+def parastell(
+    plas_eq, num_periods, build, gen_periods, num_phi = 60,
     num_theta = 100, magnets = None, source = None, export = export_def,
     logger = None):
     """Generates CadQuery workplane objects for components of a
@@ -515,14 +568,25 @@ def parametric_stellarator(
     Arguments:
         plas_eq (str): path to plasma equilibrium NetCDF file.
         num_periods (int): number of periods in stellarator geometry.
-        radial_build (dict of dicts): dictionary of component names, each with
-            a corresponding dictionary of layer thickness and optional material
-            tag to use in H5M neutronics model in the form
-            {'component': {'thickness': thickness (float, cm), 'h5m_tag':
-            h5m_tag (str)}}
-            Concentric layers will be built in the order given. If no alternate
-            material tag is supplied for the H5M file, the given component name
-            will be used.
+        build (dict): dictionary of list of toroidal and poloidal angles, as
+            well as dictionary of component names with corresponding thickness
+            matrix and optional material tag to use in H5M neutronics model.
+            The thickness matrix specifies component thickness at specified
+            (polidal angle, toroidal angle) pairs. This dictionary takes the
+            form
+            {
+                'phi_list': list of float (deg),
+                'theta_list': list of float (deg),
+                'wall_s': closed flux index extrapolation at wall (float),
+                'radial_build': {
+                    'component': {
+                        'thickness_matrix': list of list of float (cm),
+                        'h5m_tag': h5m_tag (str)
+                    }
+                }
+            }
+            If no alternate material tag is supplied for the H5M file, the
+            given component name will be used.
         gen_periods (int): number of stellarator periods to build in model.
         num_phi (int): number of phi geometric cross-sections to make for each
             period (defaults to 60).
@@ -534,6 +598,9 @@ def parametric_stellarator(
                 'cross_section': coil cross-section definition (list),
                 'start': starting line index for data in file (int),
                 'stop': stopping line index for data in file (int),
+                'sample': sampling modifier for filament points (int). For a
+                    user-supplied value of n, sample every n points in list of
+                    points in each filament,
                 'name': name to use for STEP export (str),
                 'h5m_tag': material tag to use in H5M neutronics model (str)
                 'meshing': setting for tetrahedral mesh generation (bool)
@@ -567,6 +634,10 @@ def parametric_stellarator(
                 'plas_h5m_tag': optional alternate material tag to use for
                     plasma. If none is supplied and the plasma is not excluded,
                     'plasma' will be used (str, defaults to None),
+                'sol_h5m_tag': optional alternate material tag to use for 
+                    scrape-off layer. If none is supplied and the scrape-off
+                    layer is not excluded, 'sol' will be used (str, defaults to
+                    None),
                 'facet_tol': maximum distance a facet may be from surface of
                     CAD representation for Cubit export (float, defaults to
                     None),
@@ -593,7 +664,7 @@ def parametric_stellarator(
 
     Returns:
         strengths (list): list of source strengths for each tetrahedron (1/s).
-            Only returned if source mesh is generated
+            Returned only if source mesh is generated.
     """
     # Conditionally instantiate logger
     if logger == None or not logger.hasHandlers():
@@ -621,40 +692,80 @@ def parametric_stellarator(
     # Initialize component storage dictionary
     components = {}
     
-    # Prepend plasma to radial build
-    full_radial_build = {'plasma': {'thickness': 0}, **radial_build}
+    # Extract toroidal (phi) and poloidal (theta) arrays
+    phi_list = build['phi_list']
+    theta_list = build['theta_list']
+    # Extract closed flux surface index extrapolation at wall
+    wall_s = build['wall_s']
+    # Extract radial build
+    radial_build = build['radial_build']
+
+    # Extract array dimensions
+    n_phi = len(phi_list)
+    n_theta = len(theta_list)
+
+    # Conditionally prepend scrape-off layer to radial build
+    if wall_s != 1.0:
+        sol_thickness_mat = np.zeros((n_phi, n_theta))
+        radial_build = {
+            'sol': {'thickness_matrix': sol_thickness_mat},
+            **radial_build
+        }
     
-    # Initialize offset value
-    offset = 0.0
+    # Prepend plasma to radial build
+    plas_thickness_mat = np.zeros((n_phi, n_theta))
+    radial_build = {
+            'plasma': {'thickness_matrix': plas_thickness_mat},
+            **radial_build
+        }
 
     # Initialize volume used to cut periods
     cutter = None
+
+    # Linearly interpolate angles to expand phi and theta lists
+    phi_list_exp = expand_ang(phi_list, num_phi)
+    theta_list_exp = expand_ang(theta_list, num_theta)
+
+    # Initialize offset matrix
+    offset_mat = np.zeros((n_phi, n_theta))
     
     # Generate components in radial build
-    for name, layer_data in full_radial_build.items():
-        
+    for name, layer_data in radial_build.items():
         # Notify which component is being generated
         logger.info(f'Building {name}...')
         
-        # Conditionally assign plasma h5m tag
+        # Conditionally assign plasma h5m tag and reference closed flux surface
         if name == 'plasma':
             if export_dict['plas_h5m_tag'] is not None:
                 layer_data['h5m_tag'] = export_dict['plas_h5m_tag']
+            s = 1.0
+        else:
+            s = wall_s
+
+        # Conditionally assign scrape-off layer h5m tag
+        if name == 'sol':
+            if export_dict['sol_h5m_tag'] is not None:
+                layer_data['h5m_tag'] = export_dict['sol_h5m_tag']
         
         # Conditionally populate h5m tag for layer
         if 'h5m_tag' not in layer_data:
             layer_data['h5m_tag'] = name
 
-        # Extract layer thickness
-        thickness = layer_data['thickness']
-        # Compute offset, converting from cm to m
-        offset += thickness/100
+        # Extract layer thickness matrix
+        thickness_mat = layer_data['thickness_matrix']
+        # Compute offset list, converting from cm to m
+        offset_mat += np.array(thickness_mat)/100
 
+        # Build offset interpolator
+        interp = RegularGridInterpolator((phi_list, theta_list), offset_mat)
+        
         # Generate component
         try:
             torus, cutter = stellarator_torus(
-                vmec, num_periods, offset, cutter, gen_periods, num_phi,
-                num_theta
+                vmec, num_periods, s,
+                cutter, gen_periods,
+                phi_list_exp, theta_list_exp,
+                interp
             )
         except ValueError as e:
             logger.error(e.args[0])
@@ -668,6 +779,9 @@ def parametric_stellarator(
 
     # Conditionally build graveyard volume
     if export_dict['graveyard']:
+        # Extract maximum offset
+        offset = 2*max(max(offset_mat))
+        # Build graveyard
         components = graveyard(vmec, offset, components, logger)
 
     # Conditionally initialize Cubit
diff --git a/source_mesh.py b/source_mesh.py
index c85224c5..9e5ef817 100644
--- a/source_mesh.py
+++ b/source_mesh.py
@@ -182,13 +182,14 @@ def create_tets_from_wedge(
         )
 
 
-def get_vertex_id(vert_ids, num_s, num_theta):
+def get_vertex_id(vert_ids, num_phi, num_s, num_theta):
     """Computes vertex index in row-major order as stored by MOAB from three-dimensional n x 3 matrix indices.
 
     Arguments:
         vert_ids (list of int): list of vertex indices expressed in n x 3
             matrix order. The list format is:
             [toroidal angle index, flux surface index, poloidal angle index]
+        num_phi (int): number of toroidal angles defining mesh.
         num_s (int): number of closed magnetic flux surfaces defining mesh.
         num_theta (int): number of poloidal angles defining mesh.
     
@@ -196,11 +197,23 @@ def get_vertex_id(vert_ids, num_s, num_theta):
         id (int): vertex index in row-major order as stored by MOAB
     """
     # Compute index offset from magnetic axis
-    ma_offset = vert_ids[0] * (1 + ((num_s - 1) * num_theta))
+    if vert_ids[0] == num_phi - 1:
+        # Ensure logical connectivity at final toroidal angle
+        ma_offset = 0
+    else:
+        # Otherwise, compute index of magnetic axis
+        ma_offset = vert_ids[0] * (1 + ((num_s - 1) * (num_theta - 1)))
     # Compute index offset from closed flux surface
-    s_offset = vert_ids[1] * num_theta
+    s_offset = vert_ids[1] * (num_theta - 1)
+    # Compute index offset from poloidal angle
+    if vert_ids[2] == num_theta:
+        # Ensure logical connectivity at final poloidal angle
+        theta_offset = 1
+    else:
+        # Otherwise, compute index of poloidal angle
+        theta_offset = vert_ids[2]
     # Compute vertex index
-    id = ma_offset + s_offset + vert_ids[2]
+    id = ma_offset + s_offset + theta_offset
 
     return id
 
@@ -236,10 +249,10 @@ def create_mesh(
     # Create tetrahedra, looping through vertices
     for i in range(num_phi - 1):
         # Create tetrahedra for wedges at center of plasma
-        for k in range(num_theta - 1):
+        for k in range(1, num_theta):
             # Define six wedge vertex indices in matrix format as well as index
             # offset accounting for magnetic axis vertex in the form
-            # [toroidal index, flux surface index, poloidal index, offset]
+            # [toroidal index, flux surface index, poloidal index]
             wedge_id_data = [
                 [i,     0, 0,   ],
                 [i,     0, k,   ],
@@ -252,7 +265,7 @@ def create_mesh(
             ids_wedge = []
             # Define vertex indices in row-major order for wedges
             for vertex in wedge_id_data:
-                ids_wedge += [get_vertex_id(vertex, num_s, num_theta)]
+                ids_wedge += [get_vertex_id(vertex, num_phi, num_s, num_theta)]
             # Create tetrahedra from wedge
             create_tets_from_wedge(
                 mbc, tag_handle, tet_set, mbc_verts, verts, verts_s, ids_wedge,
@@ -262,11 +275,11 @@ def create_mesh(
         # Create tetrahedra for hexahedra beyond center of plasma
         for j in range(num_s - 2):
             # Create tetrahedra for current hexahedron
-            for k in range(num_theta - 1):
+            for k in range(1, num_theta):
                 # Define eight hexahedron vertex indices in matrix format as
                 # well as index offset accounting for magnetic axis vertex in
                 # the form
-                # [toroidal index, flux surface index, poloidal index, offset]
+                # [toroidal index, flux surface index, poloidal index]
                 hex_id_data = [
                     [i,     j,     k    ],
                     [i,     j + 1, k,   ],
@@ -281,7 +294,9 @@ def create_mesh(
                 ids_hex = []
                 # Define vertex indices in row-major order for hexahedra
                 for vertex in hex_id_data:
-                    ids_hex += [get_vertex_id(vertex, num_s, num_theta)]
+                    ids_hex += [
+                        get_vertex_id(vertex, num_phi, num_s, num_theta)
+                    ]
                 # Create tetrahedra from hexahedron
                 create_tets_from_hex(
                     mbc, tag_handle, tet_set, mbc_verts, verts, verts_s,
@@ -309,8 +324,8 @@ def create_vertices(vmec, mbc, num_s, num_theta, num_phi):
             vertex.
     """
     # Compute total number of vertices
-    num_verts = num_phi * (1 + ((num_s - 1) * num_theta))
-    # Initialize n x 3 array of verticex coordinates
+    num_verts = (num_phi - 1) * (1 + ((num_s - 1) * (num_theta - 1)))
+    # Initialize n x 3 array of vertex coordinates
     verts = np.zeros((num_verts, 3))
     # Initialize array of vertex closed flux surface indices
     verts_s = np.zeros(num_verts)
@@ -331,7 +346,7 @@ def create_vertices(vmec, mbc, num_s, num_theta, num_phi):
     # Initialize vertex index
     i_vert = 0
     # Compute vertices in Cartesian space
-    for i_phi in range(num_phi):
+    for i_phi in range(num_phi - 1):
         # Compute toroidal angle for vertex
         phi = i_phi * delta_phi
         # Determine vertex at magnetic axis, convertic to cm
@@ -346,7 +361,8 @@ def create_vertices(vmec, mbc, num_s, num_theta, num_phi):
             # Compute closed flux surface index for vertex
             s = i_s * delta_s
 
-            for i_theta in range(num_theta):
+            # Exclude final poloidal angle as the final and initial are the same
+            for i_theta in range(num_theta - 1):
                 # Compute poloidal angle for vertex
                 theta = i_theta * delta_theta
                 # Determine vertex at specified toroidal angle, closed flux