Updated from last commit (master branch)

tutorials/02_EP_tissue/04_tagging/run.py

@@ -176,46 +176,36 @@ to colour them differently for each tag region.
 
 For a more detailed introduction to the region concept as implemented in openCARP see
 the :ref:`Region Tagging Section <region-definition>` of the manual.
-
-
 """
-
+import os
 EXAMPLE_DESCRIPTIVE_NAME = 'Region tagging'
 EXAMPLE_AUTHOR = 'Edward Vigmond <edward.vigmond@u-bordeaux.fr>'
+EXAMPLE_DIR = os.path.dirname(__file__)
+CALLER_DIR = os.getcwd()
 
-import os
-import sys
-from datetime import date
-
+import numpy as np
+from datetime  import date
+from carputils import mesh
 from carputils import settings
 from carputils import tools
-from carputils import mesh
-from carputils import testing
+from carputils.carpio import txt
 
 def parser():
     parser = tools.standard_parser()
     group  = parser.add_argument_group('experiment specific options')
-
     group.add_argument('--mesher',
                         default=[],
-                        choices=[ 'sphere',
-                                  'block',
-                                  'steps'],
+                        choices=['sphere', 'block', 'steps'],
                         nargs='*',
                         help='regions implemented in mesher')
-
     group.add_argument('--dynamic',
                         default=[], 
-                        choices=['cylinder',
-                                 'sphere',
-                                 'block' ], 
+                        choices=['cylinder', 'sphere', 'block'],
                         nargs='*',
                         help='regions implemented dynamically')
-
     group.add_argument('--add-bath',
                         action='store_true',
-                        help='mesher regions are bath')    
-
+                        help='mesher regions are bath')
     return parser
 
 def jobID(args):
@@ -224,96 +214,110 @@ def jobID(args):
     """
     today = date.today()
     regspec = ''
-    for r in args.mesher :
-        if r is None :
+    for reg in args.mesher:
+        if reg is None:
             break
-        regspec += r[0].upper()+r[1]
-    if len(regspec) : regspec += '_'
-    for r in args.dynamic :
-        if r is None :
+        regspec += reg[0].upper() + reg[1]
+    if len(regspec):
+        regspec += '_'
+    for reg in args.dynamic:
+        if reg is None:
             break
-        regspec += r[0]
+        regspec += reg[0]
 
-    if args.add_bath  : regspec += '_Bth'
+    if args.add_bath:
+        regspec += '_Bth'
 
-    #print (regspec)
-        
     return 'regions_{}'.format(regspec) 
 
 
-def tagregopt( reg, field, val ) :
-    return ['-tagreg['+str(reg)+'].'+field, val ]
+def tagregopt(reg, field, val):
+    return ['-tagreg['+str(reg)+'].'+field, val]
 
 
-@tools.carpexample(parser, jobID)
+@tools.carpexample(parser, jobID, clean_pattern='regions_B*|regions_C*|regions_C*|meshes')
 def run(args, job):
+    ExtraTags = [] # labels used for extracellular domain
 
-    x = 5.
-    y = 5.
-    z = 5.
+    x = 5. # slab x dim in mm
+    y = 5. # slab y dim in mm
+    z = 5. # slab z dim in mm
 
     # Block which is thin in z direction
     geom = mesh.Block(size=(x, y, z))
-    geom.set_bath( (1,1,1), both_sides=True )
+    geom.set_bath((1, 1, 1), both_sides=True)
 
     # Set fibre angle to 0, sheet angle to 0
     geom.set_fibres(-60, 60, 90, 90)
 
     # define a tag region
-    for i,r in enumerate(args.mesher) :
+    for i, reg in enumerate(args.mesher):
 
         bath = 1 if args.add_bath else 0
 
-        if r == 'block':
-          geom.add_region(mesh.BoxRegion((-x/2-0.1,-y/2.-0.1,-z/2.-0.1), (0,y/2.+0.1,z/2.+0.1), tag=10*(i+1),bath=bath))
+        if reg == 'block':
+          geom.add_region(mesh.BoxRegion((-x/2-0.1,-y/2.-0.1,-z/2.-0.1),
+                                         ( 0,       y/2.+0.1, z/2.+0.1), tag=10*(i+1), bath=bath))
     
-        elif r == 'steps':
-            nstep=4
+        elif reg == 'steps':
+            nstep = 4
             dz = z/float(nstep)
             for s in range(nstep) :
                 geom.add_region(mesh.BoxRegion((-x/2-0.1,-y/2.-0.1,-z/2.+dz*s), 
-                                        (x/2+0.1,y/2.+0.1,-z/2.+dz*(s+1)), tag=10*(1+i)+s,bath=bath))
-    
-        elif r == 'sphere':
-          reg = geom.add_region(mesh.SphereRegion(1., (0,0,0), bath=bath, tag=10*(i+1)))
-    
+                                                (x/2+0.1, y/2.+0.1,-z/2.+dz*(s+1)), tag=10*(1+i)+s, bath=bath))
+
+        elif reg == 'sphere':
+            geom.add_region(mesh.SphereRegion(1., (0,0,0), bath=bath, tag=10*(i+1)))
+
         else:
            # None, don't add region 
            pass
 
-    dyn_reg = [ '-numtagreg', len(args.dynamic) ]
-    for reg,r in enumerate(args.dynamic) :
+    # Generate and return base name
+    meshname = mesh.generate(geom)
+
+    # Query for element labels
+    _, etags,_ = txt.read(meshname + '.elem')
+    etags = np.unique(etags)
+
+    IntraTags = etags[etags != 0].tolist()
+    ExtraTags = etags.tolist().copy()
+
+
+    dyn_reg = ['-numtagreg', len(args.dynamic)]
+    for reg, r in enumerate(args.dynamic) :
 
+        dyn_tag = 100 * (reg + 1)
         if r == 'block':
-          dyn_reg.extend( tagregopt(reg, 'type',    2) )
-          dyn_reg.extend( tagregopt(reg, 'p0',  '{ -1000 -1000 -200 }') )
-          dyn_reg.extend( tagregopt(reg, 'p1',  '{ 1000 1000 100 }') )
-          dyn_reg.extend( tagregopt(reg, 'tag', 100*(reg+1)) )
-    
+            dyn_reg.extend( tagregopt(reg, 'type', 2) )
+            dyn_reg.extend( tagregopt(reg, 'p0', '{ -1000 -1000 -200 }') )
+            dyn_reg.extend( tagregopt(reg, 'p1', '{ 1000 1000 100 }') )
+            dyn_reg.extend( tagregopt(reg, 'tag', dyn_tag) )
+
         elif r == 'cylinder':
-          dyn_reg.extend( tagregopt(reg, 'type',        3) )
-          dyn_reg.extend( tagregopt(reg, 'p0',      '{ '+str(-x*1000/2.)+' 2000 1500 }') )
-          dyn_reg.extend( tagregopt(reg, 'p1',      '{ -1000 2000 1500  }') )
-          dyn_reg.extend( tagregopt(reg, 'radius',  250) )
-          dyn_reg.extend( tagregopt(reg, 'tag',     100*(reg+1)) )
-    
+            dyn_reg.extend( tagregopt(reg, 'type', 3) )
+            dyn_reg.extend( tagregopt(reg, 'p0', '{ '+str(-x*1000/2.)+' 2000 1500 }') )
+            dyn_reg.extend( tagregopt(reg, 'p1', '{ -1000 2000 1500  }') )
+            dyn_reg.extend( tagregopt(reg, 'radius', 250) )
+            dyn_reg.extend( tagregopt(reg, 'tag', dyn_tag) )
+
         elif r == 'sphere':
-          dyn_reg.extend( tagregopt(reg, 'type',        1) )
-          dyn_reg.extend( tagregopt(reg, 'p0',      '{ 0 700 0 }') )
-          dyn_reg.extend( tagregopt(reg, 'radius',  500 ) )
-          dyn_reg.extend( tagregopt(reg, 'tag',     100*(reg+1)) )
-    
+            dyn_reg.extend( tagregopt(reg, 'type', 1) )
+            dyn_reg.extend( tagregopt(reg, 'p0', '{ 0 700 0 }') )
+            dyn_reg.extend( tagregopt(reg, 'radius', 500) )
+            dyn_reg.extend( tagregopt(reg, 'tag', dyn_tag) )
+
         else:
-           # None, don't add region 
+           # None, don't add region
+           dyn_tag = None
            pass
 
+        if dyn_tag:
+            # even tags generated in a pre-processing step (within openCARP)
+            # need to be registered here!
+            IntraTags.append(dyn_tag)
+            ExtraTags.append(dyn_tag)
 
-    # Generate and return base name
-    meshname = mesh.generate(geom)
-
-    # define default view settings
-    view = 'tagged.mshz'
-  
     # Add all the non-general arguments
     cmd = tools.carp_cmd() 
 
@@ -321,37 +325,36 @@ def run(args, job):
             '-gridout_i',  2, 
             '-gridout_e',  2, 
             '-experiment', 3,
-            '-simID',     job.ID,
-            '-retagfile', 'retagged.dat' ]  
-
+            '-simID',      job.ID,
+            '-retagfile', 'retagged.dat']
+    cmd += tools.gen_physics_opts(ExtraTags=ExtraTags, IntraTags=IntraTags)
     cmd += dyn_reg
 
     # Run simulations
-    if args.dry :
-        print (cmd)
-    else :
-        job.carp(cmd)
+    job.carp(cmd, 'Building model only.')
 
     # Do visualization
-    view = 'tagging.mshz'
+    if args.visualize and not settings.platform.BATCH:
+        # prepare file paths
+        geom = os.path.join(CALLER_DIR, job.ID, os.path.basename(meshname) + '_e')
+        view = os.path.join(EXAMPLE_DIR, 'tagging.mshz')
         
-    geom = os.path.join(job.ID, os.path.basename(meshname)+'_e')
-
-    # Visualize with GUI
-    if args.webGUI:
-        folderNameExp = f'/experiments/04_tagging_{job.ID}'
-        os.mkdir(folderNameExp)
-        cmdMesh = f'meshtool collect -imsh={geom} -omsh={folderNameExp}/{job.ID} -ifmt=carp_txt -ofmt=ens_bin'
-        outMesh = os.system(cmdMesh)
-        if(outMesh == 0):
-            print('Meshtool - conversion successfully')
-        else:
-            print('Meshtool - conversion failure')
+        # Visualize with GUI
+        if args.webGUI:
+            folderNameExp = f'/experiments/04_tagging_{job.ID}'
+            os.mkdir(folderNameExp)
+            cmdMesh = f'meshtool collect -imsh={geom} -omsh={folderNameExp}/{job.ID} -ifmt=carp_txt -ofmt=ens_bin'
+            outMesh = os.system(cmdMesh)
+            if(outMesh == 0):
+                print('Meshtool - conversion successfully')
+            else:
+                print('Meshtool - conversion failure')
     
-    #Visualize with meshalyzer
+        #Visualize with meshalyzer
+        else:
+            job.meshalyzer(geom, view, compsurf=True)
     else:
-        job.meshalyzer(geom, view, compsurf=True)
-
+        print('Use --visualize option to display the results!')
 
 if __name__ == '__main__':
     run()