Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
v3d8_me.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 SUBROUTINE v3d8_me(coor,MatType,ElConnVol,Rnet,disp,dmat, &
54  s11,s22,s33,s12,s23,s13, &
55  numnp,nstart,nend,numel,nummattype, aenh,enhanced_map,mixed_map)
56 
57 !!****f* Rocfrac/Source/v3d8_me.f90
58 !!
59 !! NAME
60 !! v3d8_me
61 !!
62 !! FUNCTION
63 !!
64 !! Computes the internal force vector for an 8-node hexahedral
65 !! mixed enhanced element formulation.
66 !!
67 !! Assumes small deformation, linear elastic.
68 !!
69 !! Original author: Ertugrul Taciroglu
70 !!
71 !! Based on the CSAR internal document,
72 !! "A Linear Mixed-Enhanced Strain Element:
73 !! Formulation, Algorithms and Verification"
74 !!
75 !! Foundation of Theory:
76 !!
77 !! "A mixed-enhanced strain method Part I: Geometraically linear problems"
78 !!
79 !!
80 !! INPUTS
81 !!
82 !! NumNP -- Number of nodes
83 !! NumEL -- Number of elements
84 !! Coor -- number of coordinates
85 !! MatType -- Material id
86 !! disp -- Nodal Displacement
87 !! ElConnVol -- Nodal connectivity
88 !! Dmat -- material D matrix
89 !! S11, S22, S33, S12, S23, S13 -- Stress
90 !! nstart, nend -- element beginning and end loop counter
91 !! NumMatVol -- number of materials
92 !! Aenh -- Element enhancement parameter
93 !! enhanced_map -- enhancement map
94 !! mixed_map -- mixed map
95 !!
96 !!
97 !! OUTPUT
98 !!
99 !! Rnet -- internal force vector
100 !!
101 !!****
102 
103  IMPLICIT NONE
104 !-----Global variables
105  INTEGER :: numnp ! number of nodes
106  INTEGER :: numat_vol ! number of volumetric materials
107  INTEGER :: numel ! number of LSTets
108  integer :: nummattype
109 !-- coordinate array
110  REAL*8, DIMENSION(1:3,1:numnp) :: coor
111 !-- connectivity table for Brick
112  INTEGER, DIMENSION(1:8,1:NumEl) :: elconnvol
113  INTEGER, DIMENSION(1:NumEl) :: mattype
114 !-- elastic stiffness consts
115  REAL*8, DIMENSION(1:9,NumMatType) :: ci
116 !-- internal force
117  REAL*8, DIMENSION(1:3*numnp) :: r_in
118 !-- displacement vector
119  REAL*8, DIMENSION(1:3*numnp) :: disp
120 !---- Local variables
121 !-- node numbers
122  INTEGER :: n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
123 
124  REAL*8 :: rnet(1:3*numnp)
125 
126  INTEGER :: k1n1,k1n2,k1n3,k1n4,k1n5,k1n6,k1n7,k1n8
127  INTEGER :: k2n1,k2n2,k2n3,k2n4,k2n5,k2n6,k2n7,k2n8
128  INTEGER :: k3n1,k3n2,k3n3,k3n4,k3n5,k3n6,k3n7,k3n8
129 
130  REAL*8, DIMENSION(1:3,1:8) :: coord
131  REAL*8, DIMENSION(1:24) :: udisp
132 
133  REAL*8 :: element_volume
134  INTEGER :: ielem, imat
135  integer :: igpt
136  REAL*8, DIMENSION(1:24) :: fint
137  REAL*8 nnn(8), dn(8,3), jac(3,3), jacinv(3,3), &
138  dn20(20,3),nnn20(20),&
139  t(3,3), tinv(3,3), &
140  dmat(nummattype,9,9), bmat(9,24), bmat2(1:8,1:9,1:24), &
141  grad(9,24), tmp38(3,8), coeff(8),&
142  tkront(9,9), tkrontinv(9,9),&
143  e_mixed(9,12),e_mixedt(12,9), &
144  gsm(12,24), gbg(12,12), bmat_avg(9,24),&
145  bu(9,24), ba(9,9), bsm(12,24),& ! print edge
146  kuu(24,24), kau(9,24), kaa(9,9),atol, bavec(1:9,1:9,1:8), buvec(1:9,1:24,1:8),&
147  tmp249(24,9),tmp91(9,1),kaa_copy(9,9),enh2(9),dnorm,bmatigs(1:8,1:9,1:24),&
148  stresstmp(9,1), tmp9(9),dmatinf(nummattype,9,9),stress(1:8,1:9),sumvect(1:9),tkrontinvt(1:9,1:9),&
149  tkrontt(1:9,1:9)
150 
151 
152 ! weight for hex element
153 
154  REAL*8, DIMENSION(1:8) :: wi = &
155  (/1.000000000000000,1.000000000000000,1.000000000000000, &
156  1.000000000000000,1.000000000000000,1.000000000000000, &
157  1.000000000000000,1.000000000000000/)
158 
159 ! root3 = 1./sqrt(3.)
160 
161  REAL*8, DIMENSION(1:3,1:8) :: ri = reshape( &
162  (/-0.577350269189626,-0.577350269189626,-0.577350269189626, &
163  0.577350269189626,-0.577350269189626,-0.577350269189626, &
164  0.577350269189626, 0.577350269189626,-0.577350269189626, &
165  -0.577350269189626, 0.577350269189626,-0.577350269189626, &
166  -0.577350269189626,-0.577350269189626, 0.577350269189626, &
167  0.577350269189626,-0.577350269189626, 0.577350269189626, &
168  0.577350269189626, 0.577350269189626, 0.577350269189626, &
169  -0.577350269189626, 0.577350269189626, 0.577350269189626/),(/3,8/) )
170 
171  REAL*8, DIMENSION(1:9,1:9) :: sop = reshape( &
172  (/1.000000000000000,0.000000000000000,0.000000000000000, &
173  0.000000000000000,0.000000000000000,0.000000000000000, &
174  0.000000000000000,0.000000000000000,0.000000000000000, &
175  0.000000000000000,0.500000000000000,0.000000000000000, &
176  0.500000000000000,0.000000000000000,0.000000000000000, &
177  0.000000000000000,0.000000000000000,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.500000000000000,0.000000000000000, &
182  0.500000000000000,0.000000000000000,0.000000000000000, &
183  0.000000000000000,0.000000000000000,0.000000000000000, &
184  0.000000000000000,0.000000000000000,0.000000000000000, &
185  0.000000000000000,1.000000000000000,0.000000000000000, &
186  0.000000000000000,0.000000000000000,0.000000000000000, &
187  0.000000000000000,0.000000000000000,0.000000000000000, &
188  0.000000000000000,0.000000000000000,0.500000000000000, &
189  0.000000000000000,0.500000000000000,0.000000000000000, &
190  0.000000000000000,0.000000000000000,0.500000000000000, &
191  0.000000000000000,0.000000000000000,0.000000000000000, &
192  0.500000000000000,0.000000000000000,0.000000000000000, &
193  0.000000000000000,0.000000000000000,0.000000000000000, &
194  0.000000000000000,0.000000000000000,0.500000000000000, &
195  0.000000000000000,0.500000000000000,0.000000000000000, &
196  0.000000000000000,0.000000000000000,0.000000000000000, &
197  0.000000000000000,0.000000000000000,0.000000000000000, &
198  0.000000000000000,0.000000000000000,1.000000000000000/),(/9,9/) )
199 
200  REAL*8, DIMENSION(1:3,1:3) :: dident = reshape( &
201  (/1.000000000000000,0.000000000000000,0.000000000000000, &
202  0.000000000000000,1.000000000000000,0.000000000000000, &
203  0.000000000000000,0.000000000000000,1.000000000000000/),(/3,3/) )
204 
205  REAL*8 :: zero = 0.0, one = 1.d0, detj,det, three = 3.d0
206  LOGICAL error, debug
207  REAL*8 :: alpha2
208  REAL*8, DIMENSION(1:9) :: fenh
209  REAL*8, DIMENSION(1:9,1:9) :: e_enhanced
210  REAL*8 :: maxdisp, con1, tp, sum
211  REAL*8 :: threeeighth = 3./8.
212  REAL*8 :: denh
213  REAL*8 :: aenh(1:9,1:numel)
214  INTEGER :: igauss
215  REAL*8, DIMENSION(8) :: xi, eta, zeta
216  REAL*8, DIMENSION(8) :: xie, etae, zetae
217  REAL*8, DIMENSION(1:8,1:9,1:12) :: mixed_map
218  REAL*8, DIMENSION(1:8,1:9,1:9) :: enhanced_map
219  REAL*8, DIMENSION(1:8,1:8) :: shapefun,shapefune
220  INTEGER :: otdev,mcrd,nnode
221  INTEGER,parameter :: ngpts = 8
222  REAL*8 :: strainenh(1:ngpts,1:9,1:numel)
223  REAL*8, dimension(1:NumEl) :: s11,s22,s33,s12,s23,s13
224  INTEGER :: i,k,j
225  INTEGER :: nstart, nend
226 
227 
228 
229  DO ielem = nstart, nend ! Loop over elements
230 
231  imat = mattype(ielem)
232 
233  n1 = elconnvol(1,ielem)
234  n2 = elconnvol(2,ielem)
235  n3 = elconnvol(3,ielem)
236  n4 = elconnvol(4,ielem)
237  n5 = elconnvol(5,ielem)
238  n6 = elconnvol(6,ielem)
239  n7 = elconnvol(7,ielem)
240  n8 = elconnvol(8,ielem)
241 
242  k3n1 = 3*n1
243  k3n2 = 3*n2
244  k3n3 = 3*n3
245  k3n4 = 3*n4
246  k3n5 = 3*n5
247  k3n6 = 3*n6
248  k3n7 = 3*n7
249  k3n8 = 3*n8
250 
251  k2n1 = k3n1 - 1
252  k2n2 = k3n2 - 1
253  k2n3 = k3n3 - 1
254  k2n4 = k3n4 - 1
255  k2n5 = k3n5 - 1
256  k2n6 = k3n6 - 1
257  k2n7 = k3n7 - 1
258  k2n8 = k3n8 - 1
259 
260  k1n1 = k3n1 - 2
261  k1n2 = k3n2 - 2
262  k1n3 = k3n3 - 2
263  k1n4 = k3n4 - 2
264  k1n5 = k3n5 - 2
265  k1n6 = k3n6 - 2
266  k1n7 = k3n7 - 2
267  k1n8 = k3n8 - 2
268 
269  coord(1,1) = coor(1,n1)
270  coord(2,1) = coor(2,n1)
271  coord(3,1) = coor(3,n1)
272 
273  coord(1,2) = coor(1,n2)
274  coord(2,2) = coor(2,n2)
275  coord(3,2) = coor(3,n2)
276 
277  coord(1,3) = coor(1,n3)
278  coord(2,3) = coor(2,n3)
279  coord(3,3) = coor(3,n3)
280 
281  coord(1,4) = coor(1,n4)
282  coord(2,4) = coor(2,n4)
283  coord(3,4) = coor(3,n4)
284 
285  coord(1,5) = coor(1,n5)
286  coord(2,5) = coor(2,n5)
287  coord(3,5) = coor(3,n5)
288 
289  coord(1,6) = coor(1,n6)
290  coord(2,6) = coor(2,n6)
291  coord(3,6) = coor(3,n6)
292 
293  coord(1,7) = coor(1,n7)
294  coord(2,7) = coor(2,n7)
295  coord(3,7) = coor(3,n7)
296 
297  coord(1,8) = coor(1,n8)
298  coord(2,8) = coor(2,n8)
299  coord(3,8) = coor(3,n8)
300 
301  udisp(1) = disp(k1n1)
302  udisp(2) = disp(k2n1)
303  udisp(3) = disp(k3n1)
304 
305  udisp(4) = disp(k1n2)
306  udisp(5) = disp(k2n2)
307  udisp(6) = disp(k3n2)
308 
309  udisp(7) = disp(k1n3)
310  udisp(8) = disp(k2n3)
311  udisp(9) = disp(k3n3)
312 
313  udisp(10) = disp(k1n4)
314  udisp(11) = disp(k2n4)
315  udisp(12) = disp(k3n4)
316 
317  udisp(13) = disp(k1n5)
318  udisp(14) = disp(k2n5)
319  udisp(15) = disp(k3n5)
320 
321  udisp(16) = disp(k1n6)
322  udisp(17) = disp(k2n6)
323  udisp(18) = disp(k3n6)
324 
325  udisp(19) = disp(k1n7)
326  udisp(20) = disp(k2n7)
327  udisp(21) = disp(k3n7)
328 
329  udisp(22) = disp(k1n8)
330  udisp(23) = disp(k2n8)
331  udisp(24) = disp(k3n8)
332 
333  fint = 0.d0
334 
335 ! Compute average element quantities
336 
337  element_volume = 0. ! Initialize
338  t(1:3,1:3) = 0.
339  bmat_avg(1:9,1:24) = 0.
340 
341 !!$! Calculate the volume of each element. First calculate the Jacobian
342 !!$! at the "center" of the element.
343 !!$
344 !!$ jac(1,1) = eighth*( coord(1,1) - coord(1,2) - coord(1,3) + coord(1,4) + &
345 !!$ coord(1,5) - coord(1,6) - coord(1,7) + coord(1,8) )
346 !!$
347 !!$ jac(1,2) = eighth*( coord(2,1) - coord(2,2) - coord(2,3) + coord(2,4) + &
348 !!$ coord(2,5) - coord(2,6) - coord(2,7) + coord(2,8) )
349 !!$
350 !!$ jac(1,3) = eighth*( coord(3,1) - coord(3,2) - coord(3,3) + coord(3,4) + &
351 !!$ coord(3,5) - coord(3,6) - coord(3,7) + coord(3,8) )
352 !!$
353 !!$ jac(2,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(2,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(2,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(3,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(3,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(3,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 !!$
372 !!$ temp_det = jac(1,1)*( jac(2,2)*jac(3,3) - jac(2,3)*jac(3,2) ) - &
373 !!$ jac(1,2)*( jac(2,1)*jac(3,3) - jac(2,3)*jac(3,1) ) + &
374 !!$ jac(1,3)*( jac(2,1)*jac(3,2) - jac(2,2)*jac(3,1) )
375 !!$
376 !!$ volume = eight*temp_det
377  !
378 ! 1st LOOP OVER GUAUSS POINTS
379 
380  DO igpt = 1, 8 ! LOOP over gauss points
381 
382 !!$! Calculate the derivatives of the shape functions at the gauss points
383 !!$
384 !!$ dn(1,1) = eighth*(one + eta(igauss))*(one + zeta(igauss))
385 !!$ dn(1,2) = eighth*(one + xi(igauss) )*(one + zeta(igauss))
386 !!$ dn(1,3) = eighth*(one + xi(igauss) )*(one + eta(igauss) )
387 !!$
388 !!$ dn(2,1) = -eighth*(one + eta(igauss))*(one + zeta(igauss))
389 !!$ dn(2,2) = eighth*(one - xi(igauss) )*(one + zeta(igauss))
390 !!$ dn(2,3) = eighth*(one - xi(igauss) )*(one + eta(igauss) )
391 !!$
392 !!$ dn(3,1) = -eighth*(one - eta(igauss))*(one + zeta(igauss))
393 !!$ dn(3,2) = -eighth*(one - xi(igauss) )*(one + zeta(igauss))
394 !!$ dn(3,3) = eighth*(one - xi(igauss) )*(one - eta(igauss) )
395 !!$
396 !!$ dn(4,1) = eighth*(one - eta(igauss))*(one + zeta(igauss))
397 !!$ dn(4,2) = -eighth*(one + xi(igauss) )*(one + zeta(igauss))
398 !!$ dn(4,3) = eighth*(one + xi(igauss) )*(one - eta(igauss) )
399 !!$
400 !!$ dn(5,1) = eighth*(one + eta(igauss))*(one - zeta(igauss))
401 !!$ dn(5,2) = eighth*(one + xi(igauss) )*(one - zeta(igauss))
402 !!$ dn(5,3) = -eighth*(one + xi(igauss) )*(one + eta(igauss) )
403 !!$
404 !!$ dn(6,1) = -eighth*(one + eta(igauss))*(one - zeta(igauss))
405 !!$ dn(6,2) = eighth*(one - xi(igauss) )*(one - zeta(igauss))
406 !!$ dn(6,3) = -eighth*(one - xi(igauss) )*(one + eta(igauss) )
407 !!$
408 !!$ dn(7,1) = -eighth*(one - eta(igauss))*(one - zeta(igauss))
409 !!$ dn(7,2) = -eighth*(one - xi(igauss) )*(one - zeta(igauss))
410 !!$ dn(7,3) = -eighth*(one - xi(igauss) )*(one - eta(igauss) )
411 !!$
412 !!$ dn(8,1) = eighth*(one - eta(igauss))*(one - zeta(igauss))
413 !!$ dn(8,2) = -eighth*(one + xi(igauss) )*(one - zeta(igauss))
414 !!$ dn(8,3) = -eighth*(one + xi(igauss) )*(one - eta(igauss) )
415 !!$
416 !!$! Calculate the Jacobian, its determinant and inverse for each gauss point.
417 !!$
418 !!$ jac(1,1) = dn(1,1)*coord(1,1) + dn(2,1)*coord(1,2) &
419 !!$ + dn(3,1)*coord(1,3) + dn(4,1)*coord(1,4) &
420 !!$ + dn(5,1)*coord(1,5) + dn(6,1)*coord(1,6) &
421 !!$ + dn(7,1)*coord(1,7) + dn(8,1)*coord(1,8)
422 !!$
423 !!$ jac(2,1) = dn(1,1)*coord(2,1) + dn(2,1)*coord(2,2) &
424 !!$ + dn(3,1)*coord(2,3) + dn(4,1)*coord(2,4) &
425 !!$ + dn(5,1)*coord(2,5) + dn(6,1)*coord(2,6) &
426 !!$ + dn(7,1)*coord(2,7) + dn(8,1)*coord(2,8)
427 !!$
428 !!$ jac(3,1) = dn(1,1)*coord(3,1) + dn(2,1)*coord(3,2) &
429 !!$ + dn(3,1)*coord(3,3) + dn(4,1)*coord(3,4) &
430 !!$ + dn(5,1)*coord(3,5) + dn(6,1)*coord(3,6) &
431 !!$ + dn(7,1)*coord(3,7) + dn(8,1)*coord(3,8)
432 !!$
433 !!$ jac(1,2) = dn(1,2)*coord(1,1) + dn(2,2)*coord(1,2) &
434 !!$ + dn(3,2)*coord(1,3) + dn(4,2)*coord(1,4) &
435 !!$ + dn(5,2)*coord(1,5) + dn(6,2)*coord(1,6) &
436 !!$ + dn(7,2)*coord(1,7) + dn(8,2)*coord(1,8)
437 !!$
438 !!$ jac(2,2) = dn(1,2)*coord(2,1) + dn(2,2)*coord(2,2) &
439 !!$ + dn(3,2)*coord(2,3) + dn(4,2)*coord(2,4) &
440 !!$ + dn(5,2)*coord(2,5) + dn(6,2)*coord(2,6) &
441 !!$ + dn(7,2)*coord(2,7) + dn(8,2)*coord(2,8)
442 !!$
443 !!$ jac(3,2) = dn(1,2)*coord(3,1) + dn(2,2)*coord(3,2) &
444 !!$ + dn(3,2)*coord(3,3) + dn(4,2)*coord(3,4) &
445 !!$ + dn(5,2)*coord(3,5) + dn(6,2)*coord(3,6) &
446 !!$ + dn(7,2)*coord(3,7) + dn(8,2)*coord(3,8)
447 !!$
448 !!$ jac(1,3) = dn(1,3)*coord(1,1) + dn(2,3)*coord(1,2) &
449 !!$ + dn(3,3)*coord(1,3) + dn(4,3)*coord(1,4) &
450 !!$ + dn(5,3)*coord(1,5) + dn(6,3)*coord(1,6) &
451 !!$ + dn(7,3)*coord(1,7) + dn(8,3)*coord(1,8)
452 !!$
453 !!$ jac(2,3) = dn(1,3)*coord(2,1) + dn(2,3)*coord(2,2) &
454 !!$ + dn(3,3)*coord(2,3) + dn(4,3)*coord(2,4) &
455 !!$ + dn(5,3)*coord(2,5) + dn(6,3)*coord(2,6) &
456 !!$ + dn(7,3)*coord(2,7) + dn(8,3)*coord(2,8)
457 !!$
458 !!$ jac(3,3) = dn(1,3)*coord(3,1) + dn(2,3)*coord(3,2) &
459 !!$ + dn(3,3)*coord(3,3) + dn(4,3)*coord(3,4) &
460 !!$ + dn(5,3)*coord(3,5) + dn(6,3)*coord(3,6) &
461 !!$ + dn(7,3)*coord(3,7) + dn(8,3)*coord(3,8)
462 !!$
463 !!$ det(ielem,igauss) = jac(1,1)*( jac(2,2)*jac(3,3) - jac(2,3)*jac(3,2) ) - &
464 !!$ jac(1,2)*( jac(2,1)*jac(3,3) - jac(2,3)*jac(3,1) ) + &
465 !!$ jac(1,3)*( jac(2,1)*jac(3,2) - jac(2,2)*jac(3,1) )
466 !!$
467 !!$ jacinv(1,1) = ( jac(2,2)*jac(3,3) - jac(2,3)*jac(3,2) ) / det(ielem,igauss)
468 !!$ jacinv(1,2) = -( jac(1,2)*jac(3,3) - jac(1,3)*jac(3,2) ) / det(ielem,igauss)
469 !!$ jacinv(1,3) = ( jac(1,2)*jac(2,3) - jac(1,3)*jac(2,2) ) / det(ielem,igauss)
470 !!$ jacinv(2,1) = -( jac(2,1)*jac(3,3) - jac(2,3)*jac(3,1) ) / det(ielem,igauss)
471 !!$ jacinv(2,2) = ( jac(1,1)*jac(3,3) - jac(1,3)*jac(3,1) ) / det(ielem,igauss)
472 !!$ jacinv(2,3) = -( jac(1,1)*jac(2,3) - jac(1,3)*jac(2,1) ) / det(ielem,igauss)
473 !!$ jacinv(3,1) = ( jac(2,1)*jac(3,2) - jac(2,2)*jac(3,1) ) / det(ielem,igauss)
474 !!$ jacinv(3,2) = -( jac(1,1)*jac(3,2) - jac(1,2)*jac(3,1) ) / det(ielem,igauss)
475 !!$ jacinv(3,3) = ( jac(1,1)*jac(2,2) - jac(1,2)*jac(2,1) ) / det(ielem,igauss)
476 
477 ! Get shape functions and derivatives
478  CALL get_shape(ri,nnn,dn,igpt)
479 
480 ! Calculate the volume of each element. First calculate the Jacobian
481 ! at the "center" of the element.
482 
483 ! Compute jacobien
484  CALL get_jacobien(coord,3,8,dn,jac,jacinv,detj,error)
485 
486  coeff(igpt) = detj ! * wi(igpt) = 1.d0
487 
488  element_volume = element_volume + coeff(igpt)
489 
490  t(1:3,1:3) = t(1:3,1:3) + jac * coeff(igpt)
491 
492 ! Compute conventional strain displacement matrix <bmat> at each gauss point
493 
494 ! CALL TENSORMUL3(jacinv,3,3,'t',dn,8,3,'t',tmp38,3,8,one,zero)
495 
496  DO i = 1, 3
497  DO j = 1, 8
498  sum = 0.d0
499  DO k = 1, 3
500  sum = sum + jacinv(k,i) * dn(j,k)
501  ENDDO
502  tmp38(i,j) = sum
503  end DO
504  enddo
505 
506  CALL kronecker_product(tmp38,3,8,dident,3,3,grad,9)
507 
508 ! CALL TENSORMUL3(sop,9,9,'n',grad,9,24,'n',bmat,9,24,one,zero)
509 
510  DO i = 1, 9
511  DO j = 1, 24
512  sum = 0.d0
513  DO k = 1, 9
514  sum = sum + sop(i,k) * grad(k,j)
515  ENDDO
516  bmat(i,j) = sum
517  end DO
518  enddo
519 
520  bmat2(igpt,:,:) = bmat(:,:)
521 
522  bmat_avg(:,:) = bmat_avg(:,:) + bmat(:,:) * coeff(igpt)
523 
524  END DO ! END 1st LOOP over gauss points
525 
526  t(1:3,1:3) = t(1:3,1:3) / element_volume ! Compute average jacobien <t>
527  tinv(1:3,1:3) = t(1:3,1:3) ! and its inverse
528  CALL invert3x3(tinv,det)
529 
530 ! Compute average <bmat>
531  bmat_avg = bmat_avg / element_volume
532 
533 ! Compute product matrices
534  CALL kronecker_product(t,3,3,t,3,3,tkront,9)
535  CALL kronecker_product(tinv,3,3,tinv,3,3,tkrontinv,9)
536 
537 ! Compute the mixed strain-displacement matrix <bsm>
538 
539  gsm(1:12,1:24) = 0.
540 
541 !! 2nd LOOP OVER GAUSS POINTS
542 
543  DO igpt = 1, 8 ! LOOP over gauss points
544 
545  alpha2 = coeff(igpt)
546 
547  e_mixed(1:9,1:12) = mixed_map(igpt,1:9,1:12)
548 
549  bmat(1:9,1:24) = bmat2(igpt,1:9,1:24) - bmat_avg(1:9,1:24)
550 ! CALL TENSORMUL3(tkront,9,9,'t',bmat,9,24,'n',grad,9,24,one,zero)
551 
552  tkrontt = transpose(tkront)
553  DO i = 1, 9
554  DO j = 1, 24
555  sum = 0.d0
556  DO k = 1, 9
557  sum = sum + tkrontt(i,k) * bmat(k,j)
558  ENDDO
559  grad(i,j) = sum
560  end DO
561  enddo
562 
563 
564 ! Accumilate G small
565 ! CALL TENSORMUL3(e_mixed,9,12,'t',grad,9,24,'n',gsm,12,24,alpha2,one)
566 
567 ! e_mixedT = transpose(e_mixed)
568  DO i = 1, 12
569  DO j = 1, 24
570  sum = 0.d0
571  DO k = 1, 9
572  sum = sum + e_mixed(k,i) * grad(k,j)
573  ENDDO
574  gsm(i,j) = gsm(i,j) + alpha2*sum
575  end DO
576  enddo
577 
578 
579  END DO ! END LOOP over gauss points
580 
581 ! Compute 12x24 matrix b in (Eqn. 2.27)
582 !
583 ! (Eqn. 2.29)
584 
585  bsm(1,1:24) = gsm(1,1:24)*threeeighth
586  bsm(2,1:24) = gsm(2,1:24)*threeeighth
587  bsm(3,1:24) = gsm(3,1:24)* three*threeeighth
588  bsm(4,1:24) = gsm(4,1:24)*threeeighth
589  bsm(5,1:24) = gsm(5,1:24)*threeeighth
590  bsm(6,1:24) = gsm(6,1:24)* three*threeeighth
591  bsm(7,1:24) = gsm(7,1:24)*threeeighth
592  bsm(8,1:24) = gsm(8,1:24)*threeeighth
593  bsm(9,1:24) = gsm(9,1:24)* three*threeeighth
594  bsm(10,1:24) = gsm(10,1:24)*threeeighth*0.5
595  bsm(11,1:24) = gsm(11,1:24)*threeeighth*0.5
596  bsm(12,1:24) = gsm(12,1:24)*threeeighth*0.5
597 
598 
599 ! Compute stiffness matrices <kaa>
600 
601  kaa(1:9, 1:9) = 0.d0
602  fenh(:) = 0.d0
603 
604 ! Compute mixed/enhanced strain-displacement matrices <bu> & <ba> at
605 ! each gauss point
606 
607  DO igpt = 1, 8 ! Start 3rd loop over gauss points
608 
609  detj = coeff(igpt) ! / wi(igpt) but wi(igpt) = 1. for all gauss points
610 
611 ! CALL get_mixed_map(e_mixed,ri,igpt)
612  e_mixed(1:9,1:12) = mixed_map(igpt,1:9,1:12)
613 
614 ! CALL get_enhanced_map(e_enhanced,ri,igpt)
615  e_enhanced(1:9,1:9) = enhanced_map(igpt,1:9,1:9)
616 
617 ! Compute <bu> (9 x 24) (Eqn. 2.27)
618 
619 ! CALL TENSORMUL3(e_mixed,9,12,'n',bsm,12,24,'n',grad,9,24,one,zero)
620 
621  DO i = 1, 9
622  DO j = 1, 24
623  sum = 0.d0
624  DO k = 1, 12
625  sum = sum + e_mixed(i,k) * bsm(k,j)
626  ENDDO
627  grad(i,j) = sum
628  end DO
629  enddo
630 
631 
632  ! CALL TENSORMUL3(tkrontinv,9,9,'t',grad,9,24,'n',bu,9,24,one,zero)
633 
634  tkrontinvt = transpose(tkrontinv)
635  DO i = 1, 9
636  DO j = 1, 24
637  sum = 0.d0
638  DO k = 1, 9
639  sum = sum + tkrontinvt(i,k) * grad(k,j)
640  ENDDO
641  bu(i,j) = sum
642  end DO
643  enddo
644 
645  bu(1:9,1:24) = bu(1:9,1:24) / detj + bmat_avg(1:9,1:24)
646 
647 ! Compute <ba> (9 x 9) (Eqn. 2.31)
648 
649 ! CALL TENSORMUL3(tkrontinv,9,9,'t',e_enhanced,9,9,'n',ba,9,9,one,zero)
650 
651  DO i = 1, 9
652  DO j = 1, 9
653  sum = 0.d0
654  DO k = 1, 9
655  sum = sum + tkrontinvt(i,k) * e_enhanced(k,j)
656  ENDDO
657  ba(i,j) = sum
658  end DO
659  enddo
660 
661 ! Compute the enhanced strain field
662 ! gauss pt, component, element
663 
664  DO i = 1, 9
665  sum = 0.
666  DO k = 1, 24
667  sum = sum + bu(i,k)*udisp(k)
668  ENDDO
669  DO k = 1, 9
670  sum = sum + ba(i,k)*aenh(k,ielem)/detj ! times ba()/detj
671  ENDDO
672  strainenh(igpt,i,ielem) = sum
673  ENDDO
674  ba(1:9,1:9) = ba(1:9,1:9) /detj ! for later
675 
676 ! tmp9(1:9) = Aenh(:,ielem)
677 ! StrainEnh(igpt,:,ielem) = MATMUL(bu,Udisp) + MATMUL(ba,tmp9)
678 
679 ! Compute the enhanced stress field
680 
681  DO i = 1, 9
682  sum = 0.
683  DO k = 1,9
684  sum = sum + dmat(imat,i,k)*strainenh(igpt,k,ielem)
685  ENDDO
686  stress(igpt,i) = sum
687  ENDDO
688 
689 ! stress(igpt,:,ielem) = MATMUL(dmat,StrainEnh(igpt,:,ielem))
690 
691 ! Compute the enhanced field parameter
692 
693 ! T
694 ! F = B S dV (Accumulate Fenh from each guass point)
695 ! enh a
696 
697  DO i = 1, 9
698  sum = 0.d0
699  DO k = 1, 9
700  sum = sum + ba(k,i)*stress(igpt,k)
701  ENDDO
702  fenh(i) = fenh(i) + detj*sum
703  ENDDO
704 
705 ! CALL TENSORMUL(ba,9,9,'t',stresstmp,9,1,'n',fenh,9,detj,one)
706 
707 ! Compute the enhanced field parameter using the enhanced stress
708 
709 ! T
710 ! F = B S dV (Accumulate Fint from each guass point)
711 ! int u
712 
713  DO i = 1, 24
714  sum = 0.
715  DO k = 1, 9
716  sum = sum + bu(k,i)*stress(igpt,k)
717  ENDDO
718  fint(i) = fint(i) + detj*sum
719  ENDDO
720 
721 ! CALL TENSORMUL(bu,9,24,'t',stresstmp,9,1,'n',fint,24,detj,one)
722 !
723 ! Compute <kaa>
724 !
725 
726 ! CALL TENSORMUL3(dmat(imat,:,:),9,9,'n',ba,9,9,'n',tkront,9,9,one,zero)
727 
728  DO i = 1, 9
729  DO j = 1, 9
730  sum = 0.d0
731  DO k = 1, 9
732  sum = sum + dmat(imat,i,k) * ba(k,j)
733  ENDDO
734  tkront(i,j) = sum
735  end DO
736  enddo
737 
738 
739 ! CALL TENSORMUL3(ba,9,9,'t',tkront,9,9,'n',kaa,9,9,detj,one)
740 
741  DO i = 1, 9
742  DO j = 1, 9
743  sum = 0.d0
744  DO k = 1, 9
745  sum = sum + ba(k,i) * tkront(k,j)
746  ENDDO
747  kaa(i,j) = kaa(i,j) + detj*sum
748  end DO
749  enddo
750 
751  END DO
752 
753 
754 ! -1
755 ! [ K ]
756 ! aa
757 
758  CALL invert(kaa,9,det)
759 
760 ! Compute Eqn. (2.44) modified for explicit integration
761 !
762 ! -1 enh
763 ! -[ K ] * ( F ) = DA = Denh
764 ! aa enh
765 
766  DO i = 1, 9
767 
768  denh = 0.
769 
770 ! Matrix-Vector multiplication
771 
772  DO k = 1, 9
773  denh = denh - kaa(i,k)*fenh(k)
774  ENDDO
775 !
776 ! enh enh enh
777 ! A = A + DA
778 ! i+1 i i
779 !
780  aenh(i,ielem) = aenh(i,ielem) + denh
781 
782  ENDDO
783 
784 
785 ! ASSEMBLE THE INTERNAL + Enhanced FORCE VECTOR
786 !
787 ! local node 1
788  rnet(k1n1) = rnet(k1n1) - fint(1)
789  rnet(k2n1) = rnet(k2n1) - fint(2)
790  rnet(k3n1) = rnet(k3n1) - fint(3)
791 ! local node 2
792  rnet(k1n2) = rnet(k1n2) - fint(4)
793  rnet(k2n2) = rnet(k2n2) - fint(5)
794  rnet(k3n2) = rnet(k3n2) - fint(6)
795 ! local node 3
796  rnet(k1n3) = rnet(k1n3) - fint(7)
797  rnet(k2n3) = rnet(k2n3) - fint(8)
798  rnet(k3n3) = rnet(k3n3) - fint(9)
799 ! local node 4
800  rnet(k1n4) = rnet(k1n4) - fint(10)
801  rnet(k2n4) = rnet(k2n4) - fint(11)
802  rnet(k3n4) = rnet(k3n4) - fint(12)
803 ! local node 5
804  rnet(k1n5) = rnet(k1n5) - fint(13)
805  rnet(k2n5) = rnet(k2n5) - fint(14)
806  rnet(k3n5) = rnet(k3n5) - fint(15)
807 ! local node 6
808  rnet(k1n6) = rnet(k1n6) - fint(16)
809  rnet(k2n6) = rnet(k2n6) - fint(17)
810  rnet(k3n6) = rnet(k3n6) - fint(18)
811 ! local node 7
812  rnet(k1n7) = rnet(k1n7) - fint(19)
813  rnet(k2n7) = rnet(k2n7) - fint(20)
814  rnet(k3n7) = rnet(k3n7) - fint(21)
815 ! local node 8
816  rnet(k1n8) = rnet(k1n8) - fint(22)
817  rnet(k2n8) = rnet(k2n8) - fint(23)
818  rnet(k3n8) = rnet(k3n8) - fint(24)
819 
820 ! average stress over an element
821 
822  sumvect = 0.d0
823  DO k = 1, 9
824  DO i = 1,8
825  sumvect(k) = sumvect(k) + stress(i,k)*.125d0
826  enddo
827  enddo
828 
829  s11(ielem) = sumvect(1)
830  s22(ielem) = sumvect(5)
831  s33(ielem) = sumvect(9)
832  s12(ielem) = sumvect(4)
833  s23(ielem) = sumvect(8)
834  s13(ielem) = sumvect(7)
835 
836  ENDDO
837 
838 ! PRINT*,'Max Aenh =',MAXVAL(Aenh(:,:))
839 
840  END SUBROUTINE v3d8_me
841 
842 
843 
844 !--------------------------------------------------------------GET_SHAPE
845 
846 SUBROUTINE get_shape(r,n,dn,igpt)
847 ! *--------------------------------------------------------------------*
848 ! | |
849 ! | Returns the shape function N and its derivative matrix DN |
850 ! | |
851 ! | <r> : 3 x ngpts vector of coordinates of the isoparametric |
852 ! | point (usually an integration point) |
853 ! | <igpt> : no. of the particular gauss point |
854 ! | <n> : 8 x 1 shape function matrix evaluated at <ri> |
855 ! | <dn> : 8 x 3 shape function derivative matrix |
856 ! | evaluated at <ri> |
857 ! | |
858 ! *--------------------------------------------------------------------*
859 
860  IMPLICIT real*8 (a-h, o-z)
861 ! include 'ABA_PARAM.INC'
862 
863 ! Argument variables
864 
865  DOUBLE PRECISION r(3,*), n(*), dn(8,*)
866 
867 ! Data statements
868 
869  DATA one /1.000000000000000/
870 
871 ! Initializations
872 
873  r1 = r(1,igpt)
874  r2 = r(2,igpt)
875  r3 = r(3,igpt)
876 
877 ! Compute shape functions <n>
878 
879  n(1) = -(r1 - one) * (r2 - one) * (r3 - one)
880  n(2) = (r1 + one) * (r2 - one) * (r3 - one)
881  n(3) = -(r1 + one) * (r2 + one) * (r3 - one)
882  n(4) = (r1 - one) * (r2 + one) * (r3 - one)
883  n(5) = (r1 - one) * (r2 - one) * (r3 + one)
884  n(6) = -(r1 + one) * (r2 - one) * (r3 + one)
885  n(7) = (r1 + one) * (r2 + one) * (r3 + one)
886  n(8) = -(r1 - one) * (r2 + one) * (r3 + one)
887 
888  n(1:8) = 0.125000000000000 * n(1:8)
889 
890 ! Compute shape function derivative <dn>
891 
892  dn(1, 1) = -(r2 - one) * (r3 - one)
893  dn(2, 1) = (r2 - one) * (r3 - one)
894  dn(3, 1) = -(r2 + one) * (r3 - one)
895  dn(4, 1) = (r2 + one) * (r3 - one)
896  dn(5, 1) = (r2 - one) * (r3 + one)
897  dn(6, 1) = -(r2 - one) * (r3 + one)
898  dn(7, 1) = (r2 + one) * (r3 + one)
899  dn(8, 1) = -(r2 + one) * (r3 + one)
900 
901  dn(1, 2) = -(r1 - one) * (r3 - one)
902  dn(2, 2) = (r1 + one) * (r3 - one)
903  dn(3, 2) = -(r1 + one) * (r3 - one)
904  dn(4, 2) = (r1 - one) * (r3 - one)
905  dn(5, 2) = (r1 - one) * (r3 + one)
906  dn(6, 2) = -(r1 + one) * (r3 + one)
907  dn(7, 2) = (r1 + one) * (r3 + one)
908  dn(8, 2) = -(r1 - one) * (r3 + one)
909 
910  dn(1, 3) = -(r1 - one) * (r2 - one)
911  dn(2, 3) = (r1 + one) * (r2 - one)
912  dn(3, 3) = -(r1 + one) * (r2 + one)
913  dn(4, 3) = (r1 - one) * (r2 + one)
914  dn(5, 3) = (r1 - one) * (r2 - one)
915  dn(6, 3) = -(r1 + one) * (r2 - one)
916  dn(7, 3) = (r1 + one) * (r2 + one)
917  dn(8, 3) = -(r1 - one) * (r2 + one)
918 
919  dn(1:8, 1:3) = 0.125000000000000 * dn(1:8, 1:3)
920 
921  RETURN
922 END SUBROUTINE get_shape
923 
924 SUBROUTINE get_shape20(r,n,dn,igpt)
925 ! *--------------------------------------------------------------------*
926 ! | |
927 ! | Returns the shape function N and its derivative matrix DN |
928 ! | |
929 ! | <r> : 3 x ngpts vector of coordinates of the isoparametric |
930 ! | point (usually an integration point) |
931 ! | <igpt> : no. of the particular gauss point |
932 ! | <n> : 8 x 1 shape function matrix evaluated at <ri> |
933 ! | <dn> : 8 x 3 shape function derivative matrix |
934 ! | evaluated at <ri> |
935 ! | |
936 ! *--------------------------------------------------------------------*
937 
938  IMPLICIT real*8 (a-h, o-z)
939 ! include 'ABA_PARAM.INC'
940 
941 ! Argument variables
942 
943  DOUBLE PRECISION r(3,*), n(*), dn(20,*)
944 
945  INTEGER :: xii(20), etai(20), zetai(20)
946 
947 ! Data statements
948 
949  DATA one /1.000000000000000/
950 
951 ! Initializations
952 
953  xi = r(1,igpt)
954  eta = r(2,igpt)
955  zeta = r(3,igpt)
956 
957  xii=(/-1,-1,-1,0,1,1,1,0,-1,-1,1,1,-1,-1,-1,0,1,1,1,0/)
958  etai=(/-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,0,1,1,1,1,1,1,1,1/)
959  zetai=(/-1,0,1,1,1,0,-1,-1,-1,1,1,-1,-1,0,1,1,1,0,-1,-1/)
960 
961  DO l = 1, 20
962  xi0=xi*xii(l)
963  eta0=eta*etai(l)
964  zeta0=zeta*zetai(l)
965 
966 ! Compute shape functions <n>
967  IF(l==4.OR.l==8.OR.l==16.OR.l==20) THEN
968  n(l)=.25*(1.-xi*xi)*(1.+eta0)*(1.+zeta0)
969  ELSE IF(l>=9.AND.l<=12)THEN
970  n(l)=.25*(1.+xi0)*(1.-eta*eta)*(1.+zeta0)
971  ELSE IF(l==2.OR.l==6.OR.l==14.OR.l==18) THEN
972  n(l)=.25*(1.+xi0)*(1.+eta0)*(1.-zeta*zeta)
973  ELSE
974  n(l)=.125*(1.+xi0)*(1.+eta0)*(1.+zeta0)*(xi0+eta0+zeta0-2)
975  END IF
976 
977 ! Compute shape function derivative <dn>
978  IF(l==4.OR.l==8.OR.l==16.OR.l==20) THEN
979  dn(l,1)=-.5*xi*(1.+eta0)*(1.+zeta0)
980  dn(l,2)=.25*etai(l)*(1.-xi*xi)*(1.+zeta0)
981  dn(l,3)=.25*zetai(l)*(1.-xi*xi)*(1.+eta0)
982  ELSE IF(l>=9.AND.l<=12)THEN
983  dn(l,1)=.25*xii(l)*(1.-eta*eta)*(1.+zeta0)
984  dn(l,2)=-.5*eta*(1.+xi0)*(1.+zeta0)
985  dn(l,3)=.25*zetai(l)*(1.+xi0)*(1.-eta*eta)
986  ELSE IF(l==2.OR.l==6.OR.l==14.OR.l==18) THEN
987  dn(l,1)=.25*xii(l)*(1.+eta0)*(1.-zeta*zeta)
988  dn(l,2)=.25*etai(l)*(1.+xi0)*(1.-zeta*zeta)
989  dn(l,3)=-.5*zeta*(1.+xi0)*(1.+eta0)
990  ELSE
991  dn(l,1)=.125*xii(l)*(1.+eta0)*(1.+zeta0)*(2.*xi0+eta0+zeta0-1.)
992  dn(l,2)=.125*etai(l)*(1.+xi0)*(1.+zeta0)*(xi0+2.*eta0+zeta0-1.)
993  dn(l,3)=.125*zetai(l)*(1.+xi0)*(1.+eta0)*(xi0+eta0+2.*zeta0-1.)
994  END IF
995  ENDDO
996 
997  RETURN
998 END SUBROUTINE get_shape20
999 
1000 
1001 
1002 !-----------------------------------------------------------GET_JACOBIEN
1003 SUBROUTINE get_jacobien(coords,mcrd,nnode,dn,&
1004  jac,jacinv,detj,error)
1005 ! *--------------------------------------------------------------------*
1006 ! | |
1007 ! | Returns the jacobien, its inverse and its determinant |
1008 ! | |
1009 ! | <coords> : mcrd x nnode matrix of element nodal coordinates |
1010 ! | <dn> : 8 x 3 shape function derivative matrix |
1011 ! | <jac> : x jacobien matrix |
1012 ! | <dn> : 8 x 3 inverse jacobien matrix |
1013 ! | <detj> : determinant of the jacobien matrix |
1014 ! | <error> : logical varible; set to true if determinant is zero |
1015 ! | |
1016 ! | calls : TENSORMUL, INVERT3x3 |
1017 ! | |
1018 ! *--------------------------------------------------------------------*
1019 
1020  IMPLICIT real*8 (a-h, o-z)
1021 ! include 'ABA_PARAM.INC'
1022 
1023 ! Argument variables
1024 
1025  DOUBLE PRECISION jac(mcrd,mcrd), jacinv(mcrd,mcrd)
1026  dimension coords(1:3,1:nnode), dn(nnode,mcrd)
1027  LOGICAL error
1028 
1029 ! Data statements
1030 
1031  DATA zero,one /0.000000000000000,1.000000000000000/
1032 
1033 ! Initialize variables
1034 
1035  det = zero
1036  error = .true.
1037 
1038 ! Compute jacobien
1039 
1040  CALL tensormul(coords, mcrd, nnode,'n', &
1041  dn, nnode, mcrd, 'n', jac, mcrd, one, zero)
1042 
1043 ! call printmat(jac,mcrd,mcrd,7,' jac ')
1044 
1045 ! Compute inverse jacobien and determinant
1046 
1047  jacinv = jac
1048 
1049  CALL invert3x3(jacinv,detj)
1050 
1051 ! call printmat(jacinv,mcrd,mcrd,7,' jinv')
1052 
1053 ! PRINT*, 'detj =', detj
1054  IF(detj.NE.zero) error = .false.
1055 
1056  RETURN
1057 END SUBROUTINE get_jacobien
1058 
1059 !-----------------------------------------------------------------TENSORMUL
1060 SUBROUTINE tensormul(a,lda,n,achar,b,ldb,m,bchar,c,ldc,alpha,beta)
1061 ! *--------------------------------------------------------------------*
1062 ! | Matrix multiplication routine |
1063 ! | |
1064 ! | computes one of the following : |
1065 ! | C = alpha * A B + beta C |
1066 ! | C = alpha * Trans(A) B + beta C |
1067 ! | C = alpha * A Trans(B) + beta C |
1068 ! | C = alpha * Trans(A) Trans(B) + beta C |
1069 ! | depending the character input |
1070 ! | aform : 'n' 'N' 't' 'T' |
1071 ! | bform : 'n' 'N' 't' 'T' |
1072 ! | where 'n' stands for normal 't' stands for transpose |
1073 ! | |
1074 ! | Note that: |
1075 ! | |
1076 ! | INPUTs----lda : leading dimension of A |
1077 ! | n : the other dimension of A |
1078 ! | ldb : leading dimension of B |
1079 ! | m : the other dimension of B |
1080 ! | ldc : leading dim. of C (must be large enough) |
1081 ! | alpha : a scalar multiplier to initialize C |
1082 ! | |
1083 ! | OUTPUTs-----C : is the product matrix |
1084 ! | |
1085 ! *--------------------------------------------------------------------*
1086 
1087 
1088  IMPLICIT REAL*8 (a-h,o-z)
1089 
1090 ! Argument variables
1091 
1092  dimension a(lda,*), b(ldb,*), c(ldc,*)
1093  CHARACTER*1 achar,bchar
1094 
1095 ! Local variables
1096 
1097  parameter(max = 30)
1098  dimension anew(max,max), bnew(max,max)
1099  INTEGER crow, ccol
1100 
1101 ! Check character switches and setup indices
1102 
1103  IF(achar.EQ.'n'.OR.achar.EQ.'N') THEN
1104  ia = lda
1105  ja = n
1106  crow = lda
1107  DO i=1,ia
1108  DO j=1,ja
1109  anew(i,j) = a(i,j)
1110  END DO
1111  END DO
1112  ELSEIF(achar.EQ.'t'.OR.achar.EQ.'T') THEN
1113  ia = n
1114  ja = lda
1115  crow = n
1116  DO i=1,ia
1117  DO j=1,ja
1118  anew(i,j) = a(j,i)
1119  END DO
1120  END DO
1121  ELSE
1122  WRITE(7,*) ' >> unrecognized character, ',achar
1123  stop
1124  END IF
1125 
1126 
1127  IF(bchar.EQ.'n'.OR.bchar.EQ.'N') THEN
1128  ib = ldb
1129  jb = m
1130  ccol = m
1131  DO i=1,ib
1132  DO j=1,jb
1133  bnew(i,j) = b(i,j)
1134  END DO
1135  END DO
1136  ELSEIF(bchar.EQ.'t'.OR.bchar.EQ.'T') THEN
1137  ib = m
1138  jb = ldb
1139  ccol = ldb
1140  DO i=1,ib
1141  DO j=1,jb
1142  bnew(i,j) = b(j,i)
1143  END DO
1144  END DO
1145  ELSE
1146  WRITE(7,*) ' >> unrecognized character, ',bchar
1147  stop
1148  END IF
1149 
1150  IF (ja.NE.ib) THEN
1151  WRITE(7,*) ' >> incompatible matrices'
1152  stop
1153  END IF
1154 
1155 ! Compute C
1156 
1157  c(1:crow,1:ccol) = beta * c(1:crow,1:ccol)
1158  DO i=1, crow
1159  DO j=1,ccol
1160  DO k=1, ja
1161  c(i,j) = c(i,j) + alpha * anew(i,k)*bnew(k,j)
1162  END DO
1163  END DO
1164  END DO
1165 
1166  RETURN
1167  END SUBROUTINE tensormul
1168 
1169 !!$!-----------------------------------------------------------------TENSORMUL
1170 !!$ SUBROUTINE TENSORMUL3(a,lda,n,achar,b,ldb,m,bchar,c,ldc,nn,alpha,beta)
1171 !!$! *--------------------------------------------------------------------*
1172 !!$! | Matrix multiplication routine |
1173 !!$! | |
1174 !!$! | computes one of the following : |
1175 !!$! | C = alpha * A B + beta C |
1176 !!$! | C = alpha * Trans(A) B + beta C |
1177 !!$! | C = alpha * A Trans(B) + beta C |
1178 !!$! | C = alpha * Trans(A) Trans(B) + beta C |
1179 !!$! | depending the character input |
1180 !!$! | aform : 'n' 'N' 't' 'T' |
1181 !!$! | bform : 'n' 'N' 't' 'T' |
1182 !!$! | where 'n' stands for normal 't' stands for transpose |
1183 !!$! | |
1184 !!$! | Note that: |
1185 !!$! | |
1186 !!$! | INPUTs----lda : leading dimension of A |
1187 !!$! | n : the other dimension of A |
1188 !!$! | ldb : leading dimension of B |
1189 !!$! | m : the other dimension of B |
1190 !!$! | ldc : leading dim. of C (must be large enough) |
1191 !!$! | alpha : a scalar multiplier to initialize C |
1192 !!$! | |
1193 !!$! | OUTPUTs-----C : is the product matrix |
1194 !!$! | |
1195 !!$! *--------------------------------------------------------------------*
1196 !!$
1197 !!$
1198 !!$ IMPLICIT REAL*8(a-h,o-z)
1199 !!$
1200 !!$! Argument variables
1201 !!$
1202 !!$ DIMENSION a(lda,n), b(ldb,m), c(ldc,nn)
1203 !!$ CHARACTER*1 achar,bchar
1204 !!$
1205 !!$! Local variables
1206 !!$
1207 !!$ PARAMETER (max = 30)
1208 !!$ DIMENSION anew(max,max), bnew(max,max)
1209 !!$ INTEGER crow, ccol
1210 !!$
1211 !!$! Check character switches and setup indices
1212 !!$
1213 !!$ IF(achar.EQ.'n'.OR.achar.EQ.'N') THEN
1214 !!$ ia = lda
1215 !!$ ja = n
1216 !!$ crow = lda
1217 !!$ DO i=1,ia
1218 !!$ DO j=1,ja
1219 !!$ anew(i,j) = a(i,j)
1220 !!$ END DO
1221 !!$ END DO
1222 !!$ ELSEIF(achar.EQ.'t'.OR.achar.EQ.'T') THEN
1223 !!$ ia = n
1224 !!$ ja = lda
1225 !!$ crow = n
1226 !!$ DO i=1,ia
1227 !!$ DO j=1,ja
1228 !!$ anew(i,j) = a(j,i)
1229 !!$ END DO
1230 !!$ END DO
1231 !!$ ELSE
1232 !!$ WRITE(7,*) ' >> unrecognized character, ',achar
1233 !!$ STOP
1234 !!$ END IF
1235 !!$
1236 !!$
1237 !!$ IF(bchar.EQ.'n'.OR.bchar.EQ.'N') THEN
1238 !!$ ib = ldb
1239 !!$ jb = m
1240 !!$ ccol = m
1241 !!$ DO i=1,ib
1242 !!$ DO j=1,jb
1243 !!$ bnew(i,j) = b(i,j)
1244 !!$ END DO
1245 !!$ END DO
1246 !!$ ELSEIF(bchar.EQ.'t'.OR.bchar.EQ.'T') THEN
1247 !!$ ib = m
1248 !!$ jb = ldb
1249 !!$ ccol = ldb
1250 !!$ DO i=1,ib
1251 !!$ DO j=1,jb
1252 !!$ bnew(i,j) = b(j,i)
1253 !!$ END DO
1254 !!$ END DO
1255 !!$ ELSE
1256 !!$ WRITE(7,*) ' >> unrecognized character, ',bchar
1257 !!$ STOP
1258 !!$ END IF
1259 !!$
1260 !!$ IF (ja.NE.ib) THEN
1261 !!$ WRITE(7,*) ' >> incompatible matrices'
1262 !!$ STOP
1263 !!$ END IF
1264 !!$
1265 !!$! Compute C
1266 !!$
1267 !!$ c(1:crow,1:ccol) = beta * c(1:crow,1:ccol)
1268 !!$ DO i=1, crow
1269 !!$ DO j=1,ccol
1270 !!$ DO k=1, ja
1271 !!$ c(i,j) = c(i,j) + alpha * anew(i,k)*bnew(k,j)
1272 !!$ END DO
1273 !!$ END DO
1274 !!$ END DO
1275 !!$
1276 !!$ RETURN
1277 !!$ END SUBROUTINE TENSORMUL3
1278 !!$
1279 !!$!-----------------------------------------------------------------TENSORMUL3
1280 !!$SUBROUTINE TENSORMULSQR(a,lda,n,achar,b,ldb,m,bchar,c,ldc,alpha,beta)
1281 !!$! *--------------------------------------------------------------------*
1282 !!$! | Matrix multiplication routine |
1283 !!$! | |
1284 !!$! | computes one of the following : |
1285 !!$! | C = alpha * A B + beta C |
1286 !!$! | C = alpha * Trans(A) B + beta C |
1287 !!$! | C = alpha * A Trans(B) + beta C |
1288 !!$! | C = alpha * Trans(A) Trans(B) + beta C |
1289 !!$! | depending the character input |
1290 !!$! | aform : 'n' 'N' 't' 'T' |
1291 !!$! | bform : 'n' 'N' 't' 'T' |
1292 !!$! | where 'n' stands for normal 't' stands for transpose |
1293 !!$! | |
1294 !!$! | Note that: |
1295 !!$! | |
1296 !!$! | INPUTs----lda : leading dimension of A |
1297 !!$! | n : the other dimension of A |
1298 !!$! | ldb : leading dimension of B |
1299 !!$! | m : the other dimension of B |
1300 !!$! | ldc : leading dim. of C (must be large enough) |
1301 !!$! | alpha : a scalar multiplier to initialize C |
1302 !!$! | |
1303 !!$! | OUTPUTs-----C : is the product matrix |
1304 !!$! | |
1305 !!$! *--------------------------------------------------------------------*
1306 !!$
1307 !!$
1308 !!$ IMPLICIT REAL*8 (a-h,o-z)
1309 !!$
1310 !!$! Argument variables
1311 !!$
1312 !!$ DIMENSION a(9,9), b(9,9), c(9,9)
1313 !!$ CHARACTER*1 achar,bchar
1314 !!$
1315 !!$! Local variables
1316 !!$
1317 !!$ PARAMETER (max = 30)
1318 !!$ DIMENSION anew(max,max), bnew(max,max)
1319 !!$ INTEGER crow, ccol
1320 !!$
1321 !!$! Check character switches and setup indices
1322 !!$
1323 !!$ IF(achar.EQ.'n'.OR.achar.EQ.'N') THEN
1324 !!$ ia = lda
1325 !!$ ja = n
1326 !!$ crow = lda
1327 !!$ DO i=1,ia
1328 !!$ DO j=1,ja
1329 !!$ anew(i,j) = a(i,j)
1330 !!$ END DO
1331 !!$ END DO
1332 !!$ ELSEIF(achar.EQ.'t'.OR.achar.EQ.'T') THEN
1333 !!$ ia = n
1334 !!$ ja = lda
1335 !!$ crow = n
1336 !!$ DO i=1,ia
1337 !!$ DO j=1,ja
1338 !!$ anew(i,j) = a(j,i)
1339 !!$ END DO
1340 !!$ END DO
1341 !!$ ELSE
1342 !!$ WRITE(7,*) ' >> unrecognized character, ',achar
1343 !!$ STOP
1344 !!$ END IF
1345 !!$
1346 !!$
1347 !!$ IF(bchar.EQ.'n'.OR.bchar.EQ.'N') THEN
1348 !!$ ib = ldb
1349 !!$ jb = m
1350 !!$ ccol = m
1351 !!$ DO i=1,ib
1352 !!$ DO j=1,jb
1353 !!$ bnew(i,j) = b(i,j)
1354 !!$ END DO
1355 !!$ END DO
1356 !!$ ELSEIF(bchar.EQ.'t'.OR.bchar.EQ.'T') THEN
1357 !!$ ib = m
1358 !!$ jb = ldb
1359 !!$ ccol = ldb
1360 !!$ DO i=1,ib
1361 !!$ DO j=1,jb
1362 !!$ bnew(i,j) = b(j,i)
1363 !!$ END DO
1364 !!$ END DO
1365 !!$ ELSE
1366 !!$ WRITE(7,*) ' >> unrecognized character, ',bchar
1367 !!$ STOP
1368 !!$ END IF
1369 !!$
1370 !!$ IF (ja.NE.ib) THEN
1371 !!$ WRITE(7,*) ' >> incompatible matrices'
1372 !!$ STOP
1373 !!$ END IF
1374 !!$
1375 !!$! Compute C
1376 !!$
1377 !!$ c(1:crow,1:ccol) = beta * c(1:crow,1:ccol)
1378 !!$ DO i=1, crow
1379 !!$ DO j=1,ccol
1380 !!$ DO k=1, ja
1381 !!$ c(i,j) = c(i,j) + alpha * anew(i,k)*bnew(k,j)
1382 !!$ END DO
1383 !!$ END DO
1384 !!$ END DO
1385 !!$
1386 !!$ RETURN
1387 !!$ END SUBROUTINE TENSORMULSQR
1388 !!$
1389 
1390 !-----------------------------------------------------------------INVERT
1391 
1392 SUBROUTINE invert(a,nrow,det)
1393 
1394 ! *--------------------------------------------------------------------*
1395 ! | |
1396 ! | Inverts a matrix by gauss jordan elimination and computes the |
1397 ! | the determinant. The matrix may be symmetric or nonsymmetric. |
1398 ! | |
1399 ! | <a> : a square, double precision matrix |
1400 ! | <nrow> : dimensioned number of rows for 'a'. also the actual |
1401 ! | sizes for inversion computations |
1402 ! | <det> : determinant of the matrix ( zero if matrix is singular)|
1403 ! | |
1404 ! | NOTE : the matrix is overwritten by the inverse. |
1405 ! | |
1406 ! *--------------------------------------------------------------------*
1407  IMPLICIT NONE
1408 
1409 ! Argument variables
1410 
1411  INTEGER nrow
1412  REAL*8, DIMENSION(1:nrow,1:nrow) :: a
1413  REAL*8 :: det
1414 
1415 ! Local variables
1416 
1417  INTEGER ncol, k, j, i
1418 
1419 ! Invert matrix
1420 
1421  det = 1.
1422  ncol = nrow
1423 
1424  DO k = 1, nrow
1425 
1426  IF( a(k,k) .EQ. 0. ) THEN ! Check for singular matrix
1427  det = 0.
1428  RETURN
1429  END IF
1430 
1431  det = det * a(k, k)
1432 
1433  DO j = 1, ncol
1434  IF( j .NE. k ) a(k,j) = a(k,j) / a(k,k)
1435  END DO
1436 
1437  a(k,k) = 1. / a(k,k)
1438 
1439  DO i = 1, nrow
1440  IF( i .EQ. k) cycle
1441  DO j = 1, ncol
1442  IF( j .NE. k) a(i,j) = a(i,j) - a(i,k) * a(k,j)
1443  END DO
1444  a(i,k) = -a(i,k) * a(k,k)
1445  END DO
1446 
1447  END DO
1448 
1449  RETURN
1450 END SUBROUTINE invert
1451 
1452 !--------------------------------------------------------------INVERT3x3
1453 
1454 SUBROUTINE invert3x3(a,det)
1455 ! *--------------------------------------------------------------------*
1456 ! | |
1457 ! | Inverts 3 by 3 matrix and computes the determinant |
1458 ! | The matrix may be symmetric or nonsymmetric. |
1459 ! | |
1460 ! | <a> : 3x3 input matrix |
1461 ! | <det> : determinant of the matrix ( zero if matrix is singular)|
1462 ! | |
1463 ! | NOTE : the matrix is overwritten by the inverse. |
1464 ! | |
1465 ! *--------------------------------------------------------------------*
1466  IMPLICIT NONE
1467 
1468 ! Argument variables
1469 
1470  DOUBLE PRECISION a(3,3), det
1471 
1472 ! Local variables
1473 
1474  DOUBLE PRECISION acopy(3,3), zero
1475 
1476 ! Data statements
1477 
1478  DATA zero /0.000000000000000/
1479 
1480 ! Invert matrix
1481 
1482  acopy = a
1483 
1484  det = -acopy(1,3)*acopy(2,2)*acopy(3,1) &
1485  +acopy(1,2)*acopy(2,3)*acopy(3,1) &
1486  +acopy(1,3)*acopy(2,1)*acopy(3,2) &
1487  -acopy(1,1)*acopy(2,3)*acopy(3,2) &
1488  -acopy(1,2)*acopy(2,1)*acopy(3,3) &
1489  +acopy(1,1)*acopy(2,2)*acopy(3,3)
1490 
1491  IF(det.EQ.zero) THEN ! Check for singular matrix
1492  a(1:3,1:3) = zero
1493  RETURN
1494  END IF
1495 
1496  a(1,1) = -acopy(2,3)*acopy(3,2) + acopy(2,2)*acopy(3,3)
1497  a(2,1) = acopy(2,3)*acopy(3,1) - acopy(2,1)*acopy(3,3)
1498  a(3,1) = -acopy(2,2)*acopy(3,1) + acopy(2,1)*acopy(3,2)
1499  a(1,2) = acopy(1,3)*acopy(3,2) - acopy(1,2)*acopy(3,3)
1500  a(2,2) = -acopy(1,3)*acopy(3,1) + acopy(1,1)*acopy(3,3)
1501  a(3,2) = acopy(1,2)*acopy(3,1) - acopy(1,1)*acopy(3,2)
1502  a(1,3) = -acopy(1,3)*acopy(2,2) + acopy(1,2)*acopy(2,3)
1503  a(2,3) = acopy(1,3)*acopy(2,1) - acopy(1,1)*acopy(2,3)
1504  a(3,3) = -acopy(1,2)*acopy(2,1) + acopy(1,1)*acopy(2,2)
1505 
1506  a = a / det
1507 
1508  RETURN
1509  END SUBROUTINE invert3x3
1510 !----------------------------------------------------------GET_MIXED_MAP
1511  SUBROUTINE get_mixed_map(e_mixed,r,igpt)
1512 ! *--------------------------------------------------------------------*
1513 ! | |
1514 ! | Returns the linear mapping matrix <e_mixed> |
1515 ! | evaluated at a given isoparametric coordinate <r> |
1516 ! | |
1517 ! | <r> : 3 x ? vector of coordinates of the isoparametric |
1518 ! | point (usually an integration point) |
1519 ! | <e_mixed> : 9 x 12 matrix that interpolates the mixed field |
1520 ! | <igpt> : no. of the particular gauss point |
1521 ! | |
1522 ! | NOTE : 9 is no. of stress/strain components |
1523 ! | 12 is no. of mixed field parameters |
1524 ! | |
1525 ! *--------------------------------------------------------------------*
1526 
1527  IMPLICIT Real*8 (a-h, o-z)
1528 ! include 'ABA_PARAM.INC'
1529 
1530 ! Argument variables
1531 
1532  dimension e_mixed(9,12), r(3,*)
1533 
1534 ! Data statements
1535 
1536  DATA zero /0.000000000000000/
1537 
1538 ! Initialize variables
1539 
1540  r1 = r(1,igpt)
1541  r2 = r(2,igpt)
1542  r3 = r(3,igpt)
1543 
1544  e_mixed = zero
1545 
1546 ! Evaluate <e_mixed>
1547 
1548  e_mixed(1,1) = r2
1549  e_mixed(1,2) = r3
1550  e_mixed(1,3) = r2 * r3
1551 
1552  e_mixed(2,10) = r3
1553  e_mixed(3,12) = r2
1554  e_mixed(4,10) = r3
1555 
1556  e_mixed(5,4) = r1
1557  e_mixed(5,5) = r3
1558  e_mixed(5,6) = r1 * r3
1559 
1560  e_mixed(6,11) = r1
1561  e_mixed(7,12) = r2
1562  e_mixed(8,11) = r1
1563 
1564  e_mixed(9,7) = r1
1565  e_mixed(9,8) = r2
1566  e_mixed(9,9) = r1 * r2
1567 
1568  RETURN
1569  END SUBROUTINE get_mixed_map
1570 
1571 !-------------------------------------------------------GET_ENHANCED_MAP
1572  SUBROUTINE get_enhanced_map(e_enhanced,r,igpt)
1573 ! *--------------------------------------------------------------------*
1574 ! | |
1575 ! | Returns the linear mapping matrix <e_enhanced> |
1576 ! | evaluated at a given isoparametric coordinate <r> |
1577 ! | |
1578 ! | <r> : 3 x ? vector of coordinates of the isoparametric |
1579 ! | point (usually an integration point) |
1580 ! | <e_enhanced> : 9 x 9 matrix that interpolates the enhanced |
1581 ! | field |
1582 ! | <igpt> : no. of the particular gauss point |
1583 ! | |
1584 ! | NOTE : 9 is no. of stress/strain components |
1585 ! | 9 is no. of enhanced field parameters |
1586 ! | |
1587 ! *--------------------------------------------------------------------*
1588 
1589  IMPLICIT real*8 (a-h, o-z)
1590 ! include 'ABA_PARAM.INC'
1591 
1592 ! Argument variables
1593 
1594  dimension e_enhanced(9,9), r(3,*)
1595 
1596 ! Data statements
1597 
1598  DATA zero /0.000000000000000/
1599 
1600 ! Initialize variables
1601 
1602  r1 = r(1,igpt)
1603  r2 = r(2,igpt)
1604  r3 = r(3,igpt)
1605 
1606  e_enhanced = zero
1607 
1608 ! Evaluate <e_enhanced>
1609 
1610  e_enhanced(1,1) = r1
1611  e_enhanced(1,2) = r1 * r2
1612  e_enhanced(1,3) = r1 * r3
1613 
1614  e_enhanced(5,4) = r2
1615  e_enhanced(5,5) = r2 * r3
1616  e_enhanced(5,6) = r2 * r1
1617 
1618  e_enhanced(9,7) = r3
1619  e_enhanced(9,8) = r3 * r2
1620  e_enhanced(9,9) = r3 * r1
1621 
1622  RETURN
1623  END SUBROUTINE get_enhanced_map
1624 !------------------------------------------------------KRONECKER_PRODUCT
1625  SUBROUTINE kronecker_product(a,lda,nda,b,ldb,ndb,c,ldc)
1626 ! *--------------------------------------------------------------------*
1627 ! | This subroutine takes the Kronecker Product of the two matrices |
1628 ! | |
1629 ! | C = A <kron_prod> B |
1630 ! | |
1631 ! | INPUTs----lda : leading dimension of A |
1632 ! | nda : the other dimension of A |
1633 ! | ldb : leading dimension of B |
1634 ! | ndb : the other dimension of B |
1635 ! | ldc : leading dim. of C ( must be >= lda*ldb ) |
1636 ! | |
1637 ! | OUTPUTs-----C : is the product matrix |
1638 ! | |
1639 ! *--------------------------------------------------------------------*
1640 
1641 ! include 'ABA_PARAM.INC'
1642  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1643 
1644 ! Argument variables
1645 
1646  dimension a(lda,nda), b(ldb,ndb), c(ldc,*)
1647 
1648 ! Local variables
1649 
1650  INTEGER ia, ja, ist, iend, jst, jend
1651 
1652 ! Compute Kronecker Product
1653 
1654  DO ja=1,nda
1655  DO ia=1,lda
1656 
1657  ist = (ia - 1) * ldb + 1
1658  iend = ist + ldb - 1
1659  jst = (ja - 1) * ndb + 1
1660  jend = jst + ndb - 1
1661 
1662  c(ist:iend,jst:jend) = a(ia,ja) * b(1:ldb,1:ndb)
1663 
1664  END DO
1665  END DO
1666 
1667  RETURN
1668  END SUBROUTINE kronecker_product
1669 
Color::r
unsigned char r() const
Definition: Color.h:68
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
ncol
const int ncol
Definition: ex1.C:95
invert
subroutine invert(a, nrow, det)
Definition: v3d8_me.f90:1392
k
j indices k indices k
Definition: Indexing.h:6
Connectivity::dimension
int dimension() const
Get the dimension of the mesh.
Definition: Connectivity.h:125
max
Vector_n max(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:354
tensormul
subroutine tensormul(a, lda, n, achar, b, ldb, m, bchar, c, ldc, alpha, beta)
Definition: v3d8_me.f90:1060
Color::b
unsigned char b() const
Definition: Color.h:70
coord
int coord[NPANE][NROW *NCOL][3]
Definition: blastest.C:86
Line_2::c
RT c() const
Definition: Line_2.h:150
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
iend
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE knode iend
Definition: RFLO_CalcGradDummySymm.h:32
get_shape20
subroutine get_shape20(r, n, dn, igpt)
Definition: v3d8_me.f90:924
z
void int int int REAL REAL REAL * z
Definition: write.cpp:76
get_jacobien
subroutine get_jacobien(coords, mcrd, nnode, dn, jac, jacinv, detj, error)
Definition: v3d8_me.f90:1003
i
blockLoc i
Definition: read.cpp:79
get_enhanced_map
subroutine get_enhanced_map(e_enhanced, r, igpt)
Definition: v3d8_me.f90:1572
n
const NT & n
Definition: rational_rotation.h:72
j
j indices j
Definition: Indexing.h:6
jend
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE knode jend
Definition: RFLO_CalcGradDummySymm.h:32
nrow
const int nrow
Definition: ex1.C:95
GEM_DomainBoundary::debug
bool debug(bool s=true)
Definition: GEM.H:193
v3d8_me
subroutine v3d8_me(coor, MatType, ElConnVol, Rnet, disp, dmat, S11, S22, S33, S12, S23, S13, numnp, nstart, nend, NumEL, NumMatType, Aenh, enhanced_map, mixed_map)
Definition: v3d8_me.f90:53
Color::alpha
unsigned char alpha() const
Definition: Color.h:75
get_mixed_map
subroutine get_mixed_map(e_mixed, r, igpt)
Definition: v3d8_me.f90:1511
Line_2::a
RT a() const
Definition: Line_2.h:140
get_shape
subroutine get_shape(r, n, dn, igpt)
Definition: v3d8_me.f90:846