clear; clear all; clc; format long; E=200000; mu=0.33; dispbcelements=[161 0 0 0;162 0 0 0;163 0 0 0; 164 0 0 0; 165 0 0 0; 166 0 0 0]; fbcelements=[167 0 0 10000;168 0 0 10000;169 0 0 10000;170 0 0 10000; 171 0 0 10000;172 0 0 10000]; fbczeroelements(1:160,1:4)=0; for i=1:160 fbczeroelements(i,1)=i; end fbcelements=[fbczeroelements;fbcelements]; xyzofelements=[1 2.000000e+000 0.000000e+000 1.000000e+001; 1 1.000000e+000 0.000000e+000 1.000000e+001; 1 1.000000e+000 0.000000e+000 5.000000e+000; 2 2.000000e+000 0.000000e+000 1.000000e+001; 2 1.000000e+000 0.000000e+000 5.000000e+000; 2 2.000000e+000 0.000000e+000 5.000000e+000; 3 3.000000e+000 0.000000e+000 1.000000e+001; 3 2.000000e+000 0.000000e+000 1.000000e+001; 3 2.000000e+000 0.000000e+000 5.000000e+000; 4 3.000000e+000 0.000000e+000 1.000000e+001; 4 2.000000e+000 0.000000e+000 5.000000e+000; 4 3.000000e+000 0.000000e+000 5.000000e+000; 5 3.000000e+000 0.000000e+000 9.500000e+001; 5 2.000000e+000 0.000000e+000 9.500000e+001; 5 2.000000e+000 0.000000e+000 9.000000e+001; 6 3.000000e+000 0.000000e+000 9.500000e+001; 6 2.000000e+000 0.000000e+000 9.000000e+001; 6 3.000000e+000 0.000000e+000 9.000000e+001; 7 2.000000e+000 0.000000e+000 9.500000e+001; 7 1.000000e+000 0.000000e+000 9.500000e+001; 7 1.000000e+000 0.000000e+000 9.000000e+001; 8 2.000000e+000 0.000000e+000 9.500000e+001; 8 1.000000e+000 0.000000e+000 9.000000e+001; 8 2.000000e+000 0.000000e+000 9.000000e+001; 9 0.000000e+000 0.000000e+000 0.000000e+000; 9 2.000000e+000 0.000000e+000 0.000000e+000; 9 1.000000e+000 0.000000e+000 5.000000e+000; 10 2.000000e+000 0.000000e+000 0.000000e+000; 10 4.000000e+000 0.000000e+000 0.000000e+000; 10 3.000000e+000 0.000000e+000 5.000000e+000; 11 4.000000e+000 0.000000e+000 1.000000e+002; 11 2.000000e+000 0.000000e+000 1.000000e+002; 11 3.000000e+000 0.000000e+000 9.500000e+001; 12 2.000000e+000 0.000000e+000 1.000000e+002; 12 0.000000e+000 0.000000e+000 1.000000e+002; 12 1.000000e+000 0.000000e+000 9.500000e+001; 13 0.000000e+000 0.000000e+000 1.000000e+002; 13 0.000000e+000 0.000000e+000 8.571429e+001; 13 1.000000e+000 0.000000e+000 9.500000e+001; 14 0.000000e+000 0.000000e+000 8.571429e+001; 14 0.000000e+000 0.000000e+000 7.142857e+001; 14 4.000000e+000 0.000000e+000 8.571429e+001; 15 0.000000e+000 0.000000e+000 7.142857e+001; 15 0.000000e+000 0.000000e+000 5.714286e+001; 15 4.000000e+000 0.000000e+000 7.142857e+001; 16 0.000000e+000 0.000000e+000 5.714286e+001; 16 0.000000e+000 0.000000e+000 4.285714e+001; 16 4.000000e+000 0.000000e+000 5.714286e+001; 17 0.000000e+000 0.000000e+000 4.285714e+001; 17 0.000000e+000 0.000000e+000 2.857143e+001; 17 4.000000e+000 0.000000e+000 4.285714e+001; 18 0.000000e+000 0.000000e+000 2.857143e+001; 18 0.000000e+000 0.000000e+000 1.428571e+001; 18 4.000000e+000 0.000000e+000 1.428571e+001; 19 0.000000e+000 0.000000e+000 1.428571e+001; 19 0.000000e+000 0.000000e+000 0.000000e+000; 19 1.000000e+000 0.000000e+000 5.000000e+000; 20 4.000000e+000 0.000000e+000 0.000000e+000; 20 4.000000e+000 0.000000e+000 1.428571e+001; 20 3.000000e+000 0.000000e+000 5.000000e+000; 21 4.000000e+000 0.000000e+000 1.428571e+001; 21 4.000000e+000 0.000000e+000 2.857143e+001; 21 0.000000e+000 0.000000e+000 2.857143e+001; 22 4.000000e+000 0.000000e+000 2.857143e+001; 22 4.000000e+000 0.000000e+000 4.285714e+001; 22 0.000000e+000 0.000000e+000 2.857143e+001; 23 4.000000e+000 0.000000e+000 4.285714e+001; 23 4.000000e+000 0.000000e+000 5.714286e+001; 23 0.000000e+000 0.000000e+000 4.285714e+001; 24 4.000000e+000 0.000000e+000 5.714286e+001; 24 4.000000e+000 0.000000e+000 7.142857e+001; 24 0.000000e+000 0.000000e+000 5.714286e+001; 25 4.000000e+000 0.000000e+000 7.142857e+001; 25 4.000000e+000 0.000000e+000 8.571429e+001; 25 0.000000e+000 0.000000e+000 7.142857e+001; 26 4.000000e+000 0.000000e+000 8.571429e+001; 26 4.000000e+000 0.000000e+000 1.000000e+002; 26 3.000000e+000 0.000000e+000 9.500000e+001; 27 1.000000e+000 0.000000e+000 1.000000e+001; 27 2.000000e+000 0.000000e+000 1.000000e+001; 27 0.000000e+000 0.000000e+000 1.428571e+001; 28 1.000000e+000 0.000000e+000 5.000000e+000; 28 1.000000e+000 0.000000e+000 1.000000e+001; 28 0.000000e+000 0.000000e+000 1.428571e+001; 29 2.000000e+000 0.000000e+000 5.000000e+000; 29 1.000000e+000 0.000000e+000 5.000000e+000; 29 2.000000e+000 0.000000e+000 0.000000e+000; 30 2.000000e+000 0.000000e+000 1.000000e+001; 30 3.000000e+000 0.000000e+000 1.000000e+001; 30 4.000000e+000 0.000000e+000 1.428571e+001; 31 3.000000e+000 0.000000e+000 5.000000e+000; 31 2.000000e+000 0.000000e+000 5.000000e+000; 31 2.000000e+000 0.000000e+000 0.000000e+000; 32 3.000000e+000 0.000000e+000 1.000000e+001; 32 3.000000e+000 0.000000e+000 5.000000e+000; 32 4.000000e+000 0.000000e+000 1.428571e+001; 33 2.000000e+000 0.000000e+000 9.500000e+001; 33 3.000000e+000 0.000000e+000 9.500000e+001; 33 2.000000e+000 0.000000e+000 1.000000e+002; 34 3.000000e+000 0.000000e+000 9.000000e+001; 34 2.000000e+000 0.000000e+000 9.000000e+001; 34 4.000000e+000 0.000000e+000 8.571429e+001; 35 3.000000e+000 0.000000e+000 9.500000e+001; 35 3.000000e+000 0.000000e+000 9.000000e+001; 35 4.000000e+000 0.000000e+000 8.571429e+001; 36 1.000000e+000 0.000000e+000 9.500000e+001; 36 2.000000e+000 0.000000e+000 9.500000e+001; 36 2.000000e+000 0.000000e+000 1.000000e+002; 37 1.000000e+000 0.000000e+000 9.000000e+001; 37 1.000000e+000 0.000000e+000 9.500000e+001; 37 0.000000e+000 0.000000e+000 8.571429e+001; 38 2.000000e+000 0.000000e+000 9.000000e+001; 38 1.000000e+000 0.000000e+000 9.000000e+001; 38 0.000000e+000 0.000000e+000 8.571429e+001; 39 0.000000e+000 0.000000e+000 8.571429e+001; 39 4.000000e+000 0.000000e+000 8.571429e+001; 39 2.000000e+000 0.000000e+000 9.000000e+001; 40 4.000000e+000 0.000000e+000 1.428571e+001; 40 0.000000e+000 0.000000e+000 1.428571e+001; 40 2.000000e+000 0.000000e+000 1.000000e+001; 41 4.000000e+000 2.000000e+000 1.000000e+001; 41 4.000000e+000 1.000000e+000 1.000000e+001; 41 4.000000e+000 1.000000e+000 5.000000e+000; 42 4.000000e+000 2.000000e+000 1.000000e+001; 42 4.000000e+000 1.000000e+000 5.000000e+000; 42 4.000000e+000 2.000000e+000 5.000000e+000; 43 4.000000e+000 3.000000e+000 1.000000e+001; 43 4.000000e+000 2.000000e+000 1.000000e+001; 43 4.000000e+000 2.000000e+000 5.000000e+000; 44 4.000000e+000 3.000000e+000 1.000000e+001; 44 4.000000e+000 2.000000e+000 5.000000e+000; 44 4.000000e+000 3.000000e+000 5.000000e+000; 45 4.000000e+000 3.000000e+000 9.500000e+001; 45 4.000000e+000 2.000000e+000 9.500000e+001; 45 4.000000e+000 2.000000e+000 9.000000e+001; 46 4.000000e+000 3.000000e+000 9.500000e+001; 46 4.000000e+000 2.000000e+000 9.000000e+001; 46 4.000000e+000 3.000000e+000 9.000000e+001; 47 4.000000e+000 2.000000e+000 9.500000e+001; 47 4.000000e+000 1.000000e+000 9.500000e+001; 47 4.000000e+000 1.000000e+000 9.000000e+001; 48 4.000000e+000 2.000000e+000 9.500000e+001; 48 4.000000e+000 1.000000e+000 9.000000e+001; 48 4.000000e+000 2.000000e+000 9.000000e+001; 49 4.000000e+000 0.000000e+000 0.000000e+000; 49 4.000000e+000 2.000000e+000 0.000000e+000; 49 4.000000e+000 1.000000e+000 5.000000e+000; 50 4.000000e+000 2.000000e+000 0.000000e+000; 50 4.000000e+000 4.000000e+000 0.000000e+000; 50 4.000000e+000 3.000000e+000 5.000000e+000; 51 4.000000e+000 4.000000e+000 1.000000e+002; 51 4.000000e+000 2.000000e+000 1.000000e+002; 51 4.000000e+000 3.000000e+000 9.500000e+001; 52 4.000000e+000 2.000000e+000 1.000000e+002; 52 4.000000e+000 0.000000e+000 1.000000e+002; 52 4.000000e+000 1.000000e+000 9.500000e+001; 53 4.000000e+000 0.000000e+000 1.000000e+002; 53 4.000000e+000 0.000000e+000 8.571429e+001; 53 4.000000e+000 1.000000e+000 9.500000e+001; 54 4.000000e+000 0.000000e+000 8.571429e+001; 54 4.000000e+000 0.000000e+000 7.142857e+001; 54 4.000000e+000 4.000000e+000 8.571429e+001; 55 4.000000e+000 0.000000e+000 7.142857e+001; 55 4.000000e+000 0.000000e+000 5.714286e+001; 55 4.000000e+000 4.000000e+000 7.142857e+001; 56 4.000000e+000 0.000000e+000 5.714286e+001; 56 4.000000e+000 0.000000e+000 4.285714e+001; 56 4.000000e+000 4.000000e+000 5.714286e+001; 57 4.000000e+000 0.000000e+000 4.285714e+001; 57 4.000000e+000 0.000000e+000 2.857143e+001; 57 4.000000e+000 4.000000e+000 4.285714e+001; 58 4.000000e+000 0.000000e+000 2.857143e+001; 58 4.000000e+000 0.000000e+000 1.428571e+001; 58 4.000000e+000 4.000000e+000 1.428571e+001; 59 4.000000e+000 0.000000e+000 1.428571e+001; 59 4.000000e+000 0.000000e+000 0.000000e+000; 59 4.000000e+000 1.000000e+000 5.000000e+000; 60 4.000000e+000 4.000000e+000 0.000000e+000; 60 4.000000e+000 4.000000e+000 1.428571e+001; 60 4.000000e+000 3.000000e+000 5.000000e+000; 61 4.000000e+000 4.000000e+000 1.428571e+001; 61 4.000000e+000 4.000000e+000 2.857143e+001; 61 4.000000e+000 0.000000e+000 2.857143e+001; 62 4.000000e+000 4.000000e+000 2.857143e+001; 62 4.000000e+000 4.000000e+000 4.285714e+001; 62 4.000000e+000 0.000000e+000 2.857143e+001; 63 4.000000e+000 4.000000e+000 4.285714e+001; 63 4.000000e+000 4.000000e+000 5.714286e+001; 63 4.000000e+000 0.000000e+000 4.285714e+001; 64 4.000000e+000 4.000000e+000 5.714286e+001; 64 4.000000e+000 4.000000e+000 7.142857e+001; 64 4.000000e+000 0.000000e+000 5.714286e+001; 65 4.000000e+000 4.000000e+000 7.142857e+001; 65 4.000000e+000 4.000000e+000 8.571429e+001; 65 4.000000e+000 0.000000e+000 7.142857e+001; 66 4.000000e+000 4.000000e+000 8.571429e+001; 66 4.000000e+000 4.000000e+000 1.000000e+002; 66 4.000000e+000 3.000000e+000 9.500000e+001; 67 4.000000e+000 1.000000e+000 1.000000e+001; 67 4.000000e+000 2.000000e+000 1.000000e+001; 67 4.000000e+000 0.000000e+000 1.428571e+001; 68 4.000000e+000 1.000000e+000 5.000000e+000; 68 4.000000e+000 1.000000e+000 1.000000e+001; 68 4.000000e+000 0.000000e+000 1.428571e+001; 69 4.000000e+000 2.000000e+000 5.000000e+000; 69 4.000000e+000 1.000000e+000 5.000000e+000; 69 4.000000e+000 2.000000e+000 0.000000e+000; 70 4.000000e+000 2.000000e+000 1.000000e+001; 70 4.000000e+000 3.000000e+000 1.000000e+001; 70 4.000000e+000 4.000000e+000 1.428571e+001; 71 4.000000e+000 3.000000e+000 5.000000e+000; 71 4.000000e+000 2.000000e+000 5.000000e+000; 71 4.000000e+000 2.000000e+000 0.000000e+000; 72 4.000000e+000 3.000000e+000 1.000000e+001; 72 4.000000e+000 3.000000e+000 5.000000e+000; 72 4.000000e+000 4.000000e+000 1.428571e+001; 73 4.000000e+000 2.000000e+000 9.500000e+001; 73 4.000000e+000 3.000000e+000 9.500000e+001; 73 4.000000e+000 2.000000e+000 1.000000e+002; 74 4.000000e+000 3.000000e+000 9.000000e+001; 74 4.000000e+000 2.000000e+000 9.000000e+001; 74 4.000000e+000 4.000000e+000 8.571429e+001; 75 4.000000e+000 3.000000e+000 9.500000e+001; 75 4.000000e+000 3.000000e+000 9.000000e+001; 75 4.000000e+000 4.000000e+000 8.571429e+001; 76 4.000000e+000 1.000000e+000 9.500000e+001; 76 4.000000e+000 2.000000e+000 9.500000e+001; 76 4.000000e+000 2.000000e+000 1.000000e+002; 77 4.000000e+000 1.000000e+000 9.000000e+001; 77 4.000000e+000 1.000000e+000 9.500000e+001; 77 4.000000e+000 0.000000e+000 8.571429e+001; 78 4.000000e+000 2.000000e+000 9.000000e+001; 78 4.000000e+000 1.000000e+000 9.000000e+001; 78 4.000000e+000 0.000000e+000 8.571429e+001; 79 4.000000e+000 0.000000e+000 8.571429e+001; 79 4.000000e+000 4.000000e+000 8.571429e+001; 79 4.000000e+000 2.000000e+000 9.000000e+001; 80 4.000000e+000 4.000000e+000 1.428571e+001; 80 4.000000e+000 0.000000e+000 1.428571e+001; 80 4.000000e+000 2.000000e+000 1.000000e+001; 81 2.000000e+000 4.000000e+000 1.000000e+001; 81 3.000000e+000 4.000000e+000 1.000000e+001; 81 3.000000e+000 4.000000e+000 5.000000e+000; 82 2.000000e+000 4.000000e+000 1.000000e+001; 82 3.000000e+000 4.000000e+000 5.000000e+000; 82 2.000000e+000 4.000000e+000 5.000000e+000; 83 1.000000e+000 4.000000e+000 1.000000e+001; 83 2.000000e+000 4.000000e+000 1.000000e+001; 83 2.000000e+000 4.000000e+000 5.000000e+000; 84 1.000000e+000 4.000000e+000 1.000000e+001; 84 2.000000e+000 4.000000e+000 5.000000e+000; 84 1.000000e+000 4.000000e+000 5.000000e+000; 85 1.000000e+000 4.000000e+000 9.500000e+001; 85 2.000000e+000 4.000000e+000 9.500000e+001; 85 2.000000e+000 4.000000e+000 9.000000e+001; 86 1.000000e+000 4.000000e+000 9.500000e+001; 86 2.000000e+000 4.000000e+000 9.000000e+001; 86 1.000000e+000 4.000000e+000 9.000000e+001; 87 2.000000e+000 4.000000e+000 9.500000e+001; 87 3.000000e+000 4.000000e+000 9.500000e+001; 87 3.000000e+000 4.000000e+000 9.000000e+001; 88 2.000000e+000 4.000000e+000 9.500000e+001; 88 3.000000e+000 4.000000e+000 9.000000e+001; 88 2.000000e+000 4.000000e+000 9.000000e+001; 89 4.000000e+000 4.000000e+000 0.000000e+000; 89 2.000000e+000 4.000000e+000 0.000000e+000; 89 3.000000e+000 4.000000e+000 5.000000e+000; 90 2.000000e+000 4.000000e+000 0.000000e+000; 90 0.000000e+000 4.000000e+000 0.000000e+000; 90 1.000000e+000 4.000000e+000 5.000000e+000; 91 0.000000e+000 4.000000e+000 1.000000e+002; 91 2.000000e+000 4.000000e+000 1.000000e+002; 91 1.000000e+000 4.000000e+000 9.500000e+001; 92 2.000000e+000 4.000000e+000 1.000000e+002; 92 4.000000e+000 4.000000e+000 1.000000e+002; 92 3.000000e+000 4.000000e+000 9.500000e+001; 93 4.000000e+000 4.000000e+000 1.000000e+002; 93 4.000000e+000 4.000000e+000 8.571429e+001; 93 3.000000e+000 4.000000e+000 9.500000e+001; 94 4.000000e+000 4.000000e+000 8.571429e+001; 94 4.000000e+000 4.000000e+000 7.142857e+001; 94 0.000000e+000 4.000000e+000 8.571429e+001; 95 4.000000e+000 4.000000e+000 7.142857e+001; 95 4.000000e+000 4.000000e+000 5.714286e+001; 95 0.000000e+000 4.000000e+000 7.142857e+001; 96 4.000000e+000 4.000000e+000 5.714286e+001; 96 4.000000e+000 4.000000e+000 4.285714e+001; 96 0.000000e+000 4.000000e+000 5.714286e+001; 97 4.000000e+000 4.000000e+000 4.285714e+001; 97 4.000000e+000 4.000000e+000 2.857143e+001; 97 0.000000e+000 4.000000e+000 4.285714e+001; 98 4.000000e+000 4.000000e+000 2.857143e+001; 98 4.000000e+000 4.000000e+000 1.428571e+001; 98 0.000000e+000 4.000000e+000 1.428571e+001; 99 4.000000e+000 4.000000e+000 1.428571e+001; 99 4.000000e+000 4.000000e+000 0.000000e+000; 99 3.000000e+000 4.000000e+000 5.000000e+000; 100 0.000000e+000 4.000000e+000 0.000000e+000; 100 0.000000e+000 4.000000e+000 1.428571e+001; 100 1.000000e+000 4.000000e+000 5.000000e+000; 101 0.000000e+000 4.000000e+000 1.428571e+001; 101 0.000000e+000 4.000000e+000 2.857143e+001; 101 4.000000e+000 4.000000e+000 2.857143e+001; 102 0.000000e+000 4.000000e+000 2.857143e+001; 102 0.000000e+000 4.000000e+000 4.285714e+001; 102 4.000000e+000 4.000000e+000 2.857143e+001; 103 0.000000e+000 4.000000e+000 4.285714e+001; 103 0.000000e+000 4.000000e+000 5.714286e+001; 103 4.000000e+000 4.000000e+000 4.285714e+001; 104 0.000000e+000 4.000000e+000 5.714286e+001; 104 0.000000e+000 4.000000e+000 7.142857e+001; 104 4.000000e+000 4.000000e+000 5.714286e+001; 105 0.000000e+000 4.000000e+000 7.142857e+001; 105 0.000000e+000 4.000000e+000 8.571429e+001; 105 4.000000e+000 4.000000e+000 7.142857e+001; 106 0.000000e+000 4.000000e+000 8.571429e+001; 106 0.000000e+000 4.000000e+000 1.000000e+002; 106 1.000000e+000 4.000000e+000 9.500000e+001; 107 3.000000e+000 4.000000e+000 1.000000e+001; 107 2.000000e+000 4.000000e+000 1.000000e+001; 107 4.000000e+000 4.000000e+000 1.428571e+001; 108 3.000000e+000 4.000000e+000 5.000000e+000; 108 3.000000e+000 4.000000e+000 1.000000e+001; 108 4.000000e+000 4.000000e+000 1.428571e+001; 109 2.000000e+000 4.000000e+000 5.000000e+000; 109 3.000000e+000 4.000000e+000 5.000000e+000; 109 2.000000e+000 4.000000e+000 0.000000e+000; 110 2.000000e+000 4.000000e+000 1.000000e+001; 110 1.000000e+000 4.000000e+000 1.000000e+001; 110 0.000000e+000 4.000000e+000 1.428571e+001; 111 1.000000e+000 4.000000e+000 5.000000e+000; 111 2.000000e+000 4.000000e+000 5.000000e+000; 111 2.000000e+000 4.000000e+000 0.000000e+000; 112 1.000000e+000 4.000000e+000 1.000000e+001; 112 1.000000e+000 4.000000e+000 5.000000e+000; 112 0.000000e+000 4.000000e+000 1.428571e+001; 113 2.000000e+000 4.000000e+000 9.500000e+001; 113 1.000000e+000 4.000000e+000 9.500000e+001; 113 2.000000e+000 4.000000e+000 1.000000e+002; 114 1.000000e+000 4.000000e+000 9.000000e+001; 114 2.000000e+000 4.000000e+000 9.000000e+001; 114 0.000000e+000 4.000000e+000 8.571429e+001; 115 1.000000e+000 4.000000e+000 9.500000e+001; 115 1.000000e+000 4.000000e+000 9.000000e+001; 115 0.000000e+000 4.000000e+000 8.571429e+001; 116 3.000000e+000 4.000000e+000 9.500000e+001; 116 2.000000e+000 4.000000e+000 9.500000e+001; 116 2.000000e+000 4.000000e+000 1.000000e+002; 117 3.000000e+000 4.000000e+000 9.000000e+001; 117 3.000000e+000 4.000000e+000 9.500000e+001; 117 4.000000e+000 4.000000e+000 8.571429e+001; 118 2.000000e+000 4.000000e+000 9.000000e+001; 118 3.000000e+000 4.000000e+000 9.000000e+001; 118 4.000000e+000 4.000000e+000 8.571429e+001; 119 4.000000e+000 4.000000e+000 8.571429e+001; 119 0.000000e+000 4.000000e+000 8.571429e+001; 119 2.000000e+000 4.000000e+000 9.000000e+001; 120 0.000000e+000 4.000000e+000 1.428571e+001; 120 4.000000e+000 4.000000e+000 1.428571e+001; 120 2.000000e+000 4.000000e+000 1.000000e+001; 121 0.000000e+000 2.000000e+000 1.000000e+001; 121 0.000000e+000 3.000000e+000 1.000000e+001; 121 0.000000e+000 3.000000e+000 5.000000e+000; 122 0.000000e+000 2.000000e+000 1.000000e+001; 122 0.000000e+000 3.000000e+000 5.000000e+000; 122 0.000000e+000 2.000000e+000 5.000000e+000; 123 0.000000e+000 1.000000e+000 1.000000e+001; 123 0.000000e+000 2.000000e+000 1.000000e+001; 123 0.000000e+000 2.000000e+000 5.000000e+000; 124 0.000000e+000 1.000000e+000 1.000000e+001; 124 0.000000e+000 2.000000e+000 5.000000e+000; 124 0.000000e+000 1.000000e+000 5.000000e+000; 125 0.000000e+000 1.000000e+000 9.500000e+001; 125 0.000000e+000 2.000000e+000 9.500000e+001; 125 0.000000e+000 2.000000e+000 9.000000e+001; 126 0.000000e+000 1.000000e+000 9.500000e+001; 126 0.000000e+000 2.000000e+000 9.000000e+001; 126 0.000000e+000 1.000000e+000 9.000000e+001; 127 0.000000e+000 2.000000e+000 9.500000e+001; 127 0.000000e+000 3.000000e+000 9.500000e+001; 127 0.000000e+000 3.000000e+000 9.000000e+001; 128 0.000000e+000 2.000000e+000 9.500000e+001; 128 0.000000e+000 3.000000e+000 9.000000e+001; 128 0.000000e+000 2.000000e+000 9.000000e+001; 129 0.000000e+000 4.000000e+000 0.000000e+000; 129 0.000000e+000 2.000000e+000 0.000000e+000; 129 0.000000e+000 3.000000e+000 5.000000e+000; 130 0.000000e+000 2.000000e+000 0.000000e+000; 130 0.000000e+000 0.000000e+000 0.000000e+000; 130 0.000000e+000 1.000000e+000 5.000000e+000; 131 0.000000e+000 0.000000e+000 1.000000e+002; 131 0.000000e+000 2.000000e+000 1.000000e+002; 131 0.000000e+000 1.000000e+000 9.500000e+001; 132 0.000000e+000 2.000000e+000 1.000000e+002; 132 0.000000e+000 4.000000e+000 1.000000e+002; 132 0.000000e+000 3.000000e+000 9.500000e+001; 133 0.000000e+000 0.000000e+000 0.000000e+000; 133 0.000000e+000 0.000000e+000 1.428571e+001; 133 0.000000e+000 1.000000e+000 5.000000e+000; 134 0.000000e+000 0.000000e+000 1.428571e+001; 134 0.000000e+000 0.000000e+000 2.857143e+001; 134 0.000000e+000 4.000000e+000 2.857143e+001; 135 0.000000e+000 0.000000e+000 2.857143e+001; 135 0.000000e+000 0.000000e+000 4.285714e+001; 135 0.000000e+000 4.000000e+000 2.857143e+001; 136 0.000000e+000 0.000000e+000 4.285714e+001; 136 0.000000e+000 0.000000e+000 5.714286e+001; 136 0.000000e+000 4.000000e+000 4.285714e+001; 137 0.000000e+000 0.000000e+000 5.714286e+001; 137 0.000000e+000 0.000000e+000 7.142857e+001; 137 0.000000e+000 4.000000e+000 7.142857e+001; 138 0.000000e+000 0.000000e+000 7.142857e+001; 138 0.000000e+000 0.000000e+000 8.571429e+001; 138 0.000000e+000 4.000000e+000 8.571429e+001; 139 0.000000e+000 0.000000e+000 8.571429e+001; 139 0.000000e+000 0.000000e+000 1.000000e+002; 139 0.000000e+000 1.000000e+000 9.500000e+001; 140 0.000000e+000 4.000000e+000 1.000000e+002; 140 0.000000e+000 4.000000e+000 8.571429e+001; 140 0.000000e+000 3.000000e+000 9.500000e+001; 141 0.000000e+000 4.000000e+000 8.571429e+001; 141 0.000000e+000 4.000000e+000 7.142857e+001; 141 0.000000e+000 0.000000e+000 7.142857e+001; 142 0.000000e+000 4.000000e+000 7.142857e+001; 142 0.000000e+000 4.000000e+000 5.714286e+001; 142 0.000000e+000 0.000000e+000 5.714286e+001; 143 0.000000e+000 4.000000e+000 5.714286e+001; 143 0.000000e+000 4.000000e+000 4.285714e+001; 143 0.000000e+000 0.000000e+000 5.714286e+001; 144 0.000000e+000 4.000000e+000 4.285714e+001; 144 0.000000e+000 4.000000e+000 2.857143e+001; 144 0.000000e+000 0.000000e+000 4.285714e+001; 145 0.000000e+000 4.000000e+000 2.857143e+001; 145 0.000000e+000 4.000000e+000 1.428571e+001; 145 0.000000e+000 0.000000e+000 1.428571e+001; 146 0.000000e+000 4.000000e+000 1.428571e+001; 146 0.000000e+000 4.000000e+000 0.000000e+000; 146 0.000000e+000 3.000000e+000 5.000000e+000; 147 0.000000e+000 3.000000e+000 1.000000e+001; 147 0.000000e+000 2.000000e+000 1.000000e+001; 147 0.000000e+000 4.000000e+000 1.428571e+001; 148 0.000000e+000 3.000000e+000 5.000000e+000; 148 0.000000e+000 3.000000e+000 1.000000e+001; 148 0.000000e+000 4.000000e+000 1.428571e+001; 149 0.000000e+000 2.000000e+000 5.000000e+000; 149 0.000000e+000 3.000000e+000 5.000000e+000; 149 0.000000e+000 2.000000e+000 0.000000e+000; 150 0.000000e+000 2.000000e+000 1.000000e+001; 150 0.000000e+000 1.000000e+000 1.000000e+001; 150 0.000000e+000 0.000000e+000 1.428571e+001; 151 0.000000e+000 1.000000e+000 5.000000e+000; 151 0.000000e+000 2.000000e+000 5.000000e+000; 151 0.000000e+000 2.000000e+000 0.000000e+000; 152 0.000000e+000 1.000000e+000 1.000000e+001; 152 0.000000e+000 1.000000e+000 5.000000e+000; 152 0.000000e+000 0.000000e+000 1.428571e+001; 153 0.000000e+000 2.000000e+000 9.500000e+001; 153 0.000000e+000 1.000000e+000 9.500000e+001; 153 0.000000e+000 2.000000e+000 1.000000e+002; 154 0.000000e+000 1.000000e+000 9.000000e+001; 154 0.000000e+000 2.000000e+000 9.000000e+001; 154 0.000000e+000 0.000000e+000 8.571429e+001; 155 0.000000e+000 1.000000e+000 9.500000e+001; 155 0.000000e+000 1.000000e+000 9.000000e+001; 155 0.000000e+000 0.000000e+000 8.571429e+001; 156 0.000000e+000 3.000000e+000 9.500000e+001; 156 0.000000e+000 2.000000e+000 9.500000e+001; 156 0.000000e+000 2.000000e+000 1.000000e+002; 157 0.000000e+000 3.000000e+000 9.000000e+001; 157 0.000000e+000 3.000000e+000 9.500000e+001; 157 0.000000e+000 4.000000e+000 8.571429e+001; 158 0.000000e+000 2.000000e+000 9.000000e+001; 158 0.000000e+000 3.000000e+000 9.000000e+001; 158 0.000000e+000 4.000000e+000 8.571429e+001; 159 0.000000e+000 4.000000e+000 8.571429e+001; 159 0.000000e+000 0.000000e+000 8.571429e+001; 159 0.000000e+000 2.000000e+000 9.000000e+001; 160 0.000000e+000 0.000000e+000 1.428571e+001; 160 0.000000e+000 4.000000e+000 1.428571e+001; 160 0.000000e+000 2.000000e+000 1.000000e+001; 161 4.000000e+000 0.000000e+000 0.000000e+000; 161 2.000000e+000 0.000000e+000 0.000000e+000; 161 4.000000e+000 2.000000e+000 0.000000e+000; 162 2.000000e+000 0.000000e+000 0.000000e+000; 162 0.000000e+000 0.000000e+000 0.000000e+000; 162 0.000000e+000 2.000000e+000 0.000000e+000; 163 4.000000e+000 4.000000e+000 0.000000e+000; 163 4.000000e+000 2.000000e+000 0.000000e+000; 163 2.000000e+000 4.000000e+000 0.000000e+000; 164 0.000000e+000 4.000000e+000 0.000000e+000; 164 2.000000e+000 4.000000e+000 0.000000e+000; 164 0.000000e+000 2.000000e+000 0.000000e+000; 165 4.000000e+000 2.000000e+000 0.000000e+000; 165 2.000000e+000 0.000000e+000 0.000000e+000; 165 0.000000e+000 2.000000e+000 0.000000e+000; 166 2.000000e+000 4.000000e+000 0.000000e+000; 166 4.000000e+000 2.000000e+000 0.000000e+000; 166 0.000000e+000 2.000000e+000 0.000000e+000; 167 0.000000e+000 0.000000e+000 1.000000e+002; 167 2.000000e+000 0.000000e+000 1.000000e+002; 167 0.000000e+000 2.000000e+000 1.000000e+002; 168 2.000000e+000 0.000000e+000 1.000000e+002; 168 4.000000e+000 0.000000e+000 1.000000e+002; 168 4.000000e+000 2.000000e+000 1.000000e+002; 169 4.000000e+000 2.000000e+000 1.000000e+002; 169 4.000000e+000 4.000000e+000 1.000000e+002; 169 2.000000e+000 4.000000e+000 1.000000e+002; 170 2.000000e+000 4.000000e+000 1.000000e+002; 170 0.000000e+000 4.000000e+000 1.000000e+002; 170 0.000000e+000 2.000000e+000 1.000000e+002; 171 0.000000e+000 2.000000e+000 1.000000e+002; 171 2.000000e+000 0.000000e+000 1.000000e+002; 171 4.000000e+000 2.000000e+000 1.000000e+002; 172 4.000000e+000 2.000000e+000 1.000000e+002; 172 2.000000e+000 4.000000e+000 1.000000e+002; 172 0.000000e+000 2.000000e+000 1.000000e+002]; G=E/(2*(1+mu)); C=1/(16*pi*G*(1-mu)); C1=3-4*mu; C2=1/(8*pi*(1-mu)); C3=1-2*mu; n=2; totalsize=size(xyzofelements); nelement=totalsize(1,1)/3; sizedispbcelements=size(dispbcelements); ndispbcelements=sizedispbcelements(1,1); sizefbcelements=size(fbcelements); nfbcelements=sizefbcelements(1,1); nx(1:nelement)=0; ny(1:nelement)=0; nz(1:nelement)=0; xbar(1:nelement)=0; ybar(1:nelement)=0; zbar(1:nelement)=0; J(1:nelement)=0; for i=1:nelement x1=xyzofelements(3*i-2,2); y1=xyzofelements(3*i-2,3); z1=xyzofelements(3*i-2,4); x2=xyzofelements(3*i-1,2); y2=xyzofelements(3*i-1,3); z2=xyzofelements(3*i-1,4); x3=xyzofelements(3*i,2); y3=xyzofelements(3*i,3); z3=xyzofelements(3*i,4); d=(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1))^2+((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1))^2+((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1))^2)^0.5; nx(i)=((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1))/d; ny(i)=((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1))/d; nz(i)=((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1))/d; alpha=sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2); beta=sqrt((x2-x3)^2+(y2-y3)^2+(z2-z3)^2); gamma=sqrt((x3-x1)^2+(y3-y1)^2+(z3-z1)^2); sigma=(alpha+beta+gamma)/2; J(i)=2*sqrt(sigma*(sigma-alpha)*(sigma-beta)*(sigma-gamma)); if (abs(nz(i))>=1/sqrt(3)) xbar(i)=(x2-x1)*0.25+(x3-x1)*0.5+x1; ybar(i)=(y2-y1)*0.25+(y3-y1)*0.5+y1; zbar(i)=-(nz(i))^(-1)*(nx(i)*(xbar(i)-x1)+ny(i)*(ybar(i)-y1))+z1; elseif (abs(nz(i))<1/sqrt(3) && abs(ny(i))>=1/sqrt(3)) xbar(i)=(x2-x1)*0.25+(x3-x1)*0.5+x1; zbar(i)=(z2-z1)*0.25+(z3-z1)*0.5+z1; ybar(i)=-(ny(i))^(-1)*(nx(i)*(xbar(i)-x1)+nz(i)*(zbar(i)-z1))+y1; elseif (abs(nz(i))<1/sqrt(3) && abs(ny(i))<1/sqrt(3)) ybar(i)=(y2-y1)*0.25+(y3-y1)*0.5+y1; zbar(i)=(z2-z1)*0.25+(z3-z1)*0.5+z1; xbar(i)=-(nx(i))^(-1)*(ny(i)*(ybar(i)-y1)+nz(i)*(zbar(i)-z1))+x1; end end t=[1/4+1/(4*sqrt(3)) 1/4+1/(4*sqrt(3)) 1/4-1/(4*sqrt(3)) 1/4-1/(4*sqrt(3)) 3/4+1/(4*sqrt(3)) 3/4+1/(4*sqrt(3)) 3/4-1/(4*sqrt(3)) 3/4-1/(4*sqrt(3)) 3/4+1/(4*sqrt(3)) 3/4+1/(4*sqrt(3)) 3/4-1/(4*sqrt(3)) 3/4-1/(4*sqrt(3)) 1/4+1/(4*sqrt(3)) 1/4+1/(4*sqrt(3)) 1/4-1/(4*sqrt(3)) 1/4-1/(4*sqrt(3))]; v=[1/4+1/(4*sqrt(3)) 1/4-1/(4*sqrt(3)) 1/4+1/(4*sqrt(3)) 1/4-1/(4*sqrt(3)) 1/4+1/(4*sqrt(3)) 1/4-1/(4*sqrt(3)) 1/4+1/(4*sqrt(3)) 1/4-1/(4*sqrt(3)) 3/4+1/(4*sqrt(3)) 3/4-1/(4*sqrt(3)) 3/4+1/(4*sqrt(3)) 3/4-1/(4*sqrt(3)) 3/4+1/(4*sqrt(3)) 3/4-1/(4*sqrt(3)) 3/4+1/(4*sqrt(3)) 3/4-1/(4*sqrt(3))]; u=t.*(1-v); K=zeros(3*i); F(1:3*i,1)=0; for i=1:nelement for j=1:nfbcelements if (abs(nz(fbcelements(j,1)))>=1/sqrt(3)) x1=xyzofelements(3*fbcelements(j,1)-2,2); y1=xyzofelements(3*fbcelements(j,1)-2,3); z1=xyzofelements(3*fbcelements(j,1)-2,4); x2=xyzofelements(3*fbcelements(j,1)-1,2); y2=xyzofelements(3*fbcelements(j,1)-1,3); z2=xyzofelements(3*fbcelements(j,1)-1,4); x3=xyzofelements(3*fbcelements(j,1),2); y3=xyzofelements(3*fbcelements(j,1),3); z3=xyzofelements(3*fbcelements(j,1),4); x=(x2-x1)*u+(x3-x1)*v+x1; y=(y2-y1)*u+(y3-y1)*v+y1; z=-(nz(fbcelements(j,1)))^(-1)*(nx(fbcelements(j,1))*(x-x1)+ny(fbcelements(j,1))*(y-y1))+z1; elseif (abs(nz(fbcelements(j,1)))<1/sqrt(3) && abs(ny(fbcelements(j,1)))>=1/sqrt(3)) x1=xyzofelements(3*fbcelements(j,1)-2,2); y1=xyzofelements(3*fbcelements(j,1)-2,3); z1=xyzofelements(3*fbcelements(j,1)-2,4); x2=xyzofelements(3*fbcelements(j,1)-1,2); y2=xyzofelements(3*fbcelements(j,1)-1,3); z2=xyzofelements(3*fbcelements(j,1)-1,4); x3=xyzofelements(3*fbcelements(j,1),2); y3=xyzofelements(3*fbcelements(j,1),3); z3=xyzofelements(3*fbcelements(j,1),4); x=(x2-x1)*u+(x3-x1)*v+x1; z=(z2-z1)*u+(z3-z1)*v+z1; y=-(ny(fbcelements(j,1)))^(-1)*(nx(fbcelements(j,1))*(x-x1)+nz(fbcelements(j,1))*(z-z1))+y1; elseif (abs(nz(fbcelements(j,1)))<1/sqrt(3) && abs(ny(fbcelements(j,1)))<1/sqrt(3)) x1=xyzofelements(3*fbcelements(j,1)-2,2); y1=xyzofelements(3*fbcelements(j,1)-2,3); z1=xyzofelements(3*fbcelements(j,1)-2,4); x2=xyzofelements(3*fbcelements(j,1)-1,2); y2=xyzofelements(3*fbcelements(j,1)-1,3); z2=xyzofelements(3*fbcelements(j,1)-1,4); x3=xyzofelements(3*fbcelements(j,1),2); y3=xyzofelements(3*fbcelements(j,1),3); z3=xyzofelements(3*fbcelements(j,1),4); y=(y2-y1)*u+(y3-y1)*v+y1; z=(z2-z1)*u+(z3-z1)*v+z1; x=-(nx(fbcelements(j,1)))^(-1)*(ny(fbcelements(j,1))*(y-y1)+nz(fbcelements(j,1))*(z-z1))+x1; end ITxxE=0; ITxyE=0; ITxzE=0; ITyxE=0; ITyyE=0; ITyzE=0; ITzxE=0; ITzyE=0; ITzzE=0; IUxxE=0; IUxyE=0; IUxzE=0; IUyxE=0; IUyyE=0; IUyzE=0; IUzxE=0; IUzyE=0; IUzzE=0; for k=1:16 r=sqrt((x(k)-xbar(i))^2+(y(k)-ybar(i))^2+(z(k)-zbar(i))^2); drbydx=(x(k)-xbar(i))/r; drbydy=(y(k)-ybar(i))/r; drbydz=(z(k)-zbar(i))/r; costheta=(1/r)*((x(k)-xbar(i))*nx(fbcelements(j,1))+(y(k)-ybar(i))*ny(fbcelements(j,1))+(z(k)-zbar(i))*nz(fbcelements(j,1))); Txx=(-C2/r^2)*((C3+3*drbydx^2)*costheta); Tyy=(-C2/r^2)*((C3+3*drbydy^2)*costheta); Tzz=(-C2/r^2)*((C3+3*drbydz^2)*costheta); Txy=(-C2/r^2)*(3*drbydx*drbydy*costheta-C3*(ny(fbcelements(j,1))*drbydx-nx(fbcelements(j,1))*drbydy)); Tyx=(-C2/r^2)*(3*drbydy*drbydx*costheta-C3*(nx(fbcelements(j,1))*drbydy-ny(fbcelements(j,1))*drbydx)); Tyz=(-C2/r^2)*(3*drbydy*drbydz*costheta-C3*(nz(fbcelements(j,1))*drbydy-ny(fbcelements(j,1))*drbydz)); Tzy=(-C2/r^2)*(3*drbydz*drbydy*costheta-C3*(ny(fbcelements(j,1))*drbydz-nz(fbcelements(j,1))*drbydy)); Txz=(-C2/r^2)*(3*drbydx*drbydz*costheta-C3*(nz(fbcelements(j,1))*drbydx-nx(fbcelements(j,1))*drbydz)); Tzx=(-C2/r^2)*(3*drbydz*drbydx*costheta-C3*(nx(fbcelements(j,1))*drbydz-nz(fbcelements(j,1))*drbydx)); ITxxE=ITxxE+Txx*(1-v(k))*J(fbcelements(j,1)); ITxyE=ITxyE+Txy*(1-v(k))*J(fbcelements(j,1)); ITxzE=ITxzE+Txz*(1-v(k))*J(fbcelements(j,1)); ITyxE=ITyxE+Tyx*(1-v(k))*J(fbcelements(j,1)); ITyyE=ITyyE+Tyy*(1-v(k))*J(fbcelements(j,1)); ITyzE=ITyzE+Tyz*(1-v(k))*J(fbcelements(j,1)); ITzxE=ITzxE+Tzx*(1-v(k))*J(fbcelements(j,1)); ITzyE=ITzyE+Tzy*(1-v(k))*J(fbcelements(j,1)); ITzzE=ITzzE+Tzz*(1-v(k))*J(fbcelements(j,1)); Uxx=(C/r)*(C1+drbydx^2); Uyy=(C/r)*(C1+drbydy^2); Uzz=(C/r)*(C1+drbydz^2); Uxy=(C/r)*drbydx*drbydy; Uyx=Uxy; Uyz=(C/r)*drbydy*drbydz; Uzy=Uyz; Uxz=(C/r)*drbydx*drbydz; Uzx=Uxz; IUxxE=IUxxE+Uxx*(1-v(k))*J(fbcelements(j,1)); IUxyE=IUxyE+Uxy*(1-v(k))*J(fbcelements(j,1)); IUxzE=IUxzE+Uxz*(1-v(k))*J(fbcelements(j,1)); IUyxE=IUyxE+Uyx*(1-v(k))*J(fbcelements(j,1)); IUyyE=IUyyE+Uyy*(1-v(k))*J(fbcelements(j,1)); IUyzE=IUyzE+Uyz*(1-v(k))*J(fbcelements(j,1)); IUzxE=IUzxE+Uzx*(1-v(k))*J(fbcelements(j,1)); IUzyE=IUzyE+Uzy*(1-v(k))*J(fbcelements(j,1)); IUzzE=IUzzE+Uzz*(1-v(k))*J(fbcelements(j,1)); end ITxxE=(1/16)*ITxxE; ITxyE=(1/16)*ITxyE; ITxzE=(1/16)*ITxzE; ITyxE=(1/16)*ITyxE; ITyyE=(1/16)*ITyyE; ITyzE=(1/16)*ITyzE; ITzxE=(1/16)*ITzxE; ITzyE=(1/16)*ITzyE; ITzzE=(1/16)*ITzzE; IUxxE=(1/16)*IUxxE; IUxyE=(1/16)*IUxyE; IUxzE=(1/16)*IUxzE; IUyxE=(1/16)*IUyxE; IUyyE=(1/16)*IUyyE; IUyzE=(1/16)*IUyzE; IUzxE=(1/16)*IUzxE; IUzyE=(1/16)*IUzyE; IUzzE=(1/16)*IUzzE; K(3*i-2,3*fbcelements(j,1)-2)=ITxxE; K(3*i-2,3*fbcelements(j,1)-1)=ITxyE; K(3*i-2,3*fbcelements(j,1))=ITxzE; K(3*i-1,3*fbcelements(j,1)-2)=ITyxE; K(3*i-1,3*fbcelements(j,1)-1)=ITyyE; K(3*i-1,3*fbcelements(j,1))=ITyzE; K(3*i,3*fbcelements(j,1)-2)=ITzxE; K(3*i,3*fbcelements(j,1)-1)=ITzyE; K(3*i,3*fbcelements(j,1))=ITzzE; F(3*i-2)=F(3*i-2)+fbcelements(j,2)*IUxxE+fbcelements(j,3)*IUxyE+fbcelements(j,4)*IUxzE; F(3*i-1)=F(3*i-1)+fbcelements(j,2)*IUyxE+fbcelements(j,3)*IUyyE+fbcelements(j,4)*IUyzE; F(3*i)=F(3*i)+fbcelements(j,2)*IUzxE+fbcelements(j,3)*IUzyE+fbcelements(j,4)*IUzzE; if i==fbcelements(j,1) K(3*fbcelements(j,1)-2,3*fbcelements(j,1)-2)=K(3*fbcelements(j,1)-2,3*fbcelements(j,1)-2)+0.5; K(3*fbcelements(j,1)-1,3*fbcelements(j,1)-1)=K(3*fbcelements(j,1)-1,3*fbcelements(j,1)-1)+0.5; K(3*fbcelements(j,1),3*fbcelements(j,1))=K(3*fbcelements(j,1),3*fbcelements(j,1))+0.5; end end for j=1:ndispbcelements if (abs(nz(dispbcelements(j,1)))>=1/sqrt(3)) x1=xyzofelements(3*dispbcelements(j,1)-2,2); y1=xyzofelements(3*dispbcelements(j,1)-2,3); z1=xyzofelements(3*dispbcelements(j,1)-2,4); x2=xyzofelements(3*dispbcelements(j,1)-1,2); y2=xyzofelements(3*dispbcelements(j,1)-1,3); z2=xyzofelements(3*dispbcelements(j,1)-1,4); x3=xyzofelements(3*dispbcelements(j,1),2); y3=xyzofelements(3*dispbcelements(j,1),3); z3=xyzofelements(3*dispbcelements(j,1),4); x=(x2-x1)*u+(x3-x1)*v+x1; y=(y2-y1)*u+(y3-y1)*v+y1; z=-(nz(dispbcelements(j,1)))^(-1)*(nx(dispbcelements(j,1))*(x-x1)+ny(dispbcelements(j,1))*(y-y1))+z1; elseif (abs(nz(dispbcelements(j,1)))<1/sqrt(3) && abs(ny(dispbcelements(j,1)))>=1/sqrt(3)) x1=xyzofelements(3*dispbcelements(j,1)-2,2); y1=xyzofelements(3*dispbcelements(j,1)-2,3); z1=xyzofelements(3*dispbcelements(j,1)-2,4); x2=xyzofelements(3*dispbcelements(j,1)-1,2); y2=xyzofelements(3*dispbcelements(j,1)-1,3); z2=xyzofelements(3*dispbcelements(j,1)-1,4); x3=xyzofelements(3*dispbcelements(j,1),2); y3=xyzofelements(3*dispbcelements(j,1),3); z3=xyzofelements(3*dispbcelements(j,1),4); x=(x2-x1)*u+(x3-x1)*v+x1; z=(z2-z1)*u+(z3-z1)*v+z1; y=-(ny(dispbcelements(j,1)))^(-1)*(nx(dispbcelements(j,1))*(x-x1)+nz(dispbcelements(j,1))*(z-z1))+y1; elseif (abs(nz(dispbcelements(j,1)))<1/sqrt(3) && abs(ny(dispbcelements(j,1)))<1/sqrt(3)) x1=xyzofelements(3*dispbcelements(j,1)-2,2); y1=xyzofelements(3*dispbcelements(j,1)-2,3); z1=xyzofelements(3*dispbcelements(j,1)-2,4); x2=xyzofelements(3*dispbcelements(j,1)-1,2); y2=xyzofelements(3*dispbcelements(j,1)-1,3); z2=xyzofelements(3*dispbcelements(j,1)-1,4); x3=xyzofelements(3*dispbcelements(j,1),2); y3=xyzofelements(3*dispbcelements(j,1),3); z3=xyzofelements(3*dispbcelements(j,1),4); y=(y2-y1)*u+(y3-y1)*v+y1; z=(z2-z1)*u+(z3-z1)*v+z1; x=-(nx(dispbcelements(j,1)))^(-1)*(ny(dispbcelements(j,1))*(y-y1)+nz(dispbcelements(j,1))*(z-z1))+x1; end ITxxE=0; ITxyE=0; ITxzE=0; ITyxE=0; ITyyE=0; ITyzE=0; ITzxE=0; ITzyE=0; ITzzE=0; IUxxE=0; IUxyE=0; IUxzE=0; IUyxE=0; IUyyE=0; IUyzE=0; IUzxE=0; IUzyE=0; IUzzE=0; for k=1:16 r=sqrt((x(k)-xbar(i))^2+(y(k)-ybar(i))^2+(z(k)-zbar(i))^2); drbydx=(x(k)-xbar(i))/r; drbydy=(y(k)-ybar(i))/r; drbydz=(z(k)-zbar(i))/r; costheta=(1/r)*((x(k)-xbar(i))*nx(dispbcelements(j,1))+(y(k)-ybar(i))*ny(dispbcelements(j,1))+(z(k)-zbar(i))*nz(dispbcelements(j,1))); Txx=(-C2/r^2)*((C3+3*drbydx^2)*costheta); Tyy=(-C2/r^2)*((C3+3*drbydy^2)*costheta); Tzz=(-C2/r^2)*((C3+3*drbydz^2)*costheta); Txy=(-C2/r^2)*(3*drbydx*drbydy*costheta-C3*(ny(dispbcelements(j,1))*drbydx-nx(dispbcelements(j,1))*drbydy)); Tyx=(-C2/r^2)*(3*drbydy*drbydx*costheta-C3*(nx(dispbcelements(j,1))*drbydy-ny(dispbcelements(j,1))*drbydx)); Tyz=(-C2/r^2)*(3*drbydy*drbydz*costheta-C3*(nz(dispbcelements(j,1))*drbydy-ny(dispbcelements(j,1))*drbydz)); Tzy=(-C2/r^2)*(3*drbydz*drbydy*costheta-C3*(ny(dispbcelements(j,1))*drbydz-nz(dispbcelements(j,1))*drbydy)); Txz=(-C2/r^2)*(3*drbydx*drbydz*costheta-C3*(nz(dispbcelements(j,1))*drbydx-nx(dispbcelements(j,1))*drbydz)); Tzx=(-C2/r^2)*(3*drbydz*drbydx*costheta-C3*(nx(dispbcelements(j,1))*drbydz-nz(dispbcelements(j,1))*drbydx)); ITxxE=ITxxE+Txx*(1-v(k))*J(dispbcelements(j,1)); ITxyE=ITxyE+Txy*(1-v(k))*J(dispbcelements(j,1)); ITxzE=ITxzE+Txz*(1-v(k))*J(dispbcelements(j,1)); ITyxE=ITyxE+Tyx*(1-v(k))*J(dispbcelements(j,1)); ITyyE=ITyyE+Tyy*(1-v(k))*J(dispbcelements(j,1)); ITyzE=ITyzE+Tyz*(1-v(k))*J(dispbcelements(j,1)); ITzxE=ITzxE+Tzx*(1-v(k))*J(dispbcelements(j,1)); ITzyE=ITzyE+Tzy*(1-v(k))*J(dispbcelements(j,1)); ITzzE=ITzzE+Tzz*(1-v(k))*J(dispbcelements(j,1)); Uxx=(C/r)*(C1+drbydx^2); Uyy=(C/r)*(C1+drbydy^2); Uzz=(C/r)*(C1+drbydz^2); Uxy=(C/r)*drbydx*drbydy; Uyx=Uxy; Uyz=(C/r)*drbydy*drbydz; Uzy=Uyz; Uxz=(C/r)*drbydx*drbydz; Uzx=Uxz; IUxxE=IUxxE+Uxx*(1-v(k))*J(dispbcelements(j,1)); IUxyE=IUxyE+Uxy*(1-v(k))*J(dispbcelements(j,1)); IUxzE=IUxzE+Uxz*(1-v(k))*J(dispbcelements(j,1)); IUyxE=IUyxE+Uyx*(1-v(k))*J(dispbcelements(j,1)); IUyyE=IUyyE+Uyy*(1-v(k))*J(dispbcelements(j,1)); IUyzE=IUyzE+Uyz*(1-v(k))*J(dispbcelements(j,1)); IUzxE=IUzxE+Uzx*(1-v(k))*J(dispbcelements(j,1)); IUzyE=IUzyE+Uzy*(1-v(k))*J(dispbcelements(j,1)); IUzzE=IUzzE+Uzz*(1-v(k))*J(dispbcelements(j,1)); end ITxxE=(1/16)*ITxxE; ITxyE=(1/16)*ITxyE; ITxzE=(1/16)*ITxzE; ITyxE=(1/16)*ITyxE; ITyyE=(1/16)*ITyyE; ITyzE=(1/16)*ITyzE; ITzxE=(1/16)*ITzxE; ITzyE=(1/16)*ITzyE; ITzzE=(1/16)*ITzzE; IUxxE=(1/16)*IUxxE; IUxyE=(1/16)*IUxyE; IUxzE=(1/16)*IUxzE; IUyxE=(1/16)*IUyxE; IUyyE=(1/16)*IUyyE; IUyzE=(1/16)*IUyzE; IUzxE=(1/16)*IUzxE; IUzyE=(1/16)*IUzyE; IUzzE=(1/16)*IUzzE; K(3*i-2,3*dispbcelements(j,1)-2)=IUxxE; K(3*i-2,3*dispbcelements(j,1)-1)=IUxyE; K(3*i-2,3*dispbcelements(j,1))=IUxzE; K(3*i-1,3*dispbcelements(j,1)-2)=IUyxE; K(3*i-1,3*dispbcelements(j,1)-1)=IUyyE; K(3*i-1,3*dispbcelements(j,1))=IUyzE; K(3*i,3*dispbcelements(j,1)-2)=IUzxE; K(3*i,3*dispbcelements(j,1)-1)=IUzyE; K(3*i,3*dispbcelements(j,1))=IUzzE; F(3*i-2)=F(3*i-2)+dispbcelements(j,2)*ITxxE+dispbcelements(j,3)*ITxyE+dispbcelements(j,4)*ITxzE; F(3*i-1)=F(3*i-1)+dispbcelements(j,2)*ITyxE+dispbcelements(j,3)*ITyyE+dispbcelements(j,4)*ITyzE; F(3*i)=F(3*i)+dispbcelements(j,2)*ITzxE+dispbcelements(j,3)*ITzyE+dispbcelements(j,4)*ITzzE; if i==dispbcelements(j,1) F(3*i-2)=F(3*i-2)+0.5*dispbcelements(j,2); F(3*i-1)=F(3*i-1)+0.5*dispbcelements(j,3); F(3*i)=F(3*i)+0.5*dispbcelements(j,4); end end end %Solve U=K\F; % Display unknown dofs for each element U