% Copyright (C) 2014, Kirana Kumara P % Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: % The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. % THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. % ------------------------------------------------------------------------------------------------------------------------- % This code is authored by Kirana Kumara P, Centre for Product Design and Manufacturing, Indian Institute of Science, Bangalore - 560012, India. % Author hereby declares that the present code is the original work of the author. Further, author hereby declares that the present code, in full or in part, is not a translation or a copy of any of the existing codes written by someone else. % Author's institution (Indian Institute of Science) has informed the author in writing that the institution is not interested in claiming any copyright on the present code. % Contact email ID of the author is: kiranakumarap@gmail.com % It may be helpful to go through the following two sources to understand the present code: % [1] Kirana Kumara P, A MATLAB Code for Three Dimensional Linear Elastostatics using Constant Boundary Elements, International Journal of Advances in Engineering Sciences, Vol. 2, No. 3, pp. 9-20, 2012. % [2] Kirana Kumara P, A study of the speed and the accuracy of the Boundary Element Method as applied to the computational simulation of biological organs, Available from: http://arxiv.org/abs/1311.4533v1 % ------------------------------------------------------------------------------------------------------------------------- % This is a MATLAB program which solves any three dimensional linear elastostatic (without body forces) problem using constant boundary elements. % This is a serial code % Strongly singular integrals are handled using rigid body modes % Weakly singular integrals are not handled by any special methods. The program uses 16 by 16 (256 in total) integration points over each of the boundary elements; high number of numerical integration points take care of weak singularity. Program can be easily modified to use 4 by 4, 8 by 8, or 32 by 32 integration points (locations and weights of the integration points for these cases are listed in this code too). Test runs showed that using 4 by 4 integration does not give very accure results whereas 8 by 8 integration is fine. But use 16 by 16 integration to ensure that one does not get erroneous results because of inaccurate integration. There is not much improvement in accuracy by using more than 16 by 16 integration points. % This code has been tested using MATLAB R2011b % This program uses a geometry that is described by 96 triangular boundary elements. However, the program can be very easily modified to solve any problem that uses a geometry that is described by any number of surface triangles. % In this program, element numbers from 65 to 80 are fixed. Element numbers from 1 to 64 are subjected to zero traction, and element numbers from 81 to 96 are subjected to a traction '4' in the y direction. Hence, in total, 16 elements are subjected to displacement boundary condition, and 80 elements are subjected to traction boundary condition. When using this program to solve another problem, suitable modifications have to be made to incorporate appropriate boundary conditions. clear; clear all; clc; format long; E=200000; mu=0.33; dispbcelements=[65 0 0 0;66 0 0 0;67 0 0 0;68 0 0 0; 69 0 0 0; 70 0 0 0; 71 0 0 0; 72 0 0 0; 73 0 0 0; 74 0 0 0; 75 0 0 0;76 0 0 0;77 0 0 0;78 0 0 0; 79 0 0 0; 80 0 0 0]; fbcelements=[81 0 4 0;82 0 4 0;83 0 4 0;84 0 4 0; 85 0 4 0;86 0 4 0; 87 0 4 0; 88 0 4 0; 89 0 4 0; 90 0 4 0; 91 0 4 0;92 0 4 0; 93 0 4 0;94 0 4 0; 95 0 4 0; 96 0 4 0]; fbczeroelements(1:64,1:4)=0; for i=1:64 fbczeroelements(i,1)=i; end % fbczeroelements(_:_,1:4)=0; % for i=_:_ % fbczeroelements(i,1)=i; % end % fbczeroelements(_:_,1:4)=0; % for i=_:_ % fbczeroelements(i,1)=i; % end % fbczeroelements(_:_,1:4)=0; % for i=_:_ % fbczeroelements(i,1)=i; % end % fbczeroelements(_:_,:)=[]; % fbczeroelements(_:_,:)=[]; fbcelements=[fbczeroelements;fbcelements]; xyzofelements=[1 1.000000 4.000000 1.000000; 1 2.000000 4.000000 2.000000; 1 2.000000 4.000000 0.000000; 2 1.000000 4.000000 1.000000; 2 2.000000 4.000000 0.000000; 2 0.000000 4.000000 0.000000; 3 0.000000 1.000000 3.000000; 3 0.000000 0.000000 4.000000; 3 0.000000 2.000000 4.000000; 4 0.000000 1.000000 3.000000; 4 0.000000 2.000000 4.000000; 4 0.000000 2.000000 2.000000; 5 3.000000 4.000000 1.000000; 5 2.000000 4.000000 2.000000; 5 4.000000 4.000000 2.000000; 6 3.000000 4.000000 1.000000; 6 4.000000 4.000000 2.000000; 6 4.000000 4.000000 0.000000; 7 4.000000 3.000000 1.000000; 7 4.000000 4.000000 0.000000; 7 4.000000 4.000000 2.000000; 8 4.000000 3.000000 1.000000; 8 4.000000 4.000000 2.000000; 8 4.000000 2.000000 2.000000; 9 0.000000 3.000000 3.000000; 9 0.000000 2.000000 2.000000; 9 0.000000 2.000000 4.000000; 10 0.000000 3.000000 3.000000; 10 0.000000 2.000000 4.000000; 10 0.000000 4.000000 4.000000; 11 4.000000 1.000000 1.000000; 11 4.000000 0.000000 0.000000; 11 4.000000 2.000000 0.000000; 12 4.000000 1.000000 1.000000; 12 4.000000 2.000000 0.000000; 12 4.000000 2.000000 2.000000; 13 4.000000 1.000000 3.000000; 13 4.000000 2.000000 2.000000; 13 4.000000 2.000000 4.000000; 14 4.000000 1.000000 3.000000; 14 4.000000 2.000000 4.000000; 14 4.000000 0.000000 4.000000; 15 4.000000 3.000000 3.000000; 15 4.000000 2.000000 2.000000; 15 4.000000 4.000000 2.000000; 16 4.000000 3.000000 3.000000; 16 4.000000 4.000000 2.000000; 16 4.000000 4.000000 4.000000; 17 0.000000 1.000000 1.000000; 17 0.000000 2.000000 2.000000; 17 0.000000 2.000000 0.000000; 18 0.000000 1.000000 1.000000; 18 0.000000 2.000000 0.000000; 18 0.000000 0.000000 0.000000; 19 3.000000 0.000000 3.000000; 19 4.000000 0.000000 4.000000; 19 2.000000 0.000000 4.000000; 20 3.000000 0.000000 3.000000; 20 2.000000 0.000000 4.000000; 20 2.000000 0.000000 2.000000; 21 3.000000 0.000000 1.000000; 21 4.000000 0.000000 0.000000; 21 4.000000 0.000000 2.000000; 22 3.000000 0.000000 1.000000; 22 4.000000 0.000000 2.000000; 22 2.000000 0.000000 2.000000; 23 1.000000 0.000000 1.000000; 23 0.000000 0.000000 0.000000; 23 2.000000 0.000000 0.000000; 24 1.000000 0.000000 1.000000; 24 2.000000 0.000000 0.000000; 24 2.000000 0.000000 2.000000; 25 1.000000 0.000000 3.000000; 25 2.000000 0.000000 2.000000; 25 2.000000 0.000000 4.000000; 26 1.000000 0.000000 3.000000; 26 2.000000 0.000000 4.000000; 26 0.000000 0.000000 4.000000; 27 0.000000 3.000000 1.000000; 27 0.000000 2.000000 2.000000; 27 0.000000 4.000000 2.000000; 28 0.000000 3.000000 1.000000; 28 0.000000 4.000000 2.000000; 28 0.000000 4.000000 0.000000; 29 3.000000 4.000000 3.000000; 29 4.000000 4.000000 4.000000; 29 4.000000 4.000000 2.000000; 30 3.000000 4.000000 3.000000; 30 4.000000 4.000000 2.000000; 30 2.000000 4.000000 2.000000; 31 1.000000 4.000000 3.000000; 31 0.000000 4.000000 4.000000; 31 2.000000 4.000000 4.000000; 32 1.000000 4.000000 3.000000; 32 2.000000 4.000000 4.000000; 32 2.000000 4.000000 2.000000; 33 2.000000 4.000000 2.000000; 33 2.000000 4.000000 4.000000; 33 3.000000 4.000000 3.000000; 34 2.000000 4.000000 2.000000; 34 3.000000 4.000000 1.000000; 34 2.000000 4.000000 0.000000; 35 4.000000 2.000000 2.000000; 35 4.000000 2.000000 0.000000; 35 4.000000 3.000000 1.000000; 36 4.000000 2.000000 2.000000; 36 4.000000 3.000000 3.000000; 36 4.000000 2.000000 4.000000; 37 0.000000 2.000000 2.000000; 37 0.000000 1.000000 1.000000; 37 0.000000 0.000000 2.000000; 38 0.000000 2.000000 2.000000; 38 0.000000 0.000000 2.000000; 38 0.000000 1.000000 3.000000; 39 2.000000 0.000000 2.000000; 39 2.000000 0.000000 0.000000; 39 3.000000 0.000000 1.000000; 40 2.000000 0.000000 2.000000; 40 4.000000 0.000000 2.000000; 40 3.000000 0.000000 3.000000; 41 0.000000 3.000000 3.000000; 41 0.000000 4.000000 2.000000; 41 0.000000 2.000000 2.000000; 42 0.000000 2.000000 0.000000; 42 0.000000 2.000000 2.000000; 42 0.000000 3.000000 1.000000; 43 0.000000 4.000000 4.000000; 43 0.000000 4.000000 2.000000; 43 0.000000 3.000000 3.000000; 44 0.000000 1.000000 3.000000; 44 0.000000 0.000000 2.000000; 44 0.000000 0.000000 4.000000; 45 0.000000 0.000000 0.000000; 45 0.000000 0.000000 2.000000; 45 0.000000 1.000000 1.000000; 46 0.000000 2.000000 0.000000; 46 0.000000 3.000000 1.000000; 46 0.000000 4.000000 0.000000; 47 1.000000 4.000000 3.000000; 47 2.000000 4.000000 2.000000; 47 0.000000 4.000000 2.000000; 48 1.000000 4.000000 1.000000; 48 0.000000 4.000000 2.000000; 48 2.000000 4.000000 2.000000; 49 4.000000 4.000000 4.000000; 49 3.000000 4.000000 3.000000; 49 2.000000 4.000000 4.000000; 50 1.000000 4.000000 3.000000; 50 0.000000 4.000000 2.000000; 50 0.000000 4.000000 4.000000; 51 0.000000 4.000000 0.000000; 51 0.000000 4.000000 2.000000; 51 1.000000 4.000000 1.000000; 52 2.000000 4.000000 0.000000; 52 3.000000 4.000000 1.000000; 52 4.000000 4.000000 0.000000; 53 4.000000 1.000000 3.000000; 53 4.000000 0.000000 2.000000; 53 4.000000 2.000000 2.000000; 54 4.000000 1.000000 1.000000; 54 4.000000 2.000000 2.000000; 54 4.000000 0.000000 2.000000; 55 4.000000 0.000000 4.000000; 55 4.000000 0.000000 2.000000; 55 4.000000 1.000000 3.000000; 56 4.000000 2.000000 4.000000; 56 4.000000 3.000000 3.000000; 56 4.000000 4.000000 4.000000; 57 4.000000 4.000000 0.000000; 57 4.000000 3.000000 1.000000; 57 4.000000 2.000000 0.000000; 58 4.000000 1.000000 1.000000; 58 4.000000 0.000000 2.000000; 58 4.000000 0.000000 0.000000; 59 1.000000 0.000000 3.000000; 59 0.000000 0.000000 2.000000; 59 2.000000 0.000000 2.000000; 60 1.000000 0.000000 1.000000; 60 2.000000 0.000000 2.000000; 60 0.000000 0.000000 2.000000; 61 0.000000 0.000000 4.000000; 61 0.000000 0.000000 2.000000; 61 1.000000 0.000000 3.000000; 62 3.000000 0.000000 3.000000; 62 4.000000 0.000000 2.000000; 62 4.000000 0.000000 4.000000; 63 4.000000 0.000000 0.000000; 63 3.000000 0.000000 1.000000; 63 2.000000 0.000000 0.000000; 64 1.000000 0.000000 1.000000; 64 0.000000 0.000000 2.000000; 64 0.000000 0.000000 0.000000; 65 3.000000 1.000000 0.000000; 65 2.000000 2.000000 0.000000; 65 4.000000 2.000000 0.000000; 66 3.000000 1.000000 0.000000; 66 4.000000 2.000000 0.000000; 66 4.000000 0.000000 0.000000; 67 3.000000 3.000000 0.000000; 67 4.000000 4.000000 0.000000; 67 4.000000 2.000000 0.000000; 68 3.000000 3.000000 0.000000; 68 4.000000 2.000000 0.000000; 68 2.000000 2.000000 0.000000; 69 1.000000 3.000000 0.000000; 69 0.000000 4.000000 0.000000; 69 2.000000 4.000000 0.000000; 70 1.000000 3.000000 0.000000; 70 2.000000 4.000000 0.000000; 70 2.000000 2.000000 0.000000; 71 1.000000 1.000000 0.000000; 71 2.000000 2.000000 0.000000; 71 2.000000 0.000000 0.000000; 72 1.000000 1.000000 0.000000; 72 2.000000 0.000000 0.000000; 72 0.000000 0.000000 0.000000; 73 2.000000 2.000000 0.000000; 73 2.000000 4.000000 0.000000; 73 3.000000 3.000000 0.000000; 74 2.000000 2.000000 0.000000; 74 3.000000 1.000000 0.000000; 74 2.000000 0.000000 0.000000; 75 2.000000 2.000000 0.000000; 75 0.000000 2.000000 0.000000; 75 1.000000 3.000000 0.000000; 76 1.000000 1.000000 0.000000; 76 0.000000 2.000000 0.000000; 76 2.000000 2.000000 0.000000; 77 1.000000 3.000000 0.000000; 77 0.000000 2.000000 0.000000; 77 0.000000 4.000000 0.000000; 78 3.000000 3.000000 0.000000; 78 2.000000 4.000000 0.000000; 78 4.000000 4.000000 0.000000; 79 0.000000 0.000000 0.000000; 79 0.000000 2.000000 0.000000; 79 1.000000 1.000000 0.000000; 80 2.000000 0.000000 0.000000; 80 3.000000 1.000000 0.000000; 80 4.000000 0.000000 0.000000; 81 1.000000 3.000000 4.000000; 81 2.000000 2.000000 4.000000; 81 2.000000 4.000000 4.000000; 82 1.000000 3.000000 4.000000; 82 2.000000 4.000000 4.000000; 82 0.000000 4.000000 4.000000; 83 3.000000 3.000000 4.000000; 83 4.000000 4.000000 4.000000; 83 2.000000 4.000000 4.000000; 84 3.000000 3.000000 4.000000; 84 2.000000 4.000000 4.000000; 84 2.000000 2.000000 4.000000; 85 1.000000 1.000000 4.000000; 85 0.000000 0.000000 4.000000; 85 2.000000 0.000000 4.000000; 86 1.000000 1.000000 4.000000; 86 2.000000 0.000000 4.000000; 86 2.000000 2.000000 4.000000; 87 3.000000 1.000000 4.000000; 87 4.000000 0.000000 4.000000; 87 4.000000 2.000000 4.000000; 88 3.000000 1.000000 4.000000; 88 4.000000 2.000000 4.000000; 88 2.000000 2.000000 4.000000; 89 2.000000 2.000000 4.000000; 89 2.000000 0.000000 4.000000; 89 3.000000 1.000000 4.000000; 90 2.000000 2.000000 4.000000; 90 4.000000 2.000000 4.000000; 90 3.000000 3.000000 4.000000; 91 1.000000 3.000000 4.000000; 91 0.000000 2.000000 4.000000; 91 2.000000 2.000000 4.000000; 92 1.000000 1.000000 4.000000; 92 2.000000 2.000000 4.000000; 92 0.000000 2.000000 4.000000; 93 0.000000 4.000000 4.000000; 93 0.000000 2.000000 4.000000; 93 1.000000 3.000000 4.000000; 94 4.000000 4.000000 4.000000; 94 3.000000 3.000000 4.000000; 94 4.000000 2.000000 4.000000; 95 4.000000 0.000000 4.000000; 95 3.000000 1.000000 4.000000; 95 2.000000 0.000000 4.000000; 96 1.000000 1.000000 4.000000; 96 0.000000 2.000000 4.000000; 96 0.000000 0.000000 4.000000]; 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)); xbar(i)=(x2-x1)*0.25+(x3-x1)*0.5+x1; ybar(i)=(y2-y1)*0.25+(y3-y1)*0.5+y1; zbar(i)=(z2-z1)*0.25+(z3-z1)*0.5+z1; end % For 4 by 4 integration % tonedim=[-0.861136311594 -0.339981043585 0.339981043585 0.861136311594]; % vonedim=[-0.861136311594 -0.339981043585 0.339981043585 0.861136311594]; % weightonedim=[0.347854845137 0.652145154863 0.652145154863 0.347854845137]; % For 8 by 8 integration % tonedim=[-0.960289856498 -0.796666477414 -0.525532409916 -0.183434642496 0.183434642496 0.525532409916 0.796666477414 0.960289856498]; % vonedim=[-0.960289856498 -0.796666477414 -0.525532409916 -0.183434642496 0.183434642496 0.525532409916 0.796666477414 0.960289856498]; % weightonedim=[0.10122853629 0.222381034453 0.313706645878 0.362683783378 0.362683783378 0.313706645878 0.222381034453 0.10122853629]; % 16 by 16 integration tonedim=[-0.989400934992 -0.944575023073 -0.865631202388 -0.755404408355 -0.617876244403 -0.458016777657 -0.281603550779 -0.0950125098376 0.0950125098376 0.281603550779 0.458016777657 0.617876244403 0.755404408355 0.865631202388 0.944575023073 0.989400934992]; vonedim=[-0.989400934992 -0.944575023073 -0.865631202388 -0.755404408355 -0.617876244403 -0.458016777657 -0.281603550779 -0.0950125098376 0.0950125098376 0.281603550779 0.458016777657 0.617876244403 0.755404408355 0.865631202388 0.944575023073 0.989400934992]; weightonedim=[0.027152459411 0.0622535239372 0.0951585116838 0.124628971256 0.149595988817 0.169156519395 0.182603415045 0.189450610455 0.189450610455 0.182603415045 0.169156519395 0.149595988817 0.124628971256 0.0951585116838 0.0622535239372 0.027152459411]; % For 32 by 32 integration % tonedim=[-0.997263861849 -0.985611511545 -0.964762255588 -0.934906075938 -0.896321155766 -0.849367613733 -0.794483795968 -0.73218211874 -0.66304426693 -0.587715757241 -0.506899908932 -0.421351276131 -0.331868602282 -0.239287362252 -0.144471961583 -0.0483076656877 0.0483076656877 0.144471961583 0.239287362252 0.331868602282 0.421351276131 0.506899908932 0.587715757241 0.66304426693 0.73218211874 0.794483795968 0.849367613733 0.896321155766 0.934906075938 0.964762255588 0.985611511545 0.997263861849]; % vonedim=[-0.997263861849 -0.985611511545 -0.964762255588 -0.934906075938 -0.896321155766 -0.849367613733 -0.794483795968 -0.73218211874 -0.66304426693 -0.587715757241 -0.506899908932 -0.421351276131 -0.331868602282 -0.239287362252 -0.144471961583 -0.0483076656877 0.0483076656877 0.144471961583 0.239287362252 0.331868602282 0.421351276131 0.506899908932 0.587715757241 0.66304426693 0.73218211874 0.794483795968 0.849367613733 0.896321155766 0.934906075938 0.964762255588 0.985611511545 0.997263861849]; % weightonedim=[0.00701814576495 0.0162774265831 0.0253910098329 0.0342745478477 0.0428359896785 0.0509978738117 0.0586839394615 0.0658220603578 0.0723456094297 0.078193695762 0.083311711103 0.0876518688047 0.0911736454878 0.0938441590423 0.0956384754512 0.0965398415811 0.0965398415811 0.0956384754512 0.0938441590423 0.0911736454878 0.0876518688047 0.083311711103 0.078193695762 0.0723456094297 0.0658220603578 0.0586839394615 0.0509978738117 0.0428359896785 0.0342745478477 0.0253910098329 0.0162774265831 0.00701814576495]; index=1; t(1:256)=0; v(1:256)=0; weight(1:256)=0; for tindex=1:16 for vindex=1:16 t(index)=tonedim(tindex); v(index)=vonedim(vindex); weight(index)=weightonedim(tindex)*weightonedim(vindex); t(index)=0.5*t(index)+0.5; v(index)=0.5*v(index)+0.5; index=index+1; end end u=t.*(1-v); K=zeros(3*nelement); F(1:3*nelement,1)=0; CPVITxxE(1:nelement)=0; CPVITxyE(1:nelement)=0; CPVITxzE(1:nelement)=0; CPVITyxE(1:nelement)=0; CPVITyyE(1:nelement)=0; CPVITyzE(1:nelement)=0; CPVITzxE(1:nelement)=0; CPVITzyE(1:nelement)=0; CPVITzzE(1:nelement)=0; for i=1:nelement for j=1:nfbcelements 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=(z2-z1)*u+(z3-z1)*v+z1; 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:256 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))*weight(k); ITxyE=ITxyE+Txy*(1-v(k))*J(fbcelements(j,1))*weight(k); ITxzE=ITxzE+Txz*(1-v(k))*J(fbcelements(j,1))*weight(k); ITyxE=ITyxE+Tyx*(1-v(k))*J(fbcelements(j,1))*weight(k); ITyyE=ITyyE+Tyy*(1-v(k))*J(fbcelements(j,1))*weight(k); ITyzE=ITyzE+Tyz*(1-v(k))*J(fbcelements(j,1))*weight(k); ITzxE=ITzxE+Tzx*(1-v(k))*J(fbcelements(j,1))*weight(k); ITzyE=ITzyE+Tzy*(1-v(k))*J(fbcelements(j,1))*weight(k); ITzzE=ITzzE+Tzz*(1-v(k))*J(fbcelements(j,1))*weight(k); 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))*weight(k); IUxyE=IUxyE+Uxy*(1-v(k))*J(fbcelements(j,1))*weight(k); IUxzE=IUxzE+Uxz*(1-v(k))*J(fbcelements(j,1))*weight(k); IUyxE=IUyxE+Uyx*(1-v(k))*J(fbcelements(j,1))*weight(k); IUyyE=IUyyE+Uyy*(1-v(k))*J(fbcelements(j,1))*weight(k); IUyzE=IUyzE+Uyz*(1-v(k))*J(fbcelements(j,1))*weight(k); IUzxE=IUzxE+Uzx*(1-v(k))*J(fbcelements(j,1))*weight(k); IUzyE=IUzyE+Uzy*(1-v(k))*J(fbcelements(j,1))*weight(k); IUzzE=IUzzE+Uzz*(1-v(k))*J(fbcelements(j,1))*weight(k); end ITxxE=0.5*ITxxE; ITxyE=0.5*ITxyE; ITxzE=0.5*ITxzE; ITyxE=0.5*ITyxE; ITyyE=0.5*ITyyE; ITyzE=0.5*ITyzE; ITzxE=0.5*ITzxE; ITzyE=0.5*ITzyE; ITzzE=0.5*ITzzE; IUxxE=0.5*IUxxE; IUxyE=0.5*IUxyE; IUxzE=0.5*IUxzE; IUyxE=0.5*IUyxE; IUyyE=0.5*IUyyE; IUyzE=0.5*IUyzE; IUzxE=0.5*IUzxE; IUzyE=0.5*IUzyE; IUzzE=0.5*IUzzE; CPVITxxE(i)=CPVITxxE(i)-ITxxE; CPVITxyE(i)=CPVITxyE(i)-ITxyE; CPVITxzE(i)=CPVITxzE(i)-ITxzE; CPVITyxE(i)=CPVITyxE(i)-ITyxE; CPVITyyE(i)=CPVITyyE(i)-ITyyE; CPVITyzE(i)=CPVITyzE(i)-ITyzE; CPVITzxE(i)=CPVITzxE(i)-ITzxE; CPVITzyE(i)=CPVITzyE(i)-ITzyE; CPVITzzE(i)=CPVITzzE(i)-ITzzE; if i==fbcelements(j,1) CPVITxxE(i)=CPVITxxE(i)+ITxxE; CPVITxyE(i)=CPVITxyE(i)+ITxyE; CPVITxzE(i)=CPVITxzE(i)+ITxzE; CPVITyxE(i)=CPVITyxE(i)+ITyxE; CPVITyyE(i)=CPVITyyE(i)+ITyyE; CPVITyzE(i)=CPVITyzE(i)+ITyzE; CPVITzxE(i)=CPVITzxE(i)+ITzxE; CPVITzyE(i)=CPVITzyE(i)+ITzyE; CPVITzzE(i)=CPVITzzE(i)+ITzzE; end 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; end for j=1:ndispbcelements 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=(z2-z1)*u+(z3-z1)*v+z1; 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:256 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))*weight(k); ITxyE=ITxyE+Txy*(1-v(k))*J(dispbcelements(j,1))*weight(k); ITxzE=ITxzE+Txz*(1-v(k))*J(dispbcelements(j,1))*weight(k); ITyxE=ITyxE+Tyx*(1-v(k))*J(dispbcelements(j,1))*weight(k); ITyyE=ITyyE+Tyy*(1-v(k))*J(dispbcelements(j,1))*weight(k); ITyzE=ITyzE+Tyz*(1-v(k))*J(dispbcelements(j,1))*weight(k); ITzxE=ITzxE+Tzx*(1-v(k))*J(dispbcelements(j,1))*weight(k); ITzyE=ITzyE+Tzy*(1-v(k))*J(dispbcelements(j,1))*weight(k); ITzzE=ITzzE+Tzz*(1-v(k))*J(dispbcelements(j,1))*weight(k); 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))*weight(k); IUxyE=IUxyE+Uxy*(1-v(k))*J(dispbcelements(j,1))*weight(k); IUxzE=IUxzE+Uxz*(1-v(k))*J(dispbcelements(j,1))*weight(k); IUyxE=IUyxE+Uyx*(1-v(k))*J(dispbcelements(j,1))*weight(k); IUyyE=IUyyE+Uyy*(1-v(k))*J(dispbcelements(j,1))*weight(k); IUyzE=IUyzE+Uyz*(1-v(k))*J(dispbcelements(j,1))*weight(k); IUzxE=IUzxE+Uzx*(1-v(k))*J(dispbcelements(j,1))*weight(k); IUzyE=IUzyE+Uzy*(1-v(k))*J(dispbcelements(j,1))*weight(k); IUzzE=IUzzE+Uzz*(1-v(k))*J(dispbcelements(j,1))*weight(k); end ITxxE=0.5*ITxxE; ITxyE=0.5*ITxyE; ITxzE=0.5*ITxzE; ITyxE=0.5*ITyxE; ITyyE=0.5*ITyyE; ITyzE=0.5*ITyzE; ITzxE=0.5*ITzxE; ITzyE=0.5*ITzyE; ITzzE=0.5*ITzzE; IUxxE=0.5*IUxxE; IUxyE=0.5*IUxyE; IUxzE=0.5*IUxzE; IUyxE=0.5*IUyxE; IUyyE=0.5*IUyyE; IUyzE=0.5*IUyzE; IUzxE=0.5*IUzxE; IUzyE=0.5*IUzyE; IUzzE=0.5*IUzzE; CPVITxxE(i)=CPVITxxE(i)-ITxxE; CPVITxyE(i)=CPVITxyE(i)-ITxyE; CPVITxzE(i)=CPVITxzE(i)-ITxzE; CPVITyxE(i)=CPVITyxE(i)-ITyxE; CPVITyyE(i)=CPVITyyE(i)-ITyyE; CPVITyzE(i)=CPVITyzE(i)-ITyzE; CPVITzxE(i)=CPVITzxE(i)-ITzxE; CPVITzyE(i)=CPVITzyE(i)-ITzyE; CPVITzzE(i)=CPVITzzE(i)-ITzzE; if i==fbcelements(j,1) CPVITxxE(i)=CPVITxxE(i)+ITxxE; CPVITxyE(i)=CPVITxyE(i)+ITxyE; CPVITxzE(i)=CPVITxzE(i)+ITxzE; CPVITyxE(i)=CPVITyxE(i)+ITyxE; CPVITyyE(i)=CPVITyyE(i)+ITyyE; CPVITyzE(i)=CPVITyzE(i)+ITyzE; CPVITzxE(i)=CPVITzxE(i)+ITzxE; CPVITzyE(i)=CPVITzyE(i)+ITzyE; CPVITzzE(i)=CPVITzzE(i)+ITzzE; end 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)-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; end end end for i=1:nelement for j=1:nfbcelements if i==fbcelements(j,1) K(3*i-2,3*fbcelements(j,1)-2)=CPVITxxE(i); K(3*i-2,3*fbcelements(j,1)-1)=CPVITxyE(i); K(3*i-2,3*fbcelements(j,1))=CPVITxzE(i); K(3*i-1,3*fbcelements(j,1)-2)=CPVITyxE(i); K(3*i-1,3*fbcelements(j,1)-1)=CPVITyyE(i); K(3*i-1,3*fbcelements(j,1))=CPVITyzE(i); K(3*i,3*fbcelements(j,1)-2)=CPVITzxE(i); K(3*i,3*fbcelements(j,1)-1)=CPVITzyE(i); K(3*i,3*fbcelements(j,1))=CPVITzzE(i); end end for j=1:ndispbcelements if i==dispbcelements(j,1) F(3*i-2)=F(3*i-2)+dispbcelements(j,2)*CPVITxxE(i)+dispbcelements(j,3)*CPVITxyE(i)+dispbcelements(j,4)*CPVITxzE(i); F(3*i-1)=F(3*i-1)+dispbcelements(j,2)*CPVITyxE(i)+dispbcelements(j,3)*CPVITyyE(i)+dispbcelements(j,4)*CPVITyzE(i); F(3*i)=F(3*i)+dispbcelements(j,2)*CPVITzxE(i)+dispbcelements(j,3)*CPVITzyE(i)+dispbcelements(j,4)*CPVITzzE(i); end end end %Solve U=K\F; % Display unknown dofs for each element U