First time here? Checkout the FAQ!
x
0 votes
by (540 points)

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

1 Answer

+1 vote
by (2.6k points)
selected by
 
Best answer
This is a bug which has been fixed with commit 484dc4032.

Cheers, Aurel
Welcome to openCARP Q&A. Ask questions and receive answers from other members of the community. For best support, please use appropriate TAGS!
architecture, carputils, documentation, experiments, installation-containers-packages, limpet, slimfem, website, governance
...