Hello
In my current simulations, I generate a block using mesh with the fibre orienation as input parameter. I have read the paper on how the angles are defined in openCARP and I have scanned the source code to look for anything that I might have misunderstood. The output I get does not correspond to the thing I believe I gave as input.
Here, you can see the function I use to generate the block :
#----------------------------------------------------
# Function : construct block
#----------------------------------------------------
def construct_block(job,args):
# mesh class expects all physical sizes to be specified in mm!
# Generate the block with bath:
block_x = args.size_block[0]
block_y = args.size_block[1]
block_z = args.size_block[2] #wall thickness
geom = mesh.Block(size=(block_x,block_y,block_z), resolution=args.res, etype='tetra')
bath_thick = args.size_bath
geom.set_bath(thickness=(bath_thick,bath_thick,bath_thick), both_sides=True)
# Put the corner of the tissue block at the origin
geom.corner_at_origin()
# Transmural gradient: thickness of layer/region
endo = block_z/4
mid = block_z/2
epi = block_z/4
assert((endo+mid+epi) == block_z)
# Apico-basal gradient: thickness of layer/region
apex = block_x/3
central = block_x/3
base = block_x/3
assert((apex+central+base) == block_x)
# Tags
tag_apex_endo = 300 + 25
tag_apex_mid = 300 + 50
tag_apex_epi = 300 + 75
tag_central_endo = 200 + 25
tag_central_mid = 200 + 50
tag_central_epi = 200 + 75
tag_base_endo = 100 + 25
tag_base_mid = 100 + 50
tag_base_epi = 100 + 75
# Add regions ------> (0,0,0) is the center of the mesh, not the left lower corner!
# Contruct mesh --> add regions to geom
# 0 --> lower left corner, 1 --> upper right corner
x0_apex = -block_x/2
x0_central = -block_x/6
x0_base = block_x/6
x1_apex = x0_central
x1_central = x0_base
x1_base = block_x/2
z0_endo = -block_z/2
z0_mid = -block_z/4
z0_epi = block_z/4
z1_endo = z0_mid
z1_mid = z0_epi
z1_epi = block_z/2
apex_endo = mesh.BoxRegion((x0_apex, -block_y/2, z0_endo),(x1_apex, block_y/2, z1_endo),tag_apex_endo)
geom.add_region(apex_endo)
apex_mid = mesh.BoxRegion((x0_apex, -block_y/2, z0_mid),(x1_apex, block_y/2, z1_mid),tag_apex_mid)
geom.add_region(apex_mid)
apex_epi = mesh.BoxRegion((x0_apex, -block_y/2, z0_epi),(x1_apex, block_y/2, z1_epi),tag_apex_epi)
geom.add_region(apex_epi)
central_endo = mesh.BoxRegion((x0_central, -block_y/2, z0_endo),(x1_central, block_y/2, z1_endo),tag_central_endo)
geom.add_region(central_endo)
central_mid = mesh.BoxRegion((x0_central, -block_y/2, z0_mid),(x1_central, block_y/2, z1_mid),tag_central_mid)
geom.add_region(central_mid)
central_epi = mesh.BoxRegion((x0_central, -block_y/2, z0_epi),(x1_central, block_y/2, z1_epi),tag_central_epi)
geom.add_region(central_epi)
base_endo = mesh.BoxRegion((x0_base, -block_y/2, z0_endo),(x1_base, block_y/2, z1_endo),tag_base_endo)
geom.add_region(base_endo)
base_mid = mesh.BoxRegion((x0_base, -block_y/2, z0_mid),(x1_base, block_y/2, z1_mid),tag_base_mid)
geom.add_region(base_mid)
base_epi = mesh.BoxRegion((x0_base, -block_y/2, z0_epi),(x1_base, block_y/2, z1_epi),tag_base_epi)
geom.add_region(base_epi)
# Set fibre angles
geom.set_fibres(fibre_endo=args.fibre_angles[1],
fibre_epi=args.fibre_angles[0],
sheet_endo=0,
sheet_epi=0)
# Generate and return base name
meshname=mesh.generate(geom)
# build tag dictionary
tags = {'tag_apex_endo' : tag_apex_endo,
'tag_apex_mid' : tag_apex_mid,
'tag_apex_epi' : tag_apex_epi,
'tag_central_endo': tag_central_endo,
'tag_central_mid' : tag_central_mid,
'tag_central_epi' : tag_central_epi,
'tag_base_endo' : tag_base_endo,
'tag_base_mid' : tag_base_mid,
'tag_base_epi' : tag_base_epi }
return meshname, geom, tags
I defined different regions, but I am not yet assigning different properties to them. I used TP06, EPI cells, with these conductivities:
g_il = 0.1527
g_it = 0.0312
g_el = 0.5485
g_et = 0.3361
g_bath = 0.7
where t=n. My stimulus is a plane wave that propagates in one direction and I already checked if the ground might be an issue by both using a corner of the medium as ground and a thin slab parallel to the wave. I compared monodomain, bidomain and pseudo-bidomain and they all give the same behaviour. I check parallel fibres and that has a nice symmetric wavefront.
The issue is that the wave front using fibre orientation is sort of linear and not symmetric. When I plot the fibre orientation, one can also see that around the middle of the region, they become parallel and stop turning.
I have pictures to show the effect I am talking about, but I cannot seem to find how to include them here. This is not the behavior one would expect when setting the endo and epi fibre angle. What is happening here?
Thank you a lot for your time and assistance!
Lore
