Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
v3d10_ale.f90
Go to the documentation of this file.
1 !*********************************************************************
2 !* Illinois Open Source License *
3 !* *
4 !* University of Illinois/NCSA *
5 !* Open Source License *
6 !* *
7 !* Copyright@2008, University of Illinois. All rights reserved. *
8 !* *
9 !* Developed by: *
10 !* *
11 !* Center for Simulation of Advanced Rockets *
12 !* *
13 !* University of Illinois *
14 !* *
15 !* www.csar.uiuc.edu *
16 !* *
17 !* Permission is hereby granted, free of charge, to any person *
18 !* obtaining a copy of this software and associated documentation *
19 !* files (the "Software"), to deal with the Software without *
20 !* restriction, including without limitation the rights to use, *
21 !* copy, modify, merge, publish, distribute, sublicense, and/or *
22 !* sell copies of the Software, and to permit persons to whom the *
23 !* Software is furnished to do so, subject to the following *
24 !* conditions: *
25 !* *
26 !* *
27 !* @ Redistributions of source code must retain the above copyright *
28 !* notice, this list of conditions and the following disclaimers. *
29 !* *
30 !* @ Redistributions in binary form must reproduce the above *
31 !* copyright notice, this list of conditions and the following *
32 !* disclaimers in the documentation and/or other materials *
33 !* provided with the distribution. *
34 !* *
35 !* @ Neither the names of the Center for Simulation of Advanced *
36 !* Rockets, the University of Illinois, nor the names of its *
37 !* contributors may be used to endorse or promote products derived *
38 !* from this Software without specific prior written permission. *
39 !* *
40 !* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
41 !* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
42 !* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
43 !* NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
44 !* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
45 !* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
46 !* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
47 !* USE OR OTHER DEALINGS WITH THE SOFTWARE. *
48 !*********************************************************************
49 !* Please acknowledge The University of Illinois Center for *
50 !* Simulation of Advanced Rockets in works and publications *
51 !* resulting from this software or its derivatives. *
52 !*********************************************************************
53 !-----------------------------------------------------------------
54 ! CODE FOR 3-D ALE, ELASTODYNAMIC, SECOND-ORDER SOLID ELEMENT
55 !
56 ! ALE formulation for 10-node tetrahedron element
57 !
58 ! written by Changyu Hwang Nov. 20, 2001
59 !
60 ! reference: Finite element procedures in engineering analysis
61 ! by K. Bathe, Figure 5.27 pp 232
62 !-----------------------------------------------------------------------
63 !
64 !
65 SUBROUTINE v3d10_ale(v_bar,a_bar,d,vhalf,Rnet, &
66  e,xnu,rho,numnp,numat_vol, &
67  numlstet,matlstet,lmlstet,meshcoor, &
68  nstart,nend)
69 
70  IMPLICIT NONE
71 
72 !--- for tetrahedron element !
73 
74  INTEGER numnp ! number of nodes
75  INTEGER numat_vol ! number of materials
76  REAL*8 v_bar(3*numnp) ! nodal mesh velocities
77  REAL*8 a_bar(3*numnp) ! nodal mesh accelerations
78  REAL*8 d(3*numnp) ! nodal displacements
79  REAL*8 vhalf(3*numnp) ! nodal velocities
80  REAL*8 rnet(3*numnp) ! internal force
81  REAL*8 e(numat_vol) ! young's moduli
82  REAL*8 xnu(numat_vol) ! Poisson's ratios
83  REAL*8 rho(numat_vol) ! densities
84  REAL*8 meshcoor(3,numnp)
85  INTEGER numlstet ! number of LSTs
86  INTEGER matlstet(numlstet) ! list of all materials
87  INTEGER lmlstet(10,numlstet) ! TET conectivity table
88  INTEGER ielem , nx, ny, nz
89 
90 ! integer ndof, nnode, ndim, nmat, nsurf, nload, ldflag(nload)
91  INTEGER ndof, nnode, ndim, nmat, nsurf, nload
92 
93 !-- for 10-node tetrahedron element !
94  REAL*8 meshpos(3,10), meshvel(3,10), meshacc(3,10), mat(3), &
95  sload(4,1,3), bforce(3), elmass(30,30), elstiff(30,30), &
96  eldamp(30,30), elload(30), elstiff2(30,30)
97 
98 ! INPUT:
99 
100 ! ndof -- # of element degrees of freedom
101 ! nnode -- # of nodes on element
102 ! ndim -- # of spatial dimensions
103 ! nmat -- # of material parameters (rho, E, nu)
104 ! nsurf -- # of surfaces on element
105 ! nload -- # of load types for element (2 - body force, surface traction,
106 ! both distributions assumed uniform)
107 ! meshpos -- coordinates of position w.r.t. to global basis of nodes
108 ! meshvel -- components of mesh velocity w.r.t to global basis at nodes
109 ! meshacc -- '' '' '' accelern. '' '' '' '' '' ''
110 ! mat -- array containing material parameters (e.g. rho, E, nu)
111 ! ldflag -- array indicating if there are body forces and surface loads
112 ! (in that order) to be taken care of in this element
113 ! (values:1/0)
114 ! sload -- array containing surface loads
115 ! (face#, pos or neg if load present or absent, load vec comp)
116 ! bforce -- body force components
117 
118 ! OUTPUT:
119 
120 ! elmass -- element mass matrix
121 ! elstiff -- element stiffness matrix
122 ! eldamp -- element damping matrix
123 ! elload -- element load vector
124 
125  INTEGER ndim_v, ndim_s, nintk_v, nintk_s,knode, &
126  knodr, knodc,i,j,rindx,cindx,k,mm,ll, &
127  fp(2,6),iface,igauss, grindx,grindy,grindz
128 !
129 !--- variables for integration points
130 ! parameter(ndim_v = 3, ndim_s = 2, nintk_v = 8,nintk_s = 4) ! 8-noded brick
131 ! parameter(ndim_v = 3, ndim_s = 2, nintk_v = 4,nintk_s = 1) ! 4-noded tetrahedron
132  parameter(ndim_v = 3, ndim_s = 2, nintk_v = 4,nintk_s = 3) ! 10-noded tetrahedron
133 
134  DOUBLE PRECISION shfn(10,4), shdx(3,10,4), dv(10), dvsum, &
135  delta(3,3), mvgss(3), magss(3), mcgss(3,3), two, &
136  wsp1,wsparr(3),trmcgss,wsp2,wsp3,wsp4,one,lambda, &
137  mu,third,wsp5,vxi(10),veta(10),vzeta(10),wvar, &
138  sxi(4,6),seta(4,6),szeta(4,6),shdxi(3,10,4), &
139  wsparr1(3),shdx_av(3,10),zero, half
140  DOUBLE PRECISION disp(30), velo(30), rimsi(30), dimsi(30)
141  DOUBLE PRECISION multiplier , ajacin(3,3), surfnormal(3,4)
142  DOUBLE PRECISION wt_lstet, da, sum,factor, wt_cstet
143  DOUBLE PRECISION sto_wsp1(10),sto_wsp2(10)
144  DOUBLE PRECISION sto_wsp3(10),sto_wsp4(10)
145  INTEGER vect1(2,10), vect2(2,10), surf_elem(3,4), k2
146  REAL*8 r_in(3*numnp) ! internal force
147  REAL*8 r_dp(3*numnp) ! damping force
148 
149  INTEGER :: nstart, nend
150 
151  parameter(zero = 0.d0, one = 1.d0, two = 2.d0)
152 
153  DATA delta/1.d0,0.d0,0.d0,0.d0,1.d0,0.d0,0.d0,0.d0,1.d0/
154 
155 
156  ndim = 3
157  nnode = 10
158  ndof = ndim * nnode
159  nsurf = 4
160  multiplier = 1.d0/6.d0
161  wt_lstet = 1.d0/4.d0
162  wvar = 1.d0/3.d0
163 
164 
165 ! reference: Finite element procedures in engineering analysis
166 ! by K. Bathe, Figure 5.27 pp 232
167 
168 ! gauss point 1 alpha beta beta beta 1/4
169 ! gauss point 2 beta alpha beta beta 1/4
170 ! gauss point 3 beta beta alpha beta 1/4
171 ! gauss point 4 beta beta beta alpha 1/4
172 ! alpha = 0.58541020d0
173 ! beta = 0.13819660d0
174 
175 !--- gauss point 1
176  vxi( 1) = 0.13819660d0
177  veta( 1) = 0.13819660d0
178  vzeta(1) = 0.13819660d0
179 
180 !--- gauss point 2
181  vxi( 2) = 0.58541020d0
182  veta( 2) = 0.13819660d0
183  vzeta(2) = 0.13819660d0
184 
185 !--- gauss point 3
186  vxi( 3) = 0.13819660d0
187  veta( 3) = 0.58541020d0
188  vzeta(3) = 0.13819660d0
189 
190 !--- gauss point 4
191  vxi( 4) = 0.13819660d0
192  veta( 4) = 0.13819660d0
193  vzeta(4) = 0.58541020d0
194 
195 !-- three points on each surface: outward normal direction
196 !-- face 1 : node 1 3 2
197 !-- face 2 : node 1 2 4
198 !-- face 3 : node 2 3 4
199 !-- face 4 : node 1 4 3
200 
201  wvar = 1.d0/dsqrt(3.d0)
202  third = 1.d0/3.d0
203  half = 0.5d0
204 
205 !-- face 1 : node 1 3 2 (szeta=t=0)
206  sxi( 1,1) = half
207  seta( 1,1) = zero
208  szeta(1,1) = zero
209  sxi( 2,1) = zero
210  seta( 2,1) = half
211  szeta(2,1) = zero
212  sxi( 3,1) = half
213  seta( 3,1) = half
214  szeta(3,1) = zero
215 
216 !-- face 2 : node 1 2 4 (seta=s=0)
217  sxi( 1,2) = half
218  seta( 1,2) = zero
219  szeta(1,2) = zero
220  sxi( 2,2) = zero
221  seta( 2,2) = zero
222  szeta(2,2) = half
223  sxi( 3,2) = half
224  seta( 3,2) = zero
225  szeta(3,2) = half
226 
227 !-- face 3 : node 3 4 2
228  sxi( 1,3) = half
229  seta( 1,3) = half
230  szeta(1,3) = zero
231  sxi( 2,3) = half
232  seta( 2,3) = zero
233  szeta(2,3) = half
234  sxi( 3,3) = zero
235  seta( 3,3) = half
236  szeta(3,3) = half
237 
238 !-- face 4 : node 1 4 3 (sxi=r=0)
239  sxi( 1,4) = zero
240  seta( 1,4) = half
241  szeta(1,4) = zero
242  sxi( 2,4) = zero
243  seta( 2,4) = zero
244  szeta(2,4) = half
245  sxi( 3,4) = zero
246  seta( 3,4) = half
247  szeta(3,4) = half
248 
249 !-- surface normal vectors
250 
251  surf_elem(1,1)=1
252  surf_elem(2,1)=3
253  surf_elem(3,1)=2
254 
255  surf_elem(1,2)=1
256  surf_elem(2,2)=2
257  surf_elem(3,2)=4
258 
259  surf_elem(1,3)=3
260  surf_elem(2,3)=4
261  surf_elem(3,3)=2
262 
263  surf_elem(1,4)=1
264  surf_elem(2,4)=4
265  surf_elem(3,4)=3
266 
267 !-- initialize R_in, R_damping
268 
269  DO i = 1, numnp*3
270  r_in(i) = 0.d0
271  r_dp(i) = 0.d0
272  ENDDO
273 
274 ! loop for elements
275 
276  DO ielem = nstart, nend
277 
278  j = matlstet(ielem)
279  mat(1)= rho(j)
280  mat(2)= e(j)
281  mat(3)= xnu(j)
282 
283 !-- obtain mesh coordinate, mesh velocity, mesh accel.,
284 !-- displacement, velocity
285 
286  DO i=1, nnode
287  nx = lmlstet(i,ielem)*3-2
288  ny = lmlstet(i,ielem)*3-1
289  nz = lmlstet(i,ielem)*3
290  meshpos(1,i)=meshcoor(1,lmlstet(i,ielem))
291  meshpos(2,i)=meshcoor(2,lmlstet(i,ielem))
292  meshpos(3,i)=meshcoor(3,lmlstet(i,ielem))
293  meshvel(1,i)=v_bar(nx)
294  meshvel(2,i)=v_bar(ny)
295  meshvel(3,i)=v_bar(nz)
296  meshacc(1,i)=a_bar(nx)
297  meshacc(2,i)=a_bar(ny)
298  meshacc(3,i)=a_bar(nz)
299  disp(i*3-2) =d(nx)
300  disp(i*3-1) =d(ny)
301  disp(i*3 ) =d(nz)
302  velo(i*3-2) =vhalf(nx)
303  velo(i*3-1) =vhalf(ny)
304  velo(i*3 ) =vhalf(nz)
305  END DO
306 
307 
308 !-- evaluate shape function and its derivative at integration points
309 
310  CALL shcalc_3d10(shfn, shdx, dv, dvsum, shdx_av, meshpos, &
311  ndim,nnode,nintk_v,vxi,veta,vzeta,shdxi,ajacin)
312 
313 
314  DO igauss = 1, nintk_v
315 
316 
317 ! fill in element matrices for contributions from
318 ! volume integration points.
319 
320 ! evaluate mesh-related quantitites at vol. int. pt.
321  DO j = 1,ndim
322  mvgss(j) = zero
323  magss(j) = zero
324  DO knode = 1,nnode
325  mvgss(j) = mvgss(j) + meshvel(j,knode)* &
326  shfn(knode,igauss)
327  magss(j) = magss(j) + meshacc(j,knode)* &
328  shfn(knode,igauss)
329  END DO
330  END DO
331 
332  mcgss(1,1) = zero
333  mcgss(1,2) = zero
334  mcgss(1,3) = zero
335  mcgss(2,1) = zero
336  mcgss(2,2) = zero
337  mcgss(2,3) = zero
338  mcgss(3,1) = zero
339  mcgss(3,2) = zero
340  mcgss(3,3) = zero
341 
342  DO knode = 1,nnode
343  mcgss(1,1) = mcgss(1,1) + meshvel(1,knode)* &
344  shdx(1,knode,igauss)
345  mcgss(2,1) = mcgss(2,1) + meshvel(2,knode)* &
346  shdx(1,knode,igauss)
347  mcgss(3,1) = mcgss(3,1) + meshvel(3,knode)* &
348  shdx(1,knode,igauss)
349  mcgss(1,2) = mcgss(1,2) + meshvel(1,knode)* &
350  shdx(2,knode,igauss)
351  mcgss(2,2) = mcgss(2,2) + meshvel(2,knode)* &
352  shdx(2,knode,igauss)
353  mcgss(3,2) = mcgss(3,2) + meshvel(3,knode)* &
354  shdx(2,knode,igauss)
355  mcgss(1,3) = mcgss(1,3) + meshvel(1,knode)* &
356  shdx(3,knode,igauss)
357  mcgss(2,3) = mcgss(2,3) + meshvel(2,knode)* &
358  shdx(3,knode,igauss)
359  mcgss(3,3) = mcgss(3,3) + meshvel(3,knode)* &
360  shdx(3,knode,igauss)
361  END DO
362 
363  trmcgss = mcgss(1,1) + mcgss(2,2) + mcgss(3,3)
364 
365  wsparr(1) = mcgss(1,1)*mvgss(1) + &
366  mcgss(1,2)*mvgss(2) + &
367  mcgss(1,3)*mvgss(3)
368  wsparr(2) = mcgss(2,1)*mvgss(1) + &
369  mcgss(2,2)*mvgss(2) + &
370  mcgss(2,3)*mvgss(3)
371  wsparr(3) = mcgss(3,1)*mvgss(1) + &
372  mcgss(3,2)*mvgss(2) + &
373  mcgss(3,3)*mvgss(3)
374 
375 
376  DO knodc = 1,nnode
377  sto_wsp1(knodc) = &
378  mvgss(1)*shdx(1,knodc,igauss) + &
379  mvgss(2)*shdx(2,knodc,igauss) + &
380  mvgss(3)*shdx(3,knodc,igauss)
381  sto_wsp2(knodc) = &
382  wsparr(1)*shdx(1,knodc,igauss) + &
383  wsparr(2)*shdx(2,knodc,igauss) + &
384  wsparr(3)*shdx(3,knodc,igauss)
385  sto_wsp3(knodc) = &
386  mvgss(1)*shdx(1,knodc,igauss) + &
387  mvgss(2)*shdx(2,knodc,igauss) + &
388  mvgss(3)*shdx(3,knodc,igauss)
389  sto_wsp4(knodc) = &
390  magss(1)*shdx(1,knodc,igauss) + &
391  magss(2)*shdx(2,knodc,igauss) + &
392  magss(3)*shdx(3,knodc,igauss)
393  END DO
394 
395  factor=mat(1)*dv(igauss)*multiplier*wt_lstet
396 
397  DO knodr = 1,nnode
398  grindx = (lmlstet(knodr,ielem)-1)*3 + 1
399  grindy = (lmlstet(knodr,ielem)-1)*3 + 2
400  grindz = (lmlstet(knodr,ielem)-1)*3 + 3
401 
402 ! damping
403 ! ----------------------------------
404 
405  r_dp(grindx) = r_dp(grindx) - &
406  factor*shfn(knodr,igauss)*two* &
407  ( sto_wsp1(1)*velo((1-1)*ndim+1) + &
408  sto_wsp1(2)*velo((2-1)*ndim+1) + &
409  sto_wsp1(3)*velo((3-1)*ndim+1) + &
410  sto_wsp1(4)*velo((4-1)*ndim+1) + &
411  sto_wsp1(5)*velo((5-1)*ndim+1) + &
412  sto_wsp1(6)*velo((6-1)*ndim+1) + &
413  sto_wsp1(7)*velo((7-1)*ndim+1) + &
414  sto_wsp1(8)*velo((8-1)*ndim+1) + &
415  sto_wsp1(9)*velo((9-1)*ndim+1) + &
416  sto_wsp1(10)*velo((10-1)*ndim+1) )
417 
418  r_dp(grindy) = r_dp(grindy) - &
419  factor*shfn(knodr,igauss)*two* &
420  ( sto_wsp1(1)*velo((1-1)*ndim+2) + &
421  sto_wsp1(2)*velo((2-1)*ndim+2) + &
422  sto_wsp1(3)*velo((3-1)*ndim+2) + &
423  sto_wsp1(4)*velo((4-1)*ndim+2) + &
424  sto_wsp1(5)*velo((5-1)*ndim+2) + &
425  sto_wsp1(6)*velo((6-1)*ndim+2) + &
426  sto_wsp1(7)*velo((7-1)*ndim+2) + &
427  sto_wsp1(8)*velo((8-1)*ndim+2) + &
428  sto_wsp1(9)*velo((9-1)*ndim+2) + &
429  sto_wsp1(10)*velo((10-1)*ndim+2) )
430 
431  r_dp(grindz) = r_dp(grindz) - &
432  factor*shfn(knodr,igauss)*two* &
433  ( sto_wsp1(1)*velo((1-1)*ndim+3) + &
434  sto_wsp1(2)*velo((2-1)*ndim+3) + &
435  sto_wsp1(3)*velo((3-1)*ndim+3) + &
436  sto_wsp1(4)*velo((4-1)*ndim+3) + &
437  sto_wsp1(5)*velo((5-1)*ndim+3) + &
438  sto_wsp1(6)*velo((6-1)*ndim+3) + &
439  sto_wsp1(7)*velo((7-1)*ndim+3) + &
440  sto_wsp1(8)*velo((8-1)*ndim+3) + &
441  sto_wsp1(9)*velo((9-1)*ndim+3) + &
442  sto_wsp1(10)*velo((10-1)*ndim+3) )
443 
444 ! stiffness (vol. integration terms)
445 ! ----------------------------------
446 
447  r_in(grindx) = r_in(grindx) &
448  + factor*shfn(knodr,igauss) * ( &
449  -trmcgss*sto_wsp1(1)*disp((1-1)*ndim+1) &
450  -trmcgss*sto_wsp1(2)*disp((2-1)*ndim+1) &
451  -trmcgss*sto_wsp1(3)*disp((3-1)*ndim+1) &
452  -trmcgss*sto_wsp1(4)*disp((4-1)*ndim+1) &
453  -trmcgss*sto_wsp1(5)*disp((5-1)*ndim+1) &
454  -trmcgss*sto_wsp1(6)*disp((6-1)*ndim+1) &
455  -trmcgss*sto_wsp1(7)*disp((7-1)*ndim+1) &
456  -trmcgss*sto_wsp1(8)*disp((8-1)*ndim+1) &
457  -trmcgss*sto_wsp1(9)*disp((9-1)*ndim+1) &
458  -trmcgss*sto_wsp1(10)*disp((10-1)*ndim+1) &
459  + sto_wsp2(1)*disp((1-1)*ndim+1) &
460  + sto_wsp2(2)*disp((2-1)*ndim+1) &
461  + sto_wsp2(3)*disp((3-1)*ndim+1) &
462  + sto_wsp2(4)*disp((4-1)*ndim+1) &
463  + sto_wsp2(5)*disp((5-1)*ndim+1) &
464  + sto_wsp2(6)*disp((6-1)*ndim+1) &
465  + sto_wsp2(7)*disp((7-1)*ndim+1) &
466  + sto_wsp2(8)*disp((8-1)*ndim+1) &
467  + sto_wsp2(9)*disp((9-1)*ndim+1) &
468  + sto_wsp2(10)*disp((10-1)*ndim+1) &
469  - sto_wsp4(1)*disp((1-1)*ndim+1) &
470  - sto_wsp4(2)*disp((2-1)*ndim+1) &
471  - sto_wsp4(3)*disp((3-1)*ndim+1) &
472  - sto_wsp4(4)*disp((4-1)*ndim+1) &
473  - sto_wsp4(5)*disp((5-1)*ndim+1) &
474  - sto_wsp4(6)*disp((6-1)*ndim+1) &
475  - sto_wsp4(7)*disp((7-1)*ndim+1) &
476  - sto_wsp4(8)*disp((8-1)*ndim+1) &
477  - sto_wsp4(9)*disp((9-1)*ndim+1) &
478  - sto_wsp4(10)*disp((10-1)*ndim+1) ) &
479  + factor*sto_wsp3(knodr) * ( &
480  - sto_wsp1(1)*disp((1-1)*ndim+1) &
481  - sto_wsp1(2)*disp((2-1)*ndim+1) &
482  - sto_wsp1(3)*disp((3-1)*ndim+1) &
483  - sto_wsp1(4)*disp((4-1)*ndim+1) &
484  - sto_wsp1(5)*disp((5-1)*ndim+1) &
485  - sto_wsp1(6)*disp((6-1)*ndim+1) &
486  - sto_wsp1(7)*disp((7-1)*ndim+1) &
487  - sto_wsp1(8)*disp((8-1)*ndim+1) &
488  - sto_wsp1(9)*disp((9-1)*ndim+1) &
489  - sto_wsp1(10)*disp((10-1)*ndim+1) )
490 
491 
492  r_in(grindy) = r_in(grindy) &
493  + factor*shfn(knodr,igauss) * ( &
494  -trmcgss*sto_wsp1(1)*disp((1-1)*ndim+2) &
495  -trmcgss*sto_wsp1(2)*disp((2-1)*ndim+2) &
496  -trmcgss*sto_wsp1(3)*disp((3-1)*ndim+2) &
497  -trmcgss*sto_wsp1(4)*disp((4-1)*ndim+2) &
498  -trmcgss*sto_wsp1(5)*disp((5-1)*ndim+2) &
499  -trmcgss*sto_wsp1(6)*disp((6-1)*ndim+2) &
500  -trmcgss*sto_wsp1(7)*disp((7-1)*ndim+2) &
501  -trmcgss*sto_wsp1(8)*disp((8-1)*ndim+2) &
502  -trmcgss*sto_wsp1(9)*disp((9-1)*ndim+2) &
503  -trmcgss*sto_wsp1(10)*disp((10-1)*ndim+2) &
504  + sto_wsp2(1)*disp((1-1)*ndim+2) &
505  + sto_wsp2(2)*disp((2-1)*ndim+2) &
506  + sto_wsp2(3)*disp((3-1)*ndim+2) &
507  + sto_wsp2(4)*disp((4-1)*ndim+2) &
508  + sto_wsp2(5)*disp((5-1)*ndim+2) &
509  + sto_wsp2(6)*disp((6-1)*ndim+2) &
510  + sto_wsp2(7)*disp((7-1)*ndim+2) &
511  + sto_wsp2(8)*disp((8-1)*ndim+2) &
512  + sto_wsp2(9)*disp((9-1)*ndim+2) &
513  + sto_wsp2(10)*disp((10-1)*ndim+2) &
514  - sto_wsp4(1)*disp((1-1)*ndim+2) &
515  - sto_wsp4(2)*disp((2-1)*ndim+2) &
516  - sto_wsp4(3)*disp((3-1)*ndim+2) &
517  - sto_wsp4(4)*disp((4-1)*ndim+2) &
518  - sto_wsp4(5)*disp((5-1)*ndim+2) &
519  - sto_wsp4(6)*disp((6-1)*ndim+2) &
520  - sto_wsp4(7)*disp((7-1)*ndim+2) &
521  - sto_wsp4(8)*disp((8-1)*ndim+2) &
522  - sto_wsp4(9)*disp((9-1)*ndim+2) &
523  - sto_wsp4(10)*disp((10-1)*ndim+2) ) &
524  + factor*sto_wsp3(knodr) * ( &
525  - sto_wsp1(1)*disp((1-1)*ndim+2) &
526  - sto_wsp1(2)*disp((2-1)*ndim+2) &
527  - sto_wsp1(3)*disp((3-1)*ndim+2) &
528  - sto_wsp1(4)*disp((4-1)*ndim+2) &
529  - sto_wsp1(5)*disp((5-1)*ndim+2) &
530  - sto_wsp1(6)*disp((6-1)*ndim+2) &
531  - sto_wsp1(7)*disp((7-1)*ndim+2) &
532  - sto_wsp1(8)*disp((8-1)*ndim+2) &
533  - sto_wsp1(9)*disp((9-1)*ndim+2) &
534  - sto_wsp1(10)*disp((10-1)*ndim+2) )
535 
536  r_in(grindz) = r_in(grindz) &
537  + factor*shfn(knodr,igauss) * ( &
538  -trmcgss*sto_wsp1(1)*disp((1-1)*ndim+3) &
539  -trmcgss*sto_wsp1(2)*disp((2-1)*ndim+3) &
540  -trmcgss*sto_wsp1(3)*disp((3-1)*ndim+3) &
541  -trmcgss*sto_wsp1(4)*disp((4-1)*ndim+3) &
542  -trmcgss*sto_wsp1(5)*disp((5-1)*ndim+3) &
543  -trmcgss*sto_wsp1(6)*disp((6-1)*ndim+3) &
544  -trmcgss*sto_wsp1(7)*disp((7-1)*ndim+3) &
545  -trmcgss*sto_wsp1(8)*disp((8-1)*ndim+3) &
546  -trmcgss*sto_wsp1(9)*disp((9-1)*ndim+3) &
547  -trmcgss*sto_wsp1(10)*disp((10-1)*ndim+3) &
548  + sto_wsp2(1)*disp((1-1)*ndim+3) &
549  + sto_wsp2(2)*disp((2-1)*ndim+3) &
550  + sto_wsp2(3)*disp((3-1)*ndim+3) &
551  + sto_wsp2(4)*disp((4-1)*ndim+3) &
552  + sto_wsp2(5)*disp((5-1)*ndim+3) &
553  + sto_wsp2(6)*disp((6-1)*ndim+3) &
554  + sto_wsp2(7)*disp((7-1)*ndim+3) &
555  + sto_wsp2(8)*disp((8-1)*ndim+3) &
556  + sto_wsp2(9)*disp((9-1)*ndim+3) &
557  + sto_wsp2(10)*disp((10-1)*ndim+3) &
558  - sto_wsp4(1)*disp((1-1)*ndim+3) &
559  - sto_wsp4(2)*disp((2-1)*ndim+3) &
560  - sto_wsp4(3)*disp((3-1)*ndim+3) &
561  - sto_wsp4(4)*disp((4-1)*ndim+3) &
562  - sto_wsp4(5)*disp((5-1)*ndim+3) &
563  - sto_wsp4(6)*disp((6-1)*ndim+3) &
564  - sto_wsp4(7)*disp((7-1)*ndim+3) &
565  - sto_wsp4(8)*disp((8-1)*ndim+3) &
566  - sto_wsp4(9)*disp((9-1)*ndim+3) &
567  - sto_wsp4(10)*disp((10-1)*ndim+3) ) &
568  + factor*sto_wsp3(knodr) * ( &
569  - sto_wsp1(1)*disp((1-1)*ndim+3) &
570  - sto_wsp1(2)*disp((2-1)*ndim+3) &
571  - sto_wsp1(3)*disp((3-1)*ndim+3) &
572  - sto_wsp1(4)*disp((4-1)*ndim+3) &
573  - sto_wsp1(5)*disp((5-1)*ndim+3) &
574  - sto_wsp1(6)*disp((6-1)*ndim+3) &
575  - sto_wsp1(7)*disp((7-1)*ndim+3) &
576  - sto_wsp1(8)*disp((8-1)*ndim+3) &
577  - sto_wsp1(9)*disp((9-1)*ndim+3) &
578  - sto_wsp1(10)*disp((10-1)*ndim+3) )
579 
580 ! --------------------------------------------
581 ! all stiffness terms above are from physical
582 ! acceleration. The isotropic linear
583 ! elastic stiffness is calculated in CST
584 ! --------------------------------------------
585 
586  END DO
587 
588  ENDDO ! for igauss=1,nintk_v
589 
590 
591 ! contribution to stiffness matrix from surface terms
592 ! ---------------------------------------------------
593 
594  DO iface = 1,nsurf
595 
596 !-- shdx is constant and already calculated !!
597 
598  CALL shcalc_3d10(shfn,shdx,dv,dvsum,shdx_av,meshpos,ndim, &
599  nnode,nintk_s,sxi(1,iface),seta(1,iface), &
600  szeta(1,iface),shdxi,ajacin)
601 
602  DO igauss = 1,nintk_s
603 
604  DO j = 1,ndim
605  mvgss(j) = zero
606  DO knode = 1,nnode
607  mvgss(j) = mvgss(j) + meshvel(j,knode)* &
608  shfn(knode,igauss)
609  END DO
610  END DO
611 
612 !-- unit out-normal vector on each surface
613 
614  DO j = 1, 3
615  wsparr(j) =meshpos(j,surf_elem(2,iface)) - &
616  meshpos(j,surf_elem(1,iface))
617  wsparr1(j)=meshpos(j,surf_elem(3,iface)) - &
618  meshpos(j,surf_elem(1,iface))
619  END DO
620 
621 ! do cross product: wsparr = nk * da
622 
623  wsp1 = wsparr(2)*wsparr1(3) - wsparr(3)*wsparr1(2)
624  wsp2 = wsparr(3)*wsparr1(1) - wsparr(1)*wsparr1(3)
625  wsp3 = wsparr(1)*wsparr1(2) - wsparr(2)*wsparr1(1)
626  wsparr(1) = wsp1 * 0.5d0
627  wsparr(2) = wsp2 * 0.5d0
628  wsparr(3) = wsp3 * 0.5d0
629 
630 
631 ! evaluate terms required in stiffness
632 
633  wsp1 = mvgss(1)*wsparr(1) &
634  + mvgss(2)*wsparr(2) &
635  + mvgss(3)*wsparr(3)
636 
637 ! fill stiffness matrix
638 
639 
640  DO knodc = 1,nnode
641  sto_wsp2(knodc) = &
642  mvgss(1)*shdx(1,knodc,igauss) + &
643  mvgss(2)*shdx(2,knodc,igauss) + &
644  mvgss(3)*shdx(3,knodc,igauss)
645  END DO
646 
647  factor = third
648  DO knodr = 1, nnode
649  grindx = (lmlstet(knodr,ielem)-1)*3 + 1
650  grindy = (lmlstet(knodr,ielem)-1)*3 + 2
651  grindz = (lmlstet(knodr,ielem)-1)*3 + 3
652 
653  r_in(grindx) = r_in(grindx) &
654  + mat(1)*shfn(knodr,igauss)*wsp1*( &
655  + sto_wsp2(1)*disp((1-1)*ndim+1) &
656  + sto_wsp2(2)*disp((2-1)*ndim+1) &
657  + sto_wsp2(3)*disp((3-1)*ndim+1) &
658  + sto_wsp2(4)*disp((4-1)*ndim+1) &
659  + sto_wsp2(5)*disp((5-1)*ndim+1) &
660  + sto_wsp2(6)*disp((6-1)*ndim+1) &
661  + sto_wsp2(7)*disp((7-1)*ndim+1) &
662  + sto_wsp2(8)*disp((8-1)*ndim+1) &
663  + sto_wsp2(9)*disp((9-1)*ndim+1) &
664  + sto_wsp2(10)*disp((10-1)*ndim+1) ) &
665  * factor
666 
667  r_in(grindy) = r_in(grindy) &
668  + mat(1)*shfn(knodr,igauss)*wsp1*( &
669  + sto_wsp2(1)*disp((1-1)*ndim+2) &
670  + sto_wsp2(2)*disp((2-1)*ndim+2) &
671  + sto_wsp2(3)*disp((3-1)*ndim+2) &
672  + sto_wsp2(4)*disp((4-1)*ndim+2) &
673  + sto_wsp2(5)*disp((5-1)*ndim+2) &
674  + sto_wsp2(6)*disp((6-1)*ndim+2) &
675  + sto_wsp2(7)*disp((7-1)*ndim+2) &
676  + sto_wsp2(8)*disp((8-1)*ndim+2) &
677  + sto_wsp2(9)*disp((9-1)*ndim+2) &
678  + sto_wsp2(10)*disp((10-1)*ndim+2) ) &
679  * factor
680 
681  r_in(grindz) = r_in(grindz) &
682  + mat(1)*shfn(knodr,igauss)*wsp1*( &
683  + sto_wsp2(1)*disp((1-1)*ndim+3) &
684  + sto_wsp2(2)*disp((2-1)*ndim+3) &
685  + sto_wsp2(3)*disp((3-1)*ndim+3) &
686  + sto_wsp2(4)*disp((4-1)*ndim+3) &
687  + sto_wsp2(5)*disp((5-1)*ndim+3) &
688  + sto_wsp2(6)*disp((6-1)*ndim+3) &
689  + sto_wsp2(7)*disp((7-1)*ndim+3) &
690  + sto_wsp2(8)*disp((8-1)*ndim+3) &
691  + sto_wsp2(9)*disp((9-1)*ndim+3) &
692  + sto_wsp2(10)*disp((10-1)*ndim+3) ) &
693  * factor
694  END DO
695 
696  ENDDO ! igauss = 1, nintk_s
697  ENDDO ! iface=1,nsurf
698 
699  ENDDO ! for ielem=1,numlstet
700 
701  DO i=1,numnp*3
702  rnet(i) = rnet(i) - r_in(i) - r_dp(i)
703  ENDDO
704 
705 111 FORMAT(i5,5e15.7)
706 
707  RETURN
708 END SUBROUTINE v3d10_ale
709 
cimg_library::CImg::sum
Tfloat sum() const
Return the sum of all the pixel values in an image.
Definition: CImg.h:13022
Mesquite::Matrix3D::zero
void zero()
Sets all entries to zero (more efficient than assignement).
Definition: includeLinks/Matrix3D.hpp:167
d
const NT & d
Definition: rational_rotation.h:73
k
j indices k indices k
Definition: Indexing.h:6
shcalc_3d10
subroutine shcalc_3d10(shfn, shdx, dv, dvsum, shdx_av, meshpos, ndim, nnode, nintk, xi, eta, zeta, shdxi, wsp2)
Definition: shcalc_3d10.f90:65
i
blockLoc i
Definition: read.cpp:79
Direction_2::delta
RT delta(int i) const
Definition: Direction_2.h:159
j
j indices j
Definition: Indexing.h:6
v3d10_ale
subroutine v3d10_ale(v_bar, a_bar, d, vhalf, Rnet, E, xnu, rho, numnp, numat_vol, numlstet, matlstet, lmlstet, meshcoor, nstart, nend)
Definition: v3d10_ale.f90:65