FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
element.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 !-------------------------------------------------------------------------------
31 ! If you wish introduce new elements with new geometry or/and
32 ! new shape functions, you need do the followings
33 !!
34 !!
35 !!- Introduce new element ID corresponding to your element in this module.
36 !!- Provide corresponding shape function and shape derivative in a new
37 !! module and include this new module here.
38 !!- Add the information of your new element into functions listed above.
39 !!- If new quadrature method needed, do some modification in MODULE
40 !! Quadrature.
41 !!- If you introduce new surface geometry and do contact calculation also,
42 !! you may need do some modification in MODULE mSurfElement.
44 
45  use shape_line2n
46  use shape_line3n
47  use shape_tri3n
48  use shape_tri6n
49  use shape_quad4n
50  use shape_quad8n
51  use shape_quad9n
52  use shape_hex8n
53  use shape_hex20n
54  use shape_tet4n
55  use shape_tet10n
56  use shape_prism6n
57  use shape_prism15n
58  implicit none
59 
60  integer, parameter, private :: kreal = kind(0.0d0)
61 
62  !-------------------------------------------
63  ! Fllowing ID of element types
64  !-------------------------------------------
65  integer, parameter :: fe_unknown = -1
66 
67  integer, parameter :: fe_line2n = 111
68  integer, parameter :: fe_line3n = 112
69  integer, parameter :: fe_tri3n = 231
70  integer, parameter :: fe_tri6n = 232
71  integer, parameter :: fe_tri6nc = 2322
72  integer, parameter :: fe_quad4n = 241
73  integer, parameter :: fe_quad8n = 242
74  integer, parameter :: fe_truss = 301
75  integer, parameter :: fe_tet4n = 341
76  integer, parameter :: fe_tet4n_pipi = 3414
77  integer, parameter :: fe_tet10n = 342
78  integer, parameter :: fe_tet10nc = 3422
79  integer, parameter :: fe_prism6n = 351
80  integer, parameter :: fe_prism15n = 352
81  integer, parameter :: fe_hex8n = 361
82  integer, parameter :: fe_hex20n = 362
83  integer, parameter :: fe_hex27n = 363
84 
85  integer, parameter :: fe_beam2n = 611
86  integer, parameter :: fe_beam3n = 612
87  integer, parameter :: fe_beam341 = 641
88 
89  integer, parameter :: fe_tri6n_shell = 732
90  integer, parameter :: fe_dsg3_shell = 733
91  integer, parameter :: fe_mitc3_shell = 731
92  integer, parameter :: fe_mitc4_shell = 741
93  integer, parameter :: fe_mitc8_shell = 742
94  integer, parameter :: fe_mitc9_shell = 743
95 
96  integer, parameter :: fe_mitc3_shell361 = 761
97  integer, parameter :: fe_mitc4_shell361 = 781
98 
99  integer, parameter :: fe_tri3n_patch = 1031
100  integer, parameter :: fe_tri6n_patch = 1032
101  integer, parameter :: fe_quad4n_patch = 1041
102  integer, parameter :: fe_quad8n_patch = 1042
103  ! ---------------------------------------------
104 
105 contains
106 
107  !************************************
108  ! Following geometric information
109  !************************************
111  integer(kind=kind(2)) function getspacedimension( etype )
112  integer, intent(in) :: etype
113 
114  select case( etype)
115  case (fe_line2n, fe_line3n)
119  case default
121  end select
122  end function
123 
125  integer(kind=kind(2)) function getnumberofnodes( etype )
126  integer, intent(in) :: etype
127 
128  select case (etype)
129  case (fe_line2n, fe_beam2n)
130  getnumberofnodes = 2
131  case (fe_line3n, fe_beam3n)
132  getnumberofnodes = 3
134  getnumberofnodes = 3
136  getnumberofnodes = 6
138  getnumberofnodes = 4
140  getnumberofnodes = 8
141  case ( fe_mitc9_shell )
142  getnumberofnodes = 9
144  getnumberofnodes = 4
145  case ( fe_tet10n, fe_tet10nc )
146  getnumberofnodes = 10
147  case ( fe_prism6n )
148  getnumberofnodes = 6
149  case ( fe_prism15n )
150  getnumberofnodes = 15
151  case ( fe_hex8n )
152  getnumberofnodes = 8
153  case ( fe_hex20n )
154  getnumberofnodes = 20
155  case default
156  getnumberofnodes = -1
157  ! error message
158  end select
159  end function
160 
162  integer(kind=kind(2)) function getnumberofsubface( etype )
163  integer, intent(in) :: etype
164 
165  select case (etype)
166  case (fe_line2n, fe_line3n)
174  case ( fe_prism6n, fe_prism15n )
176  case ( fe_hex8n, fe_hex20n)
180  case default
181  getnumberofsubface = -1
182  ! error message
183  end select
184  end function
185 
187  subroutine getsubface( intype, innumber, outtype, nodes )
188  integer, intent(in) :: intype
189  integer, intent(in) :: innumber
190  integer, intent(out) :: outtype
191  integer, intent(out) :: nodes(:)
192 
193  if( innumber>getnumberofsubface( intype ) ) stop "Error in getting subface"
194  select case ( intype )
196  outtype = fe_tri3n
197  select case ( innumber )
198  case (1)
199  nodes(1)=1; nodes(2)=2; nodes(3)=3
200  case (2)
201  nodes(1)=4; nodes(2)=2; nodes(3)=1
202  case (3)
203  nodes(1)=4; nodes(2)=3; nodes(3)=2
204  case (4)
205  nodes(1)=4; nodes(2)=1; nodes(3)=3
206  end select
207  case (fe_tet10n)
208  outtype = fe_tri6n
209  select case ( innumber )
210  case (1)
211  nodes(1)=1; nodes(2)=2; nodes(3)=3
212  nodes(4)=5; nodes(5)=6; nodes(6)=7
213  case (2)
214  nodes(1)=4; nodes(2)=2; nodes(3)=1
215  nodes(4)=9; nodes(5)=5; nodes(6)=8
216  case (3)
217  nodes(1)=4; nodes(2)=3; nodes(3)=2
218  nodes(4)=10; nodes(5)=6; nodes(6)=9
219  case (4)
220  nodes(1)=4; nodes(2)=1; nodes(3)=3
221  nodes(4)=8; nodes(5)=7; nodes(6)=10
222  end select
223  case (fe_tet10nc)
224  outtype = fe_tri6nc
225  select case ( innumber )
226  case (1)
227  nodes(1)=1; nodes(2)=2; nodes(3)=3
228  nodes(4)=5; nodes(5)=6; nodes(6)=7
229  case (2)
230  nodes(1)=4; nodes(2)=2; nodes(3)=1
231  nodes(4)=9; nodes(5)=5; nodes(6)=8
232  case (3)
233  nodes(1)=4; nodes(2)=3; nodes(3)=2
234  nodes(4)=10; nodes(5)=6; nodes(6)=9
235  case (4)
236  nodes(1)=4; nodes(2)=1; nodes(3)=3
237  nodes(4)=8; nodes(5)=7; nodes(6)=10
238  end select
239  case ( fe_hex8n )
240  outtype = fe_quad4n
241  select case ( innumber )
242  case (1)
243  nodes(1)=1; nodes(2)=2; nodes(3)=3; nodes(4)=4
244  case (2)
245  nodes(1)=8; nodes(2)=7; nodes(3)=6; nodes(4)=5
246  case (3)
247  nodes(1)=5; nodes(2)=6; nodes(3)=2; nodes(4)=1
248  case (4)
249  nodes(1)=6; nodes(2)=7; nodes(3)=3; nodes(4)=2
250  case (5)
251  nodes(1)=7; nodes(2)=8; nodes(3)=4; nodes(4)=3
252  case (6)
253  nodes(1)=8; nodes(2)=5; nodes(3)=1; nodes(4)=4
254  case default
255  ! error
256  end select
257  case (fe_hex20n)
258  outtype = fe_quad8n
259  select case ( innumber )
260  case (1)
261  nodes(1)=1; nodes(2)=2; nodes(3)=3; nodes(4)=4
262  nodes(5)=9; nodes(6)=10; nodes(7)=11; nodes(8)=12
263  case (2)
264  nodes(1)=8; nodes(2)=7; nodes(3)=6; nodes(4)=5
265  nodes(5)=15; nodes(6)=14; nodes(7)=13; nodes(8)=16
266  case (3)
267  nodes(1)=5; nodes(2)=6; nodes(3)=2; nodes(4)=1
268  nodes(5)=13; nodes(6)=18; nodes(7)=9; nodes(8)=17
269  case (4)
270  nodes(1)=6; nodes(2)=7; nodes(3)=3; nodes(4)=2
271  nodes(5)=14; nodes(6)=19; nodes(7)=10; nodes(8)=18
272  case (5)
273  nodes(1)=7; nodes(2)=8; nodes(3)=4; nodes(4)=3
274  nodes(5)=15; nodes(6)=20; nodes(7)=11; nodes(8)=19
275  case (6)
276  nodes(1)=8; nodes(2)=5; nodes(3)=1; nodes(4)=4
277  nodes(5)=16; nodes(6)=17; nodes(7)=12; nodes(8)=20
278  case default
279  ! error
280  end select
281  case (fe_prism6n)
282  select case ( innumber )
283  case (1)
284  outtype = fe_tri3n
285  nodes(1)=1; nodes(2)=2; nodes(3)=3
286  case (2)
287  outtype = fe_tri3n
288  nodes(1)=6; nodes(2)=5; nodes(3)=4
289  case (3)
290  outtype = fe_quad4n
291  nodes(1)=4; nodes(2)=5; nodes(3)=2; nodes(4)=1
292  case (4)
293  outtype = fe_quad4n
294  nodes(1)=5; nodes(2)=6; nodes(3)=3; nodes(4)=2
295  case (5)
296  outtype = fe_quad4n
297  nodes(1)=6; nodes(2)=4; nodes(3)=1; nodes(4)=3
298  end select
299  case (fe_prism15n)
300  select case ( innumber )
301  case (1)
302  outtype = fe_tri6n
303  nodes(1)=1; nodes(2)=2; nodes(3)=3
304  nodes(4)=7; nodes(5)=8; nodes(6)=9
305  case (2)
306  outtype = fe_tri6n
307  nodes(1)=6; nodes(2)=5; nodes(3)=4
308  nodes(4)=11; nodes(5)=10; nodes(6)=12
309  case (3)
310  outtype = fe_quad8n
311  nodes(1)=4; nodes(2)=5; nodes(3)=2; nodes(4)=1
312  nodes(5)=10; nodes(6)=14; nodes(7)=7; nodes(8)=13
313  case (4)
314  outtype = fe_quad8n
315  nodes(1)=5; nodes(2)=6; nodes(3)=3; nodes(4)=2
316  nodes(5)=11; nodes(6)=15; nodes(7)=8; nodes(8)=14
317  case (5)
318  outtype = fe_quad8n
319  nodes(1)=6; nodes(2)=4; nodes(3)=1; nodes(4)=3
320  nodes(5)=12; nodes(6)=13; nodes(7)=9; nodes(8)=15
321  end select
322  case ( fe_tri3n, fe_mitc3_shell )
323  outtype = fe_line2n
324  select case (innumber )
325  case (1)
326  nodes(1) = 1; nodes(2)=2
327  case (2)
328  nodes(1) = 2; nodes(2)=3
329  case (3)
330  nodes(1) = 3; nodes(2)=1
331  end select
333  outtype = fe_line3n
334  select case (innumber )
335  case (1)
336  nodes(1) = 1; nodes(2)=2; nodes(3)=4
337  case (2)
338  nodes(1) = 2; nodes(2)=3; nodes(3)=5
339  case (3)
340  nodes(1) = 3; nodes(2)=1; nodes(3)=6
341  end select
342  case ( fe_quad4n, fe_mitc4_shell )
343  outtype = fe_line2n
344  select case (innumber )
345  case (1)
346  nodes(1) = 1; nodes(2)=2
347  case (2)
348  nodes(1) = 2; nodes(2)=3
349  case (3)
350  nodes(1) = 3; nodes(2)=4
351  case (4)
352  nodes(1) = 4; nodes(2)=1
353  end select
355  outtype = fe_line3n
356  select case (innumber )
357  case (1)
358  nodes(1) = 1; nodes(2)=2; nodes(3)=5
359  case (2)
360  nodes(1) = 2; nodes(2)=3; nodes(3)=6
361  case (3)
362  nodes(1) = 3; nodes(2)=4; nodes(3)=7
363  case (4)
364  nodes(1) = 4; nodes(2)=1; nodes(3)=8
365  end select
366  case (fe_mitc3_shell361)
367  select case ( innumber )
368  case (1)
369  outtype = fe_tri3n
370  nodes(1)=1; nodes(2)=2; nodes(3)=3
371  case (2)
372  outtype = fe_tri3n
373  nodes(1)=6; nodes(2)=5; nodes(3)=4
374  case (3)
375  outtype = fe_quad4n
376  nodes(1)=4; nodes(2)=5; nodes(3)=2; nodes(4)=1
377  case (4)
378  outtype = fe_quad4n
379  nodes(1)=5; nodes(2)=6; nodes(3)=3; nodes(4)=2
380  case (5)
381  outtype = fe_quad4n
382  nodes(1)=6; nodes(2)=4; nodes(3)=1; nodes(4)=3
383  end select
384  case ( fe_mitc4_shell361 )
385  outtype = fe_quad4n
386  select case ( innumber )
387  case (1)
388  nodes(1)=1; nodes(2)=2; nodes(3)=3; nodes(4)=4
389  case (2)
390  nodes(1)=8; nodes(2)=7; nodes(3)=6; nodes(4)=5
391  case (3)
392  nodes(1)=5; nodes(2)=6; nodes(3)=2; nodes(4)=1
393  case (4)
394  nodes(1)=6; nodes(2)=7; nodes(3)=3; nodes(4)=2
395  case (5)
396  nodes(1)=7; nodes(2)=8; nodes(3)=4; nodes(4)=3
397  case (6)
398  nodes(1)=8; nodes(2)=5; nodes(3)=1; nodes(4)=4
399  case default
400  ! error
401  end select
402  case ( fe_tri3n_patch )
403  outtype = fe_tri3n
404  select case ( innumber )
405  case (1)
406  nodes(1)=1; nodes(2)=2; nodes(3)=3
407  case default
408  !error
409  end select
410  case ( fe_tri6n_patch )
411  outtype = fe_tri6n
412  select case ( innumber )
413  case (1)
414  nodes(1)=1; nodes(2)=2; nodes(3)=3
415  nodes(4)=4; nodes(5)=5; nodes(6)=6
416  case default
417  !error
418  end select
419  case ( fe_quad4n_patch )
420  outtype = fe_quad4n
421  select case ( innumber )
422  case (1)
423  nodes(1)=1; nodes(2)=2; nodes(3)=3; nodes(4)=4
424  case default
425  !error
426  end select
427  case ( fe_quad8n_patch )
428  outtype = fe_quad8n
429  select case ( innumber )
430  case (1)
431  nodes(1)=1; nodes(2)=2; nodes(3)=3; nodes(4)=4
432  nodes(5)=5; nodes(6)=6; nodes(7)=7; nodes(8)=8
433  case default
434  !error
435  end select
436  case default
437  outtype = fe_unknown
438  stop "element type not defined-sbs"
439  ! error message
440  end select
441  end subroutine
442 
444  integer function numofquadpoints( fetype )
445  integer, intent(in) :: fetype
446  select case (fetype)
448  numofquadpoints = 1
449  case ( fe_tri6n )
450  numofquadpoints = 3
451  case ( fe_tri6nc )
452  numofquadpoints = 4
453  case (fe_line3n )
454  numofquadpoints = 2
456  numofquadpoints = 4
457  case ( fe_quad8n, fe_mitc9_shell )
458  numofquadpoints = 9
459  case ( fe_hex8n )
460  numofquadpoints = 8
461  case ( fe_hex20n, fe_mitc8_shell )
462  numofquadpoints = 27
463  case ( fe_prism6n )
464  numofquadpoints = 2
466  numofquadpoints = 3
467  case ( fe_prism15n, fe_tri6n_shell )
468  numofquadpoints = 9
469  case ( fe_tet10n, fe_tet4n_pipi )
470  numofquadpoints = 4
471  case ( fe_tet10nc )
472  numofquadpoints = 12
473  case default
474  numofquadpoints = -1
475  ! error message
476  stop "element type not defined-np"
477  end select
478  end function
479 
481  subroutine getquadpoint( fetype, np, pos )
482  use quadrature
483  integer, intent(in) :: fetype
484  integer, intent(in) :: np
485  real(kind=kreal), intent(out) :: pos(:)
486 
487  if( np<1 .or. np>numofquadpoints(fetype) ) then
488  ! error
489  endif
490 
491  select case (fetype)
492  case (fe_tri3n)
493  pos(1:2)=gauss2d4(:,np)
494  case ( fe_tri6n, fe_mitc3_shell )
495  pos(1:2)=gauss2d5(:,np)
496  case (fe_tri6nc )
497  pos(1:2)=gauss2d6(:,np)
498  case ( fe_quad4n, fe_mitc4_shell )
499  pos(1:2)=gauss2d2(:,np)
500  case ( fe_quad8n, fe_mitc9_shell )
501  pos(1:2)=gauss2d3(:,np)
502  case ( fe_hex8n, fe_mitc4_shell361 )
503  pos(1:3)=gauss3d2(:,np)
504  case ( fe_hex20n, fe_mitc8_shell )
505  pos(1:3)=gauss3d3(:,np)
506  case ( fe_prism6n, fe_mitc3_shell361 )
507  pos(1:3)=gauss3d7(:,np)
508  case ( fe_prism15n, fe_tri6n_shell )
509  pos(1:3)=gauss3d8(:,np)
510  case ( fe_tet4n, fe_beam341 )
511  pos(1:3)=gauss3d4(:,np)
512  case ( fe_tet10n, fe_tet4n_pipi )
513  pos(1:3)=gauss3d5(:,np)
514  case ( fe_tet10nc )
515  pos(1:3)=np
516  case ( fe_line2n )
517  pos(1:1)=gauss1d1(:,np)
518  case ( fe_line3n )
519  pos(1:1)=gauss1d2(:,np)
520  case default
521  ! error message
522  stop "element type not defined-qp"
523  end select
524  end subroutine
525 
527  real(kind=kreal) function getweight( fetype, np )
528  use quadrature
529  integer, intent(in) :: fetype
530  integer, intent(in) :: np
531  if( np<1 .or. np>numofquadpoints(fetype) ) then
532  ! error
533  endif
534 
535  select case (fetype)
536  case (fe_tri3n)
537  getweight = weight2d4(1)
538  case ( fe_tri6n, fe_mitc3_shell )
539  getweight = weight2d5(np)
540  case ( fe_quad4n, fe_mitc4_shell )
541  getweight = weight2d2(np)
542  case ( fe_quad8n, fe_mitc9_shell )
543  getweight = weight2d3(np)
544  case ( fe_hex8n, fe_mitc4_shell361 )
545  getweight = weight3d2(np)
546  case ( fe_hex20n)
547  getweight = weight3d3(np)
548  case ( fe_prism6n, fe_mitc3_shell361 )
549  getweight = weight3d7(np)
550  case ( fe_prism15n )
551  getweight = weight3d8(np)
552  case ( fe_tet4n, fe_beam341 )
553  getweight = weight3d4(1)
554  case ( fe_tet10n, fe_tet4n_pipi )
555  getweight = weight3d5(np)
556  case ( fe_line2n )
557  getweight = weight1d1(1)
558  case ( fe_line3n )
559  getweight = weight1d2(np)
560  case default
561  getweight = 0.d0
562  ! error message
563  end select
564  end function
565 
566  !************************************
567  ! Following shape function information
568  !************************************
570  subroutine getshapederiv( fetype, localcoord, shapederiv )
571  integer, intent(in) :: fetype
572  real(kind=kreal), intent(in) :: localcoord(:)
573  real(kind=kreal), intent(out) :: shapederiv(:,:)
574 
575  select case (fetype)
576  case ( fe_tri3n, fe_mitc3_shell )
577  !error check
578  call shapederiv_tri3n(shapederiv(1:3,1:2))
579  case (fe_tri6n)
580  !error check
581  call shapederiv_tri6n(localcoord,shapederiv(1:6,1:2) )
582  case ( fe_quad4n, fe_mitc4_shell )
583  !error check
584  call shapederiv_quad4n(localcoord,shapederiv(1:4,1:2))
585  case (fe_quad8n)
586  !error check
587  call shapederiv_quad8n(localcoord,shapederiv(1:8,1:2))
588  case ( fe_mitc9_shell )
589  !error check
590  call shapederiv_quad9n(localcoord,shapederiv(1:9,1:2))
592  ! error check
593  call shapederiv_hex8n(localcoord,shapederiv(1:8,1:3))
594  case (fe_hex20n)
595  ! error check
596  call shapederiv_hex20n(localcoord, shapederiv(1:20,1:3))
598  call shapederiv_prism6n(localcoord,shapederiv(1:6,1:3))
599  case (fe_prism15n)
600  call shapederiv_prism15n(localcoord,shapederiv(1:15,1:3))
602  ! error check
603  call shapederiv_tet4n(shapederiv(1:4,1:3))
604  case (fe_tet10n)
605  ! error check
606  call shapederiv_tet10n(localcoord,shapederiv(1:10,1:3))
607  case default
608  ! error message
609  stop "Element type not defined-sde"
610  end select
611  end subroutine
612 
614  subroutine getshape2ndderiv( fetype, localcoord, shapederiv )
615  integer, intent(in) :: fetype
616  real(kind=kreal), intent(in) :: localcoord(:)
617  real(kind=kreal), intent(out) :: shapederiv(:,:,:)
618 
619  select case (fetype)
620  case ( fe_tri3n, fe_mitc3_shell )
621  !error check
622  call shape2ndderiv_tri3n(shapederiv(1:3,1:2,1:2))
623  case (fe_tri6n)
624  !error check
625  call shape2ndderiv_tri6n(shapederiv(1:6,1:2,1:2))
626  case ( fe_quad4n, fe_mitc4_shell )
627  !error check
628  call shape2ndderiv_quad4n(shapederiv(1:4,1:2,1:2))
629  case (fe_quad8n)
630  !error check
631  call shape2ndderiv_quad8n(localcoord,shapederiv(1:8,1:2,1:2))
632  case default
633  ! error message
634  stop "Cannot calculate second derivatives of shape function"
635  end select
636  end subroutine
637 
639  subroutine getshapefunc( fetype, localcoord, func )
640  integer, intent(in) :: fetype
641  real(kind=kreal), intent(in) :: localcoord(:)
642  real(kind=kreal), intent(out) :: func(:)
643 
644  select case (fetype)
645  case ( fe_tri3n, fe_mitc3_shell )
646  !error check
647  call shapefunc_tri3n(localcoord,func(1:3))
648  case (fe_tri6n)
649  !error check
650  call shapefunc_tri6n(localcoord,func(1:6))
651  case ( fe_quad4n, fe_mitc4_shell )
652  !error check
653  call shapefunc_quad4n(localcoord,func(1:4))
654  case (fe_quad8n)
655  !error check
656  call shapefunc_quad8n(localcoord,func(1:8))
658  ! error check
659  call shapefunc_hex8n(localcoord,func(1:8))
660  case ( fe_mitc9_shell )
661  !error check
662  call shapefunc_quad9n(localcoord,func(1:9))
663  case (fe_hex20n)
664  ! error check
665  call shapefunc_hex20n(localcoord,func(1:20))
667  call shapefunc_prism6n(localcoord,func(1:6))
668  case (fe_prism15n)
669  call shapefunc_prism15n(localcoord,func(1:15))
671  ! error check
672  call shapefunc_tet4n(localcoord,func(1:4))
673  case (fe_tet10n)
674  ! error check
675  call shapefunc_tet10n(localcoord,func(1:10))
676  case (fe_line2n)
677  !error check
678  call shapefunc_line2n(localcoord,func(1:2))
679  case (fe_line3n)
680  !error check
681  call shapefunc_line3n(localcoord,func(1:3))
682  case default
683  stop "Element type not defined-sf"
684  ! error message
685  end select
686  end subroutine
687 
688 
689  ! (Gaku Hashimoto, The University of Tokyo, 2012/11/15) <
690  !####################################################################
691  subroutine getnodalnaturalcoord(fetype, nncoord)
692  !####################################################################
693 
694  integer, intent(in) :: fetype
695  real(kind = kreal), intent(out) :: nncoord(:, :)
696 
697  !--------------------------------------------------------------------
698 
699  select case( fetype )
701 
702  !error check
703  call nodalnaturalcoord_tri3n( nncoord(1:3, 1:2) )
704 
706 
707  !error check
708  call nodalnaturalcoord_quad4n( nncoord(1:4, 1:2) )
709 
710  case( fe_mitc9_shell )
711 
712  !error check
713  call nodalnaturalcoord_quad9n( nncoord(1:9, 1:2) )
714 
715  case default
716 
717  ! error message
718  stop "Element type not defined-sde"
719 
720  end select
721 
722  !--------------------------------------------------------------------
723 
724  return
725 
726  !####################################################################
727  end subroutine getnodalnaturalcoord
728  !####################################################################
729  ! > (Gaku Hashimoto, The University of Tokyo, 2012/11/15)
730 
731 
733  subroutine getglobalderiv( fetype, nn, localcoord, elecoord, det, gderiv )
734  integer, intent(in) :: fetype
735  integer, intent(in) :: nn
736  real(kind=kreal), intent(in) :: localcoord(:)
737  real(kind=kreal), intent(in) :: elecoord(:,:)
738  real(kind=kreal), intent(out) :: det
739  real(kind=kreal), intent(out) :: gderiv(:,:)
740 
741  real(kind=kreal) :: dum, xj(3,3), xji(3,3), deriv(nn,3)
742  integer :: nspace
743 
744  nspace = getspacedimension( fetype )
745  call getshapederiv( fetype, localcoord(:), deriv(1:nn,:) )
746 
747  if( nspace==2 ) then
748  xj(1:2,1:2)=matmul( elecoord(1:2,1:nn), deriv(1:nn,1:2) )
749  det=xj(1,1)*xj(2,2)-xj(2,1)*xj(1,2)
750  if( det==0.d0 ) stop "Math error in GetGlobalDeriv! Determinant==0.0"
751  dum=1.d0/det
752  xji(1,1)= xj(2,2)*dum
753  xji(1,2)=-xj(1,2)*dum
754  xji(2,1)=-xj(2,1)*dum
755  xji(2,2)= xj(1,1)*dum
756  else
757  ! JACOBI MATRIX
758  xj(1:3,1:3)= matmul( elecoord(1:3,1:nn), deriv(1:nn,1:3) )
759  !DETERMINANT OF JACOBIAN
760  det=xj(1,1)*xj(2,2)*xj(3,3) &
761  +xj(2,1)*xj(3,2)*xj(1,3) &
762  +xj(3,1)*xj(1,2)*xj(2,3) &
763  -xj(3,1)*xj(2,2)*xj(1,3) &
764  -xj(2,1)*xj(1,2)*xj(3,3) &
765  -xj(1,1)*xj(3,2)*xj(2,3)
766  if( det==0.d0 ) stop "Math error in GetGlobalDeriv! Determinant==0.0"
767  ! INVERSION OF JACOBIAN
768  dum=1.d0/det
769  xji(1,1)=dum*( xj(2,2)*xj(3,3)-xj(3,2)*xj(2,3) )
770  xji(1,2)=dum*(-xj(1,2)*xj(3,3)+xj(3,2)*xj(1,3) )
771  xji(1,3)=dum*( xj(1,2)*xj(2,3)-xj(2,2)*xj(1,3) )
772  xji(2,1)=dum*(-xj(2,1)*xj(3,3)+xj(3,1)*xj(2,3) )
773  xji(2,2)=dum*( xj(1,1)*xj(3,3)-xj(3,1)*xj(1,3) )
774  xji(2,3)=dum*(-xj(1,1)*xj(2,3)+xj(2,1)*xj(1,3) )
775  xji(3,1)=dum*( xj(2,1)*xj(3,2)-xj(3,1)*xj(2,2) )
776  xji(3,2)=dum*(-xj(1,1)*xj(3,2)+xj(3,1)*xj(1,2) )
777  xji(3,3)=dum*( xj(1,1)*xj(2,2)-xj(2,1)*xj(1,2) )
778  endif
779 
780  gderiv(1:nn,1:nspace)=matmul( deriv(1:nn,1:nspace), xji(1:nspace,1:nspace) )
781  end subroutine
782 
784  real(kind=kreal) function getdeterminant( fetype, nn, localcoord, elecoord )
785  integer, intent(in) :: fetype
786  integer, intent(in) :: nn
787  real(kind=kreal), intent(in) :: localcoord(:)
788  real(kind=kreal), intent(in) :: elecoord(:,:)
789 
790  real(kind=kreal) :: xj(3,3), deriv(nn,3)
791  integer :: nspace
792 
793  nspace = getspacedimension( fetype )
794  call getshapederiv( fetype, localcoord(:), deriv(1:nn,:) )
795 
796  if( nspace==2 ) then
797  xj(1:2,1:2)=matmul( elecoord(1:2,1:nn), deriv(1:nn,1:2) )
798  getdeterminant=xj(1,1)*xj(2,2)-xj(2,1)*xj(1,2)
799  else
800  xj(1:3,1:3)= matmul( elecoord(1:3,1:nn), deriv(1:nn,1:3) )
801  getdeterminant=xj(1,1)*xj(2,2)*xj(3,3) &
802  +xj(2,1)*xj(3,2)*xj(1,3) &
803  +xj(3,1)*xj(1,2)*xj(2,3) &
804  -xj(3,1)*xj(2,2)*xj(1,3) &
805  -xj(2,1)*xj(1,2)*xj(3,3) &
806  -xj(1,1)*xj(3,2)*xj(2,3)
807  endif
808 
809  end function
810 
812  subroutine getjacobian( fetype, nn, localcoord, elecoord, det, jacobian, inverse )
813  integer, intent(in) :: fetype
814  integer, intent(in) :: nn
815  real(kind=kreal), intent(in) :: localcoord(:)
816  real(kind=kreal), intent(in) :: elecoord(:,:)
817  real(kind=kreal), intent(out) :: det
818  real(kind=kreal), intent(out) :: jacobian(:,:)
819  real(kind=kreal), intent(out) :: inverse(:,:)
820 
821  real(kind=kreal) :: dum, deriv(nn,3)
822  integer :: nspace
823 
824  nspace = getspacedimension( fetype )
825  call getshapederiv( fetype, localcoord(:), deriv(1:nn,:) )
826 
827  if( nspace==2 ) then
828  jacobian(1:2,1:2)=matmul( elecoord(1:2,1:nn), deriv(1:nn,1:2) )
829  det=jacobian(1,1)*jacobian(2,2)-jacobian(2,1)*jacobian(1,2)
830  if( det==0.d0 ) stop "Math error in getJacobain! Determinant==0.0"
831  dum=1.0/det
832  inverse(1,1)= jacobian(2,2)*dum
833  inverse(1,2)=-jacobian(1,2)*dum
834  inverse(2,1)=-jacobian(2,1)*dum
835  inverse(2,2)= jacobian(1,1)*dum
836  else
837  ! JACOBI MATRIX
838  jacobian(1:3,1:3)= matmul( elecoord(1:3,1:nn), deriv(1:nn,1:3) )
839  !DETERMINANT OF JACOBIAN
840  det=jacobian(1,1)*jacobian(2,2)*jacobian(3,3) &
841  +jacobian(2,1)*jacobian(3,2)*jacobian(1,3) &
842  +jacobian(3,1)*jacobian(1,2)*jacobian(2,3) &
843  -jacobian(3,1)*jacobian(2,2)*jacobian(1,3) &
844  -jacobian(2,1)*jacobian(1,2)*jacobian(3,3) &
845  -jacobian(1,1)*jacobian(3,2)*jacobian(2,3)
846  if( det==0.d0 ) stop "Math error in getJacobain! Determinant==0.0"
847  ! INVERSION OF JACOBIAN
848  dum=1.d0/det
849  inverse(1,1)=dum*( jacobian(2,2)*jacobian(3,3)-jacobian(3,2)*jacobian(2,3) )
850  inverse(1,2)=dum*(-jacobian(1,2)*jacobian(3,3)+jacobian(3,2)*jacobian(1,3) )
851  inverse(1,3)=dum*( jacobian(1,2)*jacobian(2,3)-jacobian(2,2)*jacobian(1,3) )
852  inverse(2,1)=dum*(-jacobian(2,1)*jacobian(3,3)+jacobian(3,1)*jacobian(2,3) )
853  inverse(2,2)=dum*( jacobian(1,1)*jacobian(3,3)-jacobian(3,1)*jacobian(1,3) )
854  inverse(2,3)=dum*(-jacobian(1,1)*jacobian(2,3)+jacobian(2,1)*jacobian(1,3) )
855  inverse(3,1)=dum*( jacobian(2,1)*jacobian(3,2)-jacobian(3,1)*jacobian(2,2) )
856  inverse(3,2)=dum*(-jacobian(1,1)*jacobian(3,2)+jacobian(3,1)*jacobian(1,2) )
857  inverse(3,3)=dum*( jacobian(1,1)*jacobian(2,2)-jacobian(2,1)*jacobian(1,2) )
858  endif
859  end subroutine
860 
862  function surfacenormal( fetype, nn, localcoord, elecoord ) result( normal )
863  integer, intent(in) :: fetype
864  integer, intent(in) :: nn
865  real(kind=kreal), intent(in) :: localcoord(2)
866  real(kind=kreal), intent(in) :: elecoord(3,nn)
867  real(kind=kreal) :: normal(3)
868  real(kind=kreal) :: deriv(nn,2), gderiv(3,2)
869 
870  select case (fetype)
871  case (fe_tri3n)
872  !error check
873  call shapederiv_tri3n(deriv(1:3,1:2))
874  case (fe_tri6n)
875  !error check
876  call shapederiv_tri6n(localcoord,deriv(1:6,1:2))
877  case (fe_quad4n)
878  !error check
879  call shapederiv_quad4n(localcoord,deriv(1:4,1:2))
880  case (fe_quad8n)
881  !error check
882  call shapederiv_quad8n(localcoord,deriv(1:8,1:2))
883  case default
884  ! error message
885  normal =0.d0
886  return
887  end select
888 
889  gderiv = matmul( elecoord, deriv )
890  normal(1) = gderiv(2,1)*gderiv(3,2) - gderiv(3,1)*gderiv(2,2)
891  normal(2) = gderiv(3,1)*gderiv(1,2) - gderiv(1,1)*gderiv(3,2)
892  normal(3) = gderiv(1,1)*gderiv(2,2) - gderiv(2,1)*gderiv(1,2)
893  ! normal = normal/dsqrt(dot_product(normal, normal))
894  end function
895 
897  function edgenormal( fetype, nn, localcoord, elecoord ) result( normal )
898  integer, intent(in) :: fetype
899  integer, intent(in) :: nn
900  real(kind=kreal), intent(in) :: localcoord(1)
901  real(kind=kreal), intent(in) :: elecoord(2,nn)
902  real(kind=kreal) :: normal(2)
903  real(kind=kreal) :: deriv(nn,1), gderiv(2,1)
904 
905  select case (fetype)
906  case (fe_line2n)
907  !error check
908  call shapederiv_line2n(deriv(1:nn,:))
909  case (fe_line3n)
910  !error check
911  call shapederiv_line3n(localcoord,deriv(1:nn,:))
912  case default
913  ! error message
914  normal =0.d0
915  return
916  end select
917 
918  gderiv = matmul( elecoord, deriv )
919  normal(1) = -gderiv(2,1)
920  normal(2) = gderiv(1,1)
921  ! normal = normal/dsqrt(dot_product(normal, normal))
922  end function
923 
925  subroutine tangentbase( fetype, nn, localcoord, elecoord, tangent )
926  integer, intent(in) :: fetype
927  integer, intent(in) :: nn
928  real(kind=kreal), intent(in) :: localcoord(2)
929  real(kind=kreal), intent(in) :: elecoord(3,nn)
930  real(kind=kreal), intent(out) :: tangent(3,2)
931  real(kind=kreal) :: deriv(nn,2)
932 
933  select case (fetype)
934  case (fe_tri3n)
935  !error check
936  call shapederiv_tri3n(deriv(1:3,1:2))
937  case (fe_tri6n)
938  !error check
939  call shapederiv_tri6n(localcoord,deriv(1:6,1:2))
940  case (fe_tri6nc)
941  !error check
942  call shapederiv_tri6n(localcoord,deriv(1:6,1:2))
943  case (fe_quad4n)
944  !error check
945  call shapederiv_quad4n(localcoord,deriv(1:4,1:2))
946  case (fe_quad8n)
947  !error check
948  call shapederiv_quad8n(localcoord,deriv(1:8,1:2))
949  case default
950  ! error message
951  tangent =0.d0
952  return
953  end select
954 
955  tangent = matmul( elecoord, deriv )
956  end subroutine tangentbase
957 
959  subroutine curvature( fetype, nn, localcoord, elecoord, l2ndderiv, normal, curv )
960  integer, intent(in) :: fetype
961  integer, intent(in) :: nn
962  real(kind=kreal), intent(in) :: localcoord(2)
963  real(kind=kreal), intent(in) :: elecoord(3,nn)
964  real(kind=kreal), intent(out) :: l2ndderiv(3,2,2)
965  real(kind=kreal), intent(in), optional :: normal(3)
966  real(kind=kreal), intent(out), optional :: curv(2,2)
967  real(kind=kreal) :: deriv2(nn,2,2)
968 
969  select case (fetype)
970  case (fe_tri3n)
971  !error check
972  call shape2ndderiv_tri3n(deriv2(1:3,1:2,1:2))
973  case (fe_tri6n)
974  !error check
975  call shape2ndderiv_tri6n(deriv2(1:6,1:2,1:2))
976  case (fe_tri6nc)
977  !error check
978  call shape2ndderiv_tri6n(deriv2(1:6,1:2,1:2))
979  ! deriv2=0.d0
980  case (fe_quad4n)
981  !error check
982  call shape2ndderiv_quad4n(deriv2(1:4,1:2,1:2))
983  case (fe_quad8n)
984  !error check
985  call shape2ndderiv_quad8n(localcoord,deriv2(1:8,1:2,1:2))
986  case default
987  ! error message
988  stop "Cannot calculate second derivatives of shape function"
989  end select
990 
991  l2ndderiv(1:3,1,1) = matmul( elecoord(1:3,1:nn), deriv2(1:nn,1,1) )
992  l2ndderiv(1:3,1,2) = matmul( elecoord(1:3,1:nn), deriv2(1:nn,1,2) )
993  l2ndderiv(1:3,2,1) = matmul( elecoord(1:3,1:nn), deriv2(1:nn,2,1) )
994  l2ndderiv(1:3,2,2) = matmul( elecoord(1:3,1:nn), deriv2(1:nn,2,2) )
995  if( present(curv) ) then
996  curv(1,1) = dot_product( l2ndderiv(:,1,1), normal(:) )
997  curv(1,2) = dot_product( l2ndderiv(:,1,2), normal(:) )
998  curv(2,1) = dot_product( l2ndderiv(:,2,1), normal(:) )
999  curv(2,2) = dot_product( l2ndderiv(:,2,2), normal(:) )
1000  endif
1001  end subroutine curvature
1002 
1004  subroutine getelementcenter( fetype, localcoord )
1005  integer, intent(in) :: fetype
1006  real(kind=kreal), intent(out) :: localcoord(2)
1007 
1008  select case (fetype)
1009  case (fe_tri3n, fe_tri6n, fe_tri6nc)
1010  localcoord(:) = 1.d0/3.d0
1011  case (fe_quad4n, fe_quad8n)
1012  localcoord(:) = 0.d0
1013  case default
1014  ! error message
1015  localcoord(:) = 0.d0
1016  end select
1017  end subroutine getelementcenter
1018 
1021  integer function isinsideelement( fetype, localcoord, clearance )
1022  integer, intent(in) :: fetype
1023  real(kind=kreal), intent(inout) :: localcoord(2)
1024  real(kind=kreal), optional :: clearance
1025  real(kind=kreal) :: clr, coord3
1026 
1027  clr = 1.d-6
1028  if( present(clearance) ) clr = clearance
1029  if( dabs(localcoord(1))<clr ) localcoord(1)=0.d0
1030  if( dabs(localcoord(2))<clr ) localcoord(2)=0.d0
1031  if( dabs(dabs(localcoord(1))-1.d0)<clr ) &
1032  localcoord(1)=sign(1.d0,localcoord(1))
1033  if( dabs(dabs(localcoord(2))-1.d0)<clr ) &
1034  localcoord(2)=sign(1.d0,localcoord(2))
1035  isinsideelement = -1
1036  select case (fetype)
1037  case (fe_tri3n, fe_tri6n, fe_tri6nc)
1038  !error check
1039  coord3 = 1.d0-(localcoord(1)+localcoord(2))
1040  if( dabs(coord3)<clr ) coord3=0.d0
1041  if( localcoord(1)>=0.d0 .and. localcoord(1)<=1.d0 .and. &
1042  localcoord(2)>=0.d0 .and. localcoord(2)<=1.d0 .and. &
1043  coord3>=0.d0 .and. coord3<=1.d0 ) then
1044  isinsideelement = 0
1045  if( localcoord(1)==1.d0 ) then
1046  isinsideelement = 1
1047  elseif( localcoord(2)==1.d0 ) then
1048  isinsideelement = 2
1049  elseif( coord3==1.d0 ) then
1050  isinsideelement = 3
1051  elseif( coord3==0.d0 ) then
1052  isinsideelement = 12
1053  elseif( localcoord(1)==0.d0 ) then
1054  isinsideelement = 23
1055  elseif( localcoord(2)==0.d0 ) then
1056  isinsideelement = 31
1057  endif
1058  endif
1059  case (fe_quad4n, fe_quad8n)
1060  !error check
1061  if( all(dabs(localcoord)<=1.d0) ) then
1062  isinsideelement = 0
1063  if( localcoord(1)==-1.d0 .and. localcoord(2)==-1.d0 ) then
1064  isinsideelement = 1
1065  elseif( localcoord(1)==1.d0 .and. localcoord(2)==-1.d0 ) then
1066  isinsideelement = 2
1067  elseif( localcoord(1)==1.d0 .and. localcoord(2)==1.d0 ) then
1068  isinsideelement = 3
1069  elseif( localcoord(1)==-1.d0 .and. localcoord(2)==1.d0 ) then
1070  isinsideelement = 4
1071  elseif( localcoord(2)==-1.d0 ) then
1072  isinsideelement = 12
1073  elseif( localcoord(1)==1.d0 ) then
1074  isinsideelement = 23
1075  elseif( localcoord(2)==1.d0 ) then
1076  isinsideelement = 34
1077  elseif( localcoord(1)==-1.d0 ) then
1078  isinsideelement = 41
1079  endif
1080  endif
1081  end select
1082  end function isinsideelement
1083 
1085  subroutine getvertexcoord( fetype, cnode, localcoord )
1086  integer, intent(in) :: fetype
1087  integer, intent(in) :: cnode
1088  real(kind=kreal), intent(out) :: localcoord(2)
1089 
1090  select case (fetype)
1091  case (fe_tri3n, fe_tri6n, fe_tri6nc)
1092  if( cnode==1 ) then
1093  localcoord(1) =1.d0
1094  localcoord(2) =0.d0
1095  elseif( cnode==2 ) then
1096  localcoord(1) =0.d0
1097  localcoord(2) =1.d0
1098  else
1099  localcoord(1) =0.d0
1100  localcoord(2) =0.d0
1101  endif
1102  case (fe_quad4n, fe_quad8n)
1103  if( cnode==1 ) then
1104  localcoord(1) =-1.d0
1105  localcoord(2) =-1.d0
1106  elseif( cnode==2 ) then
1107  localcoord(1) =1.d0
1108  localcoord(2) =-1.d0
1109  elseif( cnode==3 ) then
1110  localcoord(1) =1.d0
1111  localcoord(2) =1.d0
1112  else
1113  localcoord(1) =-1.d0
1114  localcoord(2) =1.d0
1115  endif
1116  end select
1117  end subroutine
1118 
1120  subroutine extrapolatevalue( lpos, fetype, nnode, pvalue, ndvalue )
1121  real(kind=kreal), intent(in) :: lpos(:)
1122  integer, intent(in) :: fetype
1123  integer, intent(in) :: nnode
1124  real(kind=kreal), intent(in) :: pvalue(:)
1125  real(kind=kreal), intent(out) :: ndvalue(:,:)
1126 
1127  integer :: i
1128  real(kind=kreal) :: shapefunc(nnode)
1129  call getshapefunc( fetype, lpos, shapefunc )
1130  do i=1,nnode
1131  ndvalue(i,:) = shapefunc(i)*pvalue(:)
1132  enddo
1133  end subroutine
1134 
1136  subroutine interapolatevalue( lpos, fetype, nnode, pvalue, ndvalue )
1137  real(kind=kreal), intent(in) :: lpos(:)
1138  integer, intent(in) :: fetype
1139  integer, intent(in) :: nnode
1140  real(kind=kreal), intent(out) :: pvalue(:)
1141  real(kind=kreal), intent(in) :: ndvalue(:,:)
1142 
1143  integer :: i
1144  real(kind=kreal) :: shapefunc(nnode)
1145  call getshapefunc( fetype, lpos, shapefunc )
1146  pvalue(:) = 0
1147  do i=1,nnode
1148  pvalue(:) = pvalue(:)+ shapefunc(i)*ndvalue(i,:)
1149  enddo
1150  end subroutine
1151 
1153  subroutine gauss2node( fetype, gaussv, nodev )
1154  integer, intent(in) :: fetype
1155  real(kind=kreal), intent(in) :: gaussv(:,:)
1156  real(kind=kreal), intent(out) :: nodev(:,:)
1157 
1158  integer :: i, ngauss, nnode
1159  real(kind=kreal) :: localcoord(3), func(100)
1160  ngauss = numofquadpoints( fetype )
1161  nnode = getnumberofnodes( fetype )
1162  ! error checking
1163  select case (fetype)
1164  case (fe_tri3n)
1165  !error check
1166  forall(i=1:nnode)
1167  nodev(i,:) = gaussv(1,:)
1168  end forall
1169  case (fe_tri6n)
1170  !error check
1171  ! func(1:6) = ShapeFunc_tri6n(localcoord)
1172  case (fe_quad4n)
1173  !error check
1174  ! nodev(:,:) = gaussv(1,:)
1175  case (fe_quad8n)
1176  !error check
1177  call shapefunc_quad8n(localcoord,func(1:8))
1178  case (fe_hex8n, fe_mitc4_shell361)
1179  ! error check
1180  call shapefunc_hex8n(localcoord,func(1:8))
1181  case (fe_hex20n)
1182  ! error check
1183  call shapefunc_hex20n(localcoord,func(1:20))
1185  forall(i=1:3)
1186  nodev(i,:) = gaussv(1,:)
1187  end forall
1188  forall(i=1:3)
1189  nodev(i+3,:) = gaussv(2,:)
1190  end forall
1191  case (fe_prism15n)
1192  call shapefunc_prism15n(localcoord,func(1:15))
1194  ! error check
1195  forall(i=1:nnode)
1196  nodev(i,:) = gaussv(1,:)
1197  end forall
1198  case (fe_tet10n)
1199  ! error check
1200  call shapefunc_tet10n(localcoord,func(1:10))
1201  case default
1202  stop "Element type not defined"
1203  ! error message
1204  end select
1205  end subroutine
1206 
1208  real(kind=kreal) function getreferencelength( fetype, nn, localcoord, elecoord )
1209  integer, intent(in) :: fetype
1210  integer, intent(in) :: nn
1211  real(kind=kreal),intent(in) :: localcoord(2)
1212  real(kind=kreal),intent(in) :: elecoord(3,nn)
1213  real(kind=kreal) :: detjxy, detjyz, detjxz, detj
1214  detjxy = getdeterminant( fetype, nn, localcoord, elecoord(1:2,1:nn) )
1215  detjyz = getdeterminant( fetype, nn, localcoord, elecoord(2:3,1:nn) )
1216  detjxz = getdeterminant( fetype, nn, localcoord, elecoord(1:3:2,1:nn) )
1217  detj = dsqrt( detjxy **2 + detjyz **2 + detjxz **2 )
1218  getreferencelength = dsqrt( detj )
1219  end function getreferencelength
1220 
1221 
1222 end module
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
Definition: element.f90:640
subroutine getjacobian(fetype, nn, localcoord, elecoord, det, jacobian, inverse)
calculate Jacobian matrix, its determinant and inverse
Definition: element.f90:813
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
Definition: element.f90:126
integer, parameter fe_beam341
Definition: element.f90:87
integer, parameter fe_tri3n_patch
Definition: element.f90:99
integer, parameter fe_unknown
Definition: element.f90:65
real(kind=kreal) function, dimension(2) edgenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 2d-edge.
Definition: element.f90:898
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:482
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:734
integer, parameter fe_line2n
Definition: element.f90:67
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, parameter fe_prism6n
Definition: element.f90:79
integer, parameter fe_tet10nc
Definition: element.f90:78
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_dsg3_shell
Definition: element.f90:90
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate deivatives of shape fucntion in natural coordiante system.
Definition: element.f90:571
integer, parameter fe_mitc3_shell361
Definition: element.f90:96
integer, parameter fe_prism15n
Definition: element.f90:80
subroutine getvertexcoord(fetype, cnode, localcoord)
Get the natural coord of a vertex node.
Definition: element.f90:1086
integer, parameter fe_hex20n
Definition: element.f90:82
integer, parameter fe_quad8n_patch
Definition: element.f90:102
integer, parameter fe_tri3n
Definition: element.f90:69
real(kind=kreal) function getdeterminant(fetype, nn, localcoord, elecoord)
Calculate shape derivative in global coordinate system.
Definition: element.f90:785
subroutine tangentbase(fetype, nn, localcoord, elecoord, tangent)
Calculate base vector of tangent space of 3d surface.
Definition: element.f90:926
integer, parameter fe_mitc4_shell
Definition: element.f90:92
integer, parameter fe_hex27n
Definition: element.f90:83
integer, parameter fe_truss
Definition: element.f90:74
integer, parameter fe_mitc9_shell
Definition: element.f90:94
integer, parameter fe_tet4n_pipi
Definition: element.f90:76
subroutine getshape2ndderiv(fetype, localcoord, shapederiv)
Calculate the 2nd derivative of shape function in natural coodinate system.
Definition: element.f90:615
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:528
integer, parameter fe_mitc4_shell361
Definition: element.f90:97
integer, parameter fe_quad4n
Definition: element.f90:72
integer, parameter fe_mitc8_shell
Definition: element.f90:93
integer, parameter fe_hex8n
Definition: element.f90:81
integer, parameter fe_tri6n_patch
Definition: element.f90:100
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:445
integer(kind=kind(2)) function getspacedimension(etype)
Obtain the space dimension of the element.
Definition: element.f90:112
subroutine extrapolatevalue(lpos, fetype, nnode, pvalue, ndvalue)
This subroutine extrapolate a point value into elemental nodes.
Definition: element.f90:1121
integer function isinsideelement(fetype, localcoord, clearance)
if a point is inside a surface element -1: No; 0: Yes; >0: Node's (vertex) number
Definition: element.f90:1022
integer, parameter fe_mitc3_shell
Definition: element.f90:91
integer, parameter fe_tri6nc
Definition: element.f90:71
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
Definition: element.f90:863
subroutine getnodalnaturalcoord(fetype, nncoord)
Definition: element.f90:692
integer, parameter fe_beam2n
Definition: element.f90:85
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_line3n
Definition: element.f90:68
subroutine curvature(fetype, nn, localcoord, elecoord, l2ndderiv, normal, curv)
Calculate curvature tensor at a point along 3d surface.
Definition: element.f90:960
integer, parameter fe_tet10n
Definition: element.f90:77
integer, parameter fe_tri6n_shell
Definition: element.f90:89
subroutine gauss2node(fetype, gaussv, nodev)
This subroutine extroplate value in quadrature point to element nodes.
Definition: element.f90:1154
subroutine interapolatevalue(lpos, fetype, nnode, pvalue, ndvalue)
This subroutine interapolate element nodes value into a point value.
Definition: element.f90:1137
integer, parameter fe_quad4n_patch
Definition: element.f90:101
integer, parameter fe_beam3n
Definition: element.f90:86
integer, parameter fe_quad8n
Definition: element.f90:73
integer, parameter fe_tet4n
Definition: element.f90:75
This module contains Gauss point information.
Definition: quadrature.f90:28
real(kind=kreal), dimension(3, 9) gauss3d8
Definition: quadrature.f90:32
real(kind=kreal), dimension(2, 1) gauss2d4
Definition: quadrature.f90:32
real(kind=kreal), dimension(3, 4) gauss3d5
Definition: quadrature.f90:32
real(kind=kreal), dimension(1) weight2d4
Definition: quadrature.f90:32
real(kind=kreal), dimension(3, 8) gauss3d2
Definition: quadrature.f90:32
real(kind=kreal), dimension(3, 2) gauss3d7
Definition: quadrature.f90:32
real(kind=kreal), dimension(2, 9) gauss2d3
Definition: quadrature.f90:32
real(kind=kreal), dimension(2, 4) gauss2d2
Definition: quadrature.f90:32
real(kind=kreal), dimension(9) weight3d8
Definition: quadrature.f90:32
real(kind=kreal), dimension(3, 1) gauss3d4
Definition: quadrature.f90:32
real(kind=kreal), dimension(2) weight1d2
Definition: quadrature.f90:32
real(kind=kreal), dimension(4) weight3d5
Definition: quadrature.f90:32
real(kind=kreal), dimension(2, 3) gauss2d5
Definition: quadrature.f90:32
real(kind=kreal), dimension(1, 2) gauss1d2
Definition: quadrature.f90:32
real(kind=kreal), dimension(2, 4) gauss2d6
Definition: quadrature.f90:32
real(kind=kreal), dimension(27) weight3d3
Definition: quadrature.f90:32
real(kind=kreal), dimension(3, 27) gauss3d3
Definition: quadrature.f90:32
real(kind=kreal), dimension(8) weight3d2
Definition: quadrature.f90:32
real(kind=kreal), dimension(9) weight2d3
Definition: quadrature.f90:32
real(kind=kreal), dimension(4) weight2d2
Definition: quadrature.f90:32
real(kind=kreal), dimension(2) weight3d7
Definition: quadrature.f90:32
real(kind=kreal), dimension(1) weight3d4
Definition: quadrature.f90:32
real(kind=kreal), dimension(3) weight2d5
Definition: quadrature.f90:32
real(kind=kreal), dimension(1) weight1d1
Definition: quadrature.f90:32
real(kind=kreal), dimension(1, 1) gauss1d1
Definition: quadrature.f90:32
This module contains functions for interpolation in 20 node hexahedral element (Serendipity interpola...
Definition: hex20n.f90:7
subroutine shapefunc_hex20n(localcoord, func)
Definition: hex20n.f90:12
subroutine shapederiv_hex20n(localcoord, func)
Definition: hex20n.f90:41
This module contains functions for interpolation in 8 node hexahedral element (Langrange interpolatio...
Definition: hex8n.f90:7
subroutine shapederiv_hex8n(localcoord, func)
Definition: hex8n.f90:25
subroutine shapefunc_hex8n(localcoord, func)
Definition: hex8n.f90:12
This module contains functions for interpolation in 2 node line element (Langrange interpolation)
Definition: line2n.f90:7
subroutine shapefunc_line2n(lcoord, func)
Definition: line2n.f90:12
subroutine shapederiv_line2n(func)
Definition: line2n.f90:19
This module contains functions for interpolation in 3 nodes line element (Langrange interpolation)
Definition: line3n.f90:7
subroutine shapefunc_line3n(lcoord, func)
Definition: line3n.f90:12
subroutine shapederiv_line3n(lcoord, func)
Definition: line3n.f90:20
This module contains functions for interpolation in 15 node prism element (Langrange interpolation)
Definition: prism15n.f90:7
subroutine shapefunc_prism15n(ncoord, shp)
Definition: prism15n.f90:13
subroutine shapederiv_prism15n(ncoord, func)
Definition: prism15n.f90:37
This module contains functions for interpolation in 6 node prism element (Langrange interpolation)
Definition: prism6n.f90:7
subroutine shapefunc_prism6n(ncoord, func)
Definition: prism6n.f90:12
subroutine shapederiv_prism6n(ncoord, func)
Definition: prism6n.f90:27
This module contains functions for interpolation in 4 node qudrilateral element (Langrange interpolat...
Definition: quad4n.f90:7
subroutine shapederiv_quad4n(lcoord, func)
Definition: quad4n.f90:21
subroutine shape2ndderiv_quad4n(func)
Definition: quad4n.f90:35
subroutine shapefunc_quad4n(lcoord, func)
Definition: quad4n.f90:12
subroutine nodalnaturalcoord_quad4n(nncoord)
Definition: quad4n.f90:53
This module contains functions for interpolation in 8 node quadrilateral element (Serendipity interpo...
Definition: quad8n.f90:7
subroutine shape2ndderiv_quad8n(lcoord, func)
Definition: quad8n.f90:56
subroutine shapederiv_quad8n(lcoord, func)
Definition: quad8n.f90:33
subroutine shapefunc_quad8n(lcoord, func)
Definition: quad8n.f90:14
This module contains functions for interpolation in 9 node quadrilateral element.
Definition: quad9n.f90:7
subroutine shapederiv_quad9n(lcoord, func)
Definition: quad9n.f90:91
subroutine shapefunc_quad9n(lcoord, func)
Definition: quad9n.f90:23
subroutine nodalnaturalcoord_quad9n(nncoord)
Definition: quad9n.f90:174
This module contains functions for interpolation in 10 node tetrahedron element (Langrange interpolat...
Definition: tet10n.f90:7
subroutine shapefunc_tet10n(volcoord, shp)
Definition: tet10n.f90:12
subroutine shapederiv_tet10n(volcoord, shp)
Definition: tet10n.f90:30
This module contains functions for interpolation in 4 node tetrahedron element (Langrange interpolati...
Definition: tet4n.f90:7
subroutine shapefunc_tet4n(volcoord, func)
Definition: tet4n.f90:12
subroutine shapederiv_tet4n(func)
Definition: tet4n.f90:19
This module contains functions for interpolation in 3 node trianglar element (Langrange interpolation...
Definition: tri3n.f90:7
subroutine shape2ndderiv_tri3n(func)
Definition: tri3n.f90:30
subroutine nodalnaturalcoord_tri3n(nncoord)
Definition: tri3n.f90:38
subroutine shapefunc_tri3n(areacoord, func)
Definition: tri3n.f90:12
subroutine shapederiv_tri3n(func)
Definition: tri3n.f90:19
This module contains functions for interpolation in 6 node trianglar element (Langrange interpolation...
Definition: tri6n.f90:7
subroutine shapefunc_tri6n(areacoord, func)
Definition: tri6n.f90:12
subroutine shape2ndderiv_tri6n(func)
Definition: tri6n.f90:48
subroutine shapederiv_tri6n(areacoord, func)
Definition: tri6n.f90:26