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 extrapolate_se(yi,yn,ipkon,inum,kon,lakon,nfield,nk, 20 & ne,mi,ndim,orab,ielorien,co,iorienloc,cflag, 21 & vold,force,ielmat,thicke,ielprop,prop,ialeneigh, 22 & neaneigh,nebneigh,ialnneigh,naneigh,nbneigh) 23! 24! extrapolates field values at the integration points to the 25! nodes 26! 27! the C3D20RB element has 50 integration points, however, the 28! first 8 integration points coincide with those of a C3D20R 29! element. In this routine the C3D20RBR and C3D20RBC 30! elements are treated as an ordinary C3D20R element 31! 32! the number of internal state variables is limited to 999 33! (cfr. array field) 34! 35 implicit none 36! 37 logical force 38! 39 character*1 cflag 40 character*8 lakon(*),lakonl 41! 42 integer ipkon(*),inum(*),kon(*),mi(*),ne,indexe,nope, 43 & nonei20(3,12),nfield,nonei10(3,6),nk,i,j,k,l,ndim, 44 & nonei15(3,9),iorienloc,iorien,ielorien(mi(3),*),konl, 45 & mint3d,m,iflag,jj,ll,ielmat(mi(3),*),ielprop(*), 46 & nlayer,nopeexp,ilayer,kk,mint2d,nopes,kl,ki,null, 47 & itet(4),iwedge(2,9),ialeneigh(*),neaneigh,nebneigh,ij, 48 & ialnneigh(*),naneigh,nbneigh 49! 50 real*8 yi(ndim,mi(1),*),yn(nfield,*),field(999,20*mi(3)),a8(8,8), 51 & a4(4,4),a27(20,27),a9(6,9),a2(6,2),orab(7,*),co(3,*),prop(*), 52 & coords(3,mi(1)),xi,et,ze,xl(3,20),xsj,shp(4,20),weight, 53 & yiloc(6,mi(1)),a(3,3),b(3,3),c(3,3),vold(0:mi(2),*),tlayer(4), 54 & dlayer(4),xlayer(mi(3),4),thickness,xs2(3,7),xl2(3,8), 55 & xsj2(3),shp2(7,8),thicke(mi(3),*),coloc(3,8), 56 & xwedge(2,2,9),a14(8,14),a6(6,6) 57! 58! 59! 60 include "gauss.f" 61! 62 data coloc /-1.d0,-1.d0,-1.d0, 63 & 1.d0,-1.d0,-1.d0, 64 & -1.d0,1.d0,-1.d0, 65 & 1.d0,1.d0,-1.d0, 66 & -1.d0,-1.d0,1.d0, 67 & 1.d0,-1.d0,1.d0, 68 & -1.d0,1.d0,1.d0, 69 & 1.d0,1.d0,1.d0/ 70 data null /0/ 71! 72! a 10-node tet is remeshed into 10 4-node tets at contact 73! interfaces; itet contains the linear tet elements to which 74! the integration points of the parent 10-node tet belong 75! 76 data itet /1,2,3,10/ 77! 78 data iwedge /1,0,2,0,3,0,1,5,2,6,3,7,5,0,6,0,7,0/ 79! 80 data xwedge /0.975615382715435242d0,0.243846172845647580d-1, 81 & 0.d0,0.d0, 82 & 0.975615382715435242d0,0.243846172845647580d-1, 83 & 0.d0,0.d0, 84 & 0.975615382715435242d0,0.243846172845647580d-1, 85 & 0.d0,0.d0, 86 & 0.d0,.5d0,.5d0,0.d0, 87 & 0.d0,.5d0,.5d0,0.d0, 88 & 0.d0,.5d0,.5d0,0.d0, 89 & 0.243846172845647580d-1,0.975615382715435242d0, 90 & 0.d0,0.d0, 91 & 0.243846172845647580d-1,0.975615382715435242d0, 92 & 0.d0,0.d0, 93 & 0.243846172845647580d-1,0.975615382715435242d0, 94 & 0.d0,0.d0/ 95! 96 data nonei10 /5,1,2,6,2,3,7,3,1,8,1,4,9,2,4,10,3,4/ 97! 98 data nonei15 /7,1,2,8,2,3,9,3,1,10,4,5,11,5,6,12,6,4, 99 & 13,1,4,14,2,5,15,3,6/ 100! 101 data nonei20 /9,1,2,10,2,3,11,3,4,12,4,1, 102 & 13,5,6,14,6,7,15,7,8,16,8,5, 103 & 17,1,5,18,2,6,19,3,7,20,4,8/ 104! 105 data a2 /1.1455d0,-0.1455d0,1.1455d0,-0.1455d0,1.1455d0,-0.1455d0, 106 & -0.1455d0,1.1455d0,-0.1455d0,1.1455d0,-0.1455d0,1.1455d0/ 107 data a4 / 1.92705d0, -0.30902d0, -0.30902d0, -0.30902d0, 108 & -0.30902d0, 1.92705d0, -0.30902d0, -0.30902d0, 109 & -0.30902d0, -0.30902d0, 1.92705d0, -0.30902d0, 110 & -0.30902d0, -0.30902d0, -0.30902d0, 1.92705d0/ 111! 112! extrapolation from a 6 integration point scheme in a wedge to 113! the vertex nodes 114! 115 data a6 / 2.04904d0, 0.00000d0, 0.00000d0,-0.54904d0, 0.00000d0, 116 & 0.00000d0, 117 & -0.34151d0, 1.70753d0,-0.34151d0, 0.09151d0,-0.45753d0, 118 & 0.09151d0, 119 & -0.34151d0,-0.34151d0, 1.70753d0, 0.09151d0, 0.09151d0, 120 & -0.45753d0, 121 & -0.54904d0, 0.00000d0, 0.00000d0, 2.04904d0, 0.00000d0, 122 & 0.00000d0, 123 & 0.09151d0,-0.45753d0, 0.09151d0,-0.34151d0, 1.70753d0, 124 & -0.34151d0, 125 & 0.09151d0, 0.09151d0,-0.45753d0,-0.34151d0,-0.34151d0, 126 & 1.70753d0/ 127! 128 data a9 / 1.63138d0,-0.32628d0,-0.32628d0,-0.52027d0, 0.10405d0, 129 & 0.10405d0, 130 & -0.32628d0, 1.63138d0,-0.32628d0, 0.10405d0,-0.52027d0, 131 & 0.10405d0, 132 & -0.32628d0,-0.32628d0, 1.63138d0, 0.10405d0, 0.10405d0, 133 & -0.52027d0, 134 & 0.55556d0,-0.11111d0,-0.11111d0, 0.55556d0,-0.11111d0, 135 & -0.11111d0, 136 & -0.11111d0, 0.55556d0,-0.11111d0,-0.11111d0,0.55556d0, 137 & -0.11111d0, 138 & -0.11111d0,-0.11111d0, 0.55556d0,-0.11111d0,-0.11111d0, 139 & 0.55556d0, 140 & -0.52027d0, 0.10405d0, 0.10405d0, 1.63138d0,-0.32628d0, 141 & -0.32628d0, 142 & 0.10405d0,-0.52027d0, 0.10405d0,-0.32628d0, 1.63138d0, 143 & -0.32628d0, 144 & 0.10405d0, 0.10405d0,-0.52027d0,-0.32628d0,-0.32628d0, 145 & 1.63138d0/ 146! 147! extrapolation from a 2x2x2=8 integration point scheme in a hex to 148! the vertex nodes 149! 150 data a8 /2.549d0,-.683d0,.183d0,-.683d0,-.683d0,.183d0, 151 & -.04904d0,.183d0,-.683d0,2.549d0,-.683d0,.183d0, 152 & .183d0,-.683d0,.183d0,-.04904d0,-.683d0,.183d0, 153 & -.683d0,2.549d0,.183d0,-.04904d0,.183d0,-.683d0, 154 & .183d0,-.683d0,2.549d0,-.683d0,-.04904d0,.183d0, 155 & -.683d0,.183d0,-.683d0,.183d0,-.04904d0,.183d0, 156 & 2.549d0,-.683d0,.183d0,-.683d0,.183d0,-.683d0, 157 & .183d0,-.04904d0,-.683d0,2.549d0,-.683d0,.183d0, 158 & .183d0,-.04904d0,.183d0,-.683d0,-.683d0,.183d0, 159 & -.683d0,2.549d0,-.04904d0,.183d0,-.683d0,.183d0, 160 & .183d0,-.683d0,2.549d0,-.683d0/ 161! 162! extrapolation from a 14 integration point scheme in a hex to 163! the vertex nodes 164! 165 data a14 / 166 & 0.1396d+01,-0.3026d+00,0.1124d-01,-0.3026d+00, 167 & -0.3026d+00,0.1124d-01,0.4901d-01,0.1124d-01, 168 & -0.3026d+00,0.1396d+01,-0.3026d+00,0.1124d-01, 169 & 0.1124d-01,-0.3026d+00,0.1124d-01,0.4901d-01, 170 & 0.1124d-01,-0.3026d+00,0.1396d+01,-0.3026d+00, 171 & 0.4901d-01,0.1124d-01,-0.3026d+00,0.1124d-01, 172 & -0.3026d+00,0.1124d-01,-0.3026d+00,0.1396d+01, 173 & 0.1124d-01,0.4901d-01,0.1124d-01,-0.3026d+00, 174 & -0.3026d+00,0.1124d-01,0.4901d-01,0.1124d-01, 175 & 0.1396d+01,-0.3026d+00,0.1124d-01,-0.3026d+00, 176 & 0.1124d-01,-0.3026d+00,0.1124d-01,0.4901d-01, 177 & -0.3026d+00,0.1396d+01,-0.3026d+00,0.1124d-01, 178 & 0.4901d-01,0.1124d-01,-0.3026d+00,0.1124d-01, 179 & 0.1124d-01,-0.3026d+00,0.1396d+01,-0.3026d+00, 180 & 0.1124d-01,0.4901d-01,0.1124d-01,-0.3026d+00, 181 & -0.3026d+00,0.1124d-01,-0.3026d+00,0.1396d+01, 182 & 0.2069d+00,0.2069d+00,-0.6408d-01,-0.6408d-01, 183 & 0.2069d+00,0.2069d+00,-0.6408d-01,-0.6408d-01, 184 & -0.6408d-01,0.2069d+00,0.2069d+00,-0.6408d-01, 185 & -0.6408d-01,0.2069d+00,0.2069d+00,-0.6408d-01, 186 & -0.6408d-01,-0.6408d-01,0.2069d+00,0.2069d+00, 187 & -0.6408d-01,-0.6408d-01,0.2069d+00,0.2069d+00, 188 & 0.2069d+00,-0.6408d-01,-0.6408d-01,0.2069d+00, 189 & 0.2069d+00,-0.6408d-01,-0.6408d-01,0.2069d+00, 190 & 0.2069d+00,0.2069d+00,0.2069d+00,0.2069d+00, 191 & -0.6408d-01,-0.6408d-01,-0.6408d-01,-0.6408d-01, 192 & -0.6408d-01,-0.6408d-01,-0.6408d-01,-0.6408d-01, 193 & 0.2069d+00,0.2069d+00,0.2069d+00,0.2069d+00/ 194! 195! extrapolation from a 3x3x3=27 integration point scheme in a hex to 196! the all nodes in a 20-node element 197! 198 data a27 / 199 & 2.37499d0,-0.12559d0,-0.16145d0,-0.12559d0,-0.12559d0, 200 & -0.16145d0, 0.11575d0, 201 & -0.16145d0, 0.32628d0, 0.11111d0, 0.11111d0, 0.32628d0, 202 & 0.11111d0,-0.10405d0, 203 & -0.10405d0, 0.11111d0, 0.32628d0, 0.11111d0,-0.10405d0, 204 & 0.11111d0,-0.31246d0, 205 & -0.31246d0, 0.31481d0, 0.31481d0, 0.31481d0, 0.31481d0, 206 & -0.16902d0,-0.16902d0, 207 & 1.28439d0,-0.27072d0,-0.19444d0,-0.27072d0,-0.19444d0, 208 & 0.15961d0,-0.00661d0, 209 & 0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0, 0.15961d0, 210 & -0.12559d0, 2.37499d0, 211 & -0.12559d0,-0.16145d0,-0.16145d0,-0.12559d0,-0.16145d0, 212 & 0.11575d0, 0.32628d0, 213 & 0.32628d0, 0.11111d0, 0.11111d0, 0.11111d0, 0.11111d0, 214 & -0.10405d0,-0.10405d0, 215 & 0.11111d0, 0.32628d0, 0.11111d0,-0.10405d0,-0.31246d0, 216 & 0.31481d0, 0.31481d0, 217 & -0.31246d0, 0.31481d0,-0.16902d0,-0.16902d0, 0.31481d0, 218 & -0.27072d0,-0.19444d0, 219 & -0.27072d0, 1.28439d0, 0.15961d0,-0.00661d0, 0.15961d0, 220 & -0.19444d0,-0.27072d0, 221 & 0.15961d0, 0.15961d0,-0.27072d0,-0.48824d0,-0.48824d0, 222 & -0.48824d0,-0.48824d0, 223 & 0.22898d0, 0.22898d0, 0.22898d0, 0.22898d0, 0.05556d0, 224 & 0.05556d0, 0.05556d0, 225 & 0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0, 226 & -0.22222d0,-0.22222d0, 227 & -0.22222d0,-0.22222d0, 0.31481d0,-0.31246d0,-0.31246d0, 228 & 0.31481d0,-0.16902d0, 229 & 0.31481d0, 0.31481d0,-0.16902d0,-0.27072d0, 1.28439d0, 230 & -0.27072d0,-0.19444d0, 231 & 0.15961d0,-0.19444d0, 0.15961d0,-0.00661d0, 0.15961d0, 232 & -0.27072d0,-0.27072d0, 233 & 0.15961d0,-0.12559d0,-0.16145d0,-0.12559d0, 2.37499d0, 234 & -0.16145d0, 0.11575d0, 235 & -0.16145d0,-0.12559d0, 0.11111d0, 0.11111d0, 0.32628d0, 236 & 0.32628d0,-0.10405d0, 237 & -0.10405d0, 0.11111d0, 0.11111d0, 0.11111d0,-0.10405d0, 238 & 0.11111d0, 0.32628d0, 239 & 0.31481d0, 0.31481d0,-0.31246d0,-0.31246d0,-0.16902d0, 240 & -0.16902d0, 0.31481d0, 241 & 0.31481d0,-0.19444d0,-0.27072d0, 1.28439d0,-0.27072d0, 242 & -0.00661d0, 0.15961d0, 243 & -0.19444d0, 0.15961d0, 0.15961d0, 0.15961d0,-0.27072d0, 244 & -0.27072d0,-0.16145d0, 245 & -0.12559d0, 2.37499d0,-0.12559d0, 0.11575d0,-0.16145d0, 246 & -0.12559d0,-0.16145d0, 247 & 0.11111d0, 0.32628d0, 0.32628d0, 0.11111d0,-0.10405d0, 248 & 0.11111d0, 0.11111d0, 249 & -0.10405d0,-0.10405d0, 0.11111d0, 0.32628d0, 0.11111d0, 250 & -0.31246d0, 0.31481d0, 251 & -0.16902d0, 0.31481d0,-0.31246d0, 0.31481d0,-0.16902d0, 252 & 0.31481d0,-0.27072d0, 253 & 0.15961d0, 0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0, 254 & 0.15961d0,-0.27072d0, 255 & 1.28439d0,-0.19444d0,-0.00661d0,-0.19444d0,-0.48824d0, 256 & -0.48824d0, 0.22898d0, 257 & 0.22898d0,-0.48824d0,-0.48824d0, 0.22898d0, 0.22898d0, 258 & 0.05556d0,-0.22222d0, 259 & 0.05556d0,-0.22222d0, 0.05556d0,-0.22222d0, 0.05556d0, 260 & -0.22222d0, 0.05556d0, 261 & 0.05556d0, 0.05556d0, 0.05556d0, 0.31481d0,-0.31246d0, 262 & 0.31481d0,-0.16902d0, 263 & 0.31481d0,-0.31246d0, 0.31481d0,-0.16902d0,-0.27072d0, 264 & -0.27072d0, 0.15961d0, 265 & 0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0, 0.15961d0, 266 & -0.19444d0, 1.28439d0, 267 & -0.19444d0,-0.00661d0,-0.48824d0, 0.22898d0, 0.22898d0, 268 & -0.48824d0,-0.48824d0, 269 & 0.22898d0, 0.22898d0,-0.48824d0,-0.22222d0, 0.05556d0, 270 & -0.22222d0, 0.05556d0, 271 & -0.22222d0, 0.05556d0,-0.22222d0, 0.05556d0, 0.05556d0, 272 & 0.05556d0, 0.05556d0, 273 & 0.05556d0,-0.29630d0,-0.29630d0,-0.29630d0,-0.29630d0, 274 & -0.29630d0,-0.29630d0, 275 & -0.29630d0,-0.29630d0,-0.11111d0,-0.11111d0,-0.11111d0, 276 & -0.11111d0,-0.11111d0, 277 & -0.11111d0,-0.11111d0,-0.11111d0,-0.11111d0,-0.11111d0, 278 & -0.11111d0,-0.11111d0, 279 & 0.22898d0,-0.48824d0,-0.48824d0, 0.22898d0, 0.22898d0, 280 & -0.48824d0,-0.48824d0, 281 & 0.22898d0,-0.22222d0, 0.05556d0,-0.22222d0, 0.05556d0, 282 & -0.22222d0, 0.05556d0, 283 & -0.22222d0, 0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0, 284 & 0.05556d0, 0.31481d0, 285 & -0.16902d0, 0.31481d0,-0.31246d0, 0.31481d0,-0.16902d0, 286 & 0.31481d0,-0.31246d0, 287 & 0.15961d0, 0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0, 288 & 0.15961d0,-0.27072d0, 289 & -0.27072d0,-0.19444d0,-0.00661d0,-0.19444d0, 1.28439d0, 290 & 0.22898d0, 0.22898d0, 291 & -0.48824d0,-0.48824d0, 0.22898d0, 0.22898d0,-0.48824d0, 292 & -0.48824d0, 0.05556d0, 293 & -0.22222d0, 0.05556d0,-0.22222d0, 0.05556d0,-0.22222d0, 294 & 0.05556d0,-0.22222d0, 295 & 0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0,-0.16902d0, 296 & 0.31481d0,-0.31246d0, 297 & 0.31481d0,-0.16902d0, 0.31481d0,-0.31246d0, 0.31481d0, 298 & 0.15961d0,-0.27072d0, 299 & -0.27072d0, 0.15961d0, 0.15961d0,-0.27072d0,-0.27072d0 300 & , 0.15961d0,-0.00661d0, 301 & -0.19444d0, 1.28439d0,-0.19444d0,-0.12559d0,-0.16145d0, 302 & 0.11575d0,-0.16145d0, 303 & 2.37499d0,-0.12559d0,-0.16145d0,-0.12559d0, 0.11111d0, 304 & -0.10405d0,-0.10405d0, 305 & 0.11111d0, 0.32628d0, 0.11111d0, 0.11111d0, 0.32628d0, 306 & 0.32628d0, 0.11111d0, 307 & -0.10405d0, 0.11111d0, 0.31481d0, 0.31481d0,-0.16902d0, 308 & -0.16902d0,-0.31246d0, 309 & -0.31246d0, 0.31481d0, 0.31481d0,-0.19444d0, 0.15961d0, 310 & -0.00661d0, 0.15961d0, 311 & 1.28439d0,-0.27072d0,-0.19444d0,-0.27072d0,-0.27072d0, 312 & -0.27072d0, 0.15961d0, 313 & 0.15961d0,-0.16145d0,-0.12559d0,-0.16145d0, 0.11575d0, 314 & -0.12559d0, 2.37499d0, 315 & -0.12559d0,-0.16145d0, 0.11111d0, 0.11111d0,-0.10405d0, 316 & -0.10405d0, 0.32628d0, 317 & 0.32628d0, 0.11111d0, 0.11111d0, 0.11111d0, 0.32628d0, 318 & 0.11111d0,-0.10405d0, 319 & 0.31481d0,-0.16902d0,-0.16902d0, 0.31481d0,-0.31246d0, 320 & 0.31481d0, 0.31481d0, 321 & -0.31246d0, 0.15961d0,-0.00661d0, 0.15961d0,-0.19444d0, 322 & -0.27072d0,-0.19444d0, 323 & -0.27072d0, 1.28439d0,-0.27072d0, 0.15961d0, 0.15961d0, 324 & -0.27072d0, 0.22898d0, 325 & 0.22898d0, 0.22898d0, 0.22898d0,-0.48824d0,-0.48824d0, 326 & -0.48824d0,-0.48824d0, 327 & 0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0, 0.05556d0, 328 & 0.05556d0, 0.05556d0, 329 & 0.05556d0,-0.22222d0,-0.22222d0,-0.22222d0,-0.22222d0, 330 & -0.16902d0, 0.31481d0, 331 & 0.31481d0,-0.16902d0, 0.31481d0,-0.31246d0,-0.31246d0, 332 & 0.31481d0, 0.15961d0, 333 & -0.19444d0, 0.15961d0,-0.00661d0,-0.27072d0, 1.28439d0, 334 & -0.27072d0,-0.19444d0, 335 & 0.15961d0,-0.27072d0,-0.27072d0, 0.15961d0,-0.16145d0, 336 & 0.11575d0,-0.16145d0, 337 & -0.12559d0,-0.12559d0,-0.16145d0,-0.12559d0, 2.37499d0, 338 & -0.10405d0,-0.10405d0, 339 & 0.11111d0, 0.11111d0, 0.11111d0, 0.11111d0, 0.32628d0, 340 & 0.32628d0, 0.11111d0, 341 & -0.10405d0, 0.11111d0, 0.32628d0,-0.16902d0,-0.16902d0, 342 & 0.31481d0, 0.31481d0, 343 & 0.31481d0, 0.31481d0,-0.31246d0,-0.31246d0,-0.00661d0, 344 & 0.15961d0,-0.19444d0, 345 & 0.15961d0,-0.19444d0,-0.27072d0, 1.28439d0,-0.27072d0, 346 & 0.15961d0, 0.15961d0, 347 & -0.27072d0,-0.27072d0, 0.11575d0,-0.16145d0,-0.12559d0, 348 & -0.16145d0,-0.16145d0, 349 & -0.12559d0, 2.37499d0,-0.12559d0,-0.10405d0, 0.11111d0, 350 & 0.11111d0,-0.10405d0, 351 & 0.11111d0, 0.32628d0, 0.32628d0, 0.11111d0,-0.10405d0, 352 & 0.11111d0, 0.32628d0, 353 & 0.11111d0/ 354! 355 data iflag /1/ 356! 357 do ij=naneigh,nbneigh 358 i=ialnneigh(ij) 359 inum(i)=0 360 enddo 361! 362 do ij=naneigh,nbneigh 363 i=ialnneigh(ij) 364 do j=1,nfield 365 yn(j,i)=0.d0 366 enddo 367 enddo 368! 369 do ij=neaneigh,nebneigh 370! 371 i=ialeneigh(ij) 372! 373 if(ipkon(i).lt.0) cycle 374! 375 indexe=ipkon(i) 376 lakonl=lakon(i) 377! 378 if(lakonl(7:8).eq.'LC') then 379 nlayer=0 380 do j=1,mi(3) 381 if(ielmat(j,i).gt.0) then 382 nlayer=nlayer+1 383 else 384 exit 385 endif 386 enddo 387! 388 if(lakonl(4:4).eq.'2') then 389 nopeexp=28 390 elseif(lakonl(4:5).eq.'15') then 391 nopeexp=21 392 endif 393 endif 394! 395 if(lakonl(4:4).eq.'2') then 396 nope=20 397 elseif(lakonl(4:4).eq.'8') then 398 nope=8 399 elseif(lakonl(4:5).eq.'10') then 400 nope=10 401 elseif(lakonl(4:4).eq.'4') then 402 nope=4 403 elseif(lakonl(4:5).eq.'15') then 404 nope=15 405 elseif(lakonl(4:4).eq.'6') then 406 nope=6 407 elseif((lakonl(1:1).eq.'E').and. 408 & ((lakonl(7:7).eq.'A').or. 409 & (lakonl(7:7).eq.'2'))) then 410! 411 inum(kon(indexe+1))=inum(kon(indexe+1))+1 412 inum(kon(indexe+2))=inum(kon(indexe+2))+1 413 cycle 414 elseif(lakonl(1:7).eq.'ESPRNGF') then 415 read(lakonl(8:8),'(i1)') nope 416 nope=nope+1 417 inum(kon(indexe+nope))=-1 418 cycle 419 elseif(lakonl(1:1).eq.'U') then 420 call extrapolate_u(yi,yn,ipkon,inum,kon,lakon,nfield,nk, 421 & ne,mi,ndim,orab,ielorien,co,iorienloc,cflag, 422 & vold,force,ielmat,thicke,ielprop,prop,i) 423 return 424 else 425 cycle 426 endif 427! 428! storage in local coordinates 429! 430! calculation of the integration point coordinates for 431! output in the local system 432! 433 if((iorienloc.ne.0).and. 434 & ((lakonl(7:8).eq.'LC').or.(ielorien(1,i).ne.0))) then 435! 436 if(lakonl(7:8).ne.'LC') then 437 iorien=max(0,ielorien(1,i)) 438!david -start 439 elseif(lakonl(4:5).eq.'20') then 440! 441! composite materials 442! 443 mint2d=4 444 nopes=8 445! 446! determining the layer thickness and global thickness 447! at the shell integration points 448! 449! xlayer: actual layer thickness 450! tlayer: total thickness of all layers 451! dlayer: total thickness of all layers up to the actual one 452! 453 indexe=ipkon(i) 454 do kk=1,mint2d 455 xi=gauss3d2(1,kk) 456 et=gauss3d2(2,kk) 457 call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag) 458 tlayer(kk)=0.d0 459 do k=1,nlayer 460 thickness=0.d0 461 do j=1,nopes 462 thickness=thickness+thicke(k,indexe+j)*shp2(4,j) 463 enddo 464 tlayer(kk)=tlayer(kk)+thickness 465 xlayer(k,kk)=thickness 466 enddo 467 enddo 468! 469 ilayer=0 470 do k=1,mint2d 471 dlayer(k)=0.d0 472 enddo 473! 474! S6-composite shell 475! 476 elseif(lakonl(4:5).eq.'15') then 477! 478! composite materials 479! 480 mint2d=3 481 nopes=6 482! 483! determining the layer thickness and global thickness 484! at the shell integration points 485! 486! xlayer: actual layer thickness 487! tlayer: total thickness of all layers 488! dlayer: total thickness of all layers up to the actual one 489! 490 indexe=ipkon(i) 491 do kk=1,mint2d 492 xi=gauss3d10(1,kk) 493 et=gauss3d10(2,kk) 494 call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 495 tlayer(kk)=0.d0 496 do k=1,nlayer 497 thickness=0.d0 498 do j=1,nopes 499 thickness=thickness+thicke(k,indexe+j)*shp2(4,j) 500 enddo 501 tlayer(kk)=tlayer(kk)+thickness 502 xlayer(k,kk)=thickness 503 enddo 504 enddo 505! 506 ilayer=0 507 do k=1,mint2d 508 dlayer(k)=0.d0 509 enddo 510 endif 511! 512 if((lakon(i)(4:5).eq.'8R').or.(lakon(i)(1:1).eq.'F')) then 513 mint3d=1 514 elseif((lakon(i)(4:7).eq.'20RB').and. 515 & (lakon(i)(8:8).ne.'R').and. 516 & (lakon(i)(8:8).ne.'C')) then 517 call beamintscheme(lakon(i),mint3d,ielprop(i),prop, 518 & null,xi,et,ze,weight) 519 elseif((lakon(i)(4:4).eq.'8').or. 520 & (lakon(i)(4:6).eq.'20R')) then 521 if(lakonl(7:8).eq.'LC') then 522 mint3d=8*nlayer 523 else 524 mint3d=8 525 endif 526 elseif(lakon(i)(4:4).eq.'2') then 527 mint3d=27 528 elseif((lakon(i)(4:5).eq.'10')) then 529 mint3d=4 530 elseif(lakon(i)(4:4).eq.'4') then 531 mint3d=1 532 elseif(lakon(i)(4:5).eq.'15') then 533 if(lakonl(7:8).eq.'LC') then 534 mint3d=6*nlayer 535 else 536 mint3d=9 537 endif 538 elseif(lakon(i)(4:4).eq.'6') then 539 mint3d=2 540 endif 541! 542 do j=1,nope 543 konl=kon(indexe+j) 544 do k=1,3 545 xl(k,j)=co(k,konl) 546 enddo 547 enddo 548! 549 do j=1,mint3d 550 if((lakon(i)(4:5).eq.'8R').or. 551 & (lakon(i)(1:4).eq.'F3D8')) then 552 xi=gauss3d1(1,j) 553 et=gauss3d1(2,j) 554 ze=gauss3d1(3,j) 555 weight=weight3d1(j) 556 elseif((lakon(i)(4:7).eq.'20RB').and. 557 & (lakon(i)(8:8).ne.'R').and. 558 & (lakon(i)(8:8).ne.'C')) then 559 call beamintscheme(lakon(i),mint3d,ielprop(i),prop, 560 & j,xi,et,ze,weight) 561 elseif((lakon(i)(4:4).eq.'8').or. 562 & (lakon(i)(4:6).eq.'20R')) 563 & then 564 if(lakonl(7:8).ne.'LC') then 565 xi=gauss3d2(1,j) 566 et=gauss3d2(2,j) 567 ze=gauss3d2(3,j) 568 weight=weight3d2(j) 569 else 570! 571! kl: number of the integration point within the layer 572! 573 kl=mod(j,8) 574 if(kl.eq.0) kl=8 575! 576 xi=gauss3d2(1,kl) 577 et=gauss3d2(2,kl) 578 ze=gauss3d2(3,kl) 579 weight=weight3d2(kl) 580! 581! ki: position of the integration point (1...4) 582! 583 ki=mod(j,4) 584 if(ki.eq.0) ki=4 585! 586 if(kl.eq.1) then 587 ilayer=ilayer+1 588 if(ilayer.gt.1) then 589 do k=1,4 590 dlayer(k)=dlayer(k)+xlayer(ilayer-1,k) 591 enddo 592 endif 593 endif 594 ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0* 595 & xlayer(ilayer,ki))/tlayer(ki)-1.d0 596 weight=weight*xlayer(ilayer,ki)/tlayer(ki) 597! 598! material and orientation 599! 600 iorien=max(0,ielorien(ilayer,i)) 601 endif 602 elseif(lakon(i)(4:4).eq.'2') then 603 xi=gauss3d3(1,j) 604 et=gauss3d3(2,j) 605 ze=gauss3d3(3,j) 606 weight=weight3d3(j) 607 elseif(lakon(i)(4:5).eq.'10') then 608 xi=gauss3d5(1,j) 609 et=gauss3d5(2,j) 610 ze=gauss3d5(3,j) 611 weight=weight3d5(j) 612 elseif(lakon(i)(4:4).eq.'4') then 613 xi=gauss3d4(1,j) 614 et=gauss3d4(2,j) 615 ze=gauss3d4(3,j) 616 weight=weight3d4(j) 617 elseif(lakon(i)(4:5).eq.'15') then 618 if(lakonl(7:8).ne.'LC') then 619 xi=gauss3d8(1,j) 620 et=gauss3d8(2,j) 621 ze=gauss3d8(3,j) 622 weight=weight3d8(j) 623 else 624! 625! kl: number of the integration point within the layer 626! 627 kl=mod(j,6) 628 if(kl.eq.0) kl=6 629! 630 xi=gauss3d10(1,kl) 631 et=gauss3d10(2,kl) 632 ze=gauss3d10(3,kl) 633 weight=weight3d10(kl) 634! 635! ki: position of the integration point (1...3) 636! 637 ki=mod(j,3) 638 if(ki.eq.0) ki=3 639! 640 if(kl.eq.1) then 641 ilayer=ilayer+1 642 if(ilayer.gt.1) then 643 do k=1,3 644 dlayer(k)=dlayer(k)+xlayer(ilayer-1,k) 645 enddo 646 endif 647 endif 648 ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0* 649 & xlayer(ilayer,ki))/tlayer(ki)-1.d0 650 weight=weight*xlayer(ilayer,ki)/tlayer(ki) 651! 652! material and orientation 653! 654 iorien=max(0,ielorien(ilayer,i)) 655 endif 656 elseif(lakon(i)(1:4).eq.'C3D6') then 657 xi=gauss3d7(1,j) 658 et=gauss3d7(2,j) 659 ze=gauss3d7(3,j) 660 weight=weight3d7(j) 661 elseif(lakon(i)(1:4).eq.'F3D6') then 662 xi=gauss3d14(1,j) 663 et=gauss3d14(2,j) 664 ze=gauss3d14(3,j) 665 weight=weight3d14(j) 666 endif 667! 668 if(nope.eq.20) then 669 call shape20h(xi,et,ze,xl,xsj,shp,iflag) 670 elseif(nope.eq.8) then 671 call shape8h(xi,et,ze,xl,xsj,shp,iflag) 672 elseif(nope.eq.10) then 673 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 674 elseif(nope.eq.4) then 675 call shape4tet(xi,et,ze,xl,xsj,shp,iflag) 676 elseif(nope.eq.15) then 677 call shape15w(xi,et,ze,xl,xsj,shp,iflag) 678 else 679 call shape6w(xi,et,ze,xl,xsj,shp,iflag) 680 endif 681! 682! layer without orientation in a composite 683! 684 if(iorien.eq.0) then 685 if(nfield.eq.3) then 686 do k=1,3 687 yiloc(k,j)=yi(k,j,i) 688 enddo 689 elseif(nfield.eq.6) then 690 do k=1,6 691 yiloc(k,j)=yi(k,j,i) 692 enddo 693 endif 694 cycle 695 endif 696! 697! coordinates of the integration point 698! 699 do k=1,3 700 coords(k,j)=0.d0 701 do l=1,nope 702 coords(k,j)=coords(k,j)+xl(k,l)*shp(4,l) 703 enddo 704 enddo 705! 706! transforming the vector or tensor field 707! 708 if(nfield.eq.3) then 709 call transformatrix(orab(1,iorien),coords(1,j),a) 710 yiloc(1,j)=yi(1,j,i)*a(1,1)+yi(2,j,i)*a(2,1)+ 711 & yi(3,j,i)*a(3,1) 712 yiloc(2,j)=yi(1,j,i)*a(1,2)+yi(2,j,i)*a(2,2)+ 713 & yi(3,j,i)*a(3,2) 714 yiloc(3,j)=yi(1,j,i)*a(1,3)+yi(2,j,i)*a(2,3)+ 715 & yi(3,j,i)*a(3,3) 716 elseif(nfield.eq.6) then 717 call transformatrix(orab(1,iorien),coords(1,j),a) 718 b(1,1)=yi(1,j,i) 719 b(2,2)=yi(2,j,i) 720 b(3,3)=yi(3,j,i) 721 b(1,2)=yi(4,j,i) 722 b(1,3)=yi(5,j,i) 723 b(2,3)=yi(6,j,i) 724 b(2,1)=b(1,2) 725 b(3,1)=b(1,3) 726 b(3,2)=b(2,3) 727 do k=1,3 728 do l=1,3 729 c(k,l)=0.d0 730 do m=1,3 731 c(k,l)=c(k,l)+b(k,m)*a(m,l) 732 enddo 733 enddo 734 enddo 735 do k=1,3 736 do l=k,3 737 b(k,l)=0.d0 738 do m=1,3 739 b(k,l)=b(k,l)+a(m,k)*c(m,l) 740 enddo 741 enddo 742 enddo 743 yiloc(1,j)=b(1,1) 744 yiloc(2,j)=b(2,2) 745 yiloc(3,j)=b(3,3) 746 yiloc(4,j)=b(1,2) 747 yiloc(5,j)=b(1,3) 748 yiloc(6,j)=b(2,3) 749 endif 750 enddo 751! 752 if(lakonl(1:1).eq.'F') then 753 do j=1,8 754 do k=1,nfield 755 field(k,j)=yiloc(k,1) 756 enddo 757 enddo 758 elseif((lakonl(4:7).eq.'20RB').and. 759 & (lakonl(8:8).ne.'R').and. 760 & (lakonl(8:8).ne.'C')) then 761 call beamextscheme(yi(1,1,i),ndim,nfield,lakonl, 762 & ielprop(i),prop,field,mi) 763c Bernhardi start 764 elseif((lakonl(4:6).eq.'20R').or. 765 & (lakonl(4:5).eq.'8 ').or.(lakonl(4:5).eq.'8I')) then 766c Bernhardi end 767 if(lakonl(7:8).ne.'LC') then 768 do j=1,8 769 do k=1,nfield 770 field(k,j)=0.d0 771 do l=1,8 772 field(k,j)=field(k,j)+a8(j,l)*yiloc(k,l) 773 enddo 774 enddo 775 enddo 776 else 777 do m=1,nlayer 778 jj=20*(m-1) 779 ll=8*(m-1) 780 do j=1,8 781 do k=1,nfield 782 field(k,jj+j)=0.d0 783 do l=1,8 784 field(k,jj+j)= 785 & field(k,jj+j)+a8(j,l)*yiloc(k,ll+l) 786 enddo 787 enddo 788 enddo 789 enddo 790 endif 791 elseif(lakonl(4:4).eq.'8') then 792 do j=1,8 793 do k=1,nfield 794 field(k,j)=yiloc(k,1) 795 enddo 796 enddo 797 elseif(lakonl(4:5).eq.'10') then 798 do j=1,4 799 do k=1,nfield 800 field(k,j)=0.d0 801 do l=1,4 802 field(k,j)=field(k,j)+a4(j,l)*yiloc(k,l) 803 enddo 804 enddo 805 enddo 806 elseif(lakonl(4:4).eq.'2') then 807 do j=1,20 808 do k=1,nfield 809 field(k,j)=0.d0 810 do l=1,27 811 field(k,j)=field(k,j)+a27(j,l)*yiloc(k,l) 812 enddo 813 enddo 814 enddo 815 elseif(lakonl(4:4).eq.'4') then 816 do j=1,4 817 do k=1,nfield 818 field(k,j)=yiloc(k,1) 819 enddo 820 enddo 821 elseif(lakonl(4:4).eq.'1') then 822 if(lakonl(7:8).ne.'LC') then 823 do j=1,6 824 do k=1,nfield 825 field(k,j)=0.d0 826 do l=1,9 827 field(k,j)=field(k,j)+a9(j,l)*yiloc(k,l) 828 enddo 829 enddo 830 enddo 831 else 832 do m=1,nlayer 833 jj=15*(m-1) 834 ll=6*(m-1) 835 do j=1,6 836 do k=1,nfield 837 field(k,jj+j)=0.d0 838 do l=1,6 839 field(k,jj+j)= 840 & field(k,jj+j)+a6(j,l)*yiloc(k,ll+l) 841 enddo 842 enddo 843 enddo 844 enddo 845 endif 846 else 847 do j=1,6 848 do k=1,nfield 849 field(k,j)=0.d0 850 do l=1,2 851 field(k,j)=field(k,j)+a2(j,l)*yiloc(k,l) 852 enddo 853 enddo 854 enddo 855 endif 856 else 857! 858! storage in global coordinates 859! 860! determining the field values in the vertex nodes 861! for C3D20R and C3D8: trilinear extrapolation (= use of the 862! C3D8 linear interpolation functions) 863! for C3D8R: constant field value in each element 864! for C3D10: use of the C3D4 linear interpolation functions 865! for C3D4: constant field value in each element 866! for C3D15: use of the C3D6 linear interpolation functions 867! for C3D6: use of a linear interpolation function 868! 869 if(lakonl(1:1).eq.'F') then 870 do j=1,8 871 do k=1,nfield 872 field(k,j)=yi(k,1,i) 873 enddo 874 enddo 875 elseif((lakonl(4:7).eq.'20RB').and. 876 & (lakonl(8:8).ne.'R').and. 877 & (lakonl(8:8).ne.'C')) then 878 call beamextscheme(yi(1,1,i),ndim,nfield,lakonl, 879 & ielprop(i),prop,field,mi) 880! 881c Bernhardi start 882 elseif((lakonl(4:6).eq.'20R').or. 883 & (lakonl(4:5).eq.'8 ').or.(lakonl(4:5).eq.'8I')) then 884c Bernhardi end 885 if(lakonl(7:8).ne.'LC') then 886 do j=1,8 887 do k=1,nfield 888 field(k,j)=0.d0 889 do l=1,8 890 field(k,j)=field(k,j)+a8(j,l)*yi(k,l,i) 891 enddo 892 enddo 893 enddo 894 else 895 do m=1,nlayer 896 jj=20*(m-1) 897 ll=8*(m-1) 898 do j=1,8 899 do k=1,nfield 900 field(k,jj+j)=0.d0 901 do l=1,8 902 field(k,jj+j)= 903 & field(k,jj+j)+a8(j,l)*yi(k,ll+l,i) 904 enddo 905 enddo 906 enddo 907 enddo 908 endif 909 elseif(lakonl(4:4).eq.'8') then 910 do j=1,8 911 do k=1,nfield 912 field(k,j)=yi(k,1,i) 913 enddo 914 enddo 915 elseif(lakonl(4:5).eq.'10') then 916 do j=1,4 917 do k=1,nfield 918 field(k,j)=0.d0 919 do l=1,4 920 field(k,j)=field(k,j)+a4(j,l)*yi(k,l,i) 921 enddo 922 enddo 923 enddo 924 elseif(lakonl(4:4).eq.'2') then 925 do j=1,20 926 do k=1,nfield 927 field(k,j)=0.d0 928 do l=1,27 929 field(k,j)=field(k,j)+a27(j,l)*yi(k,l,i) 930 enddo 931 enddo 932 enddo 933 elseif(lakonl(4:4).eq.'4') then 934 do j=1,4 935 do k=1,nfield 936 field(k,j)=yi(k,1,i) 937 enddo 938 enddo 939 elseif(lakonl(4:4).eq.'1') then 940 if(lakonl(7:8).ne.'LC') then 941 do j=1,6 942 do k=1,nfield 943 field(k,j)=0.d0 944 do l=1,9 945 field(k,j)=field(k,j)+a9(j,l)*yi(k,l,i) 946 enddo 947 enddo 948 enddo 949 else 950 do m=1,nlayer 951 jj=15*(m-1) 952 ll=6*(m-1) 953 do j=1,6 954 do k=1,nfield 955 field(k,jj+j)=0.d0 956 do l=1,6 957 field(k,jj+j)= 958 & field(k,jj+j)+a6(j,l)*yi(k,ll+l,i) 959 enddo 960 enddo 961 enddo 962 enddo 963 endif 964 else 965 do j=1,6 966 do k=1,nfield 967 field(k,j)=0.d0 968 do l=1,2 969 field(k,j)=field(k,j)+a2(j,l)*yi(k,l,i) 970 enddo 971 enddo 972 enddo 973 endif 974 endif 975! 976! determining the field values in the midside nodes 977! 978 if(lakonl(4:6).eq.'20R') then 979 if(lakonl(7:8).ne.'LC') then 980 do j=9,20 981 do k=1,nfield 982 field(k,j)=(field(k,nonei20(2,j-8))+ 983 & field(k,nonei20(3,j-8)))/2.d0 984 enddo 985 enddo 986 else 987 do m=1,nlayer 988 jj=20*(m-1) 989 do j=9,20 990 do k=1,nfield 991 field(k,jj+j)=(field(k,jj+nonei20(2,j-8)) 992 & +field(k,jj+nonei20(3,j-8)))/2.d0 993 enddo 994 enddo 995 enddo 996 endif 997 elseif(lakonl(4:5).eq.'10') then 998 do j=5,10 999 do k=1,nfield 1000 field(k,j)=(field(k,nonei10(2,j-4))+ 1001 & field(k,nonei10(3,j-4)))/2.d0 1002 enddo 1003 enddo 1004 elseif(lakonl(4:5).eq.'15') then 1005 if(lakonl(7:8).ne.'LC') then 1006 do j=7,15 1007 do k=1,nfield 1008 field(k,j)=(field(k,nonei15(2,j-6))+ 1009 & field(k,nonei15(3,j-6)))/2.d0 1010 enddo 1011 enddo 1012 else 1013 do m=1,nlayer 1014 jj=15*(m-1) 1015 do j=7,15 1016 do k=1,nfield 1017 field(k,jj+j)=(field(k,jj+nonei15(2,j-6)) 1018 & +field(k,jj+nonei15(3,j-6)))/2.d0 1019 enddo 1020 enddo 1021 enddo 1022 endif 1023 endif 1024! 1025! transferring the field values into yn 1026! 1027 if(lakonl(7:8).ne.'LC') then 1028 do j=1,nope 1029 do k=1,nfield 1030 yn(k,kon(indexe+j))=yn(k,kon(indexe+j))+ 1031 & field(k,j) 1032 enddo 1033 inum(kon(indexe+j))=inum(kon(indexe+j))+1 1034 enddo 1035 else 1036 do j=1,nope*nlayer 1037 do k=1,nfield 1038 yn(k,kon(indexe+nopeexp+j))= 1039 & yn(k,kon(indexe+nopeexp+j))+field(k,j) 1040 enddo 1041 inum(kon(indexe+nopeexp+j))= 1042 & inum(kon(indexe+nopeexp+j))+1 1043 enddo 1044 endif 1045! 1046c Bernhardi start 1047c incompatible modes elements 1048 if(lakonl(1:5).eq.'C3D8I') then 1049 do j=1,3 1050 do k=1,nfield 1051 yn(k,kon(indexe+nope+j))=0.0d0 1052 enddo 1053 enddo 1054 endif 1055c Bernhardi end 1056! 1057 enddo 1058! 1059! taking the mean of nodal contributions coming from different 1060! elements having the node in common 1061! 1062 do ij=naneigh,nbneigh 1063 i=ialnneigh(ij) 1064 if(inum(i).gt.0) then 1065 do j=1,nfield 1066 yn(j,i)=yn(j,i)/inum(i) 1067 enddo 1068 endif 1069 enddo 1070! 1071! for 1d and 2d elements only: 1072! finding the solution in the original nodes 1073! 1074 if((cflag.ne.' ').and.(cflag.ne.'E')) then 1075 call map3dto1d2d(yn,ipkon,inum,kon,lakon,nfield,nk,ne,cflag, 1076 & co,vold,force,mi) 1077 endif 1078! 1079 return 1080 end 1081