Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
v3d10_thermal.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 v3d10_thermal(NumEl, NumNP, ElConnVol, Coor, Kappa,&
54  rnet, t, rho,cp, matcstet,numat_vol,meshvel,elemstart, elemend)
55 
56  IMPLICIT NONE
57 
58 !!****f* Rocfrac/Source/v3d10_thermal
59 !!
60 !! NAME
61 !! v3d10_thermal.f90
62 !!
63 !! FUNCTION
64 !! Calculates the internal force vector due to heat
65 !! transfer for a 3D 10-node tet element.
66 !!
67 !! ALE
68 !! {Rnet} = [K]{T} + [K] {T}.
69 !!
70 !! ALE
71 !! where [K] is the element matrices
72 !! [K] is the conductivity matrix
73 !! {T} is the temperature
74 !!
75 !!
76 !! USED BY
77 !! RocFracMain
78 !!
79 !! INPUTS
80 !! NumEl - Number of Elements
81 !! NumNP - Number of Nodes
82 !! ElConnVol - Element Connectivity Array
83 !! Coor - Current Mesh coordinates
84 !! Kappa - thermal conductivity
85 !! Rnet - Accumulated "force" vector (RnetHT)
86 !! T - Temperature
87 !! RhoCp - Density * Cp
88 !! MeshVel - Mesh motion velocity
89 !!
90 !! OUTPUT
91 !!
92 !! Accumulated "force" vector with the addition of
93 !! the internal "force" vector (RnetHT)
94 !!
95 !!***
96 
97  INTEGER :: numnp, numel,numat_vol
98  integer :: elemstart, elemend
99  REAL*8, DIMENSION(1:numat_vol) :: kappa
100  REAL*8, DIMENSION(1:NumNP) :: t
101  REAL*8, DIMENSION(1:NumNP) :: rnet
102  INTEGER , DIMENSION(1:10,1:NumEl) :: elconnvol
103  REAL*8, DIMENSION(1:3,1:NumNp) :: coor
104  REAL*8, DIMENSION(1:NumNp*3) :: meshvel
105  REAL*8, DIMENSION(1:numat_vol) :: rho
106  REAL*8, DIMENSION(1:numat_vol) :: cp
107  integer, dimension(1:NumEl) :: matcstet
108 
109  INTEGER :: i, imat, igpt,j
110  INTEGER :: n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
111 
112 
113  REAL*8, DIMENSION(1:10,1:10) :: kloc
114  REAL*8, DIMENSION(1:10) :: tloc, rinloc
115 
116 
117 !-- Coordinate subtractions
118  REAL*8 :: x14, x24, x34, y14, y24, y34, z14, z24, z34
119 !-- Coordinate subtractions: These are to speed up B calculation
120  REAL*8 :: x12, x13, y12, y13, z12, z13
121 
122  REAL*8 :: c11, c21, c31
123 !-- 6*volume and inverse of 6*volume
124  REAL*8 :: vx6, vx6inv, vol
125 
126 !-- coordinate holding variable
127  REAL*8 :: x1,x2,x3,x4,x5,x6,x7,x8,x9,x10
128  REAL*8 :: y1,y2,y3,y4,y5,y6,y7,y8,y9,y10
129  REAL*8 :: z1,z2,z3,z4,z5,z6,z7,z8,z9,z10
130  REAL*8 :: aux1,aux2,aux3,aux4,aux5,aux6,aux7,aux8,aux9,aux10,aux11,aux12
131  REAL*8 :: g1, g2, g3, g4
132  REAL*8 :: xn1, xn2, xn3, xn4
133  INTEGER :: k1n1,k1n2,k1n3,k1n4,k1n5,k1n6,k1n7,k1n8,k1n9,k1n10
134  INTEGER :: k2n1,k2n2,k2n3,k2n4,k2n5,k2n6,k2n7,k2n8,k2n9,k2n10
135  INTEGER :: k3n1,k3n2,k3n3,k3n4,k3n5,k3n6,k3n7,k3n8,k3n9,k3n10
136 !-- partial internal force
137  REAL*8 :: r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15,r16,r17,r18
138 
139  REAL*8,DIMENSION(1:4,1:4) :: gaussintpt = reshape( &
140  (/0.58541020d0,0.13819660d0,0.13819660d0,0.13819660d0, &
141  0.13819660d0,0.58541020d0,0.13819660d0,0.13819660d0, &
142  0.13819660d0,0.13819660d0,0.58541020d0,0.13819660d0, &
143  0.13819660d0,0.13819660d0,0.13819660d0,0.58541020d0/),(/4,4/) )
144 
145  real*8, dimension(1:3,1:10) :: b
146  real*8, dimension(1:10,1:3) :: bt
147  real*8, dimension(1:10) :: shapfnct
148  real*8, dimension(1:3) :: meshvelnd
149 
150  REAL*8, DIMENSION(1:3,1:3) :: kappamatrx = reshape( &
151  (/1.000000000000000,0.000000000000000,0.000000000000000, &
152  0.000000000000000,1.000000000000000,0.000000000000000, &
153  0.000000000000000,0.000000000000000,1.000000000000000/),(/3,3/) )
154 
155  DO i = elemstart, elemend
156  imat = matcstet(i)
157 
158  n1 = elconnvol(1,i)
159  n2 = elconnvol(2,i)
160  n3 = elconnvol(3,i)
161  n4 = elconnvol(4,i)
162  n5 = elconnvol(5,i)
163  n6 = elconnvol(6,i)
164  n7 = elconnvol(7,i)
165  n8 = elconnvol(8,i)
166  n9 = elconnvol(9,i)
167  n10 = elconnvol(10,i)
168 
169 
170  k3n1 = 3*n1
171  k3n2 = 3*n2
172  k3n3 = 3*n3
173  k3n4 = 3*n4
174  k3n5 = 3*n5
175  k3n6 = 3*n6
176  k3n7 = 3*n7
177  k3n8 = 3*n8
178  k3n9 = 3*n9
179  k3n10 = 3*n10
180 
181  k2n1 = k3n1 - 1
182  k2n2 = k3n2 - 1
183  k2n3 = k3n3 - 1
184  k2n4 = k3n4 - 1
185  k2n5 = k3n5 - 1
186  k2n6 = k3n6 - 1
187  k2n7 = k3n7 - 1
188  k2n8 = k3n8 - 1
189  k2n9 = k3n9 - 1
190  k2n10 = k3n10 - 1
191 
192  k1n1 = k3n1 - 2
193  k1n2 = k3n2 - 2
194  k1n3 = k3n3 - 2
195  k1n4 = k3n4 - 2
196  k1n5 = k3n5 - 2
197  k1n6 = k3n6 - 2
198  k1n7 = k3n7 - 2
199  k1n8 = k3n8 - 2
200  k1n9 = k3n9 - 2
201  k1n10 = k3n10 - 2
202 
203  x1 = coor(1,n1)
204  x2 = coor(1,n2)
205  x3 = coor(1,n3)
206  x4 = coor(1,n4)
207  x5 = coor(1,n5)
208  x6 = coor(1,n6)
209  x7 = coor(1,n7)
210  x8 = coor(1,n8)
211  x9 = coor(1,n9)
212  x10 = coor(1,n10)
213 
214  y1 = coor(2,n1)
215  y2 = coor(2,n2)
216  y3 = coor(2,n3)
217  y4 = coor(2,n4)
218  y5 = coor(2,n5)
219  y6 = coor(2,n6)
220  y7 = coor(2,n7)
221  y8 = coor(2,n8)
222  y9 = coor(2,n9)
223  y10 = coor(2,n10)
224 
225  z1 = coor(3,n1)
226  z2 = coor(3,n2)
227  z3 = coor(3,n3)
228  z4 = coor(3,n4)
229  z5 = coor(3,n5)
230  z6 = coor(3,n6)
231  z7 = coor(3,n7)
232  z8 = coor(3,n8)
233  z9 = coor(3,n9)
234  z10 = coor(3,n10)
235 
236  x14 = x1 - x4
237  x24 = x2 - x4
238  x34 = x3 - x4
239  y14 = y1 - y4
240  y24 = y2 - y4
241  y34 = y3 - y4
242  z14 = z1 - z4
243  z24 = z2 - z4
244  z34 = z3 - z4
245 
246  x12 = x1 - x2 ! not used in vol. calc
247  x13 = x1 - x3 ! not used in vol. calc
248  y12 = y1 - y2 ! not used in vol. calc
249  y13 = y1 - y3 ! not used in vol. calc
250  z12 = z1 - z2 ! not used in vol. calc
251  z13 = z1 - z3 ! not used in vol. calc
252 
253  c11 = y24*z34 - z24*y34
254  c21 = -( x24*z34 - z24*x34 )
255  c31 = x24*y34 - y24*x34
256 
257  vx6 = -( x14*c11 + y14*c21 + z14*c31 )
258 
259  vx6inv = 1.d0 / vx6
260 
261  aux1 = -(y3*z4 - y4*z3 - y2*z4 + y2*z3 + z2*y4 - z2*y3)
262  aux2 = (x3*z4 - x4*z3 - x2*z4 + x2*z3 + z2*x4 - z2*x3)
263  aux3 = -(x3*y4 - x4*y3 - x2*y4 + x2*y3 + y2*x4 - y2*x3)
264  aux4 = (y3*z4 - y4*z3 - y1*z4 + y1*z3 + z1*y4 - z1*y3)
265  aux5 = -(x3*z4 - x4*z3 - x1*z4 + x1*z3 + z1*x4 - z1*x3)
266  aux6 = (x3*y4 - x4*y3 - x1*y4 + x1*y3 + y1*x4 - y1*x3)
267  aux7 = -(y2*z4 - z2*y4 - y1*z4 + y1*z2 + z1*y4 - z1*y2)
268  aux8 = (x2*z4 - z2*x4 - x1*z4 + x1*z2 + z1*x4 - z1*x2)
269  aux9 = -(x2*y4 - y2*x4 - x1*y4 + x1*y2 + y1*x4 - y1*x2)
270  aux10 = (y2*z3 - z2*y3 - y1*z3 + y1*z2 + z1*y3 - z1*y2)
271  aux11 =-(x2*z3 - z2*x3 - x1*z3 + x1*z2 + z1*x3 - z1*x2)
272  aux12 = (x2*y3 - y2*x3 - x1*y3 + x1*y2 + y1*x3 - y1*x2)
273 
274  r1 = 0.
275  r2 = 0.
276  r3 = 0.
277  r4 = 0.
278  r5 = 0.
279  r6 = 0.
280  r7 = 0.
281  r8 = 0.
282  r9 = 0.
283  r10 = 0.
284 
285  tloc(1) = t(n1)
286  tloc(2) = t(n2)
287  tloc(3) = t(n3)
288  tloc(4) = t(n4)
289  tloc(5) = t(n5)
290  tloc(6) = t(n6)
291  tloc(7) = t(n7)
292  tloc(8) = t(n8)
293  tloc(9) = t(n9)
294  tloc(10) = t(n10)
295 
296  kloc = 0.
297 
298  DO igpt = 1, 4
299 
300  g1 = gaussintpt(igpt,1)
301  g2 = gaussintpt(igpt,2)
302  g3 = gaussintpt(igpt,3)
303  g4 = gaussintpt(igpt,4)
304 
305  shapfnct(1) = g1*(2.*g1-1.)
306  shapfnct(2) = g2*(2.*g2-1.)
307  shapfnct(3) = g3*(2.*g3-1.)
308  shapfnct(4) = g4*(2.*g4-1.)
309  shapfnct(5) = 4.*g1*g2
310  shapfnct(6) = 4.*g2*g3
311  shapfnct(7) = 4.*g3*g1
312  shapfnct(8) = 4.*g1*g4
313  shapfnct(9) = 4.*g2*g4
314  shapfnct(10) = 4.*g3*g4
315 
316 
317  xn1 = (4.d0*g1-1.d0) ! derivative of shape function
318  xn2 = (4.d0*g2-1.d0) ! dN_i/dzeta_i
319  xn3 = (4.d0*g3-1.d0)
320  xn4 = (4.d0*g4-1.d0)
321 ! xN5 = 4.d0*g1*g2
322 ! xN6 = 4.d0*g2*g3
323 ! xN7 = 4.d0*g3*g1
324 ! xN8 = 4.d0*g1*g4
325 ! xN9 = 4.d0*g2*g4
326 ! xN10= 4.d0*g3*g4
327 
328  b(1,1) = aux1*xn1
329  b(2,1) = aux2*xn1
330  b(3,1) = aux3*xn1
331  b(1,2) = aux4*xn2
332  b(2,2) = aux5*xn2
333  b(3,2) = aux6*xn2
334  b(1,3) = aux7*xn3
335  b(2,3) = aux8*xn3
336  b(3,3) = aux9*xn3
337  b(1,4) = aux10*xn4
338  b(2,4) = aux11*xn4
339  b(3,4) = aux12*xn4
340 
341  b(1,5) = 4.d0*(g2*aux1 + g1*aux4)
342  b(2,5) = 4.d0*(g2*aux2 + g1*aux5)
343  b(3,5) = 4.d0*(g2*aux3 + g1*aux6)
344 
345  b(1,6) = 4.d0*(g3*aux4 + g2*aux7)
346  b(2,6) = 4.d0*(g3*aux5 + g2*aux8)
347  b(3,6) = 4.d0*(g3*aux6 + g2*aux9)
348 
349  b(1,7) = 4.d0*(g1*aux7 + g3*aux1)
350  b(2,7) = 4.d0*(g1*aux8 + g3*aux2)
351  b(3,7) = 4.d0*(g1*aux9 + g3*aux3)
352 
353  b(1,8) = 4.d0*(g4*aux1 + g1*aux10)
354  b(2,8) = 4.d0*(g4*aux2 + g1*aux11)
355  b(3,8) = 4.d0*(g4*aux3 + g1*aux12)
356 
357  b(1,9) = 4.d0*(g4*aux4 + g2*aux10)
358  b(2,9) = 4.d0*(g4*aux5 + g2*aux11)
359  b(3,9) = 4.d0*(g4*aux6 + g2*aux12)
360 
361  b(1,10) = 4.d0*(g4*aux7 + g3*aux10)
362  b(2,10) = 4.d0*(g4*aux8 + g3*aux11)
363  b(3,10) = 4.d0*(g4*aux9 + g3*aux12)
364 
365  b = vx6inv*b
366 
367 ! BT = TRANSPOSE(B)
368 ! Compute the local stiffness matrix
369 
370 ! B = MATMUL(KappaMatrx,B)
371 
372  vol = vx6 / 6.d0
373 
374 ! print*,'***',Kloc(9,10)
375  kloc = kloc + kappa(imat)*matmul(matmul(transpose(b),kappamatrx),b) &
376  *0.25d0*vol
377 
378 ! *0.04166666666666667d0*vol ! w(i) * is this correct fix check
379 ! take the product K*T
380 
381 !!$ RinLoc(:) = MATMUL(Kloc,Tloc)
382 
383 !!$ r1 = r1 - RinLoc(1)
384 !!$ r2 = r2 - RinLoc(2)
385 !!$ r3 = r3 - RinLoc(3)
386 !!$ r4 = r4 - RinLoc(4)
387 !!$ r5 = r5 - RinLoc(5)
388 !!$ r6 = r6 - RinLoc(6)
389 !!$ r7 = r7 - RinLoc(7)
390 !!$ r8 = r8 - RinLoc(8)
391 !!$ r9 = r9 - RinLoc(9)
392 !!$ r10 = r10 - RinLoc(10)
393 
394 ! Mesh Motion term
395 
396  meshvelnd(1) = shapfnct(1)*meshvel(k1n1) + &
397  shapfnct(2)*meshvel(k1n2) + &
398  shapfnct(3)*meshvel(k1n3) + &
399  shapfnct(4)*meshvel(k1n4) + &
400  shapfnct(5)*meshvel(k1n5) + &
401  shapfnct(6)*meshvel(k1n6) + &
402  shapfnct(7)*meshvel(k1n7) + &
403  shapfnct(8)*meshvel(k1n8) + &
404  shapfnct(9)*meshvel(k1n9) + &
405  shapfnct(10)*meshvel(k1n10)
406 
407  meshvelnd(2) = shapfnct(1)*meshvel(k2n1) + &
408  shapfnct(2)*meshvel(k2n2) + &
409  shapfnct(3)*meshvel(k2n3) + &
410  shapfnct(4)*meshvel(k2n4) + &
411  shapfnct(5)*meshvel(k2n5) + &
412  shapfnct(6)*meshvel(k2n6) + &
413  shapfnct(7)*meshvel(k2n7) + &
414  shapfnct(8)*meshvel(k2n8) + &
415  shapfnct(9)*meshvel(k2n9) + &
416  shapfnct(10)*meshvel(k2n10)
417 
418  meshvelnd(3) = shapfnct(1)*meshvel(k3n1) + &
419  shapfnct(2)*meshvel(k3n2) + &
420  shapfnct(3)*meshvel(k3n3) + &
421  shapfnct(4)*meshvel(k3n4) + &
422  shapfnct(5)*meshvel(k3n5) + &
423  shapfnct(6)*meshvel(k3n6) + &
424  shapfnct(7)*meshvel(k3n7) + &
425  shapfnct(8)*meshvel(k3n8) + &
426  shapfnct(9)*meshvel(k3n9) + &
427  shapfnct(10)*meshvel(k3n10)
428 
429  bt(1,1) = shapfnct(1)*meshvelnd(1)
430  bt(1,2) = shapfnct(1)*meshvelnd(2)
431  bt(1,3) = shapfnct(1)*meshvelnd(3)
432  bt(2,1) = shapfnct(2)*meshvelnd(1)
433  bt(2,2) = shapfnct(2)*meshvelnd(2)
434  bt(2,3) = shapfnct(2)*meshvelnd(3)
435  bt(3,1) = shapfnct(3)*meshvelnd(1)
436  bt(3,2) = shapfnct(3)*meshvelnd(2)
437  bt(3,3) = shapfnct(3)*meshvelnd(3)
438  bt(4,1) = shapfnct(4)*meshvelnd(1)
439  bt(4,2) = shapfnct(4)*meshvelnd(2)
440  bt(4,3) = shapfnct(4)*meshvelnd(3)
441  bt(5,1) = shapfnct(5)*meshvelnd(1)
442  bt(5,2) = shapfnct(5)*meshvelnd(2)
443  bt(5,3) = shapfnct(5)*meshvelnd(3)
444  bt(6,1) = shapfnct(6)*meshvelnd(1)
445  bt(6,2) = shapfnct(6)*meshvelnd(2)
446  bt(6,3) = shapfnct(6)*meshvelnd(3)
447  bt(7,1) = shapfnct(7)*meshvelnd(1)
448  bt(7,2) = shapfnct(7)*meshvelnd(2)
449  bt(7,3) = shapfnct(7)*meshvelnd(3)
450  bt(8,1) = shapfnct(8)*meshvelnd(1)
451  bt(8,2) = shapfnct(8)*meshvelnd(2)
452  bt(8,3) = shapfnct(8)*meshvelnd(3)
453  bt(9,1) = shapfnct(9)*meshvelnd(1)
454  bt(9,2) = shapfnct(9)*meshvelnd(2)
455  bt(9,3) = shapfnct(9)*meshvelnd(3)
456  bt(10,1) = shapfnct(10)*meshvelnd(1)
457  bt(10,2) = shapfnct(10)*meshvelnd(2)
458  bt(10,3) = shapfnct(10)*meshvelnd(3)
459 
460  kloc = kloc - rho(imat)*cp(imat)*matmul(bt,b)*0.25d0*vol
461 
462 ! take the product K*T
463 
464 !!$! RinLoc(:) = MATMUL(Kloc,Tloc)
465 !!$
466 !!$ r1 = r1 - RinLoc(1)
467 !!$ r2 = r2 - RinLoc(2)
468 !!$ r3 = r3 - RinLoc(3)
469 !!$ r4 = r4 - RinLoc(4)
470 !!$ r5 = r5 - RinLoc(5)
471 !!$ r6 = r6 - RinLoc(6)
472 !!$ r7 = r7 - RinLoc(7)
473 !!$ r8 = r8 - RinLoc(8)
474 !!$ r9 = r9 - RinLoc(9)
475 !!$ r10 = r10 - RinLoc(10)
476  ENDDO
477 
478 
479 !!$ DO j = 1,10
480 !!$ WRITE(*,'(10(x,f10.1))')Kloc(j,1:10)
481 !!$ ENDDO
482  rinloc(:) = matmul(kloc,tloc)
483 
484 !!$ DO j = 1,10
485 !!$ WRITE(*,'(10(x,f10.1))')Kloc(j,1:10)
486 !!$ ENDDO
487 !!$ PRINT*,'R'
488 !!$ WRITE(*,'(10(x,f10.1))') Tloc(1:10)
489 
490 
491  r1 = rinloc(1)
492  r2 = rinloc(2)
493  r3 = rinloc(3)
494  r4 = rinloc(4)
495  r5 = rinloc(5)
496  r6 = rinloc(6)
497  r7 = rinloc(7)
498  r8 = rinloc(8)
499  r9 = rinloc(9)
500  r10 = rinloc(10)
501 
502  rnet(n1) = rnet(n1) - r1 !*0.04166666666666667d0
503  rnet(n2) = rnet(n2) - r2 !*0.04166666666666667d0
504  rnet(n3) = rnet(n3) - r3 !*0.04166666666666667d0
505  rnet(n4) = rnet(n4) - r4 !*0.04166666666666667d0
506  rnet(n5) = rnet(n5) - r5 !*0.04166666666666667d0
507  rnet(n6) = rnet(n6) - r6 !*0.04166666666666667d0
508  rnet(n7) = rnet(n7) - r7 !*0.04166666666666667d0
509  rnet(n8) = rnet(n8) - r8 !*0.04166666666666667d0
510  rnet(n9) = rnet(n9) - r9 !*0.04166666666666667d0
511  rnet(n10) = rnet(n10) - r10 !*0.04166666666666667d0
512 
513 
514  END DO
515 
516  RETURN
517 
518 END SUBROUTINE v3d10_thermal
519 
v3d10_thermal
subroutine v3d10_thermal(NumEl, NumNP, ElConnVol, Coor, Kappa, Rnet, T, Rho, Cp, matcstet, numat_vol, MeshVel, ElemStart, ElemEnd)
Definition: v3d10_thermal.f90:53
Color::b
unsigned char b() const
Definition: Color.h:70
i
blockLoc i
Definition: read.cpp:79
j
j indices j
Definition: Indexing.h:6