Merge pull request #127 from tawhai/unit_ventilation

Unit ventilation
LungNoodle · Aug 14, 2022 · 39065f6 · 39065f6
2 parents b5222c1 + 03046ab
commit 39065f6
Show file tree

Hide file tree

Showing 27 changed files with 1,770 additions and 812 deletions.
diff --git a/src/bindings/c/arrays.c b/src/bindings/c/arrays.c
@@ -1,9 +1,18 @@
 
 #include "arrays.h"
+#include "string.h"
 
 void set_node_field_value_c(int *row, int *col, double *value);
+void update_parameter_c(const char *parameter_name, int *parameter_name_len, double *parameter_value);
 
 void set_node_field_value(int row, int col, double value)
 {
   set_node_field_value_c(&row, &col, &value);
 }
+
+void update_parameter(const char *parameter_name, double parameter_value)
+{
+  int parameter_name_len = (int)strlen(parameter_name);
+  update_parameter_c(parameter_name, &parameter_name_len, &parameter_value);
+}
+
diff --git a/src/bindings/c/arrays.f90 b/src/bindings/c/arrays.f90
@@ -12,7 +12,7 @@ module arrays_c
   use arrays,only: dp
 
   implicit none
-  public set_node_field_value_c
+  public set_node_field_value_c,update_parameter_c
 
 contains
   subroutine set_node_field_value_c(row, col, value) bind(C, name="set_node_field_value_c")
@@ -26,5 +26,26 @@ subroutine set_node_field_value_c(row, col, value) bind(C, name="set_node_field_
 
   end subroutine set_node_field_value_c
 
+  subroutine update_parameter_c(parameter_name, parameter_name_len, parameter_value) bind(C, name="update_parameter_c")
+    use iso_c_binding, only: c_ptr
+    use utils_c, only: strncpy
+    use other_consts, only: max_filename_len
+    use arrays, only: update_parameter
+    implicit none
+
+    integer, intent(in) :: parameter_name_len
+    real(dp), intent(in) :: parameter_value
+    type(c_ptr),value, intent(in) :: parameter_name
+    character(len=max_filename_len) :: parameter_name_f
+
+    call strncpy(parameter_name_f, parameter_name, parameter_name_len)
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_update_parameter(parameter_name_f, parameter_value)
+#else
+    call update_parameter(parameter_name_f, parameter_value)
+#endif
+  end subroutine update_parameter_c
+
+
 
 end module arrays_c
diff --git a/src/bindings/c/arrays.h b/src/bindings/c/arrays.h
@@ -5,5 +5,6 @@
 #include "symbol_export.h"
 
 SHO_PUBLIC void set_node_field_value(int row, int col, double value);
+SHO_PUBLIC void update_parameter(const char *parameter_name, double parameter_value);
 
 #endif /* AETHER_ARRAYS_H */
diff --git a/src/bindings/c/exports.c b/src/bindings/c/exports.c
@@ -4,29 +4,54 @@
 
 #include <string.h>
 
+void export_cubic_lagrange_2d_c(const char *EXFILE, int *exfile_len,
+				const char *group_name, int *group_name_len);
 void export_1d_elem_field_c(int *ne_field, const char *EXELEMFILE, int *EXELEMFILE_LEN,
-                            const char *group_name, int *group_name_len, const char *field_name, int *field_name_len );
-void export_1d_elem_geometry_c(const char *EXELEMFILE, int *EXELEMFILE_LEN, const char *name, int *name_len);
-void export_elem_geometry_2d_c(const char *EXELEMFILE, int *EXELEMFILE_LEN, const char *name, int *name_len, int *offset_elem, int *offset_node);
+                            const char *group_name, int *group_name_len,
+			    const char *field_name, int *field_name_len );
+void export_1d_elem_geometry_c(const char *EXELEMFILE, int *EXELEMFILE_LEN,
+			       const char *name, int *name_len);
+void export_elem_geometry_2d_c(const char *EXELEMFILE, int *EXELEMFILE_LEN,
+			       const char *name, int *name_len, int *offset_elem, int *offset_node);
 void export_node_field_c(int *nj_field, const char *EXNODEFIELD, int *EXNODEFIELD_LEN,
                          const char *name, int *name_len, const char *field_name, int *field_name_len);
-void export_elem_geometry_2d_c(const char *EXELEMFILE, int *EXELEMFILE_LEN, const char *name, int *name_len, int *offset_elem, int *offset_node);
-void export_terminal_solution_c(const char *EXNODEFILE, int *EXNODEFILE_LEN, const char *name, int *name_len);
-void export_terminal_perfusion_c(const char *EXNODEFILE, int *EXNODEFILE_LEN, const char *name, int *name_len);
-void export_node_geometry_c(const char *EXNODEFILE, int *EXNODEFILE_LEN, const char *name, int *name_len);
-void export_node_geometry_2d_c(const char *EXNODEFILE, int *EXNODEFILE_LEN, const char *name, int *name_len, int *offset);
-void export_data_geometry_c(const char *EXDATAFILE, int *EXDATAFILE_LEN, const char *name, int *name_len, int *offset);
+void export_elem_geometry_2d_c(const char *EXELEMFILE, int *EXELEMFILE_LEN,
+			       const char *name, int *name_len, int *offset_elem, int *offset_node);
+void export_terminal_solution_c(const char *EXNODEFILE, int *EXNODEFILE_LEN,
+				const char *name, int *name_len);
+void export_terminal_perfusion_c(const char *EXNODEFILE, int *EXNODEFILE_LEN,
+				 const char *name, int *name_len);
+void export_node_geometry_c(const char *EXNODEFILE, int *EXNODEFILE_LEN,
+			    const char *name, int *name_len);
+void export_node_geometry_2d_c(const char *EXNODEFILE, int *EXNODEFILE_LEN,
+			       const char *name, int *name_len, int *offset);
+void export_data_geometry_c(const char *EXDATAFILE, int *EXDATAFILE_LEN,
+			    const char *name, int *name_len, int *offset);
 void export_elem_field_c(const char *EXELEMFIELD, int *EXELEMFIELD_LEN,
-                         const char *name, int *name_len, const char *field_name, int *field_name_len);
-void export_terminal_ssgexch_c(const char *EXNODEFILE, int *filename_len, const char *name, int *name_len);
+                         const char *name, int *name_len, const char *field_name,
+			 int *field_name_len);
+void export_terminal_ssgexch_c(const char *EXNODEFILE, int *filename_len,
+			       const char *name, int *name_len);
+
+
+
+void export_cubic_lagrange_2d(const char *EXFILE, const char *group_name)
+{
+  int filename_len = strlen(EXFILE);
+  int group_name_len = strlen(group_name);
+
+  export_cubic_lagrange_2d_c(EXFILE, &filename_len, group_name, &group_name_len);
+}
 
-void export_1d_elem_field(int ne_field, const char *EXELEMFILE, const char *group_name, const char *field_name )
+void export_1d_elem_field(int ne_field, const char *EXELEMFILE,
+			  const char *group_name, const char *field_name )
 {
   int filename_len = strlen(EXELEMFILE);
   int group_name_len = strlen(group_name);
   int field_name_len = strlen(field_name);
 
-  export_1d_elem_field_c(&ne_field, EXELEMFILE, &filename_len, group_name, &group_name_len, field_name, &field_name_len);
+  export_1d_elem_field_c(&ne_field, EXELEMFILE, &filename_len, group_name,
+			 &group_name_len, field_name, &field_name_len);
 }
 
 void export_1d_elem_geometry(const char *EXELEMFILE, const char *name)
@@ -37,21 +62,25 @@ void export_1d_elem_geometry(const char *EXELEMFILE, const char *name)
   export_1d_elem_geometry_c(EXELEMFILE, &filename_len, name, &name_len);
 }
 
-void export_elem_geometry_2d(const char *EXELEMFILE, const char *name, int offset_elem, int offset_node)
+void export_elem_geometry_2d(const char *EXELEMFILE, const char *name,
+			     int offset_elem, int offset_node)
 {
   int filename_len = strlen(EXELEMFILE);
   int name_len = strlen(name);
 
-  export_elem_geometry_2d_c(EXELEMFILE, &filename_len, name, &name_len, &offset_elem, &offset_node);
+  export_elem_geometry_2d_c(EXELEMFILE, &filename_len, name, &name_len,
+			    &offset_elem, &offset_node);
 }
 
-void export_node_field(int nj_field, const char *EXNODEFIELD, const char *name, const char *field_name)
+void export_node_field(int nj_field, const char *EXNODEFIELD,
+		       const char *name, const char *field_name)
 {
   int filename_len = strlen(EXNODEFIELD);
   int name_len = strlen(name);
   int field_name_len = strlen(field_name);
 
-  export_node_field_c(&nj_field, EXNODEFIELD, &filename_len, name, &name_len, field_name, &field_name_len);
+  export_node_field_c(&nj_field, EXNODEFIELD, &filename_len, name,
+		      &name_len, field_name, &field_name_len);
 }
 
 void export_terminal_solution(const char *EXNODEFILE, const char *name)
@@ -108,6 +137,7 @@ void export_elem_field(const char *EXELEMFIELD, const char *name, const char *fi
   int name_len = strlen(name);
   int field_name_len = strlen(field_name);
 
-  export_elem_field_c(EXELEMFIELD, &filename_len, name, &name_len, field_name, &field_name_len);
+  export_elem_field_c(EXELEMFIELD, &filename_len, name, &name_len,
+		      field_name, &field_name_len);
 }
 
diff --git a/src/bindings/c/exports.f90 b/src/bindings/c/exports.f90
@@ -4,6 +4,34 @@ module exports_c
   private
 
 contains
+
+
+!!!################################################################
+
+  subroutine export_cubic_lagrange_2d_c(EXFILE, filename_len, group_name, group_name_len) &
+    bind(C, name="export_cubic_lagrange_2d_c")
+
+    use iso_c_binding, only: c_ptr
+    use utils_c, only: strncpy
+    use exports, only: export_cubic_lagrange_2d
+    use other_consts, only: MAX_STRING_LEN, MAX_FILENAME_LEN
+    implicit none
+    integer,intent(in) :: filename_len, group_name_len
+    type(c_ptr), value, intent(in) :: EXFILE, group_name
+    character(len=MAX_FILENAME_LEN) :: filename_f
+    character(len=MAX_STRING_LEN) :: group_name_f
+
+    call strncpy(filename_f, EXFILE, filename_len)
+    call strncpy(group_name_f, group_name, group_name_len)
+
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_export_cubic_lagrange_2d(filename_f, group_name_f)
+#else
+    call export_cubic_lagrange_2d(filename_f, group_name_f)
+#endif
+
+  end subroutine export_cubic_lagrange_2d_c
+
 !!!################################################################
 
   subroutine export_1d_elem_field_c(ne_field, EXELEMFILE, filename_len, group_name, group_name_len, field_name, field_name_len) &

diff --git a/src/bindings/c/exports.h b/src/bindings/c/exports.h
@@ -4,6 +4,7 @@
 
 #include "symbol_export.h"
 
+SHO_PUBLIC void export_cubic_lagrange_2d(const char *EXFILE, const char *group_name);
 SHO_PUBLIC void export_1d_elem_field(int ne_field, const char *EXELEMFILE, const char *group_name, const char *field_name );
 SHO_PUBLIC void export_1d_elem_geometry(const char *EXELEMFILE, const char *name);
 SHO_PUBLIC void export_elem_geometry_2d(const char *EXELEMFILE, const char *name, int offset_elem, int offset_node);

diff --git a/src/bindings/c/gas_exchange.c b/src/bindings/c/gas_exchange.c
@@ -1,14 +1,18 @@
 
 #include "gas_exchange.h"
 
-void steadystate_gasexchange_c(double *c_art_o2, double *c_ven_o2,
-	double *p_art_co2, double *p_art_o2, double *p_i_o2, double *p_ven_co2, double *p_ven_o2, double *shunt_fraction,
-	double *VCO2, double *VO2);
+void initial_gasexchange_c(double *initial_concentration, double *surface_area, double *V_cap);
 
-void steadystate_gasexchange(double *c_art_o2, double *c_ven_o2,
-	double *p_art_co2, double *p_art_o2, double *p_i_o2, double *p_ven_co2, double *p_ven_o2, double *shunt_fraction,
-	double *VCO2, double *VO2)
+void steadystate_gasexchange_c(double *deadspace, double *p_i_o2, double *shunt_fraction, double *target_p_art_co2, double *target_p_ven_o2, double *VCO2, double *VO2);
+
+
+void initial_gasexchange(double initial_concentration, double surface_area, double V_cap)
+{
+  initial_gasexchange_c(&initial_concentration, &surface_area, &V_cap);
+}
+
+void steadystate_gasexchange(double deadspace, double p_i_o2, double shunt_fraction, double target_p_art_co2, double target_p_ven_o2, double VCO2, double VO2)
 {
-	steadystate_gasexchange_c(c_art_o2, c_ven_o2, p_art_co2, p_art_o2, p_i_o2, p_ven_co2, p_ven_o2, shunt_fraction, VCO2, VO2);
+  steadystate_gasexchange_c(&deadspace, &p_i_o2, &shunt_fraction, &target_p_art_co2, &target_p_ven_o2, &VCO2, &VO2);
 }
 
diff --git a/src/bindings/c/gas_exchange.f90 b/src/bindings/c/gas_exchange.f90
@@ -4,20 +4,40 @@ module gas_exchange_c
 
 contains
   !!!######################################################################
-  subroutine steadystate_gasexchange_c(c_art_o2,c_ven_o2,&
-       p_art_co2,p_art_o2,p_i_o2,p_ven_co2,p_ven_o2,shunt_fraction,&
-       VCO2,VO2) bind(C, name="steadystate_gasexchange_c")
+  subroutine initial_gasexchange_c(initial_concentration,surface_area,V_cap) &
+    bind(C, name="initial_gasexchange_c")
+    use gas_exchange,only: initial_gasexchange
+    use arrays,only: dp
+    implicit none
+    !!! Parameter List
+    real(dp),intent(in) :: initial_concentration,surface_area,V_cap
+
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_initial_gasexchange(initial_concentration,surface_area,V_cap)
+#else
+    call initial_gasexchange(initial_concentration,surface_area,V_cap)
+#endif
+
+  end subroutine initial_gasexchange_c
+
+!!!######################################################################
+  subroutine steadystate_gasexchange_c(deadspace,p_i_o2,shunt_fraction, &
+       target_p_art_co2,target_p_ven_o2,VCO2,VO2) bind(C, name="steadystate_gasexchange_c")
     use gas_exchange, only: steadystate_gasexchange
     use arrays,only: dp
     implicit none
 
     !!! Parameter List
-    real(dp),intent(in) :: p_i_o2,shunt_fraction,VCO2,VO2
-    real(dp), intent(inout) :: c_art_o2,c_ven_o2,p_art_co2,p_art_o2,p_ven_o2,p_ven_co2
-
-    call steadystate_gasexchange(c_art_o2,c_ven_o2,&
-       p_art_co2,p_art_o2,p_i_o2,p_ven_co2,p_ven_o2,shunt_fraction,&
-       VCO2,VO2)
+    real(dp),intent(in) :: deadspace,p_i_o2,shunt_fraction,target_p_art_co2, &
+         target_p_ven_o2,VCO2,VO2
+
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_steadystate_gasexchange(deadspace,p_i_o2,shunt_fraction, &
+       target_p_art_co2,target_p_ven_o2,VCO2,VO2)
+#else
+    call steadystate_gasexchange(deadspace,p_i_o2,shunt_fraction, &
+       target_p_art_co2,target_p_ven_o2,VCO2,VO2)
+#endif
 
   end subroutine steadystate_gasexchange_c
 

diff --git a/src/bindings/c/gas_exchange.h b/src/bindings/c/gas_exchange.h
@@ -3,6 +3,7 @@
 
 #include "symbol_export.h"
 
-SHO_PUBLIC void steadystate_gasexchange();
+SHO_PUBLIC void initial_gasexchange(double initial_concentration, double surface_area, double V_cap);
+SHO_PUBLIC void steadystate_gasexchange(double deadspace, double p_i_o2, double shunt_fraction, double target_p_art_co2, double target_p_ven_o2, double VCO2, double VO2);
 
 #endif /* AETHER_GAS_EXCHANGE_H */
diff --git a/src/bindings/c/geometry.c b/src/bindings/c/geometry.c
@@ -3,7 +3,7 @@
 
 #include "string.h"
 
-void add_mesh_c(const char *AIRWAY_MESHFILE, int *filename_len);
+void add_mesh_c(const char *AIRWAY_MESHFILE, int *filename_len, const char *BRANCHTYPE, int *branchtype_len, int *n_refine);
 void add_matching_mesh_c(void);
 void append_units_c(void);
 void define_1d_elements_c(const char *ELEMFILE, int *filename_len);
@@ -21,16 +21,19 @@ void define_rad_from_geom_c(const char *order_system, int *order_system_len, dou
                             const char *group_type, int *group_type_len, const char *group_options, int *group_options_len);
 void element_connectivity_1d_c(void);
 void evaluate_ordering_c(void);
-void set_initial_volume_c(int *Gdirn, double *COV, double *total_volume, double *Rmax, double *Rmin);
+void refine_1d_elements_c(int *elemlist_len, int elemlist[], int *nrefinements);
+void renumber_tree_in_order_c(void);
+void initialise_lung_volume_c(int *Gdirn, double *COV, double *total_volume, double *Rmax, double *Rmin);
 void volume_of_mesh_c(double *volume_model, double *volume_tree);
 void write_elem_geometry_2d_c(const char *ELEMFILE, int *filename_len);
 void write_geo_file_c(int *ntype, const char *GEOFILE, int *filename_len);
 void write_node_geometry_2d_c(const char *NODEFILE, int *filename_len);
 
-void add_mesh(const char *AIRWAY_MESHFILE)
+void add_mesh(const char *AIRWAY_MESHFILE, const char *BRANCHTYPE, int n_refine )
 {
   int filename_len = (int)strlen(AIRWAY_MESHFILE);
-  add_mesh_c(AIRWAY_MESHFILE, &filename_len);
+  int branchtype_len = (int)strlen(BRANCHTYPE);
+  add_mesh_c(AIRWAY_MESHFILE, &filename_len, BRANCHTYPE, &branchtype_len, &n_refine);
 }
 
 void add_matching_mesh()
@@ -127,9 +130,19 @@ void evaluate_ordering()
   evaluate_ordering_c();
 }
 
-void set_initial_volume(int Gdirn, double COV, double total_volume, double Rmax, double Rmin)
+void refine_1d_elements(int elemlist_len, int elemlist[], int nrefinements)
 {
-  set_initial_volume_c(&Gdirn, &COV, &total_volume, &Rmax, &Rmin);
+  refine_1d_elements_c(&elemlist_len, elemlist, &nrefinements);
+}
+
+void renumber_tree_in_order()
+{
+  renumber_tree_in_order_c();
+}
+
+void initialise_lung_volume(int Gdirn, double COV, double total_volume, double Rmax, double Rmin)
+{
+  initialise_lung_volume_c(&Gdirn, &COV, &total_volume, &Rmax, &Rmin);
 }
 
 void volume_of_mesh(double *volume_model, double *volume_tree)