Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
implicit_v3d8_me_K.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 !!****
55 !!
56 !! NAME
57 !! implicit_v3d8_me_K
58 !!
59 !! FUNCTION
60 !! This subroutine constructs the stiffness matrix associated
61 !! with all of the elements that are in this processors .inp file.
62 !!
63 !! INPUTS
64 !! coor -- Coordinates of the nodes
65 !! ElConnVol -- The connectivity table
66 !! dmat -- Material stiffness matrix
67 !! numnp -- Number of nodes
68 !! nstart -- First element to compute the stiffness matrix for
69 !! nend -- Last element to compute the stiffness matrix for
70 !! NumEl -- Number of elements
71 !! MatType -- Mapping of materials to elements
72 !! NumMatType -- Total number of materials
73 !! enhanced_map -- The enhanced map
74 !! mixed_map -- The mixed map
75 !!
76 !! OUTPUTS
77 !! none
78 !!
79 !! USES
80 !! none
81 !!
82 !!****
83 
84 SUBROUTINE implicit_v3d8_me_k(coor,ElConnVol,ci, &
85  numnp,nstart,nend,numel,mattype,nummattype,enhanced_map,mixed_map,nummatvol)
86 
87 
88  ! dmat is the 9x9 material matrix
89 
90  USE comp_row_global
91  USE implicit_global
92  USE precision
93 
94  IMPLICIT NONE
95 
96 !-----Global variables
97  INTEGER :: numnp ! number of nodes
98  INTEGER :: numat_vol ! number of volumetric materials
99  INTEGER :: numel ! number of LSTets
100  integer :: nummattype
101  INTEGER :: nummatvol
102 !-- coordinate array
103  REAL*8, DIMENSION(1:3,1:numnp) :: coor
104 !-- connectivity table for Brick
105  INTEGER, DIMENSION(1:8,1:NumEl) :: elconnvol
106  INTEGER, DIMENSION(1:NumEl) :: mattype
107 !-- elastic stiffness consts
108  REAL*8, DIMENSION(1:6,1:6,NumMatType) :: ci
109 !-- internal force
110  REAL*8, DIMENSION(1:3*numnp) :: r_in
111 !-- displacement vector
112  REAL*8, DIMENSION(1:3*numnp) :: disp
113 !---- Local variables
114 !-- node numbers
115  INTEGER :: n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
116 
117  INTEGER :: k1n1,k1n2,k1n3,k1n4,k1n5,k1n6,k1n7,k1n8
118  INTEGER :: k2n1,k2n2,k2n3,k2n4,k2n5,k2n6,k2n7,k2n8
119  INTEGER :: k3n1,k3n2,k3n3,k3n4,k3n5,k3n6,k3n7,k3n8
120 
121  REAL*8, DIMENSION(1:3,1:8) :: coord
122  !REAL*8, DIMENSION(1:24) :: Udisp
123 
124  REAL*8 :: element_volume
125  INTEGER :: ielem, imat
126  integer :: igpt
127  REAL*8 nnn(8), dn(8,3), jac(3,3), jacinv(3,3), &
128  dn20(20,3),nnn20(20),&
129  t(3,3), tinv(3,3), &
130  dmat(9,9), bmat(9,24), bmat2(1:8,1:9,1:24), &
131  grad(9,24), tmp38(3,8), coeff(8),&
132  tkront(9,9), tkrontinv(9,9),&
133  e_mixed(9,12), &
134  gsm(12,24), gbg(12,12), bmat_avg(9,24),&
135  bu(9,24), ba(9,9), bsm(12,24),& ! print edge
136  kuu(24,24), kau(9,24), kaa(9,9),atol, bavec(1:9,1:9,1:8), buvec(1:9,1:24,1:8),&
137  tmp249(24,9),tmp91(9,1),kaa_copy(9,9),enh2(9),dnorm,bmatigs(1:8,1:9,1:24),&
138  stresstmp(9,1), tmp9(9),dmatinf(5,9,9)
139 
140 ! weight for hex element
141 
142  REAL*8, DIMENSION(1:8) :: wi = &
143  (/1.000000000000000,1.000000000000000,1.000000000000000, &
144  1.000000000000000,1.000000000000000,1.000000000000000, &
145  1.000000000000000,1.000000000000000/)
146 
147 ! root3 = 1./sqrt(3.)
148 
149  REAL*8, DIMENSION(1:3,1:8) :: ri = reshape( &
150  (/-0.577350269189626,-0.577350269189626,-0.577350269189626, &
151  0.577350269189626,-0.577350269189626,-0.577350269189626, &
152  0.577350269189626, 0.577350269189626,-0.577350269189626, &
153  -0.577350269189626, 0.577350269189626,-0.577350269189626, &
154  -0.577350269189626,-0.577350269189626, 0.577350269189626, &
155  0.577350269189626,-0.577350269189626, 0.577350269189626, &
156  0.577350269189626, 0.577350269189626, 0.577350269189626, &
157  -0.577350269189626, 0.577350269189626, 0.577350269189626/),(/3,8/) )
158 
159  REAL*8, DIMENSION(1:9,1:9) :: sop = reshape( &
160  (/1.000000000000000,0.000000000000000,0.000000000000000, &
161  0.000000000000000,0.000000000000000,0.000000000000000, &
162  0.000000000000000,0.000000000000000,0.000000000000000, &
163  0.000000000000000,0.500000000000000,0.000000000000000, &
164  0.500000000000000,0.000000000000000,0.000000000000000, &
165  0.000000000000000,0.000000000000000,0.000000000000000, &
166  0.000000000000000,0.000000000000000,0.500000000000000, &
167  0.000000000000000,0.000000000000000,0.000000000000000, &
168  0.500000000000000,0.000000000000000,0.000000000000000, &
169  0.000000000000000,0.500000000000000,0.000000000000000, &
170  0.500000000000000,0.000000000000000,0.000000000000000, &
171  0.000000000000000,0.000000000000000,0.000000000000000, &
172  0.000000000000000,0.000000000000000,0.000000000000000, &
173  0.000000000000000,1.000000000000000,0.000000000000000, &
174  0.000000000000000,0.000000000000000,0.000000000000000, &
175  0.000000000000000,0.000000000000000,0.000000000000000, &
176  0.000000000000000,0.000000000000000,0.500000000000000, &
177  0.000000000000000,0.500000000000000,0.000000000000000, &
178  0.000000000000000,0.000000000000000,0.500000000000000, &
179  0.000000000000000,0.000000000000000,0.000000000000000, &
180  0.500000000000000,0.000000000000000,0.000000000000000, &
181  0.000000000000000,0.000000000000000,0.000000000000000, &
182  0.000000000000000,0.000000000000000,0.500000000000000, &
183  0.000000000000000,0.500000000000000,0.000000000000000, &
184  0.000000000000000,0.000000000000000,0.000000000000000, &
185  0.000000000000000,0.000000000000000,0.000000000000000, &
186  0.000000000000000,0.000000000000000,1.000000000000000/),(/9,9/) )
187 
188  REAL*8, DIMENSION(1:3,1:3) :: dident = reshape( &
189  (/1.000000000000000,0.000000000000000,0.000000000000000, &
190  0.000000000000000,1.000000000000000,0.000000000000000, &
191  0.000000000000000,0.000000000000000,1.000000000000000/),(/3,3/) )
192 
193  REAL*8 :: zero = 0.0, one = 1.0, detj,det, three = 3.
194  LOGICAL error, debug
195  REAL*8 :: alpha2, alpha
196  REAL*8, DIMENSION(1:9) :: fenh
197  REAL*8, DIMENSION(1:9,1:9) :: e_enhanced
198  REAL*8 :: maxdisp, con1, tp, sum
199  REAL*8 :: threeeighth = 3./8.
200  REAL*8 :: denh
201  INTEGER :: igauss
202  REAL*8, DIMENSION(8) :: xi, eta, zeta
203  REAL*8, DIMENSION(8) :: xie, etae, zetae
204  REAL(KIND=8), DIMENSION(1:8,1:9,1:12) :: mixed_map
205  REAL(KIND=8), DIMENSION(1:8,1:9,1:9) :: enhanced_map
206  REAL(KIND=8), DIMENSION(1:8,1:8) :: shapefun,shapefune
207  INTEGER :: otdev,mcrd,nnode
208  INTEGER,parameter :: ngpts = 8
209  REAL(KIND=8) :: strainenh(1:ngpts,1:9,1:numel)
210  !REAL*8 :: stress(1:ngpts,1:9,1:NumEl)
211  INTEGER :: i,k,j
212  INTEGER :: nstart, nend
213 
214  REAL(KIND=8), DIMENSION(1:24,1:24) :: amatrx
215 
216 
217 ! INTEGER :: GlNumNp
218 ! INTEGER, ALLOCATABLE, DIMENSION(:) :: rp_k, cval_k
219 ! REAL(KIND=wp), ALLOCATABLE, DIMENSION(:) :: aval_k
220  INTEGER :: nnz_k
221  INTEGER :: m
222  INTEGER :: dof1, dof2
223  INTEGER :: ldof1, ldof2
224  INTEGER :: gdof1, gdof2
225  INTEGER :: counter
226  REAL(kind=wp) :: tempval
227 
228 
229  ! Initialize the K matrix
230  CALL numknnz(numnp,numel,elconnvol,nnz_k)
231  ALLOCATE(rp_k(1:3*gnumnp+1))
232  ALLOCATE(cval_k(1:nnz_k))
233  ALLOCATE(aval_k(1:nnz_k))
234  CALL initk(numnp,numel,elconnvol,nnz_k,rp_k,cval_k,aval_k)
235  !print*,myid,' nnz_k = ',nnz_k
236 
237  DO ielem = nstart, nend ! Loop over elements
238 
239  imat = mattype(ielem)
240 
241  CALL map_ci2dmat(ci(1:6,1:6,imat),dmat)
242 
243  n1 = elconnvol(1,ielem)
244  n2 = elconnvol(2,ielem)
245  n3 = elconnvol(3,ielem)
246  n4 = elconnvol(4,ielem)
247  n5 = elconnvol(5,ielem)
248  n6 = elconnvol(6,ielem)
249  n7 = elconnvol(7,ielem)
250  n8 = elconnvol(8,ielem)
251 
252  k3n1 = 3*n1
253  k3n2 = 3*n2
254  k3n3 = 3*n3
255  k3n4 = 3*n4
256  k3n5 = 3*n5
257  k3n6 = 3*n6
258  k3n7 = 3*n7
259  k3n8 = 3*n8
260 
261  k2n1 = k3n1 - 1
262  k2n2 = k3n2 - 1
263  k2n3 = k3n3 - 1
264  k2n4 = k3n4 - 1
265  k2n5 = k3n5 - 1
266  k2n6 = k3n6 - 1
267  k2n7 = k3n7 - 1
268  k2n8 = k3n8 - 1
269 
270  k1n1 = k3n1 - 2
271  k1n2 = k3n2 - 2
272  k1n3 = k3n3 - 2
273  k1n4 = k3n4 - 2
274  k1n5 = k3n5 - 2
275  k1n6 = k3n6 - 2
276  k1n7 = k3n7 - 2
277  k1n8 = k3n8 - 2
278 
279  coord(1,1) = coor(1,n1)
280  coord(2,1) = coor(2,n1)
281  coord(3,1) = coor(3,n1)
282 
283  coord(1,2) = coor(1,n2)
284  coord(2,2) = coor(2,n2)
285  coord(3,2) = coor(3,n2)
286 
287  coord(1,3) = coor(1,n3)
288  coord(2,3) = coor(2,n3)
289  coord(3,3) = coor(3,n3)
290 
291  coord(1,4) = coor(1,n4)
292  coord(2,4) = coor(2,n4)
293  coord(3,4) = coor(3,n4)
294 
295  coord(1,5) = coor(1,n5)
296  coord(2,5) = coor(2,n5)
297  coord(3,5) = coor(3,n5)
298 
299  coord(1,6) = coor(1,n6)
300  coord(2,6) = coor(2,n6)
301  coord(3,6) = coor(3,n6)
302 
303  coord(1,7) = coor(1,n7)
304  coord(2,7) = coor(2,n7)
305  coord(3,7) = coor(3,n7)
306 
307  coord(1,8) = coor(1,n8)
308  coord(2,8) = coor(2,n8)
309  coord(3,8) = coor(3,n8)
310 
311 !!$ Udisp(1) = disp(k1n1)
312 !!$ Udisp(2) = disp(k2n1)
313 !!$ Udisp(3) = disp(k3n1)
314 !!$
315 !!$ Udisp(4) = disp(k1n2)
316 !!$ Udisp(5) = disp(k2n2)
317 !!$ Udisp(6) = disp(k3n2)
318 !!$
319 !!$ Udisp(7) = disp(k1n3)
320 !!$ Udisp(8) = disp(k2n3)
321 !!$ Udisp(9) = disp(k3n3)
322 !!$
323 !!$ Udisp(10) = disp(k1n4)
324 !!$ Udisp(11) = disp(k2n4)
325 !!$ Udisp(12) = disp(k3n4)
326 !!$
327 !!$ Udisp(13) = disp(k1n5)
328 !!$ Udisp(14) = disp(k2n5)
329 !!$ Udisp(15) = disp(k3n5)
330 !!$
331 !!$ Udisp(16) = disp(k1n6)
332 !!$ Udisp(17) = disp(k2n6)
333 !!$ Udisp(18) = disp(k3n6)
334 !!$
335 !!$ Udisp(19) = disp(k1n7)
336 !!$ Udisp(20) = disp(k2n7)
337 !!$ Udisp(21) = disp(k3n7)
338 !!$
339 !!$ Udisp(22) = disp(k1n8)
340 !!$ Udisp(23) = disp(k2n8)
341 !!$ Udisp(24) = disp(k3n8)
342 
343 
344 ! Compute average element quantities
345 
346  element_volume = 0. ! Initialize
347  t(1:3,1:3) = 0.
348  bmat_avg(1:9,1:24) = 0.
349 
350 !!$! Calculate the volume of each element. First calculate the Jacobian
351 !!$! at the "center" of the element.
352 !!$
353 !!$ jac(1,1) = eighth*( coord(1,1) - coord(1,2) - coord(1,3) + coord(1,4) + &
354 !!$ coord(1,5) - coord(1,6) - coord(1,7) + coord(1,8) )
355 !!$
356 !!$ jac(1,2) = eighth*( coord(2,1) - coord(2,2) - coord(2,3) + coord(2,4) + &
357 !!$ coord(2,5) - coord(2,6) - coord(2,7) + coord(2,8) )
358 !!$
359 !!$ jac(1,3) = eighth*( coord(3,1) - coord(3,2) - coord(3,3) + coord(3,4) + &
360 !!$ coord(3,5) - coord(3,6) - coord(3,7) + coord(3,8) )
361 !!$
362 !!$ jac(2,1) = eighth*( coord(1,1) + coord(1,2) - coord(1,3) - coord(1,4) + &
363 !!$ coord(1,5) + coord(1,6) - coord(1,7) - coord(1,8) )
364 !!$
365 !!$ jac(2,2) = eighth*( coord(2,1) + coord(2,2) - coord(2,3) - coord(2,4) + &
366 !!$ coord(2,5) + coord(2,6) - coord(2,7) - coord(2,8) )
367 !!$
368 !!$ jac(2,3) = eighth*( coord(3,1) + coord(3,2) - coord(3,3) - coord(3,4) + &
369 !!$ coord(3,5) + coord(3,6) - coord(3,7) - coord(3,8) )
370 !!$
371 !!$ jac(3,1) = eighth*( coord(1,1) + coord(1,2) + coord(1,3) + coord(1,4) - &
372 !!$ coord(1,5) - coord(1,6) - coord(1,7) - coord(1,8) )
373 !!$
374 !!$ jac(3,2) = eighth*( coord(2,1) + coord(2,2) + coord(2,3) + coord(2,4) - &
375 !!$ coord(2,5) - coord(2,6) - coord(2,7) - coord(2,8) )
376 !!$
377 !!$ jac(3,3) = eighth*( coord(3,1) + coord(3,2) + coord(3,3) + coord(3,4) - &
378 !!$ coord(3,5) - coord(3,6) - coord(3,7) - coord(3,8) )
379 !!$
380 !!$
381 !!$ temp_det = jac(1,1)*( jac(2,2)*jac(3,3) - jac(2,3)*jac(3,2) ) - &
382 !!$ jac(1,2)*( jac(2,1)*jac(3,3) - jac(2,3)*jac(3,1) ) + &
383 !!$ jac(1,3)*( jac(2,1)*jac(3,2) - jac(2,2)*jac(3,1) )
384 !!$
385 !!$ volume = eight*temp_det
386  !
387 ! 1st LOOP OVER GUAUSS POINTS
388 
389  DO igpt = 1, 8 ! LOOP over gauss points
390 
391 !!$! Calculate the derivatives of the shape functions at the gauss points
392 !!$
393 !!$ dn(1,1) = eighth*(one + eta(igauss))*(one + zeta(igauss))
394 !!$ dn(1,2) = eighth*(one + xi(igauss) )*(one + zeta(igauss))
395 !!$ dn(1,3) = eighth*(one + xi(igauss) )*(one + eta(igauss) )
396 !!$
397 !!$ dn(2,1) = -eighth*(one + eta(igauss))*(one + zeta(igauss))
398 !!$ dn(2,2) = eighth*(one - xi(igauss) )*(one + zeta(igauss))
399 !!$ dn(2,3) = eighth*(one - xi(igauss) )*(one + eta(igauss) )
400 !!$
401 !!$ dn(3,1) = -eighth*(one - eta(igauss))*(one + zeta(igauss))
402 !!$ dn(3,2) = -eighth*(one - xi(igauss) )*(one + zeta(igauss))
403 !!$ dn(3,3) = eighth*(one - xi(igauss) )*(one - eta(igauss) )
404 !!$
405 !!$ dn(4,1) = eighth*(one - eta(igauss))*(one + zeta(igauss))
406 !!$ dn(4,2) = -eighth*(one + xi(igauss) )*(one + zeta(igauss))
407 !!$ dn(4,3) = eighth*(one + xi(igauss) )*(one - eta(igauss) )
408 !!$
409 !!$ dn(5,1) = eighth*(one + eta(igauss))*(one - zeta(igauss))
410 !!$ dn(5,2) = eighth*(one + xi(igauss) )*(one - zeta(igauss))
411 !!$ dn(5,3) = -eighth*(one + xi(igauss) )*(one + eta(igauss) )
412 !!$
413 !!$ dn(6,1) = -eighth*(one + eta(igauss))*(one - zeta(igauss))
414 !!$ dn(6,2) = eighth*(one - xi(igauss) )*(one - zeta(igauss))
415 !!$ dn(6,3) = -eighth*(one - xi(igauss) )*(one + eta(igauss) )
416 !!$
417 !!$ dn(7,1) = -eighth*(one - eta(igauss))*(one - zeta(igauss))
418 !!$ dn(7,2) = -eighth*(one - xi(igauss) )*(one - zeta(igauss))
419 !!$ dn(7,3) = -eighth*(one - xi(igauss) )*(one - eta(igauss) )
420 !!$
421 !!$ dn(8,1) = eighth*(one - eta(igauss))*(one - zeta(igauss))
422 !!$ dn(8,2) = -eighth*(one + xi(igauss) )*(one - zeta(igauss))
423 !!$ dn(8,3) = -eighth*(one + xi(igauss) )*(one - eta(igauss) )
424 !!$
425 !!$! Calculate the Jacobian, its determinant and inverse for each gauss point.
426 !!$
427 !!$ jac(1,1) = dn(1,1)*coord(1,1) + dn(2,1)*coord(1,2) &
428 !!$ + dn(3,1)*coord(1,3) + dn(4,1)*coord(1,4) &
429 !!$ + dn(5,1)*coord(1,5) + dn(6,1)*coord(1,6) &
430 !!$ + dn(7,1)*coord(1,7) + dn(8,1)*coord(1,8)
431 !!$
432 !!$ jac(2,1) = dn(1,1)*coord(2,1) + dn(2,1)*coord(2,2) &
433 !!$ + dn(3,1)*coord(2,3) + dn(4,1)*coord(2,4) &
434 !!$ + dn(5,1)*coord(2,5) + dn(6,1)*coord(2,6) &
435 !!$ + dn(7,1)*coord(2,7) + dn(8,1)*coord(2,8)
436 !!$
437 !!$ jac(3,1) = dn(1,1)*coord(3,1) + dn(2,1)*coord(3,2) &
438 !!$ + dn(3,1)*coord(3,3) + dn(4,1)*coord(3,4) &
439 !!$ + dn(5,1)*coord(3,5) + dn(6,1)*coord(3,6) &
440 !!$ + dn(7,1)*coord(3,7) + dn(8,1)*coord(3,8)
441 !!$
442 !!$ jac(1,2) = dn(1,2)*coord(1,1) + dn(2,2)*coord(1,2) &
443 !!$ + dn(3,2)*coord(1,3) + dn(4,2)*coord(1,4) &
444 !!$ + dn(5,2)*coord(1,5) + dn(6,2)*coord(1,6) &
445 !!$ + dn(7,2)*coord(1,7) + dn(8,2)*coord(1,8)
446 !!$
447 !!$ jac(2,2) = dn(1,2)*coord(2,1) + dn(2,2)*coord(2,2) &
448 !!$ + dn(3,2)*coord(2,3) + dn(4,2)*coord(2,4) &
449 !!$ + dn(5,2)*coord(2,5) + dn(6,2)*coord(2,6) &
450 !!$ + dn(7,2)*coord(2,7) + dn(8,2)*coord(2,8)
451 !!$
452 !!$ jac(3,2) = dn(1,2)*coord(3,1) + dn(2,2)*coord(3,2) &
453 !!$ + dn(3,2)*coord(3,3) + dn(4,2)*coord(3,4) &
454 !!$ + dn(5,2)*coord(3,5) + dn(6,2)*coord(3,6) &
455 !!$ + dn(7,2)*coord(3,7) + dn(8,2)*coord(3,8)
456 !!$
457 !!$ jac(1,3) = dn(1,3)*coord(1,1) + dn(2,3)*coord(1,2) &
458 !!$ + dn(3,3)*coord(1,3) + dn(4,3)*coord(1,4) &
459 !!$ + dn(5,3)*coord(1,5) + dn(6,3)*coord(1,6) &
460 !!$ + dn(7,3)*coord(1,7) + dn(8,3)*coord(1,8)
461 !!$
462 !!$ jac(2,3) = dn(1,3)*coord(2,1) + dn(2,3)*coord(2,2) &
463 !!$ + dn(3,3)*coord(2,3) + dn(4,3)*coord(2,4) &
464 !!$ + dn(5,3)*coord(2,5) + dn(6,3)*coord(2,6) &
465 !!$ + dn(7,3)*coord(2,7) + dn(8,3)*coord(2,8)
466 !!$
467 !!$ jac(3,3) = dn(1,3)*coord(3,1) + dn(2,3)*coord(3,2) &
468 !!$ + dn(3,3)*coord(3,3) + dn(4,3)*coord(3,4) &
469 !!$ + dn(5,3)*coord(3,5) + dn(6,3)*coord(3,6) &
470 !!$ + dn(7,3)*coord(3,7) + dn(8,3)*coord(3,8)
471 !!$
472 !!$ det(ielem,igauss) = jac(1,1)*( jac(2,2)*jac(3,3) - jac(2,3)*jac(3,2) ) - &
473 !!$ jac(1,2)*( jac(2,1)*jac(3,3) - jac(2,3)*jac(3,1) ) + &
474 !!$ jac(1,3)*( jac(2,1)*jac(3,2) - jac(2,2)*jac(3,1) )
475 !!$
476 !!$ jacinv(1,1) = ( jac(2,2)*jac(3,3) - jac(2,3)*jac(3,2) ) / det(ielem,igauss)
477 !!$ jacinv(1,2) = -( jac(1,2)*jac(3,3) - jac(1,3)*jac(3,2) ) / det(ielem,igauss)
478 !!$ jacinv(1,3) = ( jac(1,2)*jac(2,3) - jac(1,3)*jac(2,2) ) / det(ielem,igauss)
479 !!$ jacinv(2,1) = -( jac(2,1)*jac(3,3) - jac(2,3)*jac(3,1) ) / det(ielem,igauss)
480 !!$ jacinv(2,2) = ( jac(1,1)*jac(3,3) - jac(1,3)*jac(3,1) ) / det(ielem,igauss)
481 !!$ jacinv(2,3) = -( jac(1,1)*jac(2,3) - jac(1,3)*jac(2,1) ) / det(ielem,igauss)
482 !!$ jacinv(3,1) = ( jac(2,1)*jac(3,2) - jac(2,2)*jac(3,1) ) / det(ielem,igauss)
483 !!$ jacinv(3,2) = -( jac(1,1)*jac(3,2) - jac(1,2)*jac(3,1) ) / det(ielem,igauss)
484 !!$ jacinv(3,3) = ( jac(1,1)*jac(2,2) - jac(1,2)*jac(2,1) ) / det(ielem,igauss)
485 
486 ! Get shape functions and derivatives
487  CALL get_shape(ri,nnn,dn,igpt)
488 
489 ! Calculate the volume of each element. First calculate the Jacobian
490 ! at the "center" of the element.
491 
492 ! Compute jacobien
493  CALL get_jacobien(coord,3,8,dn,jac,jacinv,detj,error)
494 
495  coeff(igpt) = detj ! * wi(igpt) = 1.d0
496 
497  element_volume = element_volume + coeff(igpt)
498 
499  t(1:3,1:3) = t(1:3,1:3) + jac * coeff(igpt)
500 
501 ! Compute conventional strain displacement matrix <bmat> at each gauss point
502 
503  CALL tensormul(jacinv,3,3,'t',dn,8,3,'t',tmp38,3,one,zero)
504  CALL kronecker_product(tmp38,3,8,dident,3,3,grad,9)
505  CALL tensormul(sop,9,9,'n',grad,9,24,'n',bmat,9,one,zero)
506 
507  bmat2(igpt,:,:) = bmat(:,:)
508 
509  bmat_avg(:,:) = bmat_avg(:,:) + bmat(:,:) * coeff(igpt)
510 
511  END DO ! END 1st LOOP over gauss points
512 
513  t(1:3,1:3) = t(1:3,1:3) / element_volume ! Compute average jacobien <t>
514  tinv(1:3,1:3) = t(1:3,1:3) ! and its inverse
515  CALL invert3x3(tinv,det)
516 
517 ! Compute average <bmat>
518  bmat_avg = bmat_avg / element_volume
519 
520 ! Compute product matrices
521  CALL kronecker_product(t,3,3,t,3,3,tkront,9)
522  CALL kronecker_product(tinv,3,3,tinv,3,3,tkrontinv,9)
523 
524 ! Compute the mixed strain-displacement matrix <bsm>
525 
526  gbg = 0.
527  gsm(1:12,1:24) = 0.
528 
529 !! 2nd LOOP OVER GAUSS POINTS
530 
531 
532  DO igpt = 1, 8 ! LOOP over gauss points
533 
534  alpha2 = coeff(igpt)
535 
536  e_mixed(1:9,1:12) = mixed_map(igpt,1:9,1:12)
537  CALL tensormul(e_mixed,9,12,'t',e_mixed,9,12,'n',gbg,12,one,one)
538 
539 
540  bmat(1:9,1:24) = bmat2(igpt,1:9,1:24) - bmat_avg(1:9,1:24)
541  CALL tensormul(tkront,9,9,'t',bmat,9,24,'n',grad,9,one,zero)
542 
543 
544 ! Accumilate G small
545  CALL tensormul(e_mixed,9,12,'t',grad,9,24,'n',gsm,12,alpha2,one)
546 
547  END DO ! END LOOP over gauss points
548 
549 ! Compute 12x24 matrix b in (Eqn. 2.27)
550 !
551 ! (Eqn. 2.29)
552 
553  bsm(1,1:24) = gsm(1,1:24)*threeeighth
554  bsm(2,1:24) = gsm(2,1:24)*threeeighth
555  bsm(3,1:24) = gsm(3,1:24)* three*threeeighth
556  bsm(4,1:24) = gsm(4,1:24)*threeeighth
557  bsm(5,1:24) = gsm(5,1:24)*threeeighth
558  bsm(6,1:24) = gsm(6,1:24)* three*threeeighth
559  bsm(7,1:24) = gsm(7,1:24)*threeeighth
560  bsm(8,1:24) = gsm(8,1:24)*threeeighth
561  bsm(9,1:24) = gsm(9,1:24)* three*threeeighth
562  bsm(10,1:24) = gsm(10,1:24)*threeeighth*0.5
563  bsm(11,1:24) = gsm(11,1:24)*threeeighth*0.5
564  bsm(12,1:24) = gsm(12,1:24)*threeeighth*0.5
565 
566 ! Compute stiffness matrices <kaa>
567  kuu(1:24,1:24) = 0. ! initialize
568  kau(1:9, 1:24) = 0.
569  kaa(1:9, 1:9) = 0.
570 
571 ! Compute mixed/enhanced strain-displacement matrices <bu> & <ba> at
572 ! each gauss point
573 
574 
575  DO igpt = 1, 8 ! Start 3rd loop over gauss points
576 
577  detj = coeff(igpt) ! / wi(igpt) but wi(igpt) = 1. for all gauss points
578 
579 ! CALL get_mixed_map(e_mixed,ri,igpt)
580  e_mixed(1:9,1:12) = mixed_map(igpt,1:9,1:12)
581 
582 ! CALL get_enhanced_map(e_enhanced,ri,igpt)
583  e_enhanced(1:9,1:9) = enhanced_map(igpt,1:9,1:9)
584 
585 ! Compute <bu> (9 x 24) (Eqn. 2.27)
586 
587  CALL tensormul(e_mixed,9,12,'n',bsm,12,24,'n',grad,9,one,zero)
588  CALL tensormul(tkrontinv,9,9,'t',grad,9,24,'n',bu,9,one,zero)
589  bu(1:9,1:24) = bu(1:9,1:24) / detj + bmat_avg(1:9,1:24)
590 
591 ! Compute <ba> (9 x 9) (Eqn. 2.31)
592 
593  CALL tensormul(tkrontinv,9,9,'t',e_enhanced,9,9,'n',ba,9,one,zero)
594 
595 ! Compute the enhanced strain field
596 ! gauss pt, component, element
597 
598  ba(1:9,1:9) = ba(1:9,1:9) /detj ! for later
599 
600 !
601 ! Compute <kaa>
602 !
603  alpha = detj ! coeff(igpt)
604 
605  CALL tensormul(dmat(:,:),9,9,'n',bu,9,24,'n',grad,9,one,zero)
606 ! Compute <kuu>
607  CALL tensormul(bu,9,24,'t',grad,9,24,'n',kuu,24,detj,one)
608 ! Compute <kau>
609  CALL tensormul(ba,9,9, 't',grad,9,24,'n',kau,9,detj,one)
610 
611  CALL tensormul(dmat(:,:),9,9,'n',ba,9,9,'n',tkront,9,one,zero)
612 
613  CALL tensormul(ba,9,9,'t',tkront,9,9,'n',kaa,9,detj,one)
614 
615  END DO
616 
617 ! Compute <amatrx>
618 
619  amatrx = kuu ! Static condensation
620 ! -1
621 ! [ K ]
622 ! aa
623 
624  call invert(kaa,9,det)
625 
626  CALL tensormul(kaa,9,9,'n',kau,9,24,'n',grad,9,one,zero)
627  call tensormul(kau,9,24,'t',grad,9,24,'n',amatrx,24,-one,one)
628 
629  DO i=1,24
630  DO j=1,24
631  atol = abs(amatrx(i,j)-amatrx(j,i))
632  IF(atol.GE.1.0e-03) THEN
633  WRITE(7,*) '>>> amatrx is UNSYMMETRIC, fatal error', &
634  ' i =',i,' j =',j, ' diff =', atol
635  END IF
636  END DO
637  END DO
638 
639 !**************************************************************
640 ! Assemble local K (i.e. amatrix) into Global stiffness matrix
641 !**************************************************************
642 
643  DO dof1 = 1, 3 ! dof1 loop
644  DO dof2 = 1, 3 ! dof2 loop
645  DO i = 1, 8 ! node1 loop
646  gdof1 = 3 * local2global(elconnvol(i,ielem)) - 3 + dof1
647  ldof1 = 3 * i - 3 + dof1
648  DO j = 1, 8 ! node2 loop
649  ldof2 = 3 * j - 3 + dof2
650  gdof2 = 3 * local2global(elconnvol(j,ielem)) - 3 + dof2
651 
652  ! Make sure the matrix is exactly symmetric by using only values in the upper triangle of amatrx
653  IF ( ldof1 <= ldof2 ) THEN
654  tempval = amatrx(ldof1,ldof2)
655  ELSE
656  tempval = amatrx(ldof2,ldof1)
657  ENDIF
658 
659  ! Place the nonzeros into the K matrix
660  IF (tempval /= 0.0) THEN
661  counter = 0
662  DO m = rp_k(gdof1)+1, rp_k(gdof1+1)
663  IF (cval_k(m) == gdof2-1) THEN
664  aval_k(m) = aval_k(m) + tempval
665  counter = 1
666  ENDIF
667  ENDDO
668  if(counter==0) print*,'WARNING: Unable to add value to K matrix at (',gdof1,',',gdof2, &
669  ') on processor ',myid
670  ELSE
671  print*,myid,'ZERO DETECTED IN LOCAL STIFFNESS MATRIX!',gdof1,gdof2
672  END IF
673 
674  ENDDO
675  ENDDO
676  ENDDO
677  ENDDO
678 
679 
680 ! Compute enhanced field parameters <enh>
681 
682 ! enh = zero
683 !
684 ! do i=1, maxiters
685 !
686 ! call tensormul(kau,9,24,'n',u,ndofel,1,'n',fenh,9,one,zero)
687 ! call tensormul(kaa_copy,9,9,'n',enh,9,1,'n',fenh,9,one,one)
688 !
689 ! call tensormul(kaa,9,9,'n',fenh,9,1,'n',denh,9,-one,zero)
690 !
691 ! enh = enh + denh
692 !
693 ! call tensormul(denh,9,1,'t',denh,9,1,'n',dnorm,1,one,zero)
694 !
695 ! write(7,*) ' >>> iterno = ',i,' dnorm =',dnorm
696 !
697 ! if(dnorm.le.0.0000001) then
698 ! write(7,*) '>>> CONVERGED'
699 ! call printmat(enh,9,1,otdev,' enh ')
700 ! goto 1
701 ! end if
702 !
703 ! end do
704 !
705 ! 1 continue ! Convergence achieved or max. iters exhausted
706 
707  ENDDO
708 
709  ALLOCATE(rp_temp(1:3*gnumnp+1),cval_temp(1:nnz_k),aval_temp(1:nnz_k))
710  nnz_temp = nnz_k
711  rp_temp = rp_k
712  cval_temp = cval_k
713  aval_temp = aval_k
714  DEALLOCATE(rp_k,cval_k,aval_k)
715 
716 
717  END SUBROUTINE implicit_v3d8_me_k
718 
719 
720  SUBROUTINE map_ci2dmat(ci,dmat)
721 
722 
723  REAL*8, DIMENSION(1:6,1:6) :: ci
724  REAL*8, DIMENSION(1:9,1:9) :: dmat
725 
726 
727  dmat(1,1) = ci(1,1)
728  dmat(1,2) = ci(1,4)
729  dmat(1,3) = ci(1,5)
730  dmat(1,4) = ci(1,4)
731  dmat(1,5) = ci(1,2)
732  dmat(1,6) = ci(1,6)
733  dmat(1,7) = ci(1,5)
734  dmat(1,8) = ci(1,6)
735  dmat(1,9) = ci(1,3)
736 
737  dmat(2,1) = ci(1,4)
738  dmat(2,2) = ci(4,4)
739  dmat(2,3) = ci(4,5)
740  dmat(2,4) = ci(4,4)
741  dmat(2,5) = ci(2,4)
742  dmat(2,6) = ci(4,6)
743  dmat(2,7) = ci(4,5)
744  dmat(2,8) = ci(4,6)
745  dmat(2,9) = ci(3,4)
746 
747  dmat(3,1) = ci(1,5)
748  dmat(3,2) = ci(4,5)
749  dmat(3,3) = ci(5,5)
750  dmat(3,4) = ci(4,5)
751  dmat(3,5) = ci(2,5)
752  dmat(3,6) = ci(5,6)
753  dmat(3,7) = ci(5,5)
754  dmat(3,8) = ci(5,6)
755  dmat(3,9) = ci(3,5)
756 
757  dmat(4,1) = ci(1,4)
758  dmat(4,2) = ci(4,4)
759  dmat(4,3) = ci(4,5)
760  dmat(4,4) = ci(4,4)
761  dmat(4,5) = ci(2,4)
762  dmat(4,6) = ci(4,6)
763  dmat(4,7) = ci(4,5)
764  dmat(4,8) = ci(4,6)
765  dmat(4,9) = ci(3,4)
766 
767  dmat(5,1) = ci(1,2)
768  dmat(5,2) = ci(2,4)
769  dmat(5,3) = ci(2,5)
770  dmat(5,4) = ci(2,4)
771  dmat(5,5) = ci(2,2)
772  dmat(5,6) = ci(2,6)
773  dmat(5,7) = ci(2,5)
774  dmat(5,8) = ci(2,6)
775  dmat(5,9) = ci(2,3)
776 
777  dmat(6,1) = ci(1,6)
778  dmat(6,2) = ci(4,6)
779  dmat(6,3) = ci(5,6)
780  dmat(6,4) = ci(4,6)
781  dmat(6,5) = ci(2,6)
782  dmat(6,6) = ci(6,6)
783  dmat(6,7) = ci(5,6)
784  dmat(6,8) = ci(6,6)
785  dmat(6,9) = ci(3,6)
786 
787  dmat(7,1) = ci(1,5)
788  dmat(7,2) = ci(4,5)
789  dmat(7,3) = ci(5,5)
790  dmat(7,4) = ci(4,5)
791  dmat(7,5) = ci(2,5)
792  dmat(7,6) = ci(5,6)
793  dmat(7,7) = ci(5,5)
794  dmat(7,8) = ci(5,6)
795  dmat(7,9) = ci(3,5)
796 
797 
798  dmat(8,1) = ci(1,6)
799  dmat(8,2) = ci(4,6)
800  dmat(8,3) = ci(5,6)
801  dmat(8,4) = ci(4,6)
802  dmat(8,5) = ci(2,6)
803  dmat(8,6) = ci(6,6)
804  dmat(8,7) = ci(5,6)
805  dmat(8,8) = ci(6,6)
806  dmat(8,9) = ci(3,6)
807 
808 
809  dmat(9,1) = ci(1,3)
810  dmat(9,2) = ci(3,4)
811  dmat(9,3) = ci(3,5)
812  dmat(9,4) = ci(3,4)
813  dmat(9,5) = ci(2,3)
814  dmat(9,6) = ci(3,6)
815  dmat(9,7) = ci(3,5)
816  dmat(9,8) = ci(3,6)
817  dmat(9,9) = ci(3,3)
818 
819  END SUBROUTINE map_ci2dmat
820 
numknnz
subroutine numknnz(NumNp, NumEl, ElConnVol, nnz)
Definition: initK.f90:78
Aff_transformation_2::m
FT m(int i, int j) const
Definition: Aff_transformation_2.h:228
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
comp_row_global
Definition: comp_row_global.f90:26
invert
subroutine invert(a, nrow, det)
Definition: v3d8_me.f90:1392
k
j indices k indices k
Definition: Indexing.h:6
tensormul
subroutine tensormul(a, lda, n, achar, b, ldb, m, bchar, c, ldc, alpha, beta)
Definition: v3d8_me.f90:1060
coord
int coord[NPANE][NROW *NCOL][3]
Definition: blastest.C:86
implicit_global
Definition: implicit_global.f90:76
kronecker_product
subroutine kronecker_product(a, lda, nda, b, ldb, ndb, c, ldc)
Definition: v3d8_me.f90:1625
Polyhedron_incremental_builder_3::error
bool error() const
Definition: Polyhedron_incremental_builder_3.h:178
invert3x3
subroutine invert3x3(a, det)
Definition: v3d8_me.f90:1454
get_jacobien
subroutine get_jacobien(coords, mcrd, nnode, dn, jac, jacinv, detj, error)
Definition: v3d8_me.f90:1003
initk
subroutine initk(NumNp, NumEl, ElConnVol, nnz, rp, cval, aval)
Definition: initK.f90:197
i
blockLoc i
Definition: read.cpp:79
map_ci2dmat
subroutine map_ci2dmat(ci, dmat)
Definition: implicit_v3d8_me_K.f90:720
Aff_transformation_repS2::print
virtual std::ostream & print(std::ostream &os) const
Definition: Aff_transformationS2.h:277
j
j indices j
Definition: Indexing.h:6
GEM_DomainBoundary::debug
bool debug(bool s=true)
Definition: GEM.H:193
Color::alpha
unsigned char alpha() const
Definition: Color.h:75
implicit_v3d8_me_k
subroutine implicit_v3d8_me_k(coor, ElConnVol, ci, numnp, nstart, nend, NumEL, MatType, NumMatType, enhanced_map, mixed_map, NumMatVol)
Definition: implicit_v3d8_me_K.f90:84
get_shape
subroutine get_shape(r, n, dn, igpt)
Definition: v3d8_me.f90:846