FrontISTR  5.9.0
Large-scale structural analysis program with finit element method
precheck_LIB_elements.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 module m_precheck_lib_elements
8 
9  use hecmw
10  use elementinfo
11 
12  implicit none
13 
14  private
15  public :: precheck_calc_vol_asp
16 
17 contains
18 
20  subroutine precheck_calc_vol_asp(ic_type, nn, xx, yy, zz, thick, vol, almax, almin)
21  integer(kind=kint), intent(in) :: ic_type, nn
22  real(kind=kreal), intent(in) :: xx(:), yy(:), zz(:), thick
23  real(kind=kreal), intent(out) :: vol, almax, almin
24 
25  call calc_volume(ic_type, nn, xx, yy, zz, thick, vol)
26  call calc_edgelength(ic_type, nn, xx, yy, zz, almax, almin)
27  end subroutine
28 
30  subroutine calc_volume(ic_type, nn, xx, yy, zz, thick, vol)
31  integer(kind=kint), intent(in) :: ic_type, nn
32  real(kind=kreal), intent(in) :: xx(:), yy(:), zz(:), thick
33  real(kind=kreal), intent(out) :: vol
34 
35  integer :: ndim, nqp, ip, i
36  real(kind=kreal) :: pos(3), wt, det
37  real(kind=kreal) :: elecoord(3, nn)
38 
39  if (hecmw_is_etype_shell(ic_type)) then
40  call calc_volume_shell(ic_type, nn, xx, yy, zz, thick, vol)
41  return
42  endif
43 
44  ndim = getspacedimension(ic_type)
45 
46  elecoord = 0.0d0
47  do i = 1, nn
48  elecoord(1, i) = xx(i)
49  elecoord(2, i) = yy(i)
50  if (ndim >= 3) elecoord(3, i) = zz(i)
51  enddo
52 
53  nqp = numofquadpoints(ic_type)
54  vol = 0.0d0
55  do ip = 1, nqp
56  call getquadpoint(ic_type, ip, pos)
57  wt = getweight(ic_type, ip)
58  det = getdeterminant(ic_type, nn, pos, elecoord)
59  vol = vol + wt * det
60  enddo
61 
62  if (ndim == 2) vol = vol * thick
63  end subroutine
64 
67  subroutine calc_volume_shell(ic_type, nn, xx, yy, zz, thick, vol)
68  integer(kind=kint), intent(in) :: ic_type, nn
69  real(kind=kreal), intent(in) :: xx(:), yy(:), zz(:), thick
70  real(kind=kreal), intent(out) :: vol
71 
72  integer :: nqp, ip, i
73  real(kind=kreal) :: pos(3), wt
74  real(kind=kreal) :: deriv(nn, 3)
75  real(kind=kreal) :: g1(3), g2(3), g3(3), xsum, area
76 
77  nqp = numofquadpoints(ic_type)
78  area = 0.0d0
79  do ip = 1, nqp
80  call getquadpoint(ic_type, ip, pos)
81  wt = getweight(ic_type, ip)
82  call getshapederiv(ic_type, pos, deriv)
83 
84  g1 = 0.0d0
85  g2 = 0.0d0
86  do i = 1, nn
87  g1(1) = g1(1) + deriv(i,1) * xx(i)
88  g1(2) = g1(2) + deriv(i,1) * yy(i)
89  g1(3) = g1(3) + deriv(i,1) * zz(i)
90  g2(1) = g2(1) + deriv(i,2) * xx(i)
91  g2(2) = g2(2) + deriv(i,2) * yy(i)
92  g2(3) = g2(3) + deriv(i,2) * zz(i)
93  enddo
94 
95  ! Cross product magnitude = surface area element
96  g3(1) = g1(2)*g2(3) - g1(3)*g2(2)
97  g3(2) = g1(3)*g2(1) - g1(1)*g2(3)
98  g3(3) = g1(1)*g2(2) - g1(2)*g2(1)
99  xsum = sqrt(g3(1)**2 + g3(2)**2 + g3(3)**2)
100 
101  area = area + wt * xsum
102  enddo
103 
104  vol = area * thick
105  end subroutine
106 
110  subroutine calc_edgelength(ic_type, nn, xx, yy, zz, almax, almin)
111  integer(kind=kint), intent(in) :: ic_type, nn
112  real(kind=kreal), intent(in) :: xx(:), yy(:), zz(:)
113  real(kind=kreal), intent(out) :: almax, almin
114 
115  integer :: ndim
116 
117  ndim = getspacedimension(ic_type)
118  if (hecmw_is_etype_shell(ic_type)) ndim = 3
119 
120  almax = 0.0d0
121  almin = 1.0d20
122 
123  select case(ic_type)
124  ! 3D 1st order
125  case(341)
126  call upd1(1,2); call upd1(2,3); call upd1(1,3)
127  call upd1(1,4); call upd1(2,4); call upd1(3,4)
128  case(351)
129  call upd1(1,2); call upd1(2,3); call upd1(1,3)
130  call upd1(4,5); call upd1(5,6); call upd1(4,6)
131  call upd1(1,4); call upd1(2,5); call upd1(3,6)
132  case(361)
133  call upd1(1,2); call upd1(2,3); call upd1(3,4); call upd1(1,4)
134  call upd1(5,6); call upd1(6,7); call upd1(7,8); call upd1(5,8)
135  call upd1(1,5); call upd1(2,6); call upd1(3,7); call upd1(4,8)
136  ! 3D 2nd order
137  case(342)
138  call upd2(1,5,2); call upd2(2,6,3); call upd2(3,7,1)
139  call upd2(1,8,4); call upd2(2,9,4); call upd2(3,10,4)
140  case(352)
141  call upd2(1,7,2); call upd2(2,8,3); call upd2(1,9,3)
142  call upd2(4,10,5); call upd2(5,11,6); call upd2(4,12,6)
143  call upd2(1,13,4); call upd2(2,14,5); call upd2(3,15,6)
144  case(362)
145  call upd2(1,9,2); call upd2(2,10,3); call upd2(3,11,4); call upd2(4,12,1)
146  call upd2(5,13,6); call upd2(6,14,7); call upd2(7,15,8); call upd2(8,16,5)
147  call upd2(1,17,5); call upd2(2,18,6); call upd2(3,19,7); call upd2(4,20,8)
148  ! 2D 1st order
149  case(231)
150  call upd1(1,2); call upd1(2,3); call upd1(1,3)
151  case(241)
152  call upd1(1,2); call upd1(2,3); call upd1(3,4); call upd1(1,4)
153  ! 2D 2nd order
154  case(232)
155  call upd2(1,4,2); call upd2(2,5,3); call upd2(3,6,1)
156  case(242)
157  call upd2(1,5,2); call upd2(2,6,3); call upd2(3,7,4); call upd2(4,8,1)
158  ! Shell 1st order
159  case(731)
160  call upd1(1,2); call upd1(2,3); call upd1(1,3)
161  case(741)
162  call upd1(1,2); call upd1(2,3); call upd1(3,4); call upd1(1,4)
163  case default
164  almax = 0.0d0
165  almin = 0.0d0
166  end select
167 
168  contains
169 
170  subroutine upd1(n1, n2)
171  integer, intent(in) :: n1, n2
172  real(kind=kreal) :: al
173  if (ndim == 2) then
174  al = sqrt((xx(n2)-xx(n1))**2 + (yy(n2)-yy(n1))**2)
175  else
176  al = sqrt((xx(n2)-xx(n1))**2 + (yy(n2)-yy(n1))**2 + (zz(n2)-zz(n1))**2)
177  endif
178  if (al > almax) almax = al
179  if (al < almin) almin = al
180  end subroutine
181 
182  subroutine upd2(n1, nm, n2)
183  integer, intent(in) :: n1, nm, n2
184  real(kind=kreal) :: a1, a2, al
185  if (ndim == 2) then
186  a1 = sqrt((xx(nm)-xx(n1))**2 + (yy(nm)-yy(n1))**2)
187  a2 = sqrt((xx(n2)-xx(nm))**2 + (yy(n2)-yy(nm))**2)
188  else
189  a1 = sqrt((xx(nm)-xx(n1))**2 + (yy(nm)-yy(n1))**2 + (zz(nm)-zz(n1))**2)
190  a2 = sqrt((xx(n2)-xx(nm))**2 + (yy(n2)-yy(nm))**2 + (zz(n2)-zz(nm))**2)
191  endif
192  al = a1 + a2
193  if (al > almax) almax = al
194  if (al < almin) almin = al
195  end subroutine
196 
197  end subroutine calc_edgelength
198 
199 end module m_precheck_lib_elements
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
elementinfo::getquadpoint
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:489
elementinfo::getshapederiv
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate derivatives of shape function in natural coordinate system.
Definition: element.f90:578
elementinfo::getdeterminant
real(kind=kreal) function getdeterminant(fetype, nn, localcoord, elecoord)
Calculate shape derivative in global coordinate system.
Definition: element.f90:792
elementinfo::getweight
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:535
elementinfo::numofquadpoints
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:450
elementinfo::getspacedimension
integer(kind=kind(2)) function getspacedimension(etype)
Obtain the space dimension of the element.
Definition: element.f90:117
hecmw
Definition: hecmw.f90:6
m_precheck_lib_elements
This module provides unified element quality check functions using the elementInfo common API.
Definition: precheck_LIB_elements.f90:7
m_precheck_lib_elements::precheck_calc_vol_asp
subroutine, public precheck_calc_vol_asp(ic_type, nn, xx, yy, zz, thick, vol, almax, almin)
Calculate volume and edge lengths (for aspect ratio) for any supported element type.
Definition: precheck_LIB_elements.f90:21