FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
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 !-------------------------------------------------------------------------------
7 contains
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
449 end 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 hecmw_nonzero_profile(hecMESH, hecMAT)
subroutine fstr_get_thickness(hecMESH, mid, thick)
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_242(XX, YY, thick, vol, almax, almin)
subroutine pre_231(XX, YY, thick, vol, almax, almin)
subroutine pre_232(XX, YY, thick, vol, almax, almin)
subroutine pre_241(XX, YY, thick, vol, almax, almin)
This module provides function to check input data of 3d static analysis.
subroutine pre_351(XX, YY, ZZ, vol, almax, almin)
subroutine pre_362(XX, YY, ZZ, vol, almax, almin)
subroutine pre_342(XX, YY, ZZ, vol, almax, almin)
subroutine pre_361(XX, YY, ZZ, vol, almax, almin)
subroutine pre_352(XX, YY, ZZ, vol, almax, almin)
subroutine pre_341(XX, YY, ZZ, vol, almax, almin)
This module provides function to check input data of shell elements.
subroutine pre_741(XX, YY, ZZ, thick, vol, almax, almin)
subroutine pre_731(XX, YY, ZZ, thick, vol, almax, almin)