FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
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
105contains
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)
131 case (fe_line3n, fe_beam3n)
141 case ( fe_mitc9_shell )
145 case ( fe_tet10n, fe_tet10nc )
147 case ( fe_prism6n )
149 case ( fe_prism15n )
151 case ( fe_hex8n )
153 case ( fe_hex20n )
155 case default
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
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)
449 case ( fe_tri6n )
451 case ( fe_tri6nc )
453 case (fe_line3n )
457 case ( fe_quad8n, fe_mitc9_shell )
459 case ( fe_hex8n )
461 case ( fe_hex20n, fe_mitc8_shell )
462 numofquadpoints = 27
463 case ( fe_prism6n )
469 case ( fe_tet10n, fe_tet4n_pipi )
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)
503 pos(1:3)=gauss3d2(:,np)
504 case ( fe_hex20n, fe_mitc8_shell )
505 pos(1:3)=gauss3d3(:,np)
507 pos(1:3)=gauss3d7(:,np)
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)
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)
545 getweight = weight3d2(np)
546 case ( fe_hex20n)
547 getweight = weight3d3(np)
549 getweight = weight3d7(np)
550 case ( fe_prism15n )
551 getweight = weight3d8(np)
552 case ( fe_tet4n, fe_beam341 )
554 case ( fe_tet10n, fe_tet4n_pipi )
555 getweight = weight3d5(np)
556 case ( fe_line2n )
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))
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
1222end 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