FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
34contains
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
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
308end 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
subroutine find_surface_neighbor(surf, bktdb)
Find neighboring surface elements.
Definition: surf_ele.f90:80
integer(kind=kint), parameter l_max_elem_surf
Definition: surf_ele.f90:14
integer(kind=kint), parameter n_neighbor_max_init
Definition: surf_ele.f90:16
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
subroutine update_surface_bucket_info(surf, bktdb)
Update bucket info for searching surface elements.
Definition: surf_ele.f90:272
integer(kind=kint) function next_position(surf, cpos)
Tracing next contact position.
Definition: surf_ele.f90:149
Structure for bucket search.
Structure to define surface group.
Definition: surf_ele.f90:19