7 use hecmw,
only : kint, kreal
20 subroutine stf_c2( ETYPE,NN,ecoord,gausses,PARAM1,stiff &
25 integer(kind=kint),
intent(in) :: ETYPE
26 integer(kind=kint),
intent(in) :: NN
28 real(kind=kreal),
intent(in) :: ecoord(2,nn),param1
29 real(kind=kreal),
intent(out) :: stiff(:,:)
30 integer(kind=kint),
intent(in) :: ISET
31 real(kind=kreal),
intent(in),
optional :: u(:,:)
34 integer(kind=kint) :: flag
35 integer(kind=kint) NDOF
37 real(kind=kreal) d(4,4),b(4,ndof*nn),db(4,ndof*nn)
38 real(kind=kreal) h(nn),stress(4)
39 real(kind=kreal) thick,pai
40 real(kind=kreal) det,rr,wg
41 real(kind=kreal) localcoord(2)
42 real(kind=kreal) gderiv(nn,2), cdsys(3,3)
43 integer(kind=kint) J,LX
44 real(kind=kreal) gdispderiv(2,2)
45 real(kind=kreal) b1(4,nn*ndof)
46 real(kind=kreal) smat(4,4)
47 real(kind=kreal) bn(4,nn*ndof), sbn(4,nn*ndof)
50 flag = gausses(1)%pMaterial%nlgeom_flag
60 call matlmatrix( gausses(lx), iset, d, 1.d0, 1.d0, cdsys, 0.d0 )
61 if( .not.
present(u) ) flag=infinitesimal
63 if( flag==1 .and. iset == 2 )
then
64 write(*,
'(a)')
' PROGRAM STOP : non-TL element for axisymmetric element'
73 rr=dot_product( h(1:nn), ecoord(1,1:nn) )
81 b(1,2*j-1)=gderiv(j,1)
83 b(3,2*j-1)=gderiv(j,2)
94 gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof,1:nn), gderiv(1:nn,1:ndof) )
95 b1(1:4,1:nn*ndof) = 0.d0
97 b1(1,2*j-1) = gdispderiv(1,1)*gderiv(j,1)
98 b1(2,2*j-1) = gdispderiv(1,2)*gderiv(j,2)
99 b1(3,2*j-1) = gdispderiv(1,2)*gderiv(j,1)+gdispderiv(1,1)*gderiv(j,2)
100 b1(1,2*j ) = gdispderiv(2,1)*gderiv(j,1)
101 b1(2,2*j ) = gdispderiv(2,2)*gderiv(j,2)
102 b1(3,2*j ) = gdispderiv(2,2)*gderiv(j,1)+gdispderiv(2,1)*gderiv(j,2)
107 b(:,j) = b(:,j)+b1(:,j)
111 db(1:4,1:nn*ndof) = 0.d0
112 db(1:4,1:nn*ndof) = matmul( d(1:4,1:4), b(1:4,1:nn*ndof) )
113 stiff(1:nn*ndof,1:nn*ndof) = stiff(1:nn*ndof,1:nn*ndof) + &
114 matmul( transpose( b(1:4,1:nn*ndof) ), db(1:4,1:nn*ndof) )*wg
118 stress(1:4)=gausses(lx)%stress(1:4)
119 bn(1:4,1:nn*ndof) = 0.d0
121 bn(1,2*j-1) = gderiv(j,1)
122 bn(2,2*j ) = gderiv(j,1)
123 bn(3,2*j-1) = gderiv(j,2)
124 bn(4,2*j ) = gderiv(j,2)
128 smat(j ,j ) = stress(1)
129 smat(j ,j+2) = stress(3)
130 smat(j+2,j ) = stress(3)
131 smat(j+2,j+2) = stress(2)
133 sbn(1:4,1:nn*ndof) = matmul( smat(1:4,1:4), bn(1:4,1:nn*ndof) )
134 stiff(1:nn*ndof,1:nn*ndof) = stiff(1:nn*ndof,1:nn*ndof) + &
135 matmul( transpose( bn(1:4,1:nn*ndof) ), sbn(1:4,1:nn*ndof) )*wg
144 subroutine dl_c2(ETYPE,NN,XX,YY,RHO,PARAM1,LTYPE,PARAMS,VECT,NSIZE,ISET)
158 integer(kind=kint),
intent(in) :: ETYPE,NN
159 real(kind=kreal),
intent(in) :: xx(nn),yy(nn),params(0:6)
160 real(kind=kreal),
intent(out) :: vect(nn*2)
161 real(kind=kreal) rho,param1
162 integer(kind=kint) LTYPE,NSIZE,ISET
164 integer(kind=kint),
parameter :: NDOF=2
165 real(kind=kreal) h(nn)
166 real(kind=kreal) plx(nn),ply(nn)
167 real(kind=kreal) xj(2,2),det,rr,wg
168 real(kind=kreal) pai,val,thick
169 integer(kind=kint) IVOL,ISUF
170 real(kind=kreal) ax,ay,rx,ry
171 real(kind=kreal) normal(2)
172 real(kind=kreal) coefx,coefy,xcod,ycod,hx,hy,phx,phy
173 integer(kind=kint) NOD(NN)
174 integer(kind=kint) LX,I,SURTYPE,NSUR
175 real(kind=kreal) vx,vy,localcoord(2),deriv(nn,2),elecoord(2,nn)
177 rx = 0.0d0; ry = 0.0d0
178 ay = 0.0d0; ax = 0.0d0
194 if( ltype.LT.10 )
then
196 if( ltype.EQ.5 )
then
202 else if( ltype.GE.10 )
then
204 call getsubface( etype, ltype/10, surtype, nod )
210 elecoord(1,i)=xx(nod(i))
211 elecoord(2,i)=yy(nod(i))
216 call getshapefunc( surtype, localcoord(1:1), h(1:nsur) )
217 normal=
edgenormal( surtype, nsur, localcoord(1:1), elecoord(:,1:nsur) )
222 rr=rr+h(i)*xx(nod(i))
229 vect(2*nod(i)-1)=vect(2*nod(i)-1)+val*wg*h(i)*normal(1)
230 vect(2*nod(i) )=vect(2*nod(i) )+val*wg*h(i)*normal(2)
242 xj(1,1:2)=matmul( xx(1:nn), deriv(1:nn,1:2) )
243 xj(2,1:2)=matmul( yy(1:nn), deriv(1:nn,1:2) )
245 det=xj(1,1)*xj(2,2)-xj(2,1)*xj(1,2)
249 rr=dot_product( h(1:nn),xx(1:nn) )
250 wg=wg*det*rr*2.d0*pai
259 xcod=dot_product( h(1:nn),xx(1:nn) )
260 ycod=dot_product( h(1:nn),yy(1:nn) )
261 hx=ax + ( (xcod-ax)*rx+(ycod-ay)*ry )/(rx**2+ry**2)*rx
262 hy=ay + ( (xcod-ax)*rx+(ycod-ay)*ry )/(rx**2+ry**2)*ry
265 coefx=rho*val*val*phx
266 coefy=rho*val*val*phy
269 plx(i)=plx(i)+h(i)*wg*coefx
270 ply(i)=ply(i)+h(i)*wg*coefy
275 vect(2*i-1)=val*plx(i)
278 else if( ltype.EQ.2 )
then
281 vect(2*i )=val*ply(i)
283 else if( ltype.EQ.4 )
then
286 vx=vx/sqrt( params(1)**2 + params(2)**2 )
287 vy=vy/sqrt( params(1)**2 + params(2)**2 )
289 vect(2*i-1)=val*plx(i)*rho*vx
290 vect(2*i )=val*ply(i)*rho*vy
292 else if( ltype==5 )
then
302 subroutine tload_c2(ETYPE,NN,XX,YY,TT,T0,gausses,PARAM1,ISET,VECT)
311 integer(kind=kint),
intent(in) :: ETYPE,NN
312 real(kind=kreal),
intent(in) :: xx(nn),yy(nn),tt(nn),t0(nn),param1
313 type(tgaussstatus),
intent(in) :: gausses(:)
314 real(kind=kreal),
intent(out) :: vect(nn*2)
315 integer(kind=kint) iset
317 integer(kind=kint) ndof
319 real(kind=kreal) alp,pp,d(4,4),b(4,ndof*nn)
320 real(kind=kreal) h(nn)
321 real(kind=kreal) eps(4),sgm(4),localcoord(2)
322 real(kind=kreal) deriv(nn,2),dnde(2,nn)
323 real(kind=kreal) xj(2,2),det,rr,dum
324 real(kind=kreal) xji(2,2)
325 real(kind=kreal) pai,thick,wg
326 integer(kind=kint) j,lx
327 real(kind=kreal) tempc,temp0,thermal_eps
337 if(iset==2) thick=1.d0
341 pp = gausses(1)%pMaterial%variables(
m_poisson)
347 xj(1,1:2)=matmul( xx(1:nn), deriv(1:nn,1:2) )
348 xj(2,1:2)=matmul( yy(1:nn), deriv(1:nn,1:2) )
350 det=xj(1,1)*xj(2,2)-xj(2,1)*xj(1,2)
352 tempc=dot_product(h(1:nn),tt(1:nn))
353 temp0=dot_product(h(1:nn),t0(1:nn))
355 rr=dot_product(h(1:nn),xx(1:nn))
356 wg=
getweight( etype, lx )*det*rr*2.d0*pai
362 xji(1,1)= xj(2,2)*dum
363 xji(1,2)=-xj(2,1)*dum
364 xji(2,1)=-xj(1,2)*dum
365 xji(2,2)= xj(1,1)*dum
367 dnde=matmul( xji, transpose(deriv) )
384 thermal_eps=alp*(tempc-temp0)
385 if( iset .EQ. 2 )
then
390 else if( iset.EQ.0 )
then
391 eps(1)=thermal_eps*(1.d0+pp)
392 eps(2)=thermal_eps*(1.d0+pp)
395 else if( iset.EQ.1 )
then
404 sgm = matmul( eps(1:4), d(1:4,1:4) )
408 vect(1:nn*ndof)=vect(1:nn*ndof)+matmul( sgm(1:4), b(1:4,1:nn*ndof) )*wg
416 subroutine update_c2( etype,nn,ecoord,gausses,PARAM1,iset, &
417 u, ddu, qf, TT, T0, TN )
423 integer(kind=kint),
intent(in) :: etype, nn
424 real(kind=kreal),
intent(in) :: ecoord(3,nn),param1
425 integer(kind=kint),
intent(in) :: iset
427 real(kind=kreal),
intent(in) :: u(2,nn)
428 real(kind=kreal),
intent(in) :: ddu(2,nn)
429 real(kind=kreal),
intent(out) :: qf(:)
430 real(kind=kreal),
intent(in) :: tt(nn)
431 real(kind=kreal),
intent(in) :: t0(nn)
432 real(kind=kreal),
intent(in) :: tn(nn)
435 integer(kind=kint),
parameter :: ndof=2
436 real(kind=kreal) :: d(4,4), b(4,ndof*nn)
437 real(kind=kreal) :: h(nn)
438 real(kind=kreal) :: thick,pai,rr
439 real(kind=kreal) :: det, wg
440 real(kind=kreal) :: localcoord(2)
441 real(kind=kreal) :: gderiv(nn,2), ttc, tt0, ttn
442 real(kind=kreal) :: cdsys(3,3)
443 integer(kind=kint) :: j, LX
444 real(kind=kreal) :: totaldisp(2*nn)
445 real(kind=kreal) :: epsth(4), alp, alp0, ina(1), outa(1)
450 totaldisp((j-1)*2+1) = u(1,j) + ddu(1,j)
451 totaldisp(j*2) = u(2,j) + ddu(2,j)
463 call matlmatrix( gausses(lx), iset, d, 1.d0, 1.d0, cdsys, 0.d0 )
470 ttc = dot_product(tt(:), h(:))
471 tt0 = dot_product(t0(:), h(:))
472 ttn = dot_product(tn(:), h(:))
473 if( dabs(ttc-tt0) > 1.d-14 )
then
475 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
476 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
479 call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
480 if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
488 rr=dot_product( h(1:nn), ecoord(1,1:nn) )
494 b(1,2*j-1)=gderiv(j,1)
496 b(3,2*j-1)=gderiv(j,2)
498 b(2,2*j )=gderiv(j,2)
499 b(3,2*j )=gderiv(j,1)
504 gausses(lx)%strain(1:4) = matmul( b(:,:), totaldisp )
505 gausses(lx)%stress(1:4) = matmul( d(1:4, 1:4), gausses(lx)%strain(1:4)-epsth(1:4) )
508 gausses(lx)%strain_out(1:4) = gausses(lx)%strain(1:4)
509 gausses(lx)%stress_out(1:4) = gausses(lx)%stress(1:4)
513 qf(1:nn*ndof) = qf(1:nn*ndof) + &
514 matmul( gausses(lx)%stress(1:4), b(1:4,1:nn*ndof) )*wg
517 gausses(lx)%strain_energy = 0.5d0*dot_product(gausses(lx)%stress(1:6),gausses(lx)%strain(1:6))*wg
530 integer(kind=kint),
intent(in) :: ETYPE,NN
532 real(kind=kreal),
intent(out) :: ndstrain(nn,4)
533 real(kind=kreal),
intent(out) :: ndstress(nn,4)
536 real(kind=kreal) :: temp(8)
541 temp(1:4) = temp(1:4) + gausses(i)%strain_out(1:4)
542 temp(5:8) = temp(5:8) + gausses(i)%stress_out(1:4)
544 temp(1:8) = temp(1:8)/ic
546 ndstrain(i,1:4) = temp(1:4)
547 ndstress(i,1:4) = temp(5:8)
559 integer(kind=kint),
intent(in) :: ETYPE
561 real(kind=kreal),
intent(out) :: strain(4)
562 real(kind=kreal),
intent(out) :: stress(4)
566 strain(:)=0.d0; stress(:)=0.d0
569 strain(:) = strain(:) + gausses(i)%strain_out(1:4)
570 stress(:) = stress(:) + gausses(i)%stress_out(1:4)
572 strain(:) = strain(:)/ic
573 stress(:) = stress(:)/ic
This module encapsulate the basic functions of all elements provide by this software.
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coordinate system.
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
real(kind=kreal) function, dimension(2) edgenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 2d-edge.
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
subroutine getsubface(intype, innumber, outtype, nodes)
Find the definition of surface of the element.
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate derivatives of shape function in natural coordinate system.
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
This module provides functions for elastic material.
subroutine calelasticmatrix(matl, sectType, D, temp, hdflag)
Calculate isotropic elastic matrix.
This module defines common data and basic structures for analysis.
real(kind=kreal), pointer ref_temp
REFTEMP.
This module manages calculation relates with materials.
subroutine matlmatrix(gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp, hdflag)
Calculate constituive matrix.
This module provides common functions of Plane deformation elements.
subroutine dl_c2(ETYPE, NN, XX, YY, RHO, PARAM1, LTYPE, PARAMS, VECT, NSIZE, ISET)
subroutine tload_c2(ETYPE, NN, XX, YY, TT, T0, gausses, PARAM1, ISET, VECT)
subroutine stf_c2(ETYPE, NN, ecoord, gausses, PARAM1, stiff, ISET, u)
subroutine update_c2(etype, nn, ecoord, gausses, PARAM1, iset, u, ddu, qf, TT, T0, TN)
Update strain and stress inside element.
subroutine elementstress_c2(ETYPE, gausses, strain, stress)
subroutine nodalstress_c2(ETYPE, NN, gausses, ndstrain, ndstress)
This module summarizes all information of material properties.
integer(kind=kint), parameter m_exapnsion
integer(kind=kint), parameter m_poisson
This modules defines a structure to record history dependent parameter in static analysis.
All data should be recorded in every quadrature points.