FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
fstr_precheck.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!-------------------------------------------------------------------------------
7contains
8
9 subroutine fstr_get_thickness(hecMESH,mid,thick)
10 use hecmw
11 use m_fstr
12 implicit none
13 type (hecmwST_local_mesh) :: hecMESH
14 integer(kind=kint) :: mid, ihead
15 real(kind=kreal) :: thick
16
17 ihead = hecmesh%section%sect_R_index(mid-1)
18 thick = hecmesh%section%sect_R_item(ihead+1)
19 ! if(thick.LE.0.0) then
20 ! write(*,*) "Zero thickness <= 0 is illegal"
21 ! call hecmw_abort( hecmw_comm_get_comm())
22 ! endif
23 end subroutine fstr_get_thickness
24 !C
25 !C***
26 !C*** Pre Check for FSTR solver
27 !C***
28 !C
29 subroutine fstr_precheck ( hecMESH, hecMAT, soltype )
30 use m_fstr
31 implicit none
32 type (hecmwST_local_mesh) :: hecMESH
33 type (hecmwST_matrix ) :: hecMAT
34 integer(kind=kint) :: soltype
35
36 if(myrank .EQ. 0) then
37 write(imsg,*)
38 write(imsg,*) ' **** STAGE PreCheck **'
39 endif
40
41 call fstr_precheck_elem ( hecmesh, hecmat )
42 write(idbg,*) 'fstr_precheck_elem: OK'
43
44 !C
45 !C Output sparsity pattern
46 !C
47 if( soltype == kstnzprof )then
48 call hecmw_nonzero_profile(hecmesh, hecmat)
49 endif
50
51 end subroutine fstr_precheck
52 !C
53 !C
54 subroutine fstr_precheck_elem ( hecMESH, hecMAT )
55 use m_fstr
59
60 implicit none
61
62 type (hecmwST_matrix) :: hecMAT
63 type (hecmwST_local_mesh) :: hecMESH
64 type (fstr_solid) :: fstrSOLID
65
66 !** Local variables
67 integer(kind=kint) :: nelem, mid, j, isect, nline, tline, icel, iiS
68 integer(kind=kint) :: ndof2, nelem_wo_mpc
69 integer(kind=kint) :: ie, ia, jelem, ic_type, nn, jS, jE, itype
70 integer(kind=kint) :: nodLOCAL(20),NTOTsum(1)
71 integer(kind=kint) :: nonzero
72 real(kind=kreal) :: ntdof2
73 real(kind=kreal) :: al, almin, almax, aa, thick, vol, avvol
74 real(kind=kreal) :: tvol, tvmax, tvmin, tlmax, tlmin, asp, aspmax
75 real(kind=kreal) :: xx(20),yy(20),zz(20)
76 real(kind=kreal) :: totsum(1),totmax(3),totmin(2)
77 !C
78 !C INIT
79 !C
80 nelem = 0
81 tvol = 0.0
82 tvmax = 0.0
83 tvmin = 1.0e+20
84 tlmax = 0.0
85 tlmin = 1.0e+20
86 aspmax = 0.0
87
88 !C
89 !C Mesh Summary
90 !C
91 write(ilog,"(a)") '### Mesh Summary ###'
92 write(ilog,"(a,i12)") ' Num of node:',hecmesh%n_node
93 write(* ,"(a,i12)") ' Num of node:',hecmesh%n_node
94 write(ilog,"(a,i12)") ' Num of DOF :',hecmesh%n_node*hecmesh%n_dof
95 write(* ,"(a,i12)") ' Num of DOF :',hecmesh%n_node*hecmesh%n_dof
96 ndof2 = hecmesh%n_dof**2
97 ntdof2 = dble(hecmesh%n_node*hecmesh%n_dof)**2
98 write(ilog,"(a,i12)") ' Num of elem:',hecmesh%n_elem
99 write(* ,"(a,i12)") ' Num of elem:',hecmesh%n_elem
100 nelem_wo_mpc = 0
101 do itype = 1, hecmesh%n_elem_type
102 js = hecmesh%elem_type_index(itype-1)
103 je = hecmesh%elem_type_index(itype )
104 ic_type = hecmesh%elem_type_item(itype)
105 if(.not. (hecmw_is_etype_link(ic_type) .or. hecmw_is_etype_patch(ic_type))) &
106 nelem_wo_mpc = nelem_wo_mpc + je-js
107 enddo
108 write(ilog,"(a,i12)") ' ** w/o MPC/Patch:',nelem_wo_mpc
109 write(* ,"(a,i12)") ' ** w/o MPC/Patch:',nelem_wo_mpc
110 do itype = 1, hecmesh%n_elem_type
111 js = hecmesh%elem_type_index(itype-1)
112 je = hecmesh%elem_type_index(itype )
113 ic_type = hecmesh%elem_type_item(itype)
114 write(ilog,"(a,i4,a,i12)") ' Num of ',ic_type,':',je-js
115 write(* ,"(a,i4,a,i12)") ' Num of ',ic_type,':',je-js
116 enddo
117 nonzero = ndof2*(hecmat%NP + hecmat%NPU + hecmat%NPL)
118 write(ilog,"(a,i12)") ' Num of NZ :',nonzero
119 write(* ,"(a,i12)") ' Num of NZ :',nonzero
120 write(ilog,"(a,1pe12.5,a)") ' Sparsity :',100.0d0*dble(nonzero)/dble(ntdof2),"[%]"
121 write(* ,"(a,1pe12.5,a)") ' Sparsity :',100.0d0*dble(nonzero)/dble(ntdof2),"[%]"
122
123 !C
124 !C 3D
125 !C
126 if( hecmesh%n_dof .eq. 3 ) then
127 do ie = 1, hecmesh%n_elem
128 ia = hecmesh%elem_ID(ie*2)
129 if( ia.ne.hecmesh%my_rank ) cycle
130 ! je = hecMESH%elem_ID(ie*2-1)
131 jelem = hecmesh%global_elem_ID(ie)
132 !C
133 ic_type = hecmesh%elem_type(ie)
134 !C
135 if (.not. (hecmw_is_etype_rod(ic_type) .or. hecmw_is_etype_solid(ic_type) &
136 & .or. hecmw_is_etype_beam(ic_type) .or. hecmw_is_etype_shell(ic_type))) then
137 write(ilog,*) jelem, ' This Element cannot be checked. Type=',ic_type
138 cycle
139 endif
140 nn = hecmw_get_max_node(ic_type)
141 !C
142 js = hecmesh%elem_node_index(ie-1)
143 je = hecmesh%elem_node_index(ie)
144
145 do j = 1, nn
146 nodlocal(j) = hecmesh%elem_node_item (js+j)
147 xx(j) = hecmesh%node(3*nodlocal(j)-2)
148 yy(j) = hecmesh%node(3*nodlocal(j)-1)
149 zz(j) = hecmesh%node(3*nodlocal(j) )
150 enddo
151 !C
152 if ( ic_type.eq.111 ) then
153 isect = hecmesh%section_ID(ie)
154 mid = hecmesh%section%sect_mat_ID_item(isect)
155 call fstr_get_thickness( hecmesh,mid,aa )
156 al = sqrt( (xx(2)-xx(1))**2+(yy(2)-yy(1))**2+(zz(2)-zz(1))**2 )
157 nline = 1
158 tline = al
159 vol = aa*al
160 almax = al
161 almin = al
162 elseif( ic_type.eq.341 ) then
163 call pre_341 ( xx,yy,zz,vol,almax,almin )
164 elseif( ic_type.eq.351 ) then
165 call pre_351 ( xx,yy,zz,vol,almax,almin )
166 elseif( ic_type.eq.361 ) then
167 call pre_361 ( xx,yy,zz,vol,almax,almin )
168 elseif( ic_type.eq.342 ) then
169 call pre_342 ( xx,yy,zz,vol,almax,almin )
170 elseif( ic_type.eq.352 ) then
171 call pre_352 ( xx,yy,zz,vol,almax,almin )
172 elseif( ic_type.eq.362 ) then
173 call pre_362 ( xx,yy,zz,vol,almax,almin )
174 elseif( ic_type.eq.641 ) then
175 vol = 1.0d-12
176 elseif( ic_type.eq.761 ) then
177 vol = 1.0d-12
178 elseif( ic_type.eq.781 ) then
179 vol = 1.0d-12
180 endif
181 !C
182 if( vol.le.0.0 ) then
183 write(ilog,*) ' %%% ERROR %%% Volume of Element no.=',jelem,' is zero or negative.'
184 endif
185 nelem = nelem + 1
186 tvol = tvol + vol
187 if( vol.gt.tvmax ) tvmax = vol
188 if( vol.lt.tvmin ) tvmin = vol
189 if( almax.gt.tlmax ) tlmax = almax
190 if( almin.lt.tlmin ) tlmin = almin
191 asp = almax/almin
192 if( asp.gt.aspmax ) aspmax = asp
193 if( asp.gt.50 ) then
194 write(ilog,*) ' %%% WARNIG %%% Aspect ratio of Element no.=',jelem,' exceeds 50.'
195 write(ilog,*) ' Maximum length =',almax
196 write(ilog,*) ' Minimum length =',almin
197 endif
198 enddo
199 !C
200 !C 2D
201 !C
202 elseif( hecmesh%n_dof .eq. 2 ) then
203 do ie = 1, hecmesh%n_elem
204 ia = hecmesh%elem_ID(ie*2)
205 if( ia.ne.hecmesh%my_rank ) cycle
206 ! je = hecMESH%elem_ID(ie*2-1)
207 jelem = hecmesh%global_elem_ID(ie)
208 !C
209 ic_type = hecmesh%elem_type(ie)
210 !C
211 if (.not. (hecmw_is_etype_rod(ic_type) .or. hecmw_is_etype_surface(ic_type))) then
212 write(ilog,*) jelem, ' This Element cannot be checked. Type=',ic_type
213 cycle
214 endif
215 nn = hecmw_get_max_node(ic_type)
216 !C
217 js = hecmesh%elem_node_index(ie-1)
218 do j = 1, nn
219 nodlocal(j) = hecmesh%elem_node_item (js+j)
220 xx(j) = hecmesh%node(3*nodlocal(j)-2)
221 yy(j) = hecmesh%node(3*nodlocal(j)-1)
222 enddo
223 !C
224 isect = hecmesh%section_ID(ie)
225 mid = hecmesh%section%sect_mat_ID_item(isect)
226 call fstr_get_thickness( hecmesh,mid,aa )
227 !C
228 if ( ic_type.eq.111 ) then
229 al = sqrt( (xx(2)-xx(1))**2+(yy(2)-yy(1))**2 )
230 vol = aa*al
231 if( al.gt.tlmax ) tlmax = al
232 if( al.lt.tlmin ) tlmin = al
233 aspmax = 1.0
234 elseif( ic_type.eq.231 ) then
235 call pre_231 ( xx,yy,aa,vol,almax,almin )
236 elseif( ic_type.eq.241 ) then
237 call pre_241 ( xx,yy,aa,vol,almax,almin )
238 elseif( ic_type.eq.232 ) then
239 call pre_232 ( xx,yy,aa,vol,almax,almin )
240 elseif( ic_type.eq.242 ) then
241 call pre_242 ( xx,yy,aa,vol,almax,almin )
242 else
243 vol = 0.0
244 endif
245 !C
246 if( vol.le.0.0 ) then
247 write(ilog,*) ' %%% ERROR %%% Volume of Element no.=',jelem,' is zero or negative.'
248 endif
249 nelem = nelem + 1
250 tvol = tvol + vol
251 if( vol.gt.tvmax ) tvmax = vol
252 if( vol.lt.tvmin ) tvmin = vol
253 if( almax.gt.tlmax ) tlmax = almax
254 if( almin.lt.tlmin ) tlmin = almin
255 asp = almax/almin
256 if( asp.gt.aspmax ) aspmax = asp
257 if( asp.gt.50 ) then
258 write(ilog,*) ' %%% WARNIG %%% Aspect ratio of Element no.=',jelem,' exceeds 50.'
259 write(ilog,*) ' Maximum length =',almax
260 write(ilog,*) ' Minimum length =',almin
261 endif
262 enddo
263 !C
264 !C SHELL
265 !C
266 elseif( hecmesh%n_dof .eq. 6 ) then
267 do ie = 1, hecmesh%n_elem
268 ia = hecmesh%elem_ID(ie*2)
269 if( ia.ne.hecmesh%my_rank ) cycle
270 ! je = hecMESH%elem_ID(ie*2-1)
271 jelem = hecmesh%global_elem_ID(ie)
272 !C
273 ic_type = hecmesh%elem_type(ie)
274 !C
275 if (.not. (hecmw_is_etype_beam(ic_type) .or. hecmw_is_etype_shell(ic_type))) then
276 write(ilog,*) jelem, ' This Element cannot be checked. Type=',ic_type
277 cycle
278 endif
279 nn = hecmw_get_max_node(ic_type)
280 !C
281 js = hecmesh%elem_node_index(ie-1)
282 do j = 1, nn
283 nodlocal(j) = hecmesh%elem_node_item (js+j)
284 xx(j) = hecmesh%node(3*nodlocal(j)-2)
285 yy(j) = hecmesh%node(3*nodlocal(j)-1)
286 zz(j) = hecmesh%node(3*nodlocal(j) )
287 enddo
288 !C
289 isect = hecmesh%section_ID(ie)
290 mid = hecmesh%section%sect_mat_ID_item(isect)
291 call fstr_get_thickness( hecmesh,mid,aa )
292 !C
293 if ( ic_type.eq.111 ) then
294 al = sqrt( (xx(2)-xx(1))**2+(yy(2)-yy(1))**2+(zz(2)-zz(1))**2 )
295 nline = nline + 1
296 tline = tline + al
297 vol = aa*al
298 if( al.gt.tlmax ) tlmax = al
299 if( al.lt.tlmin ) tlmin = al
300 aspmax = 1.0
301 elseif( ic_type.eq.731 ) then
302 call pre_731 ( xx,yy,zz,aa,vol,almax,almin )
303 elseif( ic_type.eq.741 ) then
304 call pre_741 ( xx,yy,zz,aa,vol,almax,almin )
305 endif
306 !C
307 if( vol.le.0.0 ) then
308 write(ilog,*) ' %%% ERROR %%% Volume of Element no.=',jelem,' is zero or negative.'
309 endif
310 nelem = nelem + 1
311 tvol = tvol + vol
312 if( vol.gt.tvmax ) tvmax = vol
313 if( vol.lt.tvmin ) tvmin = vol
314 if( almax.gt.tlmax ) tlmax = almax
315 if( almin.lt.tlmin ) tlmin = almin
316 asp = almax/almin
317 if( asp.gt.aspmax ) aspmax = asp
318 if( asp.gt.50 ) then
319 write(ilog,*) ' %%% WARNIG %%% Aspect ratio of Element no.=',jelem,' exceeds 50.'
320 write(ilog,*) ' Maximum length =',almax
321 write(ilog,*) ' Minimum length =',almin
322 endif
323 enddo
324 endif
325 !C
326 avvol = tvol / nelem
327 write(ilog,*) '### Sub Summary ###'
328 write(ilog,*) ' Total Volumes in this region = ',tvol
329 write(ilog,*) ' Average Volume of elements = ',avvol
330 write(ilog,*) ' Maximum Volume of elements = ',tvmax
331 write(ilog,*) ' Minimum Volume of elements = ',tvmin
332 write(ilog,*) ' Maximum length of element edges = ',tlmax
333 write(ilog,*) ' Minimum length of element edges = ',tlmin
334
335 write(ilog,*) ' Maximum aspect ratio in this region = ',aspmax
336 totsum(1) = tvol
337 call hecmw_allreduce_r(hecmesh,totsum,1,hecmw_sum)
338 ntotsum(1) = nelem
339 call hecmw_allreduce_i(hecmesh,ntotsum,1,hecmw_sum)
340 totmax(1) = tvmax
341 totmax(2) = tlmax
342 totmax(3) = aspmax
343 call hecmw_allreduce_r(hecmesh,totmax,3,hecmw_max)
344 totmin(1) = tvmin
345 totmin(2) = tlmin
346 call hecmw_allreduce_r(hecmesh,totmin,2,hecmw_min)
347 if( hecmesh%my_rank .eq. 0 ) then
348 avvol = totsum(1) / ntotsum(1)
349 write(ilog,*) '### Global Summary ###'
350 write(ilog,*) ' TOTAL VOLUME = ',totsum(1)
351 write(*,*) ' TOTAL VOLUME = ',totsum(1)
352 write(ilog,*) ' AVERAGE VOLUME OF ELEMENTS = ',avvol
353 write(*,*) ' AVERAGE VOLUME OF ELEMENTS = ',avvol
354 write(ilog,*) ' MAXIMUM VOLUME OF ELEMENTS = ',totmax(1)
355 write(*,*) ' MAXIMUM VOLUME OF ELEMENTS = ',totmax(1)
356 write(ilog,*) ' MINIMUM VOLUME OF ELEMENTS = ',totmin(1)
357 write(*,*) ' MINIMUM VOLUME OF ELEMENTS = ',totmin(1)
358 write(ilog,*) ' MAXIMUM LENGTH OF ELEMENT EDGES = ',totmax(2)
359 write(*,*) ' MAXIMUM LENGTH OF ELEMENT EDGES = ',totmax(2)
360 write(ilog,*) ' MINIMUM LENGTH OF ELEMENT EDGES = ',totmin(2)
361 write(*,*) ' MINIMUM LENGTH OF ELEMENT EDGES = ',totmin(2)
362 write(ilog,*) ' MAXIMUM ASPECT RATIO = ',totmax(3)
363 write(*,*) ' MAXIMUM ASPECT RATIO = ',totmax(3)
364 endif
365 !C
366 end subroutine fstr_precheck_elem
367
368 subroutine hecmw_nonzero_profile(hecMESH, hecMAT)
369 use hecmw_util
370 implicit none
371 type (hecmwST_local_mesh) :: hecMESH
372 type (hecmwST_matrix) :: hecMAT
373
374 integer(kind=kint) :: i, j, in, jS, jE, ftype, n, ndof, nnz, fio
375 real(kind=kreal) :: rnum, dens, cond
376 character :: fileid*3
377
378 fio = 70 + hecmesh%my_rank
379 write(fileid,"(i3.3)")hecmesh%my_rank
380
381 !ftype = 2: eps
382 !ftype = 4: png
383 ftype = 4
384
385 n = hecmat%N
386 ndof = 3*n
387 nnz = 9*n + 9*2*hecmat%indexL(n)
388 dens = 100*dble(nnz)/dble(9*n*n)
389 rnum = (7.21d0+0.01*dlog10(dble(hecmat%N)))*10.0d0/dble(hecmat%N)
390 cond = 1.0d0
391 !rnum = (7.25d0)*10.0d0/dble(hecMAT%N)
392
393 open(fio,file='nonzero.dat.'//fileid, status='replace')
394 !write(fio,"(a,f12.5,i0)")"##magic number 10 : 7.2, ",rnum,hecMAT%N
395 do i= 1, n
396 js= hecmat%indexL(i-1) + 1
397 je= hecmat%indexL(i )
398 write(fio,"(i0,a,i0)")i," ",i
399 do j= js, je
400 in = hecmat%itemL(j)
401 write(fio,"(i0,a,i0)")i, " ",in
402 write(fio,"(i0,a,i0)")in," ",i
403 enddo
404 enddo
405 close(fio)
406
407 open(fio,file='nonzero.plt.'//fileid, status='replace')
408 if(ftype == 4)then
409 write(fio,"(a)")'set terminal png size 1500,1500'
410 else
411 write(fio,"(a)")'set terminal postscript eps enhanced color solid "TimesNewRomanPSMT" 20'
412 endif
413
414 write(fio,"(a)")'unset key'
415 write(fio,"(a)")'unset xtics'
416 write(fio,"(a)")'unset ytics'
417 write(fio,"(a)")'set size ratio 1.0'
418 write(fio,"(a)")'set border lw 1.0'
419 write(fio,"(a,i0,a)")'set xrange[0.5:',n,'.5]'
420 write(fio,"(a,i0,a)")'set yrange[0.5:',n,'.5] reverse '
421
422 if(ftype == 4)then
423 write(fio,"(a)")'set out "image.'//fileid//'.png"'
424 else
425 write(fio,"(a)")'set out "image.'//fileid//'.eps"'
426 write(fio,"(a)" )'set label 1 "Name" at graph 1.1,0.9'
427 write(fio,"(a)")'set label 2 "N" at graph 1.1,0.85'
428 write(fio,"(a)")'set label 3 "Non-Zero Elem." at graph 1.1,0.8'
429 write(fio,"(a)")'set label 4 "Density [%]" at graph 1.1,0.75'
430 write(fio,"(a)")'set label 9 "Condition Num." at graph 1.1,0.7'
431 write(fio,"(a)" )'set label 5 ": matrix" at graph 1.4,0.9'
432 write(fio,"(a,i0,a)")'set label 6 ": ',ndof,'" at graph 1.4,0.85'
433 write(fio,"(a,i0,a)")'set label 7 ": ',nnz,'" at graph 1.4,0.8'
434 write(fio,"(a,1pe9.2,a)")'set label 8 ": ',dens,'" at graph 1.4,0.75'
435 write(fio,"(a,1pe9.2,a)")'set label 10 ": ',cond,'" at graph 1.4,0.7'
436 endif
437
438 write(fio,"(a,f12.5,a)")'plot "nonzero.dat.'//fileid//'" pointtype 5 pointsize ',rnum,' linecolor rgb "#F96566"'
439 close(fio)
440
441 write(*,*)''
442 write(*,*)' ### Command recommendation'
443 write(*,*)' gnuplot -persist "nonzero.plt"'
444
445 !call system('gnuplot -persist "nonzero.plt"')
446 !open(fio,file='nonzero.dat',status='old')
447 !close(fio,status="delete")
448 end subroutine hecmw_nonzero_profile
449end module m_fstr_precheck
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal
Definition: hecmw.f90:6
This module provides function to check input data of IFSTR solver.
subroutine fstr_get_thickness(hecmesh, mid, thick)
subroutine hecmw_nonzero_profile(hecmesh, hecmat)
subroutine fstr_precheck_elem(hecmesh, hecmat)
subroutine fstr_precheck(hecmesh, hecmat, soltype)
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint) myrank
PARALLEL EXECUTION.
Definition: m_fstr.f90:80
integer(kind=kint), parameter imsg
Definition: m_fstr.f90:94
integer(kind=kint), parameter idbg
Definition: m_fstr.f90:95
integer(kind=kint), parameter ilog
FILE HANDLER.
Definition: m_fstr.f90:91
integer(kind=kint), parameter kstnzprof
Definition: m_fstr.f90:43
This module provides function to check input data of 2d static analysis.
subroutine pre_241(xx, yy, thick, vol, almax, almin)
subroutine pre_232(xx, yy, thick, vol, almax, almin)
subroutine pre_242(xx, yy, thick, vol, almax, almin)
subroutine pre_231(xx, yy, thick, vol, almax, almin)
This module provides function to check input data of 3d static analysis.
subroutine pre_352(xx, yy, zz, vol, almax, almin)
subroutine pre_342(xx, yy, zz, vol, almax, almin)
subroutine pre_362(xx, yy, zz, vol, almax, almin)
subroutine pre_361(xx, yy, zz, vol, almax, almin)
subroutine pre_341(xx, yy, zz, vol, almax, almin)
subroutine pre_351(xx, yy, zz, vol, almax, almin)
This module provides function to check input data of shell elements.
subroutine pre_731(xx, yy, zz, thick, vol, almax, almin)
subroutine pre_741(xx, yy, zz, thick, vol, almax, almin)