1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19      subroutine wcoef(v,vo,al,um)
20!
21!     computation of the coefficients of w in the derivation of the
22!     second order element stiffness matrix
23!
24      implicit none
25!
26      real*8 v(3,3,3,3),vo(3,3)
27!
28      real*8 a2u,al,um,au,p1,p2,p3
29!
30!
31!
32      a2u=al+2.d0*um
33      au=al+um
34!
35      p1=vo(1,1)+1.d0
36      p2=vo(2,2)+1.d0
37      p3=vo(3,3)+1.d0
38!
39      v(1,1,1,1)=a2u*p1*p1+um*(vo(1,2)**2+vo(1,3)**2)
40      v(2,1,1,1)=au*vo(1,2)*p1
41      v(3,1,1,1)=au*vo(1,3)*p1
42      v(1,2,1,1)=v(2,1,1,1)
43      v(2,2,1,1)=a2u*vo(1,2)**2+um*(p1*p1+vo(1,3)**2)
44      v(3,2,1,1)=au*vo(1,2)*vo(1,3)
45      v(1,3,1,1)=v(3,1,1,1)
46      v(2,3,1,1)=v(3,2,1,1)
47      v(3,3,1,1)=a2u*vo(1,3)**2+um*(p1*p1+vo(1,2)**2)
48!
49      v(1,1,2,1)=al*vo(2,1)*p1+
50     &  um*(2.d0*vo(2,1)*p1+vo(1,2)*p2+vo(2,3)*vo(1,3))
51      v(2,1,2,1)=al*p1*p2+um*vo(2,1)*vo(1,2)
52      v(3,1,2,1)=al*vo(2,3)*p1+um*vo(2,1)*vo(1,3)
53      v(1,2,2,1)=al*vo(2,1)*vo(1,2)+um*p1*p2
54      v(2,2,2,1)=al*vo(1,2)*p2+
55     &  um*(vo(2,1)*p1+2.d0*vo(1,2)*p2+vo(2,3)*vo(1,3))
56      v(3,2,2,1)=al*vo(2,3)*vo(1,2)+um*vo(1,3)*p2
57      v(1,3,2,1)=al*vo(2,1)*vo(1,3)+um*vo(2,3)*p1
58      v(2,3,2,1)=al*vo(1,3)*p2+um*vo(2,3)*vo(1,2)
59      v(3,3,2,1)=a2u*vo(2,3)*vo(1,3)+
60     &  um*(vo(2,1)*p1+vo(1,2)*p2)
61!
62      v(1,1,3,1)=al*vo(3,1)*p1+
63     &  um*(vo(1,3)*p3+2.d0*vo(3,1)*p1+vo(3,2)*vo(1,2))
64      v(2,1,3,1)=al*vo(3,2)*p1+um*vo(3,1)*vo(1,2)
65      v(3,1,3,1)=al*p1*p3+um*vo(3,1)*vo(1,3)
66      v(1,2,3,1)=al*vo(3,1)*vo(1,2)+um*vo(3,2)*p1
67      v(2,2,3,1)=a2u*vo(3,2)*vo(1,2)+
68     &  um*(vo(1,3)*p3+vo(3,1)*p1)
69      v(3,2,3,1)=al*vo(1,2)*p3+um*vo(3,2)*vo(1,3)
70      v(1,3,3,1)=al*vo(3,1)*vo(1,3)+um*p1*p3
71      v(2,3,3,1)=al*vo(3,2)*vo(1,3)+um*vo(1,2)*p3
72      v(3,3,3,1)=al*vo(1,3)*p3+
73     &  um*(2.d0*vo(1,3)*p3+vo(3,1)*p1+vo(3,2)*vo(1,2))
74!
75      v(1,1,1,2)=al*vo(2,1)*p1+
76     &  um*(vo(1,2)*p2+2.d0*vo(2,1)*p1+vo(1,3)*vo(2,3))
77      v(2,1,1,2)=al*vo(1,2)*vo(2,1)+um*p1*p2
78      v(3,1,1,2)=al*vo(1,3)*vo(2,1)+um*vo(2,3)*p1
79      v(1,2,1,2)=al*p1*p2+um*vo(1,2)*vo(2,1)
80      v(2,2,1,2)=al*vo(1,2)*p2+
81     &  um*(2.d0*vo(1,2)*p2+vo(2,1)*p1+vo(1,3)*vo(2,3))
82      v(3,2,1,2)=al*vo(1,3)*p2+um*vo(1,2)*vo(2,3)
83      v(1,3,1,2)=al*vo(2,3)*p1+um*vo(1,3)*vo(2,1)
84      v(2,3,1,2)=al*vo(1,2)*vo(2,3)+um*vo(1,3)*p2
85      v(3,3,1,2)=a2u*vo(1,3)*vo(2,3)+
86     &  um*(vo(1,2)*p2+vo(2,1)*p1)
87!
88      v(1,1,2,2)=a2u*vo(2,1)**2+um*(p2*p2+vo(2,3)**2)
89      v(2,1,2,2)=au*vo(2,1)*p2
90      v(3,1,2,2)=au*vo(2,3)*vo(2,1)
91      v(1,2,2,2)=v(2,1,2,2)
92      v(2,2,2,2)=a2u*p2*p2+um*(vo(2,1)**2+vo(2,3)**2)
93      v(3,2,2,2)=au*vo(2,3)*p2
94      v(1,3,2,2)=v(3,1,2,2)
95      v(2,3,2,2)=v(3,2,2,2)
96      v(3,3,2,2)=a2u*vo(2,3)**2+um*(p2*p2+vo(2,1)**2)
97!
98      v(1,1,3,2)=a2u*vo(3,1)*vo(2,1)+
99     &  um*(vo(3,2)*p2+vo(2,3)*p3)
100      v(2,1,3,2)=al*vo(3,2)*vo(2,1)+um*vo(3,1)*p2
101      v(3,1,3,2)=al*vo(2,1)*p3+um*vo(3,1)*vo(2,3)
102      v(1,2,3,2)=al*vo(3,1)*p2+um*vo(3,2)*vo(2,1)
103      v(2,2,3,2)=al*vo(3,2)*p2+
104     &  um*(2.d0*vo(3,2)*p2+vo(2,3)*p3+vo(3,1)*vo(2,1))
105      v(3,2,3,2)=al*p2*p3+um*vo(3,2)*vo(2,3)
106      v(1,3,3,2)=al*vo(3,1)*vo(2,3)+um*vo(2,1)*p3
107      v(2,3,3,2)=al*vo(3,2)*vo(2,3)+um*p2*p3
108      v(3,3,3,2)=al*vo(2,3)*p3+
109     &  um*(vo(3,2)*p2+2.d0*vo(2,3)*p3+vo(3,1)*vo(2,1))
110!
111      v(1,1,1,3)=al*vo(3,1)*p1+
112     &  um*(vo(1,3)*p3+2.d0*vo(3,1)*p1+vo(1,2)*vo(3,2))
113      v(2,1,1,3)=al*vo(1,2)*vo(3,1)+um*vo(3,2)*p1
114      v(3,1,1,3)=al*vo(1,3)*vo(3,1)+um*p1*p3
115      v(1,2,1,3)=al*vo(3,2)*p1+um*vo(1,2)*vo(3,1)
116      v(2,2,1,3)=a2u*vo(1,2)*vo(3,2)+
117     &  um*(vo(1,3)*p3+vo(3,1)*p1)
118      v(3,2,1,3)=al*vo(1,3)*vo(3,2)+um*vo(1,2)*p3
119      v(1,3,1,3)=al*p1*p3+um*vo(1,3)*vo(3,1)
120      v(2,3,1,3)=al*vo(1,2)*p3+um*vo(1,3)*vo(3,2)
121      v(3,3,1,3)=al*vo(1,3)*p3+
122     &  um*(2.d0*vo(1,3)*p3+vo(3,1)*p1+vo(1,2)*vo(3,2))
123!
124      v(1,1,2,3)=a2u*vo(2,1)*vo(3,1)+
125     &  um*(vo(2,3)*p3+vo(3,2)*p2)
126      v(2,1,2,3)=al*vo(3,1)*p2+um*vo(2,1)*vo(3,2)
127      v(3,1,2,3)=al*vo(2,3)*vo(3,1)+um*vo(2,1)*p3
128      v(1,2,2,3)=al*vo(2,1)*vo(3,2)+um*vo(3,1)*p2
129      v(2,2,2,3)=al*vo(3,2)*p2+
130     &  um*(vo(2,3)*p3+2.d0*vo(3,2)*p2+vo(2,1)*vo(3,1))
131      v(3,2,2,3)=al*vo(2,3)*vo(3,2)+um*p2*p3
132      v(1,3,2,3)=al*vo(2,1)*p3+um*vo(2,3)*vo(3,1)
133      v(2,3,2,3)=al*p2*p3+um*vo(2,3)*vo(3,2)
134      v(3,3,2,3)=al*vo(2,3)*p3+
135     &  um*(2.d0*vo(2,3)*p3+vo(3,2)*p2+vo(2,1)*vo(3,1))
136!
137      v(1,1,3,3)=a2u*vo(3,1)**2+um*(p3*p3+vo(3,2)**2)
138      v(2,1,3,3)=au*vo(3,2)*vo(3,1)
139      v(3,1,3,3)=au*vo(3,1)*p3
140      v(1,2,3,3)=v(2,1,3,3)
141      v(2,2,3,3)=a2u*vo(3,2)**2+um*(p3*p3+vo(3,1)**2)
142      v(3,2,3,3)=au*vo(3,2)*p3
143      v(1,3,3,3)=v(3,1,3,3)
144      v(2,3,3,3)=v(3,2,3,3)
145      v(3,3,3,3)=a2u*p3*p3+um*(vo(3,1)**2+vo(3,2)**2)
146!
147      return
148      end
149