CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
newton.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine newton (icalccg, ne, ipkon, lakon, kon, t0, co, rhcon, nrhcon, ntmat_, physcon, nelem, cgr, bodyf, ielmat, ithermal, vold, mi)
 

Function/Subroutine Documentation

◆ newton()

subroutine newton ( integer, intent(inout)  icalccg,
integer, intent(in)  ne,
integer, dimension(*), intent(in)  ipkon,
character*8, dimension(*), intent(in)  lakon,
integer, dimension(*), intent(in)  kon,
real*8, dimension(*), intent(in)  t0,
real*8, dimension(3,*), intent(in)  co,
real*8, dimension(0:1,ntmat_,*), intent(in)  rhcon,
integer, dimension(*), intent(in)  nrhcon,
integer, intent(in)  ntmat_,
real*8, dimension(*), intent(in)  physcon,
integer, intent(in)  nelem,
real*8, dimension(4,*), intent(inout)  cgr,
real*8, dimension(3), intent(inout)  bodyf,
integer, dimension(mi(3),*), intent(in)  ielmat,
integer, intent(in)  ithermal,
real*8, dimension(0:mi(2),*), intent(in)  vold,
integer, dimension(*), intent(in)  mi 
)
22 !
23 ! assigns the body forces to the elements by use of field ipobody
24 !
25  implicit none
26 !
27  character*8 lakon(*)
28 !
29  integer i,j,k,ne,icalccg,ipkon(*),nope,konl(20),kon(*),two,id,
30  & nrhcon(*),ntmat_,nelem,indexe,imat,mi(*),ielmat(mi(3),*),
31  & ithermal,iflag
32 !
33  real*8 xi,et,ze,weight,xl(3,20),shp(4,20),xsj,rho,cgr(4,*),
34  & t0l,t0(*),rhcon(0:1,ntmat_,*),physcon(*),co(3,*),dd,bodyf(3),
35  & vold(0:mi(2),*)
36 !
37  intent(in) ne,ipkon,lakon,kon,t0,co,rhcon,
38  & nrhcon,ntmat_,physcon,nelem,ielmat,ithermal,
39  & vold,mi
40 !
41  intent(inout) bodyf,cgr,icalccg
42 !
43 c data two /2/
44  two=2
45 !
46  if(icalccg.eq.0) then
47 !
48 ! first call: calculate the center of gravity of all elements
49 !
50  icalccg=1
51  do i=1,ne
52  if(ipkon(i).lt.0) cycle
53  if(lakon(i)(4:4).eq.'2') then
54  nope=20
55  xi=0.d0
56  et=0.d0
57  ze=0.d0
58  weight=8.d0
59  elseif(lakon(i)(4:4).eq.'8') then
60  nope=8
61  xi=0.d0
62  et=0.d0
63  ze=0.d0
64  weight=8.d0
65  elseif(lakon(i)(4:5).eq.'10') then
66  nope=10
67  xi=0.25d0
68  et=0.25d0
69  ze=0.25d0
70  weight=1.d0/6.d0
71  elseif(lakon(i)(4:4).eq.'4') then
72  nope=4
73  xi=0.25d0
74  et=0.25d0
75  ze=0.25d0
76  weight=1.d0/6.d0
77  elseif(lakon(i)(4:5).eq.'15') then
78  nope=15
79  xi=1.d0/3.d0
80  et=1.d0/3.d0
81  ze=0.d0
82  weight=1.d0
83  elseif(lakon(i)(4:4).eq.'6') then
84  nope=6
85  xi=1.d0/3.d0
86  et=1.d0/3.d0
87  ze=0.d0
88  weight=1.d0
89  else
90  cycle
91  endif
92 !
93 ! coordinates of the nodes in the deformed configuration
94 !
95  indexe=ipkon(i)
96  do j=1,nope
97  konl(j)=kon(indexe+j)
98  do k=1,3
99  xl(k,j)=co(k,konl(j))+vold(k,konl(j))
100  enddo
101  enddo
102 !
103 ! calculation of the shape functions
104 ! in the gauss point
105 !
106  iflag=1
107  if(nope.eq.20) then
108  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
109  elseif(nope.eq.8) then
110  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
111  elseif(nope.eq.10) then
112  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
113  elseif(nope.eq.4) then
114  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
115  elseif(nope.eq.15) then
116  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
117  else
118  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
119  endif
120 !
121 ! calculation of the center of gravity
122 !
123  do k=1,3
124  cgr(k,i)=0.d0
125  do j=1,nope
126  cgr(k,i)=cgr(k,i)+shp(4,j)*xl(k,j)
127  enddo
128  enddo
129 !
130 ! determining the density
131 !
132  imat=ielmat(1,nelem)
133  if(ithermal.eq.0) then
134  rho=rhcon(1,1,imat)
135  else
136 !
137 ! calculation of the initial temperature
138 !
139  t0l=0.d0
140  do j=1,nope
141  t0l=t0l+shp(4,j)*t0(konl(j))
142  enddo
143  call ident2(rhcon(0,1,imat),t0l,nrhcon(imat),two,id)
144  if(nrhcon(imat).eq.0) then
145  continue
146  elseif(nrhcon(imat).eq.1) then
147  rho=rhcon(1,1,imat)
148  elseif(id.eq.0) then
149  rho=rhcon(1,1,imat)
150  elseif(id.eq.nrhcon(imat)) then
151  rho=rhcon(1,id,imat)
152  else
153  rho=rhcon(1,id,imat)+
154  & (rhcon(1,id+1,imat)-rhcon(1,id,imat))*
155  & (t0l-rhcon(0,id,imat))/
156  & (rhcon(0,id+1,imat)-rhcon(0,id,imat))
157  endif
158  endif
159 !
160 ! coordinates of the nodes in the undeformed configuration
161 !
162  indexe=ipkon(i)
163  do j=1,nope
164  konl(j)=kon(indexe+j)
165  do k=1,3
166  xl(k,j)=co(k,konl(j))+vold(k,konl(j))
167  enddo
168  enddo
169 !
170 ! calculation of the Jacobian determinant
171 ! in the gauss point
172 !
173  iflag=2
174  if(nope.eq.20) then
175  if(lakon(i)(7:7).eq.'A') then
176  call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag)
177  elseif((lakon(i)(7:7).eq.'E').or.
178  & (lakon(i)(7:7).eq.'S')) then
179  call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag)
180  else
181  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
182  endif
183  elseif(nope.eq.8) then
184  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
185  elseif(nope.eq.10) then
186  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
187  elseif(nope.eq.4) then
188  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
189  elseif(nope.eq.15) then
190  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
191  else
192  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
193  endif
194 !
195 ! calculating the "constant" term in the gravity force
196 !
197  cgr(4,i)=physcon(3)*rho*xsj*weight
198 !
199  enddo
200  endif
201 !
202 ! calculating the force per unit mass
203 !
204  do j=1,3
205  bodyf(j)=0.d0
206  enddo
207 !
208  do i=1,ne
209  if(ipkon(i).lt.0) cycle
210  if(i.eq.nelem) cycle
211  dd=(cgr(1,i)-cgr(1,nelem))**2+(cgr(2,i)-cgr(2,nelem))**2+
212  & (cgr(3,i)-cgr(3,nelem))**2
213  do j=1,3
214  bodyf(j)=bodyf(j)+(cgr(j,i)-cgr(j,nelem))*cgr(4,i)/
215  & (dd*dsqrt(dd))
216  enddo
217  enddo
218 c write(*,*) 'newton',nelem,bodyf(1),bodyf(2),bodyf(3)
219 !
220  return
shape20h_pl
subroutine shape20h_pl(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h_pl.f:20
shape6w
subroutine shape6w(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape6w.f:20
shape10tet
subroutine shape10tet(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape10tet.f:20
shape8h
subroutine shape8h(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape8h.f:20
ident2
subroutine ident2(x, px, n, ninc, id)
Definition: ident2.f:27
shape15w
subroutine shape15w(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape15w.f:20
shape20h_ax
subroutine shape20h_ax(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h_ax.f:20
shape20h
subroutine shape20h(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape20h.f:20
shape4tet
subroutine shape4tet(xi, et, ze, xl, xsj, shp, iflag)
Definition: shape4tet.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)