Added options for boundary conditions; relax_ic_in_csf; relax_ic_coeff;

Better management of boundary conditions with dirichlet_at_walls vs
dirichlet_at_skull. Some exception handling when setting up the
AdLemModel. BUG REMAINING: the null space when using dirichlet_at_skull
and no relax_ic_in_csf is not supported yet! Currently just provides a
warning for this case.

scripts/runAdLemModel.py

@@ -1,7 +1,6 @@
 #!/usr/bin/env python
 import os
 import os.path as op
-import sys
 import argparse as ag
 import bish_utils as bu
 
@@ -16,36 +15,46 @@ def get_input_options():
         ' written here.')
     parser.add_argument(
         'res_prefix', help='prefix to be added to all the output filenames')
-    parser.add_argument('lame_paras', help='input lame parameters: '
-                        'Format: "muTissue muCsf lambdaTissue lambdaCsf" '
-                        'If both --mu_image and --in_dti used, lame_paras wont'
-                        ' be used and those images will be used instead. ')
+    parser.add_argument(
+        'lame_paras', help='input lame parameters: Format: "muTissue muCsf '
+        'lambdaTissue lambdaCsf" If both --mu_image and --in_dti used, '
+        'lame_paras wont be used and those images will be used instead.')
+    parser.add_argument('boundary_condition', help='Possible values: '
+                        'dirichlet_at_walls dirichlet_at_skull')
     parser.add_argument('atrophy', help='input atrophy file.')
     parser.add_argument('in_seg', help='input brain seg ')
     parser.add_argument(
         'in_img', help='valid input image that will be warped by the obtained'
         ' displacement fields from the model.')
+    parser.add_argument('petsc_op_file', help='file with petsc options')
     parser.add_argument('time_steps', help='(integer) '
                         'number of time steps you want to solve the system.')
+    parser.add_argument(
+        '--relax_ic_in_csf', action='store_true', help='If given, relaxes IC.')
+    parser.add_argument(
+        '-k', '--relax_ic_coeff', help='relaxation coeff. Relevant only when '
+        'using --relax_ic_in_csf. If --relax_ic_in_csf given and this option is'
+        ' not provided, uses reciprocal of lamda as coefficient.')
     parser.add_argument('-m', '--mu_file', help='input mu image file: ')
-    parser.add_argument('--use_dti', dest='use_dti', action='store_true',
-                        help='If the option is provided, input DTI must be '
-                        'given with the option -d')
-    parser.add_argument('-d', '--in_dti', help='valid tensor image file present'
-                        ' in patientID directory.')
-    parser.add_argument('-p', '--petsc_op_file', help='file with petsc options')
-    parser.add_argument('--wrt_press', action='store_true',
-                        help='Write pressure map on the disk.')
-    parser.add_argument('--wrt_force', action='store_true',
-                        help='Write force on the the disk.')
-    parser.add_argument('--wrt_residual', action='store_true',
-                        help='Write residual on the disk.')
-    parser.add_argument('--in_cluster', action='store_true',
-                        help='If provided, launches as a job on the cluster.')
-    # parser.add_argument('--no_petsc_summary', action='store_true',
-    #                     help='If provided, does not detail the petsc outputs.')
+    parser.add_argument(
+        '--use_dti', dest='use_dti', action='store_true', help='If the option '
+        'is provided, input DTI must be given with the option -d')
+    parser.add_argument(
+        '-d', '--in_dti', help='tensor image file in patientID directory.')
+    parser.add_argument(
+        '--wrt_press', action='store_true', help='Write pressure file.')
+    parser.add_argument(
+        '--wrt_force', action='store_true', help='Write force file.')
+    parser.add_argument(
+        '--wrt_residual', action='store_true', help='Write residual file.')
+    parser.add_argument(
+        '--in_cluster', action='store_true', help='launch as a job in cluster.')
 
     ops = parser.parse_args()
+    if not ops.relax_ic_in_csf:
+        if ops.relax_ic_coeff is not None:
+            raise ValueError(
+                'Cannot use relax_ic_coeff when not using --relax_ic_in_csf')
     if ops.use_dti:
         if ops.in_dti is None:
             ops.in_dti = raw_input('Enter a DTI file (since you are using '
@@ -57,53 +66,48 @@ def main():
     """
     Run the model adlem for a given patient, segmentation mask, atrophy map and
     [optional DTI]."""
-    ops = get_input_options()
     work_dir = os.getenv('ADLEM_DIR')
     if work_dir is None:
         raise ValueError('environment variable ADLEM_DIR not set.')
     target = op.join(work_dir, 'build/src/AdLemMain')
+
     ops = get_input_options()
     res_dir = op.join(work_dir, 'results/patients', ops.patient)
     in_img = op.join(res_dir, ops.in_img)
     atrophy = op.join(res_dir, ops.atrophy)
     in_seg = op.join(res_dir, ops.in_seg)
     res_path = res_dir+'/' #due to stupid limitation in my AdLem code
     # where it concatenates simply the -resPath argument with the file prefix.
-    if ops.use_dti is True:
-        use_dti = 'true'
-        in_dti = op.join(res_dir, ops.in_dti)
-    else:
-        use_dti = 'false'
-        in_dti = "dummyFileName"
-    if ops.wrt_press is True:
-        wrt_press = 'true'
-    else:
-        wrt_press = 'false'
-    if ops.wrt_force is True:
-        wrt_force = 'true'
-    else:
-        wrt_force = 'false'
-    if ops.wrt_residual is True:
-        wrt_residual = 'true'
-    else:
-        wrt_residual = 'false'
-    if ops.mu_file:
-        mu_file = op.join(res_dir, ops.mu_file)
-        mu_ops = '-parameters "%s" -muFile %s' % (ops.lame_paras, mu_file)
-    else:
-        mu_ops = '-parameters "%s"' % (ops.lame_paras)
     with open(ops.petsc_op_file, 'r') as fil:
-        # Join all lines in the file that does not start with '#'
-        # with a whitespace separator.
-        petsc_ops = ' '.join(
-            line.strip() for line in fil if not line.startswith('#'))
-    cmd = ('%s %s -atrophyFile %s -maskFile %s -imageFile %s -useTensorLambda '
-           '%s -lambdaFile %s -numOfTimeSteps %s -resPath %s '
-           '-resultsFilenamesPrefix %s -writePressure %s -writeForce %s '
-           '-writeResidual %s %s'
-           % (target, mu_ops, atrophy, in_seg, in_img, use_dti, in_dti,
-              ops.time_steps, res_path, ops.res_prefix, wrt_press, wrt_force,
-              wrt_residual, petsc_ops))
+        # Get all lines in the file that does not start with '#'
+        # into a list
+        petsc_ops = [
+            line.strip() for line in fil if not line.startswith('#')]
+
+    optional_args, bool_args = [], []
+    if ops.relax_ic_in_csf:
+        bool_args.append('--relax_ic_in_csf')
+        if ops.relax_ic_coeff is not None:
+            optional_args.append('-relax_ic_coeff ' + ops.relax_ic_coeff)
+    if ops.mu_file is not None:
+        optional_args.append('-muFile ' + ops.mu_file)
+    if ops.use_dti:
+        bool_args.append('--useTensorLambda')
+        optional_args.append('-lambdaFile ' + op.join(res_dir, ops.in_dti))
+    if ops.wrt_press:
+        bool_args.append('--writePressure')
+    if ops.wrt_force:
+        bool_args.append('--writeForce')
+    if ops.wrt_residual:
+        bool_args.append('--writeResidual')
+
+    cmd = ('%s -parameters "%s" -boundary_condition %s -atrophyFile %s '
+           '-maskFile %s -imageFile %s -numOfTimeSteps %s -resPath %s '
+           '-resultsFilenamesPrefix %s %s'
+           % (target, ops.lame_paras, ops.boundary_condition, atrophy,
+              in_seg, in_img, ops.time_steps, res_path, ops.res_prefix,
+              ' '.join(bool_args + optional_args + petsc_ops)))
+
     if ops.in_cluster is True:
         cluster_mpi = '/opt/openmpi-gcc/current/bin/mpiexec '
         cmd = bu.sophia_nef_pbs_setting() + cluster_mpi + cmd