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