diff --git a/.github/workflows/build_test.yml b/.github/workflows/build_test.yml index e2349d50..2e4eb2d9 100644 --- a/.github/workflows/build_test.yml +++ b/.github/workflows/build_test.yml @@ -24,15 +24,15 @@ jobs: os: windows-2022, } - { - name: "Ubuntu 20.04", + name: "Ubuntu 24.04", build_type: "Release", - os: ubuntu-20.04, + os: ubuntu-24.04, } - { - name: "Ubuntu 20.04 with SuperLU", + name: "Ubuntu 24.04 with SuperLU", build_type: "Release", super_lu: true, - os: ubuntu-20.04, + os: ubuntu-24.04, } - { name: "macOS Ventura", diff --git a/src/bindings/c/geometry.c b/src/bindings/c/geometry.c index e117ffb4..96d34a54 100644 --- a/src/bindings/c/geometry.c +++ b/src/bindings/c/geometry.c @@ -23,6 +23,7 @@ int get_local_node_f_c(const char *ndimension, int *dimension_len, const char *n void define_rad_from_geom_c(const char *order_system, int *order_system_len, double *control_param, const char *start_from, int *start_from_len, double *start_rad, const char *group_type, int *group_type_len, const char *group_options, int *group_options_len); +void occlude_vessel_c(int *VESSEL_NUMBER, double *RATIO); void element_connectivity_1d_c(void); void evaluate_ordering_c(void); void volume_of_mesh_c(double *volume_model, double *volume_tree); @@ -141,6 +142,13 @@ void define_rad_from_geom(const char *order_system, double control_param, const } +void occlude_vessel(int VESSEL_NUMBER, double RATIO) +{ + + occlude_vessel_c(&VESSEL_NUMBER, &RATIO); + +} + void element_connectivity_1d() { element_connectivity_1d_c(); diff --git a/src/bindings/c/geometry.f90 b/src/bindings/c/geometry.f90 index e6cd6bcd..7a44b732 100644 --- a/src/bindings/c/geometry.f90 +++ b/src/bindings/c/geometry.f90 @@ -4,7 +4,7 @@ module geometry_c use indices !use mesh_functions !use precision ! sets dp for precision - !use math_constants !pi + !use math_constants !pi implicit none private @@ -184,7 +184,7 @@ subroutine list_tree_statistics_c(filename, filename_len) & integer,intent(in) :: filename_len type(c_ptr), value, intent(in) :: filename character(len=MAX_FILENAME_LEN) :: filename_f - + call strncpy(filename_f, filename, filename_len) call list_tree_statistics(filename_f) @@ -194,34 +194,34 @@ end subroutine list_tree_statistics_c !################################################################################### ! subroutine internal_mesh_reorder_c() bind(C, name="internal_mesh_reorder_c") - + use geometry,only: internal_mesh_reorder implicit none call internal_mesh_reorder() end subroutine internal_mesh_reorder_c - + ! !################################################################################### ! subroutine make_data_grid_c(surface_elems_len, surface_elems, num_target, offset, spacing)& bind(C, name="make_data_grid_c") - + use arrays,only: dp use iso_c_binding, only: c_ptr use utils_c, only: strncpy use other_consts, only: MAX_FILENAME_LEN, MAX_STRING_LEN use geometry, only: make_data_grid implicit none - + integer,intent(in) :: surface_elems_len integer,intent(in) :: surface_elems(surface_elems_len) integer,intent(in) :: num_target real(dp),intent(in) :: offset, spacing - + call make_data_grid(surface_elems, num_target, offset, spacing) - + end subroutine make_data_grid_c !!!################################################################################### @@ -236,7 +236,7 @@ subroutine make_2d_vessel_from_1d_c(elemlist, elemlist_len) bind(C, name="make_2 call make_2d_vessel_from_1d(elemlist) end subroutine make_2d_vessel_from_1d_c - + ! !################################################################################### ! @@ -329,6 +329,26 @@ subroutine define_rad_from_geom_c(order_system, order_system_len, control_param, call define_rad_from_geom(order_system_f, control_param, start_from_f, start_rad, group_type_f, group_options_f) end subroutine define_rad_from_geom_c + ! + !################################################################################## + ! + !*define_rad_from_geom:* Defines vessel or airway radius based on their geometric structure + subroutine occlude_vessel_c(VESSEL_NUMBER, RATIO) bind(C, name="occlude_vessel_c") + + !use iso_c_binding, only: c_ptr + !use utils_c, only: strncpy + !use other_consts, only: MAX_STRING_LEN + use arrays, only: dp + use geometry, only: occlude_vessel + implicit none + + integer,intent(in) :: VESSEL_NUMBER + real(dp),intent(in) :: RATIO + !character(len=MAX_STRING_LEN) :: order_system_f, start_from_f, group_type_f, group_options_f + + call occlude_vessel(VESSEL_NUMBER, RATIO) + + end subroutine occlude_vessel_c ! !########################################################################### ! @@ -373,10 +393,10 @@ function get_local_node_f_c(ndimension,np_global) result(get_local_node) bind(C, use arrays, only: dp use geometry, only: get_local_node_f implicit none - + integer :: ndimension,np_global integer :: get_local_node - + get_local_node=get_local_node_f(ndimension,np_global) end function get_local_node_f_c @@ -443,5 +463,3 @@ end subroutine write_node_geometry_2d_c ! end module geometry_c - - diff --git a/src/bindings/c/geometry.h b/src/bindings/c/geometry.h index 255e0db4..ee87950d 100644 --- a/src/bindings/c/geometry.h +++ b/src/bindings/c/geometry.h @@ -23,6 +23,7 @@ SHO_PUBLIC void define_rad_from_file(const char *FIELDFILE, const char *radius_t SHO_PUBLIC int get_local_node_f(const char *ndimenstion, const char *np_global); SHO_PUBLIC void define_rad_from_geom(const char *ORDER_SYSTEM, double CONTROL_PARAM, const char *START_FROM, double START_RAD, const char *GROUP_TYPE, const char *GROUP_OPTIONS); +SHO_PUBLIC void occlude_vessel(int VESSEL_NUMBER, double RATIO); SHO_PUBLIC void element_connectivity_1d(); SHO_PUBLIC void evaluate_ordering(); SHO_PUBLIC void volume_of_mesh(double *volume_model, double *volume_tree); diff --git a/src/bindings/c/pressure_resistance_flow.c b/src/bindings/c/pressure_resistance_flow.c index 2a287b20..f44af941 100644 --- a/src/bindings/c/pressure_resistance_flow.c +++ b/src/bindings/c/pressure_resistance_flow.c @@ -1,8 +1,21 @@ #include "pressure_resistance_flow.h" #include +extern void compliance_list_c(int *elemlist2_len, int elemlist2[]); +extern void occlusion_list_c(int *elemlist_len, int elemlist[]); + void evaluate_prq_c(const char *mesh_type, int *mesh_type_len, const char *vessel_type, int *vessel_type_len,int *grav_dirn, double *grav_factor, const char *bc_type, int *bc_type_len, double *inlet_bc, double *outlet_bc, double *remodeling_grade); +void compliance_list(int elemlist2_len, int elemlist2[]) +{ + compliance_list_c(&elemlist2_len, elemlist2); +} + +void occlusion_list(int elemlist_len, int elemlist[]) +{ + occlusion_list_c(&elemlist_len, elemlist); +} + void evaluate_prq(const char *mesh_type, const char *vessel_type,int grav_dirn, double grav_factor, const char *bc_type, double inlet_bc, double outlet_bc, double remodeling_grade) { int mesh_type_len = strlen(mesh_type); diff --git a/src/bindings/c/pressure_resistance_flow.f90 b/src/bindings/c/pressure_resistance_flow.f90 index 91e0e1b9..f9c5a920 100644 --- a/src/bindings/c/pressure_resistance_flow.f90 +++ b/src/bindings/c/pressure_resistance_flow.f90 @@ -4,6 +4,47 @@ module pressure_resistance_flow_c contains + + + ! + !################################################################################### + ! + ! the main growing subroutine. Generates a volume-filling tree into a closed surface. + subroutine occlusion_list_c(surface_elems_len, surface_elems) bind(C, name="occlusion_list_c") + + !use arrays,only: dp + !use iso_c_binding, only: c_ptr + !use utils_c, only: strncpy + !use other_consts, only: MAX_FILENAME_LEN + use pressure_resistance_flow,only: occlusion_list + implicit none + + integer,intent(in) :: surface_elems_len + integer,intent(in) :: surface_elems(surface_elems_len) + + call occlusion_list(surface_elems) + + end subroutine occlusion_list_c + ! + !################################################################################### + ! + ! makes the vessels given rigid. + subroutine compliance_list_c(surface_elems_len, surface_elems) bind(C, name="compliance_list_c") + + !use arrays,only: dp + !use iso_c_binding, only: c_ptr + !use utils_c, only: strncpy + !use other_consts, only: MAX_FILENAME_LEN + use pressure_resistance_flow,only: compliance_list + implicit none + + integer,intent(in) :: surface_elems_len + integer,intent(in) :: surface_elems(surface_elems_len) + + call compliance_list(surface_elems) + + end subroutine compliance_list_c + !!!################################################################################### subroutine evaluate_prq_c(mesh_type,mesh_type_len,vessel_type,vessel_type_len,grav_dirn,grav_factor,bc_type,bc_type_len,inlet_bc, & diff --git a/src/bindings/c/pressure_resistance_flow.h b/src/bindings/c/pressure_resistance_flow.h index 729750d0..2e9864ca 100644 --- a/src/bindings/c/pressure_resistance_flow.h +++ b/src/bindings/c/pressure_resistance_flow.h @@ -3,7 +3,8 @@ #include "symbol_export.h" - +SHO_PUBLIC void occlusion_list(int elemlist_len, int elemlist[]); +SHO_PUBLIC void compliance_list(int elemlist2_len, int elemlist2[]); SHO_PUBLIC void evaluate_prq(const char *mesh_type, const char *vessel_type, int grav_dirn, double grav_factor, const char *bc_type, double inlet_bc, double outlet_bc, double remodeling_grade); diff --git a/src/bindings/interface/pressure_resistance_flow.i b/src/bindings/interface/pressure_resistance_flow.i index faa2dbf4..88736c88 100644 --- a/src/bindings/interface/pressure_resistance_flow.i +++ b/src/bindings/interface/pressure_resistance_flow.i @@ -1,8 +1,55 @@ %module(package="aether") pressure_resistance_flow %include symbol_export.h -%include pressure_resistance_flow.h + +%typemap(in) (int elemlist_len, int elemlist[]) { +int i; +if (!PyList_Check($input)) { + PyErr_SetString(PyExc_ValueError, "Expecting a list"); + SWIG_fail; +} +$1 = PyList_Size($input); +$2 = (int *) malloc(($1)*sizeof(int)); +for (i = 0; i < $1; i++) { + PyObject *o = PyList_GetItem($input, i); + if (!PyInt_Check(o)) { + free($2); + PyErr_SetString(PyExc_ValueError, "List items must be integers"); + SWIG_fail; + } + $2[i] = PyInt_AsLong(o); +} +} + +%typemap(in) (int elemlist2_len, int elemlist2[]) { +int i; +if (!PyList_Check($input)) { + PyErr_SetString(PyExc_ValueError, "Expecting a list"); + SWIG_fail; +} +$1 = PyList_Size($input); +$2 = (int *) malloc(($1)*sizeof(int)); +for (i = 0; i < $1; i++) { + PyObject *o = PyList_GetItem($input, i); + if (!PyInt_Check(o)) { + free($2); + PyErr_SetString(PyExc_ValueError, "List items must be integers"); + SWIG_fail; + } + $2[i] = PyInt_AsLong(o); +} +} + +%typemap(freearg) (int elemlist_len, int elemlist[]) { +if ($2) free($2); +} + +%typemap(freearg) (int elemlist2_len, int elemlist2[]) { +if ($2) free($2); +} %{ #include "pressure_resistance_flow.h" %} + +%include pressure_resistance_flow.h diff --git a/src/lib/geometry.f90 b/src/lib/geometry.f90 index ad8ed8c6..14870b97 100644 --- a/src/lib/geometry.f90 +++ b/src/lib/geometry.f90 @@ -37,6 +37,7 @@ module geometry public define_data_geometry public define_rad_from_file public define_rad_from_geom + public occlude_vessel public element_connectivity_1d public element_connectivity_2d public evaluate_ordering @@ -297,7 +298,7 @@ subroutine add_matching_mesh() do noelem=1,num_elems ne_m=elems(noelem) ne=ne_global+ne_m ! element num ordering for veins matchs artery - elem_field(ne_group,ne)=2.0_dp!VEIN + elem_field(ne_group,ne)=2.0_dp!VEIN elem_field(ne_group,ne_m)=0.0_dp!ARTERY elems(ne0+noelem)=ne if(.NOT.REVERSE)then @@ -323,7 +324,7 @@ subroutine add_matching_mesh() elem_cnct(-1,n,ne)=elem_cnct(1,n,ne_m)+ne0 enddo endif - + !if worrying about regions and versions do it here elems_at_node(elem_nodes(1,ne),0)=elems_at_node(elem_nodes(1,ne),0)+1 elems_at_node(elem_nodes(1,ne),elems_at_node(elem_nodes(1,ne),0))=ne @@ -336,7 +337,7 @@ subroutine add_matching_mesh() nindex=no_hord elem_ordrs(nindex,ne)=elem_ordrs(nindex,ne_m) enddo - + !update current no of nodes and elements to determine connectivity np0=np !current highest node ne1=maxval(elems) !current highest element @@ -569,7 +570,7 @@ subroutine define_1d_elements(ELEMFILE) call evaluate_ordering elem_ordrs(no_type,:) = 1 ! 0 for respiratory, 1 for conducting endif - + call enter_exit(sub_name,2) end subroutine define_1d_elements @@ -1257,7 +1258,7 @@ subroutine import_ply_triangles(ply_file) end subroutine import_ply_triangles !!!############################################################################# - + subroutine list_tree_statistics(filename) character(len=*),intent(in) :: filename @@ -1274,9 +1275,9 @@ subroutine list_tree_statistics(filename) logical :: add,writefile character(len=300) :: treefile character(len=60) :: sub_name - + ! -------------------------------------------------------------------------- - + sub_name = 'list_tree_statistics' call enter_exit(sub_name,1) @@ -1320,7 +1321,7 @@ subroutine list_tree_statistics(filename) num_ddp = 0 num_llp = 0 N = 0 - + do ne = 1,num_elems ne0 = elem_cnct(-1,1,ne) ! parent element number forall(i=1:3) ind(i) = elem_ordrs(i,ne) ! the gen, Hord, Sord for i=1,2,3 @@ -1332,7 +1333,7 @@ subroutine list_tree_statistics(filename) else if(ne0.ne.0.and.elem_ordrs(1,ne0).ne.ind(1))then ! gen not same as parent so add add = .true. endif - + if(add)then N = N + 1 nbranches(:,N) = ind(:) ! generation, H order, S order, D-D S order @@ -1341,7 +1342,7 @@ subroutine list_tree_statistics(filename) !!!...... Add length of all segments along branch, calculate their mean diameter n_segments=1 mean_diam = diameters(ne) - branches(1,N) = elem_field(ne_length,ne) + branches(1,N) = elem_field(ne_length,ne) ne_next = ne do while(elem_cnct(1,0,ne_next).eq.1) ! while a line of elements ne_next = elem_cnct(1,1,ne_next) !next element @@ -1388,7 +1389,7 @@ subroutine list_tree_statistics(filename) !!!... count the terminal branches in each generation if(elem_cnct(1,0,ne).eq.0) & ! this is a terminal element n_terminal(ind(1)) = n_terminal(ind(1)) + 1 ! ind(1) is element generation - + !!!... Calculate geometric properties of tree branches(4,N) = -1.0_dp !initialise to no rotation angle if(elem_cnct(-1,0,ne).gt.0.and.elem_cnct(1,0,ne).gt.1)then @@ -1407,11 +1408,11 @@ subroutine list_tree_statistics(filename) xp5(:) = node_xyz(:,np5) call make_plane_from_3points(norm_1,2,xp1,xp2,xp3) ! calculate unit normal and plane call make_plane_from_3points(norm_2,2,xp2,xp4,xp5) ! calculate unit normal and plane - branches(4,N) = angle_btwn_vectors(norm_1,norm_2)*180.0_dp/pi ! rotation angle + branches(4,N) = angle_btwn_vectors(norm_1,norm_2)*180.0_dp/pi ! rotation angle endif endif enddo ! ne - + n_br = N do ne = 1,num_elems ne0 = elem_cnct(-1,1,ne) ! parent element @@ -1427,7 +1428,7 @@ subroutine list_tree_statistics(filename) if(elem_cnct(1,0,ne).ge.2)then !'bi'furcations only ne1 = elem_cnct(1,1,ne) !first child ne2 = elem_cnct(1,2,ne) !second child - + !!!.... Summary statistics if(stats(6,ne1).ge.0.0_dp.and.stats(6,ne2).ge.0.0_dp)then if(stats(6,ne1).ge.stats(6,ne2))then !diameter classification @@ -1441,7 +1442,7 @@ subroutine list_tree_statistics(filename) stats(11,ne) = stats(2,ne_minor)*180.0_dp/pi stats(12,ne) = stats(2,ne_major)*180.0_dp/pi endif - + if(diameters(ne_minor).gt.0.0_dp.and.diameters(ne_major).gt.0.0_dp)then stats(13,ne) = stats(5,ne_minor)/diameters(ne_minor) !L/D minor stats(14,ne) = stats(5,ne_major)/diameters(ne_major) !L/D major @@ -1453,37 +1454,37 @@ subroutine list_tree_statistics(filename) endif endif ! elem_cnct enddo ! ne - + !!! Calculate mean branching statistics from values in 'branches' (not elements!) do N = 1,N_BR do i = 1,3 !for generations, Horsfield orders, Strahler orders ind(i) = nbranches(i,N) ! the branch gen, Hord, Sord -!!!...... length and diameter +!!!...... length and diameter do j = 1,2 sum_mean(i,j,ind(i)) = sum_mean(i,j,ind(i)) + branches(j,N) if(i.eq.3.and.j.eq.1)then - if(ind(i).ne.nbranches(5,N))then !not same as parent + if(ind(i).ne.nbranches(5,N))then !not same as parent ntally(i,j,ind(i)) = ntally(i,j,ind(i)) + 1 endif else ntally(i,j,ind(i)) = ntally(i,j,ind(i)) + 1 endif enddo !j -!!!...... branching angle and rotation angle +!!!...... branching angle and rotation angle do j = 3,4 if(branches(j,N).ge.0.0_dp)then sum_mean(i,j,ind(i)) = sum_mean(i,j,ind(i)) + branches(j,N) ntally(i,j,ind(i)) = ntally(i,j,ind(i)) + 1 endif enddo !j -!!!...... ratio of L:D +!!!...... ratio of L:D j = 5 if(branches(j,N).ge.0.0_dp)then sum_mean(i,j,ind(i)) = sum_mean(i,j,ind(i)) + branches(j,N) ntally(i,j,ind(i)) = ntally(i,j,ind(i)) + 1 endif enddo !i - + !!!... Summary statistics from branches do j = 3,5 !branching angle, rotation angle, L/D if(branches(j,N).ge.0.0_dp)then @@ -1511,7 +1512,7 @@ subroutine list_tree_statistics(filename) bins(N) = bins(N)/dble(nbins(N))*180.0_dp/PI endif enddo ! N - + !!!...... Summary statistics from branches do j = 3,5 !branching angle, rotation angle, L/D if(ntotaln(j-2).ne.0)then @@ -1520,7 +1521,7 @@ subroutine list_tree_statistics(filename) means(j-2) = 0.0_dp endif enddo !j - + i = 2 !Horsfield orders j = 6 !Nw/Nw-1 do N = 1,genm-1 @@ -1530,7 +1531,7 @@ subroutine list_tree_statistics(filename) sum_mean(i,j,N)=0.0_dp endif enddo !N - + !!! Summary statistics from CE do ne = 1,num_elems do j = 11,21 @@ -1540,14 +1541,14 @@ subroutine list_tree_statistics(filename) endif enddo !j enddo ! ne - + do j = 11,21 if(ntotaln(j-7).gt.0)then means(j-7) = means(j-7)/dble(ntotaln(j-7)) endif enddo !j !!! End of mean calculations - + !!! Calculate the standard deviations...... sum of (value-mean)^2 SD = 0.0_dp do N = 1,n_br @@ -1572,7 +1573,7 @@ subroutine list_tree_statistics(filename) endif enddo !j enddo !noelem - + !!! SD = sqrt(sum/(n-1)) do N = 1,genm do i = 1,3 !for generations, Horsfield orders, Strahler orders @@ -1588,8 +1589,8 @@ subroutine list_tree_statistics(filename) SDT(j) = sqrt(SDT(j)/dble(ntotaln(j)-1)) endif enddo !j -!!! End of standard deviation calculation - +!!! End of standard deviation calculation + !!! Output tree statistics average_term_gen = 0.0_dp sum_term = 0 @@ -1605,7 +1606,7 @@ subroutine list_tree_statistics(filename) & angle(deg)'')') write(10,'(115(''-''))') endif - + i = 1 do N = 1,nmax_gen(i) write(*,'(3(i10),5(f8.2,'' ('',f6.2,'')''))') N,ntally(i,1,N),n_terminal(N), & @@ -1624,7 +1625,7 @@ subroutine list_tree_statistics(filename) else average_term_gen = 0.0_dp endif - + write(*,'(/'' Horsfield #branches Length'',11x,''Diameter& & Branching Rotation ratio L:D Nw/Nw-1'')') write(*,'(4x,''order'',20x,''(mm)'',14x,''(mm)'',9x,''angle(deg)& @@ -1637,7 +1638,7 @@ subroutine list_tree_statistics(filename) & angle(deg)'')') write(10,'(115(''-''))') endif - + i = 2 do N = 1,nmax_gen(i) write(*,'(2(i10),5(f8.2,'' ('',f6.2,'')''),f8.2)') N,ntally(2,1,N),sum_mean(i,1,N), & @@ -1649,7 +1650,7 @@ subroutine list_tree_statistics(filename) SD(i,4,N),sum_mean(i,5,N),SD(i,5,N),sum_mean(i,6,N) endif enddo - + write(*,'(/'' Strahler #branches Length'',10x,''Diameter'',8x,''Branching& & Rotation'',10x,''ratio L:D'')') write(*,'('' order'',18x,''(mm)'',13x,''(mm)'',10x,''angle(deg) angle(deg)'')') @@ -1671,7 +1672,7 @@ subroutine list_tree_statistics(filename) sum_mean(i,5,N),SD(i,5,N) endif enddo - + do i = 2,3 !Horsfield and Strahler orders do N = 1,nmax_gen(i) X(N) = N @@ -1702,20 +1703,20 @@ subroutine list_tree_statistics(filename) write(*,'('' L/Lp = '',f7.3,'' ('',f6.3,'')'')') means(12),SDT(12) write(*,'('' %L/Lp < 1 = '',f7.3)') dble(num_llp)/dble(num_elems-1)*100.0_dp write(*,'('' L1/L2 (L1 < L2) = '',f7.3,'' ('',f6.3,'')'')') means(13),SDT(13) - + write(*,'('' Rb Strahler = '',f7.3,'' Rsq = '',f6.3)') ratios(3,1),r_sq(3,1) write(*,'('' Rl Strahler = '',f7.3,'' Rsq = '',f6.3)') ratios(3,2),r_sq(3,2) write(*,'('' Rd Strahler = '',f7.3,'' Rsq = '',f6.3)') ratios(3,3),r_sq(3,3) write(*,'('' Rb Horsfield = '',f7.3,'' Rsq = '',f6.3)') ratios(2,1),r_sq(2,1) write(*,'('' Rl Horsfield = '',f7.3,'' Rsq = '',f6.3)') ratios(2,2),r_sq(2,2) write(*,'('' Rd Horsfield = '',f7.3,'' Rsq = '',f6.3)') ratios(2,3),r_sq(2,3) - + write(*,'('' mean angle Dp 4.0+ = '',f7.3)') bins(1) write(*,'('' mean angle Dp 3.0+ = '',f7.3)') bins(2) write(*,'('' mean angle Dp 2.0+ = '',f7.3)') bins(3) write(*,'('' mean angle Dp 1.0+ = '',f7.3)') bins(4) write(*,'('' mean angle Dp 0.7+ = '',f7.3)') bins(5) - + if(writefile)then write(10,'(/''SUMMARY OF MEAN GEOMETRY STATISTICS'')') write(10,'(60(''-''))') @@ -1735,32 +1736,32 @@ subroutine list_tree_statistics(filename) write(10,'('' L/Lp = '',f7.3,'' ('',f6.3,'')'')') means(12),SDT(12) write(10,'('' %L/Lp < 1 = '',f7.3)') dble(num_llp)/dble(num_elems-1)*100.0_dp write(10,'('' L1/L2 (L1 < L2) = '',f7.3,'' ('',f6.3,'')'')') means(13),SDT(13) - + write(10,'('' Rb Strahler = '',f7.3,'' Rsq = '',f6.3)') ratios(3,1),r_sq(3,1) write(10,'('' Rl Strahler = '',f7.3,'' Rsq = '',f6.3)') ratios(3,2),r_sq(3,2) write(10,'('' Rd Strahler = '',f7.3,'' Rsq = '',f6.3)') ratios(3,3),r_sq(3,3) write(10,'('' Rb Horsfield = '',f7.3,'' Rsq = '',f6.3)') ratios(2,1),r_sq(2,1) write(10,'('' Rl Horsfield = '',f7.3,'' Rsq = '',f6.3)') ratios(2,2),r_sq(2,2) write(10,'('' Rd Horsfield = '',f7.3,'' Rsq = '',f6.3)') ratios(2,3),r_sq(2,3) - + write(10,'('' mean angle Dp 4.0+ = '',f7.3)') bins(1) write(10,'('' mean angle Dp 3.0+ = '',f7.3)') bins(2) write(10,'('' mean angle Dp 2.0+ = '',f7.3)') bins(3) write(10,'('' mean angle Dp 1.0+ = '',f7.3)') bins(4) write(10,'('' mean angle Dp 0.7+ = '',f7.3)') bins(5) endif - + deallocate(diameters) deallocate(stats) deallocate(branches) deallocate(nbranches) - + close(10) - + call enter_exit(sub_name,2) - + end subroutine list_tree_statistics - + !!!############################################################################# subroutine linregress(n,r_squared,slope,x,y) @@ -1782,10 +1783,10 @@ subroutine linregress(n,r_squared,slope,x,y) xxsum = xxsum + x(i) * x(i) enddo !!! calculate least squares estimate of straight line thru solution - slope = (xysum-xsum*ysum/n)/(xxsum-xsum*xsum/n) - intercept = ysum/n - slope*xsum/n + slope = (xysum-xsum*ysum/n)/(xxsum-xsum*xsum/n) + intercept = ysum/n - slope*xsum/n !!! calculate r-squared correlation coefficient -!!! see Numerical Recipes, Fortran 77, 2nd edition, page 632. +!!! see Numerical Recipes, Fortran 77, 2nd edition, page 632. ax = xsum/N !mean of X ay = ysum/N !mean of Y sxx = 0.0_dp @@ -1800,9 +1801,9 @@ subroutine linregress(n,r_squared,slope,x,y) enddo r = sxy/sqrt(sxx*syy) r_squared = r**2.0_dp - + end subroutine linregress - + !!!############################################################################# subroutine triangles_from_surface(surface_elems) @@ -3060,7 +3061,7 @@ subroutine merge_trifurcations(short_elements) end subroutine merge_trifurcations !!!############################################################################# - + subroutine define_rad_from_file(FIELDFILE, radius_type_in) !*define_rad_from_file:* reads in a radius field associated with an ! airway tree and assigns radius information to each element, also @@ -3351,6 +3352,38 @@ subroutine define_rad_from_geom(ORDER_SYSTEM, CONTROL_PARAM, START_FROM, & end subroutine define_rad_from_geom + !!!############################################################################# + + subroutine occlude_vessel(VESSEL_NUMBER, RATIO) + !*occlude_vessel:* Occludes/modifies vessel or airway radius based on + ! the ratio provided by user. This subroutine is made for partial occlusions + ! where the vessel (artery or vein) element number and the ratio is provided by user + ! and this will be applied on unstrained radius of the vessels. This subroutine should + ! be called after define_rad_from_geom/file so that the tree radii are identified. + ! This subroutine is useful for running Pulmonary hypertension or pulmonary embolism cases. + + integer, intent(in) :: VESSEL_NUMBER ! Element number that you want to apply occlusion on + real(dp), intent(in) :: RATIO ! partial/or full occlsion ratio (100 means full occlusions + ! and any other number between 0 and 100 is partial occlusion) + + character(len=60) :: sub_name + + ! -------------------------------------------------------------------------- + + sub_name = 'occlude_vessel' + call enter_exit(sub_name,1) + ! write(*,*) "before:", elem_field(ne_radius_in, VESSEL_NUMBER) + ! write(*,*) "ne_radius_in:", elem_field(ne_radius_in,21) + + elem_field(ne_radius, VESSEL_NUMBER) = RATIO * elem_field(ne_radius, VESSEL_NUMBER) + elem_field(ne_radius_in, VESSEL_NUMBER) = RATIO * elem_field(ne_radius_in, VESSEL_NUMBER) + elem_field(ne_radius_out, VESSEL_NUMBER) = RATIO * elem_field(ne_radius_out, VESSEL_NUMBER) + ! write(*,*) "after:", elem_field(ne_radius_in, VESSEL_NUMBER) + ! pause + call enter_exit(sub_name,2) + + end subroutine occlude_vessel + !!!############################################################################# subroutine element_connectivity_1d() @@ -3600,9 +3633,9 @@ subroutine internal_mesh_reorder real(dp),allocatable :: temp_elem_direction(:,:),temp_elem_field(:,:) logical,allocatable :: temp_expansile(:) character(len=60) :: sub_name - + ! -------------------------------------------------------------------------- - + sub_name = 'internal_mesh_reorder' call enter_exit(sub_name,1) @@ -3610,12 +3643,12 @@ subroutine internal_mesh_reorder allocate(list_term_branches(num_elems)) allocate(map_to_new(num_elems)) allocate(map_to_old(num_elems)) - + ! work through each successive generation, incrementing elements one by one count_elems = 0 num_term_branches = 1 ngen = 0 - list_term_branches(1) = 1 ! this assumes that the first element is the stem! + list_term_branches(1) = 1 ! this assumes that the first element is the stem! do while(num_term_branches.ne.0) num_branches = num_term_branches ! temporary, to loop over num_term_branches = 0 ! reset to zero and count for this generation @@ -3651,7 +3684,7 @@ subroutine internal_mesh_reorder allocate(temp_elem_field(num_ne,num_elems)) allocate(temp_elem_direction(3,num_elems)) if(model_type.eq.'gas_mix') allocate(temp_expansile(num_elems)) - + do ne = 1,num_elems ! for the ordered elements ne_old = map_to_old(ne) ! the unordered element number temp_elems(ne) = elems(ne_old) ! mapping to global @@ -3686,13 +3719,13 @@ subroutine internal_mesh_reorder call element_connectivity_1d call evaluate_ordering - + elem_ordrs(no_type,:) = 1 ! all conducting - + call enter_exit(sub_name,2) end subroutine internal_mesh_reorder - + !!!############################################################################# subroutine evaluate_ordering() diff --git a/src/lib/pressure_resistance_flow.f90 b/src/lib/pressure_resistance_flow.f90 index 583fc487..ad301d5e 100644 --- a/src/lib/pressure_resistance_flow.f90 +++ b/src/lib/pressure_resistance_flow.f90 @@ -17,6 +17,7 @@ module pressure_resistance_flow use solve, only: BICGSTAB_LinSolv,pmgmres_ilu_cr implicit none + integer,allocatable :: occ_list(:),compl_list(:) !Module parameters @@ -26,7 +27,7 @@ module pressure_resistance_flow !Interfaces private - public evaluate_prq,calculate_ppl + public evaluate_prq,calculate_ppl,occlusion_list,compliance_list contains !################################################################################### ! @@ -450,6 +451,65 @@ subroutine boundary_conditions(ADD,FIX,bc_type,grav_vect,density,inletbc,outletb endif call enter_exit(sub_name,2) end subroutine boundary_conditions + +! +!################################################################################### +! + !* occlusion_list gets an element list from user which includes the element that were + ! occluded for CTEPH and now the occlusions are removed and we want to make sure there + ! is no remodelling downstream of those elements + subroutine occlusion_list(elem_list) + + integer,intent(in) :: elem_list(:) ! list of surface elements defining the host region + !local variables + integer :: i + character(len=60) :: sub_name + + sub_name = 'occlusion_list' + call enter_exit(sub_name,1) + + if(count(elem_list.ne.0).gt.0)then ! a surface element list is given for converting to + ! create a list of occluded elements + allocate(occ_list(count(elem_list.ne.0))) + !!! get the list of occlusion list from the given list + do i = 1,count(elem_list.ne.0) + occ_list(i) = elem_list(i) + enddo + endif + ! + ! write(*,*) 'occlusion_list:', occ_list + ! pause + + + call enter_exit(sub_name,2) + end subroutine occlusion_list +! +!################################################################################### +! + !* compliance_list gets an element list from user which includes the element that were + ! occluded for CTEPH the compliance of these vessels is reduced to zero - basically make these + ! vessels rigid + subroutine compliance_list(elem_list2) + + integer,intent(in) :: elem_list2(:) ! list of surface elements defining the host region + !local variables + integer :: i + character(len=60) :: sub_name + + sub_name = 'compliance_list' + call enter_exit(sub_name,1) + + if(count(elem_list2.ne.0).gt.0)then ! a surface element list is given for converting to + ! create a list of occluded elements + allocate(compl_list(count(elem_list2.ne.0))) + !!! get the list of occlusion list from the given list + do i = 1,count(elem_list2.ne.0) + compl_list(i) = elem_list2(i) + enddo + endif + + call enter_exit(sub_name,2) + end subroutine compliance_list ! !################################################################################### ! @@ -931,10 +991,11 @@ subroutine calc_press_area(grav_vect,KOUNT,depvar_at_node,prq_solution,& real(dp),intent(in) :: elasticity_parameters(3),mechanics_parameters(2) !local variables - integer :: nj,np,ne,ny,nn + integer :: nj,np,ne,ny,nn,k,no_compl_ne real(dp) :: h,Ptm,R0,Pblood,Ppl,counter,cc1,cc2,cc3 real(dp) :: alt_hyp,alt_fib,prox_fib,narrow_rad_one,narrow_rad_two,narrow_factor,prune_rad,prune_fraction,counter1,counter2 integer,allocatable :: templss(:) + logical:: check, FOUND character(len=60) :: sub_name sub_name = 'calc_press_area' @@ -957,45 +1018,45 @@ subroutine calc_press_area(grav_vect,KOUNT,depvar_at_node,prq_solution,& Ptm=Pblood+Ppl ! Pa if(nn.eq.1)R0=elem_field(ne_radius_in0,ne) if(nn.eq.2)R0=elem_field(ne_radius_out0,ne) - if(vessel_type.eq.'elastic_g0_beta')then - if(Ptm.lt.elasticity_parameters(3).and.elasticity_parameters(1).gt.0.0_dp)then - if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((Ptm/elasticity_parameters(1))+1.d0)**(1.d0/elasticity_parameters(2)) - if(nn.eq.2) elem_field(ne_radius_out,ne)=R0*((Ptm/elasticity_parameters(1))+1.d0)**(1.d0/elasticity_parameters(2)) - elseif(Ptm.lt.0.0_dp.or.elasticity_parameters(1).LT.TOLERANCE)then - if(Ptm.lt.0)write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl - if(nn.eq.1) elem_field(ne_radius_in,ne)=R0 - if(nn.eq.2) elem_field(ne_radius_out,ne)=R0 - else !ptm>ptmmax - if(nn.eq.1)then - elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(3)/elasticity_parameters(1))+1.d0) & - **(1.d0/elasticity_parameters(2)) - endif - if(nn.eq.2)then - elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(3)/elasticity_parameters(1))+1.d0) & - **(1.d0/elasticity_parameters(2)) - endif - endif - elseif(vessel_type.eq.'elastic_alpha')then - if(Ptm.LT.elasticity_parameters(2))then - if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((Ptm*elasticity_parameters(1))+1.d0) - if(nn.eq.2) elem_field(ne_radius_out,ne)=R0*((Ptm*elasticity_parameters(1))+1.d0) - elseif(Ptm.lt.0.0_dp)then - if(Ptm.lt.0)write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl - if(nn.eq.1) elem_field(ne_radius_in,ne)=R0 - if(nn.eq.2) elem_field(ne_radius_out,ne)=R0 - else !ptm>ptmmax - if(nn.eq.1)then - elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(2)*elasticity_parameters(1))+1.d0) + if(vessel_type.eq.'elastic_g0_beta')then + if(Ptm.lt.elasticity_parameters(3).and.elasticity_parameters(1).gt.0.0_dp)then + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((Ptm/elasticity_parameters(1))+1.d0)**(1.d0/elasticity_parameters(2)) + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0*((Ptm/elasticity_parameters(1))+1.d0)**(1.d0/elasticity_parameters(2)) + elseif(Ptm.lt.0.0_dp.or.elasticity_parameters(1).LT.TOLERANCE)then + if(Ptm.lt.0)write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0 + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0 + else !ptm>ptmmax + if(nn.eq.1)then + elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(3)/elasticity_parameters(1))+1.d0) & + **(1.d0/elasticity_parameters(2)) + endif + if(nn.eq.2)then + elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(3)/elasticity_parameters(1))+1.d0) & + **(1.d0/elasticity_parameters(2)) + endif endif - if(nn.eq.2)then - elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(2)*elasticity_parameters(1))+1.d0) + elseif(vessel_type.eq.'elastic_alpha')then + if(Ptm.LT.elasticity_parameters(2))then + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((Ptm*elasticity_parameters(1))+1.d0) + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0*((Ptm*elasticity_parameters(1))+1.d0) + elseif(Ptm.lt.0.0_dp)then + if(Ptm.lt.0)write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0 + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0 + else !ptm>ptmmax + if(nn.eq.1)then + elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(2)*elasticity_parameters(1))+1.d0) + endif + if(nn.eq.2)then + elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(2)*elasticity_parameters(1))+1.d0) + endif endif - endif - elseif(vessel_type.eq.'elastic_hooke')then - h=elasticity_parameters(2)*R0 - if(nn.eq.1) elem_field(ne_radius_in,ne)=R0+3.0_dp*R0**2*Ptm/(4.0_dp*elasticity_parameters(1)*h) - if(nn.eq.2) elem_field(ne_radius_out,ne)=R0+3.0_dp*R0**2*Ptm/(4.0_dp*elasticity_parameters(1)*h) - else + elseif(vessel_type.eq.'elastic_hooke')then + h=elasticity_parameters(2)*R0 + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0+3.0_dp*R0**2*Ptm/(4.0_dp*elasticity_parameters(1)*h) + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0+3.0_dp*R0**2*Ptm/(4.0_dp*elasticity_parameters(1)*h) + else print *, 'no vessel type defined, assuming rigid' if(nn.eq.1) elem_field(ne_radius_in,ne)=R0 if(nn.eq.2) elem_field(ne_radius_out,ne)=R0 @@ -1056,211 +1117,272 @@ subroutine calc_press_area(grav_vect,KOUNT,depvar_at_node,prq_solution,& if(nn.eq.1) R0=elem_field(ne_radius_in0,ne) if(nn.eq.2) R0=elem_field(ne_radius_out0,ne) if(elem_field(ne_group,ne).eq.0.0_dp) then !only applying on arteries - if(nn.eq.1) then - if(R0.lt.prune_rad.and.elem_ordrs(no_sord,ne).eq.1) then - if(counter1/100.le.prune_fraction) then ! pruning the right percentage based on the fraction defined - cc1 = cc1+1.0_dp - R0=0.005_dp ! Setting the radius to a small value - else ! the remaining of the canditates that are not pruned because of the fraction - R0=elem_field(ne_radius_in0,ne) + FOUND = .False. + if (allocated(occ_list))then + do k=1,count(occ_list.ne.0) + check = is_downstream(ne,occ_list(k)) + if(check)then + FOUND = .True. endif - counter1 = counter1 + 1.0_dp ! since a canditate was found, one is added to the counter1 - if(counter1.ge.101.0_dp) counter1=1.0_dp ! now that the fraction out of hundred was blocked set the counter back to start - else ! R0 greater than prune_rad - cc3=cc3+1.0_dp - R0=elem_field(ne_radius_in0,ne) ! treating the artery as normal unstrained radius (no constraints) - endif + end do endif - if(nn.eq.2) then ! same thing as nn=1 - if(R0.lt.prune_rad.and.elem_ordrs(no_sord,ne).eq.1) then - if(counter2/100.le.prune_fraction) then - cc2=cc2+1.0_dp - R0=0.005_dp - else - R0=elem_field(ne_radius_out0,ne) - endif - counter2 = counter2 + 1.0_dp - if(counter2.ge.101.0_dp) counter2=1.0_dp - else - R0=elem_field(ne_radius_out0,ne) - endif - endif - endif - if(vessel_type.eq.'elastic_g0_beta') then - if(elem_field(ne_group,ne).eq.0.0_dp) then !only applying on arteries - if(Ptm.LT.elasticity_parameters(3).and.elasticity_parameters(1).gt.0.0_dp)then - if(nn.eq.1) then - if((R0.gt.0.015).and.(R0.lt.0.15)) then - elem_field(ne_radius_in,ne)=0.55_dp*R0*((Ptm/(0.16_dp*elasticity_parameters(1)))+1.d0) & - **(1.d0/elasticity_parameters(2)) - else - elem_field(ne_radius_in,ne)=R0*((Ptm/elasticity_parameters(1))+1.d0)**(1.d0/elasticity_parameters(2)) + if(.NOT.FOUND)then ! FOUND + if(nn.eq.1) then ! nn.eq.1 + if(R0.lt.prune_rad.and.elem_ordrs(no_sord,ne).eq.1) then + if(counter1/100.le.prune_fraction) then ! pruning the right percentage based on the fraction defined + cc1 = cc1+1.0_dp + R0=0.005_dp ! Setting the radius to a small value + else ! the remaining of the canditates that are not pruned because of the fraction + R0=elem_field(ne_radius_in0,ne) + endif + counter1 = counter1 + 1.0_dp ! since a canditate was found, one is added to the counter1 + if(counter1.ge.101.0_dp) counter1=1.0_dp ! now that the fraction out of hundred was blocked set the counter back to start + else ! R0 greater than prune_rad + cc3=cc3+1.0_dp + R0=elem_field(ne_radius_in0,ne) ! treating the artery as normal unstrained radius (no constraints) endif - endif - if(nn.eq.2) then - if((R0.gt.0.015).and.(R0.lt.0.15)) then - elem_field(ne_radius_out,ne)=0.55_dp*R0*((Ptm/(0.16_dp*elasticity_parameters(1)))+1.d0) & - **(1.d0/elasticity_parameters(2)) + endif ! nn.eq.1 + if(nn.eq.2) then ! same thing as nn=1 + if(R0.lt.prune_rad.and.elem_ordrs(no_sord,ne).eq.1) then + if(counter2/100.le.prune_fraction) then + cc2=cc2+1.0_dp + R0=0.005_dp + else + R0=elem_field(ne_radius_out0,ne) + endif + counter2 = counter2 + 1.0_dp + if(counter2.ge.101.0_dp) counter2=1.0_dp else - elem_field(ne_radius_out,ne)=R0*((Ptm/elasticity_parameters(1))+1.d0)**(1.d0/elasticity_parameters(2)) + R0=elem_field(ne_radius_out0,ne) endif - endif - elseif(Ptm.lt.0.0_dp.or.elasticity_parameters(1).LT.TOLERANCE)then - if(Ptm.lt.0)write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl - if(nn.eq.1) elem_field(ne_radius_in,ne)=R0 - if(nn.eq.2) elem_field(ne_radius_out,ne)=R0 - else !ptm>ptmmax - if(nn.eq.1) then - if((R0.gt.0.015).and.(R0.lt.0.15)) then - elem_field(ne_radius_in,ne)=0.55_dp*R0*((elasticity_parameters(3)/(0.16_dp*elasticity_parameters(1)))+1.d0) & + endif ! same thing as nn=1 + else ! treat as healthy + if(nn.eq.1) R0=elem_field(ne_radius_in0,ne) + if(nn.eq.2) R0=elem_field(ne_radius_out0,ne) + endif ! FOUND + if(vessel_type.eq.'elastic_g0_beta') then ! elastic_g0_beta + if(elem_field(ne_group,ne).eq.0.0_dp) then !only applying on arteries + if(Ptm.LT.elasticity_parameters(3).and.elasticity_parameters(1).gt.0.0_dp)then + if(nn.eq.1) then + if((R0.gt.0.015).and.(R0.lt.0.15)) then + elem_field(ne_radius_in,ne)=0.55_dp*R0*((Ptm/(0.16_dp*elasticity_parameters(1)))+1.d0) & + **(1.d0/elasticity_parameters(2)) + else + elem_field(ne_radius_in,ne)=R0*((Ptm/elasticity_parameters(1))+1.d0)**(1.d0/elasticity_parameters(2)) + endif + endif + if(nn.eq.2) then + if((R0.gt.0.015).and.(R0.lt.0.15)) then + elem_field(ne_radius_out,ne)=0.55_dp*R0*((Ptm/(0.16_dp*elasticity_parameters(1)))+1.d0) & + **(1.d0/elasticity_parameters(2)) + else + elem_field(ne_radius_out,ne)=R0*((Ptm/elasticity_parameters(1))+1.d0)**(1.d0/elasticity_parameters(2)) + endif + endif + elseif(Ptm.lt.0.0_dp.or.elasticity_parameters(1).LT.TOLERANCE)then + if(Ptm.lt.0)write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0 + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0 + else !ptm>ptmmax + if(nn.eq.1) then + if((R0.gt.0.015).and.(R0.lt.0.15)) then + elem_field(ne_radius_in,ne)=0.55_dp*R0*((elasticity_parameters(3)/(0.16_dp*elasticity_parameters(1)))+1.d0) & **(1.d0/elasticity_parameters(2)) - else - elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(3)/elasticity_parameters(1))+1.d0) & + else + elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(3)/elasticity_parameters(1))+1.d0) & **(1.d0/elasticity_parameters(2)) + endif + endif + if(nn.eq.2) then + if((R0.gt.0.015).and.(R0.lt.0.15)) then + elem_field(ne_radius_out,ne)=0.55_dp*R0*((elasticity_parameters(3)/(0.16_dp*elasticity_parameters(1)))+1.d0) & + **(1.d0/elasticity_parameters(2)) + else + elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(3)/elasticity_parameters(1))+1.d0) & + **(1.d0/elasticity_parameters(2)) + endif + endif endif - endif - if(nn.eq.2) then - if((R0.gt.0.015).and.(R0.lt.0.15)) then - elem_field(ne_radius_out,ne)=0.55_dp*R0*((elasticity_parameters(3)/(0.16_dp*elasticity_parameters(1)))+1.d0) & - **(1.d0/elasticity_parameters(2)) - else - elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(3)/elasticity_parameters(1))+1.d0) & + else !other than arteries + if(Ptm.LT.elasticity_parameters(3).and.elasticity_parameters(1).gt.0.0_dp)then + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((Ptm/elasticity_parameters(1))+1.d0)**(1.d0/elasticity_parameters(2)) + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0*((Ptm/elasticity_parameters(1))+1.d0)**(1.d0/elasticity_parameters(2)) + elseif(Ptm.lt.0.0_dp.or.elasticity_parameters(1).LT.TOLERANCE)then + if(Ptm.lt.0)write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0 + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0 + else !ptm>ptmmax + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(3)/elasticity_parameters(1))+1.d0) & + **(1.d0/elasticity_parameters(2)) + if(nn.eq.2)then + elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(3)/elasticity_parameters(1))+1.d0) & **(1.d0/elasticity_parameters(2)) + endif endif - endif - endif + endif !artery or vein + elseif(vessel_type.eq.'elastic_alpha') then ! elastic alpha + !if(elem_field(ne_group,ne).eq.0.0_dp) then !only applying on arteries BEN WAS HERE + FOUND = .False. + if (allocated(occ_list))then ! occ allocated + do k=1,count(occ_list.ne.0) + check = is_downstream(ne,occ_list(k)) + if(check)then + FOUND = .True. + endif + end do + endif ! occ allocated + if(.NOT.FOUND)then !FOUND + if(Ptm.lt.elasticity_parameters(2))then + if(nn.eq.1) then ! nn.eq.1 + if((R0.gt.narrow_rad_one).and.(R0.lt.0.5_dp)) then ! Hypertorphy+narrow factor effect + if(R0.lt.0.05_dp) then ! only Narrow_factor + elem_field(ne_radius_in,ne) = narrow_factor*R0*((Ptm*elasticity_parameters(1))+1.d0) + elseif(R0.gt.narrow_rad_two) then ! only Hypertophy + elem_field(ne_radius_in,ne) = R0*((Ptm*alt_hyp*elasticity_parameters(1))+1.d0) + else ! Both hypertophy and narrowing + elem_field(ne_radius_in,ne) = narrow_factor*R0*((Ptm*alt_hyp*alt_fib*elasticity_parameters(1))+1.d0) + endif + elseif(R0.gt.0.5_dp) then !Not within the range of our target radii, hence, no remodeling for this element + elem_field(ne_radius_in,ne) = R0*((Ptm*elasticity_parameters(1))+1.d0) + else !Pruning + elem_field(ne_radius_in,ne) = R0 + endif + endif ! nn.eq.1 + if(nn.eq.2) then ! nn.eq.2 + if((R0.gt.narrow_rad_one).and.(R0.lt.0.5_dp)) then ! Hypertorphy+narrow factor effect + if(R0.lt.0.05_dp) then ! only Narrow_factor + elem_field(ne_radius_out,ne) = narrow_factor*R0*((Ptm*elasticity_parameters(1))+1.d0) + elseif(R0.gt.narrow_rad_two) then ! hypertophy only + elem_field(ne_radius_out,ne) = R0*((Ptm*alt_hyp*elasticity_parameters(1))+1.d0) + else ! Both hypertophy and narrowing + elem_field(ne_radius_out,ne) = narrow_factor*R0*((Ptm*alt_hyp*alt_fib*elasticity_parameters(1))+1.d0) + endif + elseif(R0.gt.0.5_dp) then !Not within the range of our target radii, hence, no remodeling + elem_field(ne_radius_out,ne) = R0*((Ptm*elasticity_parameters(1))+1.d0) + else !Pruning + elem_field(ne_radius_out,ne) = R0 + endif + endif ! nn.eq.2 + if (allocated(compl_list)) then ! compl allocated + do no_compl_ne=1,count(compl_list.ne.0) + !pause + if(ne.eq.compl_list(no_compl_ne))then + if(nn.eq.1) elem_field(ne_radius_in,compl_list(no_compl_ne))=R0!*((Ptm*elasticity_parameters(1)*0.0)+1.d0) + if(nn.eq.2) elem_field(ne_radius_out,compl_list(no_compl_ne))=R0!*((Ptm*elasticity_parameters(1)*0.0)+1.d0) + endif + end do + endif ! compl allocated + elseif(Ptm.lt.0.0_dp)then !Ptm + if(Ptm.lt.0) write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl + if(nn.eq.1) then ! nn.eq.1 + if((R0.gt.narrow_rad_one).and.(R0.lt.0.5)) then ! Hypertophy+narrowing effect + if(R0.lt.narrow_rad_two) then !only narrowing + elem_field(ne_radius_in,ne)=narrow_factor*R0 + else + elem_field(ne_radius_in,ne)=R0 + endif + else ! Not within the target range, hence, no remodeling + elem_field(ne_radius_in,ne)=R0 + endif + endif ! nn.eq.1 + if(nn.eq.2) then ! nn.eq.2 + if((R0.gt.narrow_rad_one).and.(R0.lt.0.5_dp)) then + if(R0.lt.narrow_rad_two) then + elem_field(ne_radius_out,ne)=narrow_factor*R0 + else + elem_field(ne_radius_out,ne)=R0 + endif + else ! Not within the target range, hence, no remodeling + elem_field(ne_radius_out,ne)=R0 + endif + endif ! nn.eq.2 + else !ptm>ptmmax + if(nn.eq.1) then ! nn.eq.1 + if((R0.gt.narrow_rad_one).and.(R0.lt.0.5_dp)) then ! Hypertophy+narrowing effect + if(R0.lt.0.05_dp) then ! only Narrow_factor + elem_field(ne_radius_in,ne) = narrow_factor*R0*((Ptm*alt_hyp*alt_fib*elasticity_parameters(1))+1.d0) + elseif(R0.gt.narrow_rad_two) then ! hypertophy only + elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(2)*alt_hyp*elasticity_parameters(1))+1.d0) + else ! Both hypertophy and narrowing + elem_field(ne_radius_in,ne)=narrow_factor*R0*(elasticity_parameters(2)* & + alt_hyp*alt_fib*elasticity_parameters(1)+1.d0) + endif + elseif(R0.gt.0.5_dp) then ! Not within the target range, hence, no remodeling + elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(2)*(elasticity_parameters(1)))+1.d0) + else + elem_field(ne_radius_in,ne)=R0 + endif + endif ! nn.eq.1 + if(nn.eq.2) then ! nn.eq.2 + if((R0.gt.narrow_rad_one).and.(R0.lt.0.5_dp)) then ! Hypertrophy+narrowing effect + if(R0.lt.0.05_dp) then ! only Narrow_factor + elem_field(ne_radius_out,ne)=narrow_factor*R0*((elasticity_parameters(2)*elasticity_parameters(1))+1.d0) + elseif(R0.gt.narrow_rad_two) then ! hypertophy only + elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(2)*alt_hyp*elasticity_parameters(1))+1.d0) + else ! Both hypertophy and narrowing + elem_field(ne_radius_out,ne)=narrow_factor*R0*((elasticity_parameters(2)* & + alt_hyp*alt_fib*elasticity_parameters(1))+1.d0) + endif + elseif(R0.gt.0.5_dp) then ! Not within the target range, hence, no remodeling + elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(2)*(elasticity_parameters(1)))+1.d0) + else + elem_field(ne_radius_out,ne)=R0 + endif ! radius criteria + endif ! nn.eq.2 + if (allocated(compl_list)) then ! compl allocated + do no_compl_ne=1,count(compl_list.ne.0) + !pause + if(ne.eq.compl_list(no_compl_ne))then + if(nn.eq.1) elem_field(ne_radius_in,compl_list(no_compl_ne))=R0!*((Ptm*elasticity_parameters(1)*0.0)+1.d0) + if(nn.eq.2) elem_field(ne_radius_out,compl_list(no_compl_ne))=R0!*((Ptm*elasticity_parameters(1)*0.0)+1.d0) + endif + end do + endif ! compl allocated + endif !ptm + else ! FOUND + if(Ptm.lt.elasticity_parameters(2))then + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((Ptm*elasticity_parameters(1))+1.d0) + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0*((Ptm*elasticity_parameters(1))+1.d0) + elseif(Ptm.lt.0.0_dp)then !Ptm + if(Ptm.lt.0)write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0 + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0 + else !ptm>ptmmax + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(2)*elasticity_parameters(1))+1.d0) + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(2)*elasticity_parameters(1))+1.d0) + endif + endif ! FOUND + if (allocated(compl_list)) then ! compl allocated + do no_compl_ne=1,count(compl_list.ne.0) + ! pause + if(ne.eq.compl_list(no_compl_ne))then + if(nn.eq.1) elem_field(ne_radius_in,compl_list(no_compl_ne))=R0*((Ptm*elasticity_parameters(1)*0.0)+1.d0) + if(nn.eq.2) elem_field(ne_radius_out,compl_list(no_compl_ne))=R0*((Ptm*elasticity_parameters(1)*0.0)+1.d0) + endif + end do + endif ! compl allocated + elseif(vessel_type.eq.'elastic_hooke')then + h=elasticity_parameters(2)*R0 + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0+3.0_dp*R0**2*Ptm/(4.0_dp*elasticity_parameters(1)*h) + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0+3.0_dp*R0**2*Ptm/(4.0_dp*elasticity_parameters(1)*h) !vessel type + else + print *, 'no vessel type defined, assuming rigid' + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0 + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0 !vessel type + endif !vessel type else !other than arteries - if(Ptm.LT.elasticity_parameters(3).and.elasticity_parameters(1).gt.0.0_dp)then - if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((Ptm/elasticity_parameters(1))+1.d0)**(1.d0/elasticity_parameters(2)) - if(nn.eq.2) elem_field(ne_radius_out,ne)=R0*((Ptm/elasticity_parameters(1))+1.d0)**(1.d0/elasticity_parameters(2)) - elseif(Ptm.lt.0.0_dp.or.elasticity_parameters(1).LT.TOLERANCE)then + if(Ptm.lt.elasticity_parameters(2))then + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((Ptm*elasticity_parameters(1))+1.d0) + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0*((Ptm*elasticity_parameters(1))+1.d0) + elseif(Ptm.lt.0.0_dp)then !Ptm if(Ptm.lt.0)write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl if(nn.eq.1) elem_field(ne_radius_in,ne)=R0 if(nn.eq.2) elem_field(ne_radius_out,ne)=R0 else !ptm>ptmmax - if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(3)/elasticity_parameters(1))+1.d0) & - **(1.d0/elasticity_parameters(2)) - if(nn.eq.2)then - elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(3)/elasticity_parameters(1))+1.d0) & - **(1.d0/elasticity_parameters(2)) - endif + if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(2)*elasticity_parameters(1))+1.d0) + if(nn.eq.2) elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(2)*elasticity_parameters(1))+1.d0) endif - endif - elseif(vessel_type.eq.'elastic_alpha') then - if(elem_field(ne_group,ne).eq.0.0_dp) then !only applying on arteries - if(Ptm.lt.elasticity_parameters(2))then - if(nn.eq.1) then - if((R0.gt.narrow_rad_one).and.(R0.lt.0.5_dp)) then ! Hypertorphy+narrow factor effect - if(R0.lt.0.05_dp) then ! only Narrow_factor - elem_field(ne_radius_in,ne) = narrow_factor*R0*((Ptm*elasticity_parameters(1))+1.d0) - elseif(R0.gt.narrow_rad_two) then ! only Hypertophy - elem_field(ne_radius_in,ne) = R0*((Ptm*alt_hyp*elasticity_parameters(1))+1.d0) - else ! Both hypertophy and narrowing - elem_field(ne_radius_in,ne) = narrow_factor*R0*((Ptm*alt_hyp*alt_fib*elasticity_parameters(1))+1.d0) - endif - elseif(R0.gt.0.5_dp) then !Not within the range of our target radii, hence, no remodeling for this element - elem_field(ne_radius_in,ne) = R0*((Ptm*elasticity_parameters(1))+1.d0) - else !Pruning - elem_field(ne_radius_in,ne) = R0 - endif - endif - if(nn.eq.2) then - if((R0.gt.narrow_rad_one).and.(R0.lt.0.5_dp)) then ! Hypertorphy+narrow factor effect - if(R0.lt.0.05_dp) then ! only Narrow_factor - elem_field(ne_radius_out,ne) = narrow_factor*R0*((Ptm*elasticity_parameters(1))+1.d0) - elseif(R0.gt.narrow_rad_two) then ! hypertophy only - elem_field(ne_radius_out,ne) = R0*((Ptm*alt_hyp*elasticity_parameters(1))+1.d0) - else ! Both hypertophy and narrowing - elem_field(ne_radius_out,ne) = narrow_factor*R0*((Ptm*alt_hyp*alt_fib*elasticity_parameters(1))+1.d0) - endif - elseif(R0.gt.0.5_dp) then !Not within the range of our target radii, hence, no remodeling - elem_field(ne_radius_out,ne) = R0*((Ptm*elasticity_parameters(1))+1.d0) - else !Pruning - elem_field(ne_radius_out,ne) = R0 - endif - endif - elseif(Ptm.lt.0.0_dp)then !Ptm - if(Ptm.lt.0) write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl - if(nn.eq.1) then - if((R0.gt.narrow_rad_one).and.(R0.lt.0.5)) then ! Hypertophy+narrowing effect - if(R0.lt.narrow_rad_two) then !only narrowing - elem_field(ne_radius_in,ne)=narrow_factor*R0 - else - elem_field(ne_radius_in,ne)=R0 - endif - else ! Not within the target range, hence, no remodeling - elem_field(ne_radius_in,ne)=R0 - endif - endif - if(nn.eq.2) then - if((R0.gt.narrow_rad_one).and.(R0.lt.0.5_dp)) then - if(R0.lt.narrow_rad_two) then - elem_field(ne_radius_out,ne)=narrow_factor*R0 - else - elem_field(ne_radius_out,ne)=R0 - endif - else ! Not within the target range, hence, no remodeling - elem_field(ne_radius_out,ne)=R0 - endif - endif - else !ptm>ptmmax - if(nn.eq.1) then - if((R0.gt.narrow_rad_one).and.(R0.lt.0.5_dp)) then ! Hypertophy+narrowing effect - if(R0.lt.0.05_dp) then ! only Narrow_factor - elem_field(ne_radius_in,ne) = narrow_factor*R0*((Ptm*alt_hyp*alt_fib*elasticity_parameters(1))+1.d0) - elseif(R0.gt.narrow_rad_two) then ! hypertophy only - elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(2)*alt_hyp*elasticity_parameters(1))+1.d0) - else ! Both hypertophy and narrowing - elem_field(ne_radius_in,ne)=narrow_factor*R0*(elasticity_parameters(2)* & - alt_hyp*alt_fib*elasticity_parameters(1)+1.d0) - endif - elseif(R0.gt.0.5_dp) then ! Not within the target range, hence, no remodeling - elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(2)*(elasticity_parameters(1)))+1.d0) - else - elem_field(ne_radius_in,ne)=R0 - endif - endif - if(nn.eq.2) then - if((R0.gt.narrow_rad_one).and.(R0.lt.0.5_dp)) then ! Hypertrophy+narrowing effect - if(R0.lt.0.05_dp) then ! only Narrow_factor - elem_field(ne_radius_out,ne)=narrow_factor*R0*((elasticity_parameters(2)*elasticity_parameters(1))+1.d0) - elseif(R0.gt.narrow_rad_two) then ! hypertophy only - elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(2)*alt_hyp*elasticity_parameters(1))+1.d0) - else ! Both hypertophy and narrowing - - elem_field(ne_radius_out,ne)=narrow_factor*R0*((elasticity_parameters(2)* & - alt_hyp*alt_fib*elasticity_parameters(1))+1.d0) - - endif - elseif(R0.gt.0.5_dp) then ! Not within the target range, hence, no remodeling - elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(2)*(elasticity_parameters(1)))+1.d0) - else - elem_field(ne_radius_out,ne)=R0 - endif - endif - endif - else !other than arteries - if(Ptm.lt.elasticity_parameters(2))then - if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((Ptm*elasticity_parameters(1))+1.d0) - if(nn.eq.2) elem_field(ne_radius_out,ne)=R0*((Ptm*elasticity_parameters(1))+1.d0) - elseif(Ptm.lt.0.0_dp)then !Ptm - if(Ptm.lt.0)write(*,*) 'Transmural pressure < zero',ne,Ptm,Pblood,Ppl - if(nn.eq.1) elem_field(ne_radius_in,ne)=R0 - if(nn.eq.2) elem_field(ne_radius_out,ne)=R0 - else !ptm>ptmmax - if(nn.eq.1) elem_field(ne_radius_in,ne)=R0*((elasticity_parameters(2)*elasticity_parameters(1))+1.d0) - if(nn.eq.2) elem_field(ne_radius_out,ne)=R0*((elasticity_parameters(2)*elasticity_parameters(1))+1.d0) - endif - endif - elseif(vessel_type.eq.'elastic_hooke')then - h=elasticity_parameters(2)*R0 - if(nn.eq.1) elem_field(ne_radius_in,ne)=R0+3.0_dp*R0**2*Ptm/(4.0_dp*elasticity_parameters(1)*h) - if(nn.eq.2) elem_field(ne_radius_out,ne)=R0+3.0_dp*R0**2*Ptm/(4.0_dp*elasticity_parameters(1)*h) - else - print *, 'no vessel type defined, assuming rigid' - if(nn.eq.1) elem_field(ne_radius_in,ne)=R0 - if(nn.eq.2) elem_field(ne_radius_out,ne)=R0 - endif + endif ! if artery enddo!nn enddo!ne endif @@ -1295,6 +1417,36 @@ subroutine map_solution_to_mesh(prq_solution,depvar_at_elem,depvar_at_node,mesh_ call enter_exit(sub_name,2) end subroutine map_solution_to_mesh +!##############################################################################!############################################################################## +! +!*is_downstream* checks to see if a certain element is downstream of another element +!!!!!!! +! This function is useful for doing regional remodeling where you do not want to +! remodel vessels downstream of an occlusion (needed for post-PEA modelling) +!!!!!!! +recursive function is_downstream(ne,occ_ne) result(FOUND) + + integer, intent(in) :: ne,occ_ne + !local variables + logical :: FOUND + integer :: ne_temp,np + character(len=60) :: sub_name + + FOUND = .False. + np=elem_nodes(1,ne) + + if (ne == occ_ne) then + FOUND = .True. + else if (ne == 1) then + FOUND = .false. + else + ne_temp = elems_at_node(np,1) + FOUND = is_downstream(ne_temp, occ_ne) + end if + + call enter_exit(sub_name,2) +end function is_downstream + !############################################################################## ! !*map_flow_to_terminals* maps the solution array to appropriate nodal and element fields diff --git a/src/lib/wave_transmission.f90 b/src/lib/wave_transmission.f90 index 131d9cdd..40332ae4 100644 --- a/src/lib/wave_transmission.f90 +++ b/src/lib/wave_transmission.f90 @@ -2852,9 +2852,12 @@ subroutine pressure_flow_factor(no_freq,p_factor,q_factor,reflect,prop_const,cha !(1+reflect(nf,ne)*exp(-2.0_dp*prop_const(nf,ne)*elem_field(ne_length,ne))) else ne_up=elem_cnct(-1,1,ne) - q_factor(nf,ne)=q_factor(nf,ne_up)*char_admit(nf,ne)*(1+reflect(nf,ne_up))* & - exp(-1.0_dp*elem_field(ne_length,ne_up)*prop_const(nf,ne_up))/& - (char_admit(nf,ne_up)*(1+reflect(nf,ne)*exp(-2.0_dp*elem_field(ne_length,ne)*prop_const(nf,ne)))) + p_factor(nf,ne)=p_factor(nf,ne_up)*(char_admit(nf,ne_up)*(1-reflect(nf,ne_up))* & + exp(-1.0_dp*elem_field(ne_length,ne_up)*prop_const(nf,ne_up)))/& + (2*char_admit(nf,ne)*(1-reflect(nf,ne)*exp(-2.0_dp*elem_field(ne_length,ne)*prop_const(nf,ne)))) + q_factor(nf,ne)=p_factor(nf,ne_up)*char_admit(nf,ne)!(1+reflect(nf,ne_up))* & + !exp(-1.0_dp*elem_field(ne_length,ne_up)*prop_const(nf,ne_up))/& + !(1+reflect(nf,ne)*exp(-2.0_dp*elem_field(ne_length,ne)*prop_const(nf,ne))) endif!neup enddo!ne enddo!nf