1! FindAdjacency: Interacting with iMesh
2!
3! This program shows how to get more information about a mesh, by
4! getting connectivity two different ways (as connectivity and as
5! adjacent 0-dimensional entities).
6
7! Usage: FindAdjacency
8#define CHECK(a) \
9  if (ierr .ne. 0) print *, a
10
11program findadjacency
12  implicit none
13#include "iMesh_f.h"
14
15  ! declarations
16  iMesh_Instance mesh
17  IBASE_HANDLE_T rpents, rpverts, rpallverts, ipoffsets
18  integer ioffsets
19  IBASE_HANDLE_T ents, verts, allverts, verths
20  pointer (rpents, ents(0:*))
21  pointer (rpverts, verts(0:*))
22  pointer (rpallverts, allverts(0:*))
23  pointer (ipoffsets, ioffsets(0:*))
24!  for all vertices in one call
25  iBase_EntityHandle verth
26  pointer (verths, verth(0:*))
27  integer ierr, ents_alloc, ents_size
28  integer iverts_alloc, iverts_size
29  integer allverts_alloc, allverts_size
30  integer offsets_alloc, offsets_size, coords_alloc, coords_size
31  iBase_EntitySetHandle root_set
32  integer vert_uses, i, num_ents
33  real*8 coords
34  pointer (pcoord, coords(0:*))
35!
36  iBase_EntitySetHandle :: sethand
37
38  ! create the Mesh instance
39  call iMesh_newMesh("", mesh, ierr)
40  CHECK("Problems instantiating interface.")
41
42  ! load the mesh
43  call iMesh_getRootSet(%VAL(mesh), root_set, ierr)
44  CHECK("Problems getting root set")
45
46  call iMesh_load(%VAL(mesh), %VAL(root_set), &
47       "../../MeshFiles/unittest/125hex.g", "", ierr)
48  CHECK("Load failed")
49
50  ! get all 3d elements
51  ents_alloc = 0
52  call iMesh_getEntities(%VAL(mesh), %VAL(root_set), %VAL(iBase_REGION), &
53       %VAL(iMesh_ALL_TOPOLOGIES), rpents, ents_alloc, ents_size, &
54       ierr)
55  CHECK("Couldn't get entities")
56
57  vert_uses = 0
58
59  ! iterate through them;
60  do i = 0, ents_size-1
61     ! get connectivity
62     iverts_alloc = 0
63     call iMesh_getEntAdj(%VAL(mesh), %VAL(ents(i)), &
64          %VAL(iBase_VERTEX), &
65          rpverts, iverts_alloc, iverts_size, &
66          ierr)
67     CHECK("Failure in getEntAdj")
68     ! sum number of vertex uses
69
70     vert_uses = vert_uses + iverts_size
71
72     if (iverts_size .ne. 0) call iMesh_freeMemory(%VAL(mesh), rpverts)
73  end do
74
75  ! now get adjacencies in one big block
76  allverts_alloc = 0
77  offsets_alloc = 0
78  call iMesh_getEntArrAdj(%VAL(mesh), %VAL(rpents), &
79       %VAL(ents_size), %VAL(iBase_VERTEX), rpallverts,  &
80       allverts_alloc, allverts_size, ipoffsets, offsets_alloc, &
81       offsets_size, ierr)
82  CHECK("Failure in getEntArrAdj")
83
84  if (allverts_size .ne. 0) call iMesh_freeMemory(%VAL(mesh), rpallverts);
85  if (offsets_size .ne. 0) call iMesh_freeMemory(%VAL(mesh), ipoffsets);
86  if (ents_size .ne. 0) call iMesh_freeMemory(%VAL(mesh), rpents);
87
88  ! compare results of two calling methods
89  if (allverts_size .ne. vert_uses) then
90     write(*,'("Sizes didn''t agree!")')
91  else
92     write(*, *)"Sizes did agree: ", vert_uses
93  endif
94! get all vertices , and then their coordinates
95  ents_alloc = 0
96  call iMesh_getEntities(%VAL(mesh), %VAL(root_set), %VAL(iBase_VERTEX), &
97       %VAL(iMesh_ALL_TOPOLOGIES), verths, ents_alloc, ents_size, &
98       ierr)
99  write (*, *) "number of vertices: " , ents_size
100  print *, "few vertex handles: ", (verth(i), i=0,ents_size/10)
101
102! set creation
103  call iMesh_createEntSet(%VAL(mesh), %VAL(1), &
104        sethand,ierr)
105  write(0,*) "createset",ierr,sethand
106
107! we should have
108  call iMesh_addEntArrToSet(%VAL(mesh),verth,%VAL(ents_size), &
109      %VAL(sethand),ierr)
110  write(0,*) "add Ent Arr to Set",ierr,sethand
111
112  call iMesh_getNumOfType(%VAL(mesh), %VAL(sethand),  &
113   %VAL(iBase_VERTEX), num_ents, ierr)
114  write(0,*) "num verts retrieved from set", num_ents
115
116  ents_alloc = 0;
117  call iMesh_getVtxArrCoords(%VAL(mesh), verth, %VAL(ents_size), &
118    %VAL(iBase_INTERLEAVED), pcoord, ents_alloc , ents_size, ierr)
119!
120  write(*, *) "num coords: ", ents_size, " few coords: ", (coords(i), i=0, ents_size/100)
121
122  call iMesh_freeMemory(%VAL(mesh), verths);
123  call iMesh_freeMemory(%VAL(mesh), pcoord);
124
125  call iMesh_dtor(%VAL(mesh), ierr)
126  CHECK("Failed to destroy interface")
127
128end program findadjacency
129