% 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