xref: /phasta/phSolver/compressible/e3massl.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)
1        subroutine e3massl (bcool,    shape,  WdetJ,   A0,
2     &                      EGmass )
3c
4c----------------------------------------------------------------------
5c This routine calculates the time contribution to the LHS tangent
6c mass matrix.
7c
8c input:
9c  shape   (npro,nshl)        : element shape functions
10c  WdetJ   (npro)               : weighted Jacobian determinant
11c  A0      (npro)               : weighted Jacobian matrix
12c  EGmass  (npro,nedof,nedof)   : partial LHS matrix
13c
14c output:
15c  EGmass  (npro,nedof,nedof)   : partial LHS matrix
16c
17c----------------------------------------------------------------------
18c
19        include "common.h"
20c
21        dimension  A0(npro,nflow,nflow),   shape(npro,nshl),
22     &             WdetJ(npro),          EGmass(npro,nedof,nedof)
23c
24        dimension  fact(npro,5),           temp(npro),
25     &             shpij(npro),          bcool(npro)
26c
27c.... ---------------------------->  LHS  <----------------------------
28c
29        temp = WdetJ * almi/gami/alfi*Dtgl
30        bcool = WdetJ * bcool  ! note that bcool is not used after this
31                               ! (if we got in here at least).
32c
33c.... loop through columns (nodes j) and rows (nodes i)
34c
35        do j = 1, nshl
36           j0 = nflow * (j - 1)
37
38           do i = 1, nshl
39              i0 = nflow * (i - 1)
40c
41c.... get the factor
42c
43              shpij = shape(:,i) * shape(:,j)
44              fact(:,1)  = shpij * (temp + bcool*spongeContinuity)
45              fact(:,2)  = shpij * (temp + bcool*spongeMomentum1)
46              fact(:,3)  = shpij * (temp + bcool*spongeMomentum2)
47              fact(:,4)  = shpij * (temp + bcool*spongeMomentum3)
48              fact(:,5)  = shpij * (temp + bcool*spongeEnergy)
49c
50c.... loop through d.o.f.'s
51c
52              do jdof = 1, nflow
53                 EGmass(:,i0+1,j0+jdof) = EGmass(:,i0+1,j0+jdof)
54     &                                  + fact(:,1) * A0(:,1,jdof)
55                 EGmass(:,i0+2,j0+jdof) = EGmass(:,i0+2,j0+jdof)
56     &                                  + fact(:,2) * A0(:,2,jdof)
57                 EGmass(:,i0+3,j0+jdof) = EGmass(:,i0+3,j0+jdof)
58     &                                  + fact(:,3) * A0(:,3,jdof)
59                 EGmass(:,i0+4,j0+jdof) = EGmass(:,i0+4,j0+jdof)
60     &                                  + fact(:,4) * A0(:,4,jdof)
61                 EGmass(:,i0+5,j0+jdof) = EGmass(:,i0+5,j0+jdof)
62     &                                  + fact(:,5) * A0(:,5,jdof)
63              enddo
64c
65c.... end loop on row nodes
66c
67           enddo
68c
69c.... end loop on column nodes
70c
71        enddo
72c
73c.... return
74c
75        return
76        end
77c
78c
79c
80        subroutine e3masslSclr (shape,  WdetJ,   A0t,
81     &                      EGmasst,srcp)
82c
83c----------------------------------------------------------------------
84c This routine calculates the time contribution to the LHS tangent
85c mass matrix.
86c
87c input:
88c  shape   (npro,nshl)          : element shape functions
89c  WdetJ   (npro)               : weighted Jacobian determinant
90c  A0t     (npro)               : weighted Jacobian matrix
91c  EGmasst (npro,nshape, nshape): partial LHS matrix
92c
93c output:
94c  EGmasst (npro,nshape,nshape) : partial LHS matrix
95c
96c----------------------------------------------------------------------
97c
98        include "common.h"
99c
100        dimension  A0t(npro),         shape(npro,nshl),
101     &             WdetJ(npro),       EGmasst(npro,nshape,nshape)
102c
103        dimension  fact(npro),        temp(npro),
104     &             shpij(npro),       srcp(npro)
105c
106c.... ---------------------------->  LHS  <----------------------------
107c
108        temp = WdetJ * almi/gami/alfi*Dtgl
109c
110c.... loop through columns (nodes j) and rows (nodes i)
111c
112        do j = 1, nshl
113           do i = 1, nshl
114c
115c.... get the factor
116c
117              shpij = shape(:,i) * shape(:,j)
118              fact  = shpij * temp
119              EGmasst(:,i,j) = EGmasst(:,i,j) + fact * A0t(:)
120              EGmasst(:,i,j) = EGmasst(:,i,j)
121     &                       - shpij*WdetJ(:) * srcp(:)
122c
123c.... end loop on row nodes
124c
125           enddo
126c
127c.... end loop on column nodes
128c
129        enddo
130c
131c.... return
132c
133        return
134        end
135
136