Merge branch 'main' of https://github.com/simoneponcioni/HFE

.github/ISSUE_TEMPLATE/bug_report.md

@@ -0,0 +1,38 @@
+---
+name: Bug report
+about: Create a report to help us improve
+title: ''
+labels: ''
+assignees: ''
+
+---
+
+**Describe the bug**
+A clear and concise description of what the bug is.
+
+**To Reproduce**
+Steps to reproduce the behavior:
+1. Go to '...'
+2. Click on '....'
+3. Scroll down to '....'
+4. See error
+
+**Expected behavior**
+A clear and concise description of what you expected to happen.
+
+**Screenshots**
+If applicable, add screenshots to help explain your problem.
+
+**Desktop (please complete the following information):**
+ - OS: [e.g. iOS]
+ - Browser [e.g. chrome, safari]
+ - Version [e.g. 22]
+
+**Smartphone (please complete the following information):**
+ - Device: [e.g. iPhone6]
+ - OS: [e.g. iOS8.1]
+ - Browser [e.g. stock browser, safari]
+ - Version [e.g. 22]
+
+**Additional context**
+Add any other context about the problem here.

.gitignore

@@ -21,6 +21,7 @@ outputs/
 !04_SIMULATIONS/simulations_directory.text
 99_TEMP/
 QMSKI/DSC_PAPER/
+out/
 
 multirun/
 

02_CODE/cfg/hfe-repro.yaml

@@ -51,11 +51,11 @@ loadcase:
   time_for_displacement: 1
   min_step_size: 0.0000000001
   max_step_size: 0.3
-  load_displacement: -0.35           # (mm)
+  load_displacement: -0.35          # (mm)
 
 abaqus:
   nlgeom: on                        # affects NLGEOM parameter in abaqus simulation (on/off)
-  abaqus_nprocs: 8                  # (-)
+  abaqus_nprocs: 4                  # (-)
   abaqus_memory: 6000               # (MB)
   delete_odb: False                 # delete odb after retrieving data
   max_increments: 1000              # int number [1-1000]

02_CODE/cfg/hfe-tibia.yaml

@@ -21,7 +21,7 @@ mesher:
   air_elements: False               # True #if you want keep the elements to have a fullbock mesh in x-y direction 
 
 image_processing:
-  origaim_separate: False           # For Hosseini Dataset, Image parameters are read from original aim, not from processed BMD file, as there they were deleted by medtool pre-processing
+  origaim_separate: False           # For Hosseini Dataset, Image parameters are read from original aim, not from processed BMD file, there they were deleted by medtool pre-processing
   mask_separate: True               # Standard evaluation gives two separate mask file (CORTMASK and TRABMASK), False if one maskfile with two different labels
   imtype: NATIVE                    # NATIVE or BMD
   bvtv_scaling: 1                   # 0: no scaling, 1: Scaling of BVTV 61um to BVTV 11.4um [Hosseini2017]
@@ -51,7 +51,7 @@ loadcase:
   time_for_displacement: 1
   min_step_size: 0.0000000001
   max_step_size: 0.3
-  load_displacement: -1.0           # (mm)
+  load_displacement: -0.35          # (mm)
 
 abaqus:
   nlgeom: on                        # affects NLGEOM parameter in abaqus simulation (on/off)

02_CODE/cfg/mesh.yaml

@@ -12,7 +12,7 @@ meshing_settings:
   outside_val: 1                      # threshold value for the outside of the mask
   lower_thresh: 0                     # lower threshold for the mask
   upper_thresh: 0.9                   # upper threshold for the mask
-  s: 500                              # smoothing factor of the spline
+  s: 400                              # smoothing factor of the spline
   k: 3                                # degree of the spline
   interp_points: 200                  # number of points to interpolate the spline
   dp_simplification_outer: 3          # Ramer-Douglas-Peucker simplification factor for the periosteal contour
@@ -23,11 +23,12 @@ meshing_settings:
   center_square_length_factor: 0.4    # size ratio of the refinement square: 0 < l_f < 1
   mesh_order: 1                       # set order of the mesh (1: linear, 2: quadratic)
   sweep_factor: 1                     # factor for the sweep used in hydra for the sensitivity analysis
-  n_elms_longitudinal: 3
-  n_elms_transverse_trab: 15
-  n_elms_transverse_cort: 3
+  n_elms_longitudinal: 3              
+  n_elms_transverse_trab: 15          
+  n_elms_transverse_cort: 3           
   n_elms_radial: 20                   # ! Should be 10 if trab_refinement is True
   ellipsoid_fitting: True
+
   show_plots: False                   # show plots during construction
   show_gmsh: False                    # show gmsh GUI
   write_mesh: True                    # write mesh to file

02_CODE/cfg/paths-repro.yaml

@@ -1,8 +1,8 @@
 paths:
   origaimdir: 00_ORIGAIM/REPRO/
   aimdir:  01_DATA/REPRO/
-  feadir:  04_SIMULATIONS/REPRO
-  sumdir: 05_SUMMARIES/REPRO
+  feadir:  04_SIMULATIONS/REPRO_PAPER
+  sumdir: 05_SUMMARIES/REPRO_PAPER
   commondir: 01_DATA/BATCH/FEA_noReg/IMAGES/                        # if registration is true -> path to registration files
   folder_bc_psl_loadcases: 02_CODE/abq/BC_PSL/
   boundary_conditions: 02_CODE/abq/BC_PSL/boundary_conditions.inp

02_CODE/cfg/paths-tibia.yaml

@@ -1,7 +1,7 @@
 paths:
   origaimdir: 00_ORIGAIM/TIBIA
   aimdir:  01_DATA/TIBIA/
-  feadir:  04_SIMULATIONS/TIBIA
+  feadir:  04_SIMULATIONS/TIBIA_P
   sumdir: 05_SUMMARIES/TIBIA
   commondir: 01_DATA/BATCH/FEA_noReg/IMAGES/                        # if registration is true -> path to registration files
   folder_bc_psl_loadcases: 02_CODE/abq/BC_PSL/
@@ -20,6 +20,6 @@ filenames:
 
 version:
   verification_files: 1
-  current_version: 04_strain
+  current_version: 03_testmpi
   site_bone: Tibia                                                  # Radius / Tibia
 

02_CODE/docker_apptainer_hpc/build_container.sh

@@ -0,0 +1,24 @@
+#!/bin/bash
+
+# User info
+#SBATCH --mail-user=simone.poncioni@unibe.ch
+#SBATCH --mail-type=begin,end,fail
+
+# Job name
+#SBATCH --job-name="build_container"
+
+# Runtime and memory
+#SBATCH --time=01:00:00
+#SBATCH --cpus-per-task=10
+#SBATCH --mem-per-cpu=4G
+#SBATCH --tmp=64G
+
+# Workdir
+#SBATCH --chdir=/storage/workspaces/artorg_msb/hpc_abaqus/poncioni/apptainer
+#SBATCH --out=%x.out
+#SBATCH --err=%x.err
+
+
+# Run command
+srun apptainer build --force hfe_development_ifort.sif docker://simoneponcioni/hfe_development:latest
+

02_CODE/docker_apptainer_hpc/submit_hfe.sh

@@ -0,0 +1,43 @@
+#!/bin/bash
+
+# User info
+#SBATCH --mail-user=simone.poncioni@unibe.ch
+#SBATCH --mail-type=begin,end,fail
+
+# Job name
+#SBATCH --job-name="hfe_pipeline"
+
+# Runtime and memory
+#SBATCH --time=03:00:00
+#SBATCH --cpus-per-task=16
+#SBATCH --ntasks-per-node=1
+#SBATCH --mem-per-cpu=8G
+#SBATCH --tmp=100G
+#SBATCH --array=1-5%5
+
+# Workdir
+#SBATCH --chdir=/storage/workspaces/artorg_msb/hpc_abaqus/poncioni/HFE
+#SBATCH --output=out/hfe_%A_%a.out
+#SBATCH --error=out/hfe_%A_%a.err
+
+##############################################################################################################
+### Load modules
+HPC_WORKSPACE=hpc_abaqus module load Workspace
+
+unset SLURM_GTIDS
+
+### greyscale_filenames.txt contains lines with 1 greyscale_filename per line.
+greyscale_filenames=/storage/workspaces/artorg_msb/hpc_abaqus/poncioni/HFE/00_ORIGAIM/filenames.txt
+
+###    Line <i> contains greyscale_filename for run <i>
+# Get greyscale_filename                                                                                                                              
+greyscale_filename=$(cat $greyscale_filenames | awk -v var=$SLURM_ARRAY_TASK_ID 'NR==var {print $1}')
+
+
+### Zero pad the task ID to match the numbering of the input files
+n=$(printf "%04d" $SLURM_ARRAY_TASK_ID)
+
+# Run command
+srun apptainer exec /storage/workspaces/artorg_msb/hpc_abaqus/poncioni/apptainer/hfe_development_ifort.sif \
+/bin/bash -c "source /opt/intel/oneapi/setvars.sh && source /opt/miniconda/etc/profile.d/conda.sh && conda activate hfe-essentials && python 02_CODE/src/pipeline_runner.py simulations.grayscale_filenames=$greyscale_filename"
+

02_CODE/src/hfe_abq/simulation.py

@@ -1,8 +1,10 @@
 import logging
 import os
+import subprocess
 import sys
 import traceback
 from pathlib import Path
+import time
 
 LOGGING_NAME = "HFE-ACCURATE"
 logger = logging.getLogger(LOGGING_NAME)
@@ -11,7 +13,7 @@
 
 def simulate_loadcase(cfg, sample: str, inputfile: str, umat: str, loadcase: str):
     """
-    Run abaqus simulation from os subprocess.
+    Run Abaqus simulation from os subprocess.
 
     Args:
         config: configuration dictionary
@@ -23,6 +25,9 @@ def simulate_loadcase(cfg, sample: str, inputfile: str, umat: str, loadcase: str
     Returns:
         None
     """
+
+    setup_env_cmd = 'module load intel-compilers/2021.2.0'
+
     ABAQUS = cfg.solver.abaqus
     NPROCS = cfg.abaqus.abaqus_nprocs
     RAM = cfg.abaqus.abaqus_memory
@@ -40,19 +45,30 @@ def simulate_loadcase(cfg, sample: str, inputfile: str, umat: str, loadcase: str
     basepath = os.getcwd()
 
     os.chdir(simdir)
-    command = str(
-        "%s job=%s inp=%s cpus=%d memory=%d user=%s scratch=%s ask_delete=OFF -interactive"
-        % (ABAQUS, job, inputfile, NPROCS, RAM, umat, SCRATCH)
+    # command = (
+    #     f"{setup_env_cmd} && "
+    #     f"{ABAQUS} job={job} inp={inputfile} cpus={NPROCS} memory={RAM} "
+    #     f"user={umat} scratch={SCRATCH} ask_delete=OFF verbose=3 -interactive"
+    # )
+    command = (
+        f"{ABAQUS} job={job} inp={inputfile} cpus={NPROCS} memory={RAM} "
+        f"user={umat} scratch={SCRATCH} ask_delete=OFF verbose=3 -interactive"
     )
     logger.info(command)
     try:
-        os.system(command)
+        result = subprocess.run(command, capture_output=True, text=True, shell=True)
+        logger.info(result.stdout)
+        if result.returncode == 0:
+            logger.info("Abaqus simulation completed successfully")
+        else:
+            logger.error("Abaqus simulation failed with return code %d", result.returncode)
+            logger.error(result.stderr)
     except Exception as e:
         logger.error("Simulation of FZ_MAX loadcase resulted in error")
         logger.error(e)
         logger.error(traceback.format_exc())
-        logger.error(sys.stderr)
-        pass
-    os.chdir(basepath)
+    finally:
+        os.chdir(basepath)
+
     odb_path = simdir / (job + ".odb")
     return odb_path

02_CODE/src/hfe_abq/write_abaqus.py

@@ -557,7 +557,8 @@ def _write_output_history(self, REF_NODE: str):
             RF
             *NODE PRINT, NSET=REF_NODE, FREQUENCY=1, SUMMARY=NO
             U,
-            RF
+            RF,
+            CF,
             *End Step
         """
         return textwrap.dedent(out_hist)