Merge pull request #120 from tawhai/generalised_tree_growing

Generalised tree growing

CMakeFiles/cmake.check_cache

@@ -0,0 +1 @@
+# This file is generated by cmake for dependency checking of the CMakeCache.txt file

src/bindings/c/geometry.c

@@ -12,8 +12,7 @@ void define_mesh_geometry_test_c(void);
 void define_node_geometry_c(const char *NODEFILE, int *filename_len);
 void define_node_geometry_2d_c(const char *NODEFILE, int *filename_len);
 void define_data_geometry_c(const char *DATAFILE, int *filename_len);
-void group_elem_parent_term_c(int *ne_parent);
-void make_data_grid_c(int *surface_elems, double *spacing, int *to_export, const char *filename, int *filename_len, const char *groupname, int *groupname_len);
+extern void make_data_grid_c(int *elemlist_len, int elemlist[], double *offset, double *spacing, const char *filename, int *filename_len, const char *groupname, int *groupname_len);
 extern void make_2d_vessel_from_1d_c(int *elemlist_len, int elemlist[]);
 void define_rad_from_file_c(const char *FIELDFILE, int *filename_len, const char *radius_type, int *radius_type_len);
 int get_local_node_f_c(const char *ndimension, int *dimension_len, const char *np_global, int *np_global_len);
@@ -80,16 +79,11 @@ void define_data_geometry(const char *DATAFILE)
   define_data_geometry_c(DATAFILE, &filename_len);
 }
 
-void group_elem_parent_term(int ne_parent)
-{
-  group_elem_parent_term_c(&ne_parent);
-}
-
-void make_data_grid(int surface_elems, double spacing, int to_export, const char *filename, const char *groupname)
+void make_data_grid(int elemlist_len, int elemlist[], double offset, double spacing, const char *filename, const char *groupname)
 {
   int filename_len = (int)strlen(filename);
   int groupname_len = (int)strlen(groupname);
-  make_data_grid_c(&surface_elems, &spacing, &to_export, filename, &filename_len, groupname, &groupname_len);
+  make_data_grid_c(&elemlist_len, elemlist, &offset, &spacing, filename, &filename_len, groupname, &groupname_len);
 }
 
 void make_2d_vessel_from_1d(int elemlist_len, int elemlist[])

src/bindings/c/geometry.f90

@@ -130,7 +130,8 @@ end subroutine define_node_geometry_c
 !
 !###################################################################################
 !
-  subroutine make_data_grid_c(surface_elems, spacing, to_export, filename, filename_len, groupname, groupname_len)&
+  subroutine make_data_grid_c(surface_elems_len, surface_elems, offset, spacing, &
+       filename, filename_len, groupname, groupname_len)&
  bind(C, name="make_data_grid_c")
 
     use arrays,only: dp
@@ -140,9 +141,9 @@ subroutine make_data_grid_c(surface_elems, spacing, to_export, filename, filenam
     use geometry, only: make_data_grid
     implicit none
 
-    integer,intent(in) :: surface_elems(:)
-    real(dp),intent(in) :: spacing
-    logical,intent(in) :: to_export
+    integer,intent(in) :: surface_elems_len
+    integer,intent(in) :: surface_elems(surface_elems_len)
+    real(dp),intent(in) :: offset, spacing
     integer,intent(in) :: filename_len, groupname_len
     type(c_ptr), value, intent(in) :: filename, groupname
     character(len=MAX_FILENAME_LEN) :: filename_f
@@ -151,7 +152,11 @@ subroutine make_data_grid_c(surface_elems, spacing, to_export, filename, filenam
     call strncpy(filename_f, filename, filename_len)
     call strncpy(groupname_f, groupname, groupname_len)
 
-    call make_data_grid(surface_elems, spacing, to_export, filename_f, groupname_f)
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_make_data_grid(surface_elems, offset, spacing, filename_f, groupname_f)
+#else
+    call make_data_grid(surface_elems, offset, spacing, filename_f, groupname_f)
+#endif
 
   end subroutine make_data_grid_c
 
@@ -168,21 +173,6 @@ subroutine make_2d_vessel_from_1d_c(elemlist, elemlist_len) bind(C, name="make_2
 
   end subroutine make_2d_vessel_from_1d_c
 
-!
-!###################################################################################
-!
-  subroutine group_elem_parent_term_c(ne_parent) bind(C, name="group_elem_parent_term_c")
-
-    use iso_c_binding, only: c_ptr
-    use geometry, only: group_elem_parent_term
-    implicit none
-
-    integer,intent(in) :: ne_parent
-
-    call group_elem_parent_term(ne_parent)
-
-  end subroutine group_elem_parent_term_c
-
 !
 !###################################################################################
 !

src/bindings/c/geometry.h

@@ -13,8 +13,7 @@ SHO_PUBLIC void define_mesh_geometry_test();
 SHO_PUBLIC void define_node_geometry(const char *NODEFILE);
 SHO_PUBLIC void define_node_geometry_2d(const char *NODEFILE);
 SHO_PUBLIC void define_data_geometry(const char *DATAFILE);
-SHO_PUBLIC void group_elem_parent_term(int ne_parent);
-SHO_PUBLIC void make_data_grid(int surface_elems, double spacing, int to_export, const char *filename, const char *groupname);
+SHO_PUBLIC void make_data_grid(int elemlist_len, int elemlist[], double offset, double spacing, const char *filename, const char *groupname);
 SHO_PUBLIC void make_2d_vessel_from_1d(int elemlist_len, int elemlist[]);
 SHO_PUBLIC void define_rad_from_file(const char *FIELDFILE, const char *radius_type);
 SHO_PUBLIC int get_local_node_f(const char *ndimenstion, const char *np_global);

src/bindings/c/growtree.c

@@ -3,10 +3,12 @@
 
 #include "string.h"
 
-void grow_tree_c(int *parent_ne, int *surface_elems, double *angle_max, double *angle_min, double *branch_fraction, double *length_limit, double *shortest_length, double *rotation_limit, int *to_export, const char *filename, int *filename_len);
+extern void grow_tree_c(int *elemlist_len, int elemlist[], int *parent_ne, double *angle_max, double *angle_min, double *branch_fraction, double *length_limit, double *shortest_length, double *rotation_limit);
 
-void grow_tree(int parent_ne, int surface_elems, double angle_max, double angle_min, double branch_fraction, double length_limit, double shortest_length, double rotation_limit,int to_export, const char *filename)
+
+void grow_tree(int elemlist_len, int elemlist[], int parent_ne, double angle_max, double angle_min, double branch_fraction, double length_limit, double shortest_length, double rotation_limit)
 {
-  int filename_len = strlen(filename);
-  grow_tree_c(&parent_ne, &surface_elems, &angle_max, &angle_min, &branch_fraction, &length_limit, &shortest_length, &rotation_limit, &to_export, filename, &filename_len);
+  grow_tree_c(&elemlist_len, elemlist, &parent_ne, &angle_max, &angle_min, &branch_fraction, &length_limit, &shortest_length, &rotation_limit);
 }
+
+

src/bindings/c/growtree.f90

@@ -3,37 +3,39 @@ module growtree_c
   private
 
 contains
-! 
-!#########################################################################
-! 
-!*growtree:* the main growing subroutine. Generates a volume-filling tree into a closed surface.
-  subroutine grow_tree_c(parent_ne, surface_elems, angle_max, angle_min, branch_fraction, length_limit,&
-shortest_length, rotation_limit, to_export, filename, filename_len) bind(C, name="grow_tree_c")
 
+  !
+  !###################################################################################
+  !
+  ! the main growing subroutine. Generates a volume-filling tree into a closed surface.
+  subroutine grow_tree_c(surface_elems_len, surface_elems, parent_ne, &
+       angle_max, angle_min, branch_fraction, length_limit, &
+       shortest_length, rotation_limit) bind(C, name="grow_tree_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 growtree,only: grow_tree
     implicit none
-
+
+    integer,intent(in) :: surface_elems_len
+    integer,intent(in) :: surface_elems(surface_elems_len)
     integer,intent(in) :: parent_ne
-    integer,intent(in) :: surface_elems(:)
     real(dp),intent(in) :: angle_max
     real(dp),intent(in) :: angle_min
     real(dp),intent(in) :: branch_fraction
     real(dp),intent(in) :: length_limit
     real(dp),intent(in) :: shortest_length
     real(dp),intent(in) :: rotation_limit
-    logical,intent(in) :: to_export
-    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 grow_tree(parent_ne, surface_elems, angle_max, angle_min, branch_fraction, length_limit,&
-      shortest_length, rotation_limit, to_export, filename_f)
+#if defined _WIN32 && defined __INTEL_COMPILER
+    call so_grow_tree(surface_elems, parent_ne, angle_max, angle_min, branch_fraction, length_limit,&
+shortest_length, rotation_limit)
+#else
+    call grow_tree(surface_elems, parent_ne, angle_max, angle_min, branch_fraction, length_limit,&
+shortest_length, rotation_limit)
+#endif
 
   end subroutine grow_tree_c
 

src/bindings/c/growtree.h

@@ -4,6 +4,7 @@
 
 #include "symbol_export.h"
 
-SHO_PUBLIC void grow_tree(int parent_ne, int surface_elems, double angle_max, double angle_min, double branch_fraction, double length_limit, double shortest_length, double rotation_limit, int to_export, const char *filename);
+
+SHO_PUBLIC void grow_tree(int elemlist_len, int elemlist[], int parent_ne, double angle_max, double angle_min, double branch_fraction, double length_limit, double shortest_length, double rotation_limit);
 
 #endif /* AETHER_GROWTREE_H */

src/bindings/interface/growtree.i

@@ -1,9 +1,31 @@
-
 %module(package="aether") growtree
-%include symbol_export.h
-%include growtree.h
+  %include symbol_export.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(freearg) (int elemlist_len, int elemlist[]) {
+  if ($2) free($2);
+ }
 
 %{
 #include "growtree.h"
-%}
+  %}
 
+%include growtree.h

src/lib/exports.f90

@@ -29,11 +29,155 @@ module exports
        export_terminal_solution, &
        export_terminal_perfusion,&
        export_terminal_ssgexch, &
+       export_triangle_elements, &
+       export_triangle_nodes, &
        export_1d_elem_field, &
        export_data_geometry
 
 contains
-!!!################################################################
+
+!
+!##############################################################################
+!
+
+  subroutine export_triangle_elements(num_triangles,triangle,EXELEMFILE,groupname)
+    !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EXPORT_TRIANGLE_ELEMENTS" :: EXPORT_TRIANGLE_ELEMENTS
+
+!!! Parameters
+    integer :: num_triangles, triangle(:,:)
+    character(len=MAX_FILENAME_LEN),intent(in) :: EXELEMFILE
+    character(len=MAX_STRING_LEN),intent(in) :: groupname
+
+!!! Local Variables
+    integer :: ne,nj,nline,nn
+    character(len=1) :: char1
+    character(len=100) :: writefile
+    character(len=60) :: sub_name
+
+    ! --------------------------------------------------------------------------
+
+    sub_name = 'export_triangle_elements'
+    call enter_exit(sub_name,1)
+
+    if(index(EXELEMFILE, ".exelem")> 0) then !full filename is given
+       writefile = EXELEMFILE(1:100)
+    else ! need to append the correct filename extension
+       writefile = trim(EXELEMFILE)//'.exelem'
+    endif
+    open(10, file=writefile, status='replace')
+    !**     write the group name
+    write(10,'( '' Group name: '',a10)') groupname
+    !**     write the lines
+    write(10,'( '' Shape. Dimension=1'' )')
+    nline = 0
+    do ne = 1,num_triangles
+       write(10,'( '' Element: 0 0 '',I5)') nline+1
+       write(10,'( '' Element: 0 0 '',I5)') nline+2
+       write(10,'( '' Element: 0 0 '',I5)') nline+3
+       nline = nline+3
+    enddo !ne
+
+    !**        write the elements
+    write(10,'( '' Shape. Dimension=2, line*line'' )')
+    write(10,'( '' #Scale factor sets=1'' )')
+    write(10,'( '' l.Lagrange*l.Lagrange, #Scale factors=4'' )')
+    write(10,'( '' #Nodes= '',I2 )') 4
+    write(10,'( '' #Fields=1'' )')
+    write(10,'( '' 1) coordinates, coordinate, rectangular cartesian, #Components=3'')')
+
+    do nj = 1,3
+       if(nj==1) char1='x'; if(nj==2) char1='y'; if(nj==3) char1='z';
+       write(10,'(''  '',A2,''. l.Lagrange*l.Lagrange, no modify, standard node based.'')') char1
+       write(10,'( ''   #Nodes=4'')')
+       do nn = 1,4
+          write(10,'(''   '',I1,''. #Values=1'')') nn
+          write(10,'(''     Value indices: '',I4)') 1
+          write(10,'(''     Scale factor indices: '',I4)') nn
+       enddo !nn
+    enddo !nj
+
+    nline = 0
+    do ne = 1,num_triangles
+       !**         write the element
+       write(10,'(1X,''Element: '',I12,'' 0 0'' )') ne
+       !**          write the faces
+       WRITE(10,'(3X,''Faces: '' )')
+       WRITE(10,'(5X,''0 0'',I6)') nline+1
+       WRITE(10,'(5X,''0 0'',I6)') nline+2
+       WRITE(10,'(5X,''0 0'',I6)') nline+3
+       WRITE(10,'(5X,''0 0'',I6)') 0
+       nline = nline+3
+       !**          write the nodes
+       write(10,'(3X,''Nodes:'' )')
+       write(10,'(4X,4(1X,I12))') triangle(:,ne),triangle(3,ne)
+       !**          write the scale factors
+       write(10,'(3X,''Scale factors:'' )')
+       write(10,'(4X,4(1X,E12.5))') 1.0_dp,1.0_dp,1.0_dp,1.0_dp
+    enddo
+    close(10)
+
+    call enter_exit(sub_name,2)
+
+  end subroutine export_triangle_elements
+
+!
+!##############################################################################
+!
+
+  subroutine export_triangle_nodes(num_vertices, vertex_xyz, EXNODEFILE, groupname)
+  !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EXPORT_TRIANGLE_NODES" :: EXPORT_TRIANGLE_NODES
+
+!!! Parameters
+    integer :: num_vertices
+    real(dp) :: vertex_xyz(:,:)
+    character(len=MAX_FILENAME_LEN),intent(in) :: EXNODEFILE
+    character(len=MAX_STRING_LEN),intent(in) :: groupname
+
+!!! Local Variables
+    integer :: i,nj
+    character(len=100) :: writefile
+    character(len=60) :: sub_name
+
+    ! --------------------------------------------------------------------------
+
+    sub_name = 'export_triangle_nodes'
+    call enter_exit(sub_name,1)
+
+    if(index(EXNODEFILE, ".exnode")> 0) then !full filename is given
+       writefile = EXNODEFILE(1:100)
+    else ! need to append the correct filename extension
+       writefile = trim(EXNODEFILE)//'.exnode'
+    endif
+    open(10, file = writefile, status='replace')
+    !**    write the group name
+    write(10,'( '' Group name: '',A)') trim(groupname)
+    !*** Exporting Geometry
+    !*** Write the field information
+    write(10,'( '' #Fields=1'' )')
+    write(10,'('' 1) coordinates, coordinate, rectangular cartesian, '',&
+         &''#Components=3'')')
+    do nj = 1,3
+       if(nj.eq.1) write(10,'(2X,''x.  '')',advance="no")
+       if(nj.eq.2) write(10,'(2X,''y.  '')',advance="no")
+       if(nj.eq.3) write(10,'(2X,''z.  '')',advance="no")
+       write(10,'(''Value index='',I2,'', #Derivatives='',I1)', &
+            advance="no") nj,0
+       write(10,'()')
+    enddo
+    do i = 1,num_vertices
+       !***    write the node
+       write(10,'(1X,''Node: '',I12)') i
+       write(10,'(2x,3(f12.6))') vertex_xyz(:,i)
+    enddo
+    close(10)
+
+    call enter_exit(sub_name,2)
+
+  end subroutine export_triangle_nodes
+
+!
+!##############################################################################
+!
 
   subroutine export_1d_elem_field(ne_field, EXELEMFILE, group_name, field_name )
 
@@ -259,7 +403,7 @@ subroutine export_elem_geometry_2d(EXELEMFILE, name, offset_elem, offset_node)
           enddo !nj
        endif
        !**               write the element
-       WRITE(10,'(1X,''Element: '',I12,'' 0 0'' )') ne+offset_elem
+       WRITE(10,'(1X,''Element: '',I12,'' 0 0'' )') elems_2d(ne)+offset_elem
        !**                 write the faces
        WRITE(10,'(3X,''Faces: '' )')