FindAdjacencyF90.F90

Go to the documentation of this file.

! FindAdjacency: Interacting with iMesh
! 
! This program shows how to get more information about a mesh, by
! getting connectivity two different ways (as connectivity and as
! adjacent 0-dimensional entities).
 
! Usage: FindAdjacency
#define CHECK(a) \
 if (ierr .ne. 0) print *, a
 
program findadjacency
 implicit none
#include "iMesh_f.h"
 
 ! declarations
 imesh_instance mesh
 ibase_handle_t rpents, rpverts, rpallverts, ipoffsets
 integer ioffsets
 ibase_handle_t ents, verts, allverts, verths
 pointer(rpents, ents(0:*))
 pointer(rpverts, verts(0:*))
 pointer(rpallverts, allverts(0:*))
 pointer(ipoffsets, ioffsets(0:*))
! for all vertices in one call
 ibase_entityhandle verth
 pointer(verths, verth(0:*))
 integer ierr, ents_alloc, ents_size
 integer iverts_alloc, iverts_size
 integer allverts_alloc, allverts_size
 integer offsets_alloc, offsets_size, coords_alloc, coords_size
 ibase_entitysethandle root_set
 integer vert_uses, i, num_ents
 real*8 coords
 pointer(pcoord, coords(0:*))
!
 ibase_entitysethandle :: sethand
 
 ! create the Mesh instance
 call imesh_newmesh("", mesh, ierr)
 check("Problems instantiating interface.")
 
 ! load the mesh
 call imesh_getrootset(%VAL(mesh), root_set, ierr)
 check("Problems getting root set")
 
 call imesh_load(%VAL(mesh), %VAL(root_set), &
 "../../MeshFiles/unittest/125hex.g", "", ierr)
 check("Load failed")
 
 ! get all 3d elements
 ents_alloc = 0
 call imesh_getentities(%VAL(mesh), %VAL(root_set), %VAL(ibase_region), &
 %VAL(imesh_all_topologies), rpents, ents_alloc, ents_size, &
 ierr)
 check("Couldn't get entities")
 
 vert_uses = 0
 
 ! iterate through them; 
 do i = 0, ents_size-1
 ! get connectivity
 iverts_alloc = 0
 call imesh_getentadj(%VAL(mesh), %VAL(ents(i)), &
 %VAL(ibase_vertex), &
 rpverts, iverts_alloc, iverts_size, &
 ierr)
 check("Failure in getEntAdj")
 ! sum number of vertex uses
 
 vert_uses = vert_uses + iverts_size
 
 if (iverts_size .ne. 0) call imesh_freememory(%VAL(mesh), rpverts)
 end do
 
 ! now get adjacencies in one big block
 allverts_alloc = 0
 offsets_alloc = 0
 call imesh_getentarradj(%VAL(mesh), %VAL(rpents), &
 %VAL(ents_size), %VAL(ibase_vertex), rpallverts, &
 allverts_alloc, allverts_size, ipoffsets, offsets_alloc, &
 offsets_size, ierr)
 check("Failure in getEntArrAdj")
 
 if (allverts_size .ne. 0) call imesh_freememory(%VAL(mesh), rpallverts);
 if (offsets_size .ne. 0) call imesh_freememory(%VAL(mesh), ipoffsets);
 if (ents_size .ne. 0) call imesh_freememory(%VAL(mesh), rpents);
 
 ! compare results of two calling methods
 if (allverts_size .ne. vert_uses) then
 write(*,'("Sizes didn''t agree!")')
 else 
 write(*, *)"Sizes did agree: ", vert_uses
 endif
! get all vertices , and then their coordinates
 ents_alloc = 0
 call imesh_getentities(%VAL(mesh), %VAL(root_set), %VAL(ibase_vertex), &
 %VAL(imesh_all_topologies), verths, ents_alloc, ents_size, &
 ierr)
 write (*, *) "number of vertices: " , ents_size
 print *, "few vertex handles: ", (verth(i), i=0,ents_size/10)
 
! set creation
 call imesh_createentset(%VAL(mesh), %VAL(1), &
 sethand,ierr)
 write(0,*) "createset",ierr,sethand 
 
! we should have 
 call imesh_addentarrtoset(%VAL(mesh),verth,%VAL(ents_size), &
 %VAL(sethand),ierr)
 write(0,*) "add Ent Arr to Set",ierr,sethand 
 
 call imesh_getnumoftype(%VAL(mesh), %VAL(sethand), &
 %VAL(ibase_vertex), num_ents, ierr)
 write(0,*) "num verts retrieved from set", num_ents
 
 ents_alloc = 0;
 call imesh_getvtxarrcoords(%VAL(mesh), verth, %VAL(ents_size), &
 %VAL(ibase_interleaved), pcoord, ents_alloc , ents_size, ierr)
!
 write(*, *) "num coords: ", ents_size, " few coords: ", (coords(i), i=0, ents_size/100) 
 
 call imesh_freememory(%VAL(mesh), verths);
 call imesh_freememory(%VAL(mesh), pcoord);
 
 call imesh_dtor(%VAL(mesh), ierr)
 check("Failed to destroy interface")
 
end program findadjacency