I made some modifications to my code, but I'm still encountering an issue when trying to run the simulation using my spherical mesh. Below is the updated code, and the error I get is:
Modified code:
#!/usr/bin/env python
import os
from datetime import date
from carputils import settings
from carputils import tools
from carputils.carpio import txt
import numpy as np
def parser():
parser = tools.standard_parser()
group = parser.add_argument_group('experiment specific options')
group.add_argument('--duration',
type=float,
default=1000.,
help='Duration of simulation in [ms] (default: 1000.)')
group.add_argument('--S1-strength',
type=float,
default=400.,
help='Stimulus strength in [uA/cm^2] (default: 400.)')
group.add_argument('--S1-dur',
type=float,
default=2.,
help='Stimulus duration in [ms] (default: 2.)')
return parser
def jobID(args):
today = date.today()
return '{}_monodomain_{}'.format(today.isoformat(), args.duration)
def sphere_boundary_condition(nodes, center, radius, label):
"""
Generates a boundary condition for a spherical surface.
Args:
nodes: np.array with the coordinates of the mesh nodes.
center: Center of the sphere [x, y, z].
radius: Radius of the sphere.
label: Label for the boundary condition.
Returns:
List of commands to define the boundary condition.
"""
boundary_nodes = []
for i, node in enumerate(nodes):
dist = np.linalg.norm(node - np.array(center))
if np.isclose(dist, radius, atol=1e-3): # Tolerance to find the surface
boundary_nodes.append(i)
condition = []
for i, node_id in enumerate(boundary_nodes):
condition += [
f'-boundary[{i}].name', label,
f'-boundary[{i}].type', 1, # Type 1: Stimulus
f'-boundary[{i}].nodes', node_id
]
return condition
@tools.carpexample(parser, jobID)
def run(args, job):
# Mesh setup
meshname = "sphere" # Mesh name
nodes = np.loadtxt(meshname + ".pts", skiprows=1) # Read mesh nodes ignoring the first row
_, etags, _ = txt.read(meshname + '.elem')
etags = np.unique(etags)
IntraTags = etags[etags != 0].tolist()
ExtraTags = etags.tolist().copy()
# Sphere parameters
center = [0, 0, 0] # Center of the sphere
radius = 1.0 # Radius of the sphere
# Generate the boundary condition for the spherical surface
electrode = sphere_boundary_condition(nodes, center, radius, 'stimulus')
# Calculate number of stimuli
num_stim = int((args.duration - 20) / 500) + 1 # Includes the initial stimulus at 20 ms
# Create stimulus configuration
stim = ['-num_stim', num_stim]
for i in range(num_stim):
stim += [
f'-stimulus[{i}].name', 'S1',
f'-stimulus[{i}].stimtype', 0,
f'-stimulus[{i}].strength', args.S1_strength,
f'-stimulus[{i}].duration', args.S1_dur,
f'-stimulus[{i}].start', 20 + i * 500 # Stimulus every 500 ms starting from 20 ms
]
# Generate CARP command
cmd = tools.carp_cmd(os.path.join(os.path.dirname(__file__), 'monodomain.par'))
cmd += ['-simID', job.ID,
'-meshname', meshname,
'-dt', 5,
'-tend', args.duration]
cmd += tools.gen_physics_opts(ExtraTags=ExtraTags, IntraTags=IntraTags)
cmd += stim + electrode
if args.visualize:
cmd += ['-gridout_i', 3, '-gridout_e', 3, '-spacedt', 0.1]
# Run simulation
job.carp(cmd)
if __name__ == '__main__':
run()
