diff --git a/src/bindings/c/fetal.c b/src/bindings/c/fetal.c index 76953c0..de20104 100644 --- a/src/bindings/c/fetal.c +++ b/src/bindings/c/fetal.c @@ -3,13 +3,13 @@ #include -void fetal_model_c(const char *OUTDIR, int *filename_len, double *dt,int *num_heart_beats,double *T_beat,double *T_vs,double *T_as,double *T_v_delay,double *U0RV,double *EsysRV,double *EdiaRV,double *RvRv,double *U0LV,double *EsysLV,double *EdiaLV,double *RvLV,double *U0A, double *V0V, double *V0A); +void fetal_model_c(const char *OUTDIR, int *filename_len, double *dt,int *num_heart_beats,double *T_beat,double *T_vs,double *T_as,double *T_v_delay,double *U0RV,double *EsysRV,double *EdiaRV,double *RvRv,double *U0LV,double *EsysLV,double *EdiaLV,double *RvLV,double *U0A, double *V0V, double *V0A, int *use_plac_model); void assign_fetal_arrays_c(); -void fetal_model(const char *OUTDIR, double dt,int num_heart_beats,double T_beat,double T_vs,double T_as,double T_v_delay,double U0RV,double EsysRV,double EdiaRV,double RvRv,double U0LV,double EsysLV,double EdiaLV,double RvLV,double U0A, double V0V, double V0A) +void fetal_model(const char *OUTDIR, double dt,int num_heart_beats,double T_beat,double T_vs,double T_as,double T_v_delay,double U0RV,double EsysRV,double EdiaRV,double RvRv,double U0LV,double EsysLV,double EdiaLV,double RvLV,double U0A, double V0V, double V0A, int use_plac_model) { int filename_len = strlen(OUTDIR); - fetal_model_c(OUTDIR, &filename_len,&dt,&num_heart_beats,&T_beat,&T_vs,&T_as,&T_v_delay,&U0RV,&EsysRV,&EdiaRV,&RvRv,&U0LV,&EsysLV,&EdiaLV,&RvLV,&U0A,&V0V,&V0A); + fetal_model_c(OUTDIR, &filename_len,&dt,&num_heart_beats,&T_beat,&T_vs,&T_as,&T_v_delay,&U0RV,&EsysRV,&EdiaRV,&RvRv,&U0LV,&EsysLV,&EdiaLV,&RvLV,&U0A,&V0V,&V0A,&use_plac_model); } void assign_fetal_arrays() diff --git a/src/bindings/c/fetal.f90 b/src/bindings/c/fetal.f90 index 3adcff6..ce225c1 100644 --- a/src/bindings/c/fetal.f90 +++ b/src/bindings/c/fetal.f90 @@ -8,7 +8,7 @@ module fetal_c ! !> Perfusion fetal subroutine fetal_model_c(OUTDIR,filename_len,dt,num_heart_beats,T_beat,T_vs,T_as,T_v_delay,U0RV,EsysRV,EdiaRV,RvRv,& - U0LV,EsysLV,EdiaLV,RvLV,U0A,V0V,V0A) bind(C, name="fetal_model_c") + U0LV,EsysLV,EdiaLV,RvLV,U0A,V0V,V0A,use_plac_model) bind(C, name="fetal_model_c") use iso_c_binding, only: c_ptr use utils_c, only: strncpy use arrays, only: dp @@ -34,6 +34,7 @@ subroutine fetal_model_c(OUTDIR,filename_len,dt,num_heart_beats,T_beat,T_vs,T_as real(dp), intent(in) :: U0A real(dp), intent(in) :: V0V real(dp), intent(in) :: V0A + logical, intent(in) :: use_plac_model integer,intent(in) :: filename_len character(len=MAX_FILENAME_LEN) :: filename_f @@ -41,10 +42,10 @@ subroutine fetal_model_c(OUTDIR,filename_len,dt,num_heart_beats,T_beat,T_vs,T_as #if defined _WIN32 && defined __INTEL_COMPILER call so_fetal_model(filename_f,dt,num_heart_beats,T_beat,T_vs,T_as,T_v_delay,U0RV,EsysRV,EdiaRV,& - RvRv,U0LV,EsysLV,EdiaLV,RvLV,U0A,V0V,V0A) + RvRv,U0LV,EsysLV,EdiaLV,RvLV,U0A,V0V,V0A, use_plac_model) #else call fetal_model(filename_f,dt,num_heart_beats,T_beat,T_vs,T_as,T_v_delay,U0RV,EsysRV,EdiaRV,RvRv,& - U0LV,EsysLV,EdiaLV,RvLV,U0A,V0V,V0A) + U0LV,EsysLV,EdiaLV,RvLV,U0A,V0V,V0A, use_plac_model) #endif end subroutine fetal_model_c diff --git a/src/bindings/c/fetal.h b/src/bindings/c/fetal.h index 3df689a..32d156e 100644 --- a/src/bindings/c/fetal.h +++ b/src/bindings/c/fetal.h @@ -2,7 +2,7 @@ #define REPROSIM_FETAL_H #include "symbol_export.h" -SHO_PUBLIC void fetal_model(const char *OUTDIR, double dt, int num_heart_beats,double T_beat,double T_vs,double T_as,double T_v_delay,double U0RV,double EsysRV,double EdiaRV,double RvRv,double U0LV,double EsysLV,double EdiaLV,double RvLV,double U0A,double V0V, double V0A); +SHO_PUBLIC void fetal_model(const char *OUTDIR, double dt, int num_heart_beats,double T_beat,double T_vs,double T_as,double T_v_delay,double U0RV,double EsysRV,double EdiaRV,double RvRv,double U0LV,double EsysLV,double EdiaLV,double RvLV,double U0A,double V0V, double V0A, int use_plac_model); SHO_PUBLIC void assign_fetal_arrays(); #endif /* REPROSIM_FETAL_H */ \ No newline at end of file diff --git a/src/lib/fetal.f90 b/src/lib/fetal.f90 index dc3b721..79e16aa 100644 --- a/src/lib/fetal.f90 +++ b/src/lib/fetal.f90 @@ -38,7 +38,7 @@ module fetal contains subroutine fetal_model(OUTDIR,dt,num_heart_beats,T_beat,T_vs,T_as,T_v_delay,U0RV,EsysRV,EdiaRV,RvRv,U0LV,EsysLV,EdiaLV,& - RvLV,U0A,V0V,V0A) + RvLV,U0A,V0V,V0A,use_plac_model) use diagnostics, only: enter_exit,get_diagnostics_level use other_consts, only: MAX_FILENAME_LEN, MAX_STRING_LEN @@ -75,6 +75,10 @@ subroutine fetal_model(OUTDIR,dt,num_heart_beats,T_beat,T_vs,T_as,T_v_delay,U0RV real(dp) :: Aatria !Atrial activation (no units) real(dp) :: dpress,Pgrad,Qnod,dQ,Vnod,press real(dp) :: art_resistance,ven_resistance,total_volume + + logical, intent(in) :: use_plac_model !Boolean variable to control whether or not an anatomic model is used for the + !fetal compartment + character(len=60) :: mesh_type logical :: continue character(len=60) :: sub_name @@ -117,20 +121,31 @@ subroutine fetal_model(OUTDIR,dt,num_heart_beats,T_beat,T_vs,T_as,T_v_delay,U0RV write(*,*) 'Total blood volume (ml):',total_volume/1000. - write(*,*) 'Calculating placental resistance' - mesh_type = 'simple_tree' - elem_field(ne_viscfact,:) = 1.0_dp !initialise viscosity factor - call calculate_resistance(0.33600e-02_dp,mesh_type) - call tree_resistance(art_resistance,ven_resistance) - write(*,*) 'Arterial resistance (Pa.s/mm3)= ', art_resistance - write(*,*) 'Venous resistance (Pa.s/mm3)= ', ven_resistance - do ne =1,num_elems_fetal - if (abs(elem_field_fetal(ne_group,ne)-9.0_dp).lt.loose_tol)then!Umbilical artery - elem_field_fetal(ne_resist,ne) = art_resistance ! Pa s /mm3 - elseif(abs(elem_field_fetal(ne_group,ne)-10.0_dp).lt.loose_tol)then!Umbilical vein - elem_field_fetal(ne_resist,ne) = ven_resistance ! Pa s /mm3 - end if - end do + if (use_plac_model.eqv..TRUE.) then + write(*,*) 'Calculating placental resistance' + mesh_type = 'simple_tree' + elem_field(ne_viscfact,:) = 1.0_dp !initialise viscosity factor + call calculate_resistance(0.33600e-02_dp,mesh_type) + call tree_resistance(art_resistance,ven_resistance) + write(*,*) 'Arterial resistance (Pa.s/mm3)= ', art_resistance + write(*,*) 'Venous resistance (Pa.s/mm3)= ', ven_resistance + do ne =1,num_elems_fetal + if (abs(elem_field_fetal(ne_group,ne)-9.0_dp).lt.loose_tol)then!Umbilical artery + elem_field_fetal(ne_resist,ne) = art_resistance ! Pa s /mm3 + elseif(abs(elem_field_fetal(ne_group,ne)-10.0_dp).lt.loose_tol)then!Umbilical vein + elem_field_fetal(ne_resist,ne) = ven_resistance ! Pa s /mm3 + end if + end do + end if + + ! print*, ne_group + ! print*, ne_resist + ! print*, nef_K + ! print*, nef_L + + ! call print_matrix(elem_field_fetal) + ! write(*,*) 'This is where Toby wants this subroutine to stop' + ! stop 0 Write(*,*) 'Initialising flows' !Initialise flows @@ -210,7 +225,7 @@ subroutine fetal_model(OUTDIR,dt,num_heart_beats,T_beat,T_vs,T_as,T_v_delay,U0RV - continue = .true. + continue = .TRUE. n = 0 do while (continue) n = n + 1 ! increment the heart beat number @@ -344,8 +359,8 @@ subroutine fetal_model(OUTDIR,dt,num_heart_beats,T_beat,T_vs,T_as,T_v_delay,U0RV node_field_fetal(njf_netQ,23) end do if (n.eq.num_heart_beats) then - continue = .false. - endif + continue = .FALSE. + end if end do !do while close(10) close(20) @@ -828,4 +843,16 @@ subroutine tree_resistance(art_resistance,ven_resistance) call enter_exit(sub_name,2) end subroutine tree_resistance +subroutine print_matrix(A) + real(dp), intent(in) :: A(:,:) ! An assumed-shape dummy argument + + integer :: i + + do i = 1, size(A,2) + print'(F6.3,$)', A(:,i) + print*, '' + end do + + end subroutine print_matrix + end module fetal \ No newline at end of file