FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
surf_ele.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
9  use hecmw, only: kint, kreal
10  implicit none
11 
12  integer(kind=kint), parameter :: l_max_surface_node =20
13  integer(kind=kint), parameter :: l_max_elem_node = 100
14  integer(kind=kint), parameter :: l_max_elem_surf = 6
15 
16  integer(kind=kint), parameter :: n_neighbor_max_init = 8
17 
20  integer(kind=kint) :: eid
21  integer(kind=kint) :: etype
22  integer(kind=kint), pointer :: nodes(:)=>null()
23  integer(kind=kint) :: n_neighbor
24  integer(kind=kint) :: n_neighbor_max
25  integer(kind=kint), pointer :: neighbor(:)=>null()
26  real(kind=kreal) :: reflen
27  real(kind=kreal) :: xavg(3)
28  real(kind=kreal) :: dmax
29  integer(kind=kint) :: bktid
30  end type tsurfelement
31 
32  integer(kind=kint), parameter, private :: debug = 0
33 
34 contains
35 
37  subroutine initialize_surf( eid, etype, nsurf, surf )
38  use elementinfo
39  integer(kind=kint), intent(in) :: eid
40  integer(kind=kint), intent(in) :: etype
41  integer(kind=kint), intent(in) :: nsurf
42  type(tsurfelement), intent(inout) :: surf
43  integer(kind=kint) :: n, outtype, nodes(100)
44  surf%eid = eid
45  call getsubface( etype, nsurf, outtype, nodes )
46  surf%etype = outtype
47  n=getnumberofnodes( outtype )
48  allocate( surf%nodes(n) )
49  surf%nodes(1:n)=nodes(1:n)
50  n=getnumberofsubface( outtype )
51  surf%n_neighbor = 0
52  surf%n_neighbor_max = 0
53  surf%reflen = -1.d0
54  surf%xavg(:) = 0.d0
55  surf%dmax = -1.d0
56  surf%bktID = -1
57  end subroutine
58 
60  subroutine finalize_surf( surf )
61  type(tsurfelement), intent(inout) :: surf
62  if( associated(surf%nodes) ) deallocate( surf%nodes )
63  if( associated(surf%neighbor) ) deallocate( surf%neighbor )
64  end subroutine
65 
67  subroutine write_surf( file, surf )
68  integer(kind=kint), intent(in) :: file
69  type(tsurfelement), intent(in) :: surf
70  integer(kind=kint) :: i
71  write(file,*) "Element:",surf%eid,"Surface type=",surf%etype
72  if( associated( surf%nodes ) ) &
73  write(file,*) ( surf%nodes(i),i=1,size(surf%nodes) )
74  if( associated( surf%neighbor ) ) &
75  write(file,*) ( surf%neighbor(i),i=1,surf%n_neighbor )
76  end subroutine
77 
79  subroutine find_surface_neighbor( surf, bktDB )
80  use m_utilities
81  use hecmw_util
82  use bucket_search
83  type(tsurfelement), intent(inout) :: surf(:)
84  type(bucketdb), intent(in) :: bktDB
85  integer(kind=kint) :: i, j, ii,jj, nd1, nd2, nsurf
86  integer(kind=kint) :: k, oldsize, newsize
87  integer(kind=kint), pointer :: dumarray(:) => null()
88  integer(kind=kint) :: bktID, ncand, js
89  integer(kind=kint), allocatable :: indexSurf(:)
90  if (debug >= 1) write(0,*) 'DEBUG: find_surface_neighbor: start'
91 
92  nsurf = size(surf)
93 
94  !$omp parallel do default(none), &
95  !$omp&private(i,ii,nd1,j,jj,nd2,oldsize,newsize,dumarray,k,bktID,ncand,indexSurf,js), &
96  !$omp&shared(nsurf,surf,bktDB)
97  do i=1,nsurf
98  bktid = bucketdb_getbucketid(bktdb, surf(i)%xavg)
99  ncand = bucketdb_getnumcand(bktdb, bktid)
100  if (ncand == 0) cycle
101  allocate(indexsurf(ncand))
102  call bucketdb_getcand(bktdb, bktid, ncand, indexsurf)
103  jloop: do js=1,ncand
104  j = indexsurf(js)
105  if( i==j ) cycle
106  if( associated(surf(i)%neighbor) ) then
107  if ( any( surf(i)%neighbor(1:surf(i)%n_neighbor)==j ) ) cycle
108  endif
109  do ii=1, size(surf(i)%nodes)
110  nd1 = surf(i)%nodes(ii)
111  do jj=1, size(surf(j)%nodes)
112  nd2 = surf(j)%nodes(jj)
113  if( nd1==nd2 ) then
114  surf(i)%n_neighbor = surf(i)%n_neighbor+1
115 
116  if( .not. associated(surf(i)%neighbor) ) then
117  allocate( surf(i)%neighbor(n_neighbor_max_init) )
118  surf(i)%n_neighbor_max = n_neighbor_max_init
119  surf(i)%neighbor(1) = j
120  else if( surf(i)%n_neighbor > surf(i)%n_neighbor_max ) then
121  oldsize = surf(i)%n_neighbor_max
122  newsize = oldsize * 2
123  dumarray => surf(i)%neighbor
124  allocate( surf(i)%neighbor(newsize) )
125  surf(i)%n_neighbor_max = newsize
126  do k=1,oldsize
127  surf(i)%neighbor(k) = dumarray(k)
128  enddo
129  surf(i)%neighbor(oldsize+1) = j
130  deallocate( dumarray )
131  else
132  surf(i)%neighbor(surf(i)%n_neighbor) = j
133  endif
134 
135  cycle jloop
136  endif
137  enddo
138  enddo
139  enddo jloop
140  deallocate(indexsurf)
141  enddo
142  !$omp end parallel do
143 
144  if (debug >= 1) write(0,*) 'DEBUG: find_surface_neighbor: end'
145  end subroutine
146 
148  integer(kind=kint) function next_position( surf, cpos )
149  use elementinfo
150  type(tsurfelement), intent(in) :: surf
151  real(kind=kreal), intent(in) :: cpos(2)
152 
153  integer(kind=kint) :: i
154  real(kind=kreal) :: maxv(3)
155  next_position = surf%eid
156  if( .not. associated(surf%neighbor) ) return ! do nothing when not initialized
157  maxv(:) = 0.d0
158  select case(surf%etype)
159  case( fe_tri3n, fe_tri6n )
160  if( all(cpos>0.d0) .and. all(cpos<1.d0) ) return
161  if( size(surf%neighbor)<3 ) return
162  if( all(cpos(:)>0.d0) ) return
163  do i=1,3
164  if( cpos(i)< 0.d0 ) maxv(i) = -cpos(i)
165  enddo
166  next_position = maxloc(maxv,1)
167  case( fe_quad4n, fe_quad8n )
168  if( all(cpos>-1.d0) .and. all(cpos<1.d0) ) return
169  if( size(surf%neighbor)<4 ) return
170  if( cpos(1)<-1.d0 .and. dabs(cpos(2))<1.d0 ) then
171  next_position = 1
172  elseif( cpos(1)>1.d0 .and. dabs(cpos(2))<1.d0 ) then
173  next_position = 3
174  elseif( dabs(cpos(1))<1.d0 .and. cpos(2)<-1.d0 ) then
175  next_position = 2
176  elseif( dabs(cpos(1))<1.d0 .and. cpos(2)>1.d0 ) then
177  next_position = 4
178  elseif( cpos(1)<-1.d0 .and. cpos(2)<-1.d0 ) then
179  if( cpos(1)>cpos(2) ) then
180  next_position = 2
181  else
182  next_position = 1
183  endif
184  elseif( cpos(1)<-1.d0 .and. cpos(2)>1.d0 ) then
185  if( dabs(cpos(1))>cpos(2) ) then
186  next_position = 1
187  else
188  next_position = 4
189  endif
190  elseif( cpos(1)>1.d0 .and. cpos(2)<-1.d0 ) then
191  if( cpos(1)>dabs(cpos(2)) ) then
192  next_position = 3
193  else
194  next_position = 2
195  endif
196  elseif( cpos(1)>1.d0 .and. cpos(2)>1.d0 ) then
197  if( cpos(1)>cpos(2) ) then
198  next_position = 3
199  else
200  next_position = 4
201  endif
202  endif
203  case default
204  stop "type of surface element not defined"
205  end select
206  next_position = surf%neighbor( next_position )
207 
208  end function next_position
209 
211  subroutine update_surface_reflen( surf, coord )
212  use elementinfo
213  type(tsurfelement), intent(inout) :: surf(:)
214  real(kind=kreal), intent(in) :: coord(:)
215  real(kind=kreal) :: elem(3, l_max_surface_node), r0(2)
216  integer(kind=kint) :: nn, i, j, iss
217  !$omp parallel do default(none) private(j,nn,i,iss,elem,r0) shared(surf,coord)
218  do j=1,size(surf)
219  nn = size(surf(j)%nodes)
220  do i=1,nn
221  iss = surf(j)%nodes(i)
222  elem(1:3,i) = coord(3*iss-2:3*iss)
223  enddo
224  call getelementcenter( surf(j)%etype, r0 )
225  surf(j)%reflen = getreferencelength( surf(j)%etype, nn, r0, elem )
226  enddo
227  !$omp end parallel do
228  end subroutine update_surface_reflen
229 
231  subroutine update_surface_box_info( surf, currpos )
232  type(tsurfelement), intent(inout) :: surf(:)
233  real(kind=kreal), intent(in) :: currpos(:)
234  integer(kind=kint) :: nn, i, j, iss
235  real(kind=kreal) :: elem(3,l_max_surface_node),xmin(3), xmax(3), xsum(3), lx(3)
236  !$omp parallel do default(none) private(i,nn,iss,elem,xmin,xmax,xsum,lx) shared(surf,currpos)
237  do j=1, size(surf)
238  nn = size(surf(j)%nodes)
239  do i=1,nn
240  iss = surf(j)%nodes(i)
241  elem(1:3,i) = currpos(3*iss-2:3*iss)
242  enddo
243  do i=1,3
244  xmin(i) = minval(elem(i,1:nn))
245  xmax(i) = maxval(elem(i,1:nn))
246  xsum(i) = sum(elem(i,1:nn))
247  enddo
248  surf(j)%xavg(:) = xsum(:) / nn
249  lx(:) = xmax(:) - xmin(:)
250  surf(j)%dmax = maxval(lx) * 0.5d0
251  enddo
252  !$omp end parallel do
253  end subroutine update_surface_box_info
254 
256  logical function is_in_surface_box(surf, x0, exp_rate)
257  type(tsurfelement), intent(inout) :: surf
258  real(kind=kreal), intent(in) :: x0(3)
259  real(kind=kreal), intent(in) :: exp_rate
260  real(kind=kreal) :: dist(3), er
261  er = max(exp_rate, 1.d0)
262  dist(:) = abs(x0(:) - surf%xavg(:))
263  if ( maxval(dist(:)) < surf%dmax * er ) then
264  is_in_surface_box = .true.
265  else
266  is_in_surface_box = .false.
267  endif
268  end function is_in_surface_box
269 
271  subroutine update_surface_bucket_info(surf, bktDB)
272  use bucket_search
273  type(tsurfelement), intent(inout) :: surf(:)
274  type(bucketdb), intent(inout) :: bktDB
275  real(kind=kreal) :: x_min(3), x_max(3), d_max
276  integer(kind=kint) :: nsurf, i, j
277  if (debug >= 1) write(0,*) 'DEBUG: update_surface_bucket_info: start'
278  nsurf = size(surf)
279  if (nsurf == 0) return
280  x_min(:) = surf(1)%xavg(:)
281  x_max(:) = surf(1)%xavg(:)
282  d_max = surf(1)%dmax
283  do i = 2, nsurf
284  do j = 1, 3
285  if (surf(i)%xavg(j) < x_min(j)) x_min(j) = surf(i)%xavg(j)
286  if (surf(i)%xavg(j) > x_max(j)) x_max(j) = surf(i)%xavg(j)
287  if (surf(i)%dmax > d_max) d_max = surf(i)%dmax
288  enddo
289  enddo
290  do j = 1, 3
291  x_min(j) = x_min(j) - d_max
292  x_max(j) = x_max(j) + d_max
293  enddo
294  call bucketdb_setup(bktdb, x_min, x_max, d_max*2.d0, nsurf)
295  !$omp parallel do default(none) private(i) shared(nsurf,surf,bktDB)
296  do i = 1, nsurf
297  surf(i)%bktID = bucketdb_getbucketid(bktdb, surf(i)%xavg)
298  call bucketdb_registerpre(bktdb, surf(i)%bktID)
299  enddo
300  !$omp end parallel do
301  call bucketdb_allocate(bktdb)
302  do i = 1, nsurf
303  call bucketdb_register(bktdb, surf(i)%bktID, i)
304  enddo
305  if (debug >= 1) write(0,*) 'DEBUG: update_surface_bucket_info: end'
306  end subroutine update_surface_bucket_info
307 
308 end module msurfelement
This module provides bucket-search functionality It provides definition of bucket info and its access...
subroutine, public bucketdb_allocate(bktdb)
Allocate memory before actually registering members Before allocating memory, bucketDB_registerPre ha...
subroutine, public bucketdb_registerpre(bktdb, bid)
Pre-register for just counting members to be actually registered Bucket ID has to be obtained with bu...
integer(kind=kint) function, public bucketdb_getbucketid(bktdb, x)
Get bucket ID that includes given point.
subroutine, public bucketdb_register(bktdb, bid, sid)
Register member Before actually register, bucketDB_allocate has to be called.
integer(kind=kint) function, public bucketdb_getnumcand(bktdb, bid)
Get number of candidates within neighboring buckets of a given bucket Bucket ID has to be obtained wi...
subroutine, public bucketdb_getcand(bktdb, bid, ncand, cand)
Get candidates within neighboring buckets of a given bucket Number of candidates has to be obtained w...
subroutine, public bucketdb_setup(bktdb, x_min, x_max, dmin, n_tot)
Setup basic info of buckets.
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
Definition: element.f90:126
subroutine getelementcenter(fetype, localcoord)
Return natural coordinate of the center of surface element.
Definition: element.f90:1005
integer, parameter fe_tri6n
Definition: element.f90:70
integer(kind=kind(2)) function getnumberofsubface(etype)
Obtain number of sub-surface.
Definition: element.f90:163
subroutine getsubface(intype, innumber, outtype, nodes)
Find the definition of surface of the element.
Definition: element.f90:188
integer, parameter fe_tri3n
Definition: element.f90:69
integer, parameter fe_quad4n
Definition: element.f90:72
real(kind=kreal) function getreferencelength(fetype, nn, localcoord, elecoord)
This function calculates reference length at a point in surface.
Definition: element.f90:1209
integer, parameter fe_quad8n
Definition: element.f90:73
I/O and Utility.
Definition: hecmw_util_f.F90:7
Definition: hecmw.f90:6
This module provides aux functions.
Definition: utilities.f90:6
This module manage surface elements in 3D It provides basic definition of surface elements (trianglar...
Definition: surf_ele.f90:8
subroutine initialize_surf(eid, etype, nsurf, surf)
Initializer.
Definition: surf_ele.f90:38
integer(kind=kint), parameter l_max_elem_node
Definition: surf_ele.f90:13
subroutine update_surface_reflen(surf, coord)
Compute reference length of surface elements.
Definition: surf_ele.f90:212
subroutine write_surf(file, surf)
Write out elemental surface.
Definition: surf_ele.f90:68
subroutine update_surface_box_info(surf, currpos)
Update info of cubic box including surface elements.
Definition: surf_ele.f90:232
integer(kind=kint), parameter l_max_elem_surf
Definition: surf_ele.f90:14
subroutine find_surface_neighbor(surf, bktDB)
Find neighboring surface elements.
Definition: surf_ele.f90:80
integer(kind=kint), parameter n_neighbor_max_init
Definition: surf_ele.f90:16
subroutine update_surface_bucket_info(surf, bktDB)
Update bucket info for searching surface elements.
Definition: surf_ele.f90:272
integer(kind=kint), parameter l_max_surface_node
Definition: surf_ele.f90:12
subroutine finalize_surf(surf)
Memeory management subroutine.
Definition: surf_ele.f90:61
logical function is_in_surface_box(surf, x0, exp_rate)
Check if given point is within cubic box including surface element.
Definition: surf_ele.f90:257
integer(kind=kint) function next_position(surf, cpos)
Tracing next contact position.
Definition: surf_ele.f90:149
Structure to define surface group.
Definition: surf_ele.f90:19