I have copied the code below after all, this should be easier for everyone. I have tried to speed it up myself using parfor, with some improvement, but not nearly enough.
%Warning: this script will clear and therefore erase your stored variables. If you do not
%desire this, delete the next two lines of code, but beware that the
%script may not function properly.
clear;
clc;
%Define geometry
L=84e3; %length
W=24e3; %width
%Define the elasticity tensor
C1=[303973359.614198 89929777.3924441 0 0 0 0 0 0;89929777.3924441 303973359.614198 0 0 0 0 0 0;0 0 107021791.110877 0 0 0 0 0;0 0 0 497.738940473395 80.7354743011537 15.9810696297077 0 0;0 0 0 80.7354743011537 391.198476275344 15.9810696297077 0 0;0 0 0 15.9810696297077 15.9810696297077 103.524825925731 0 0;0 0 0 0 0 0 20440000 0;0 0 0 0 0 0 0 20440000];
%Specify the degree of the approximating basis polynomials in x/y
m=3;
n=3;
%Construct basis polynomials
syms x y z
B11 = sym(zeros(m+1, 1));
B21 = sym(zeros(m+1, 1));
B12 = sym(zeros(n+1, 1));
B22 = sym(zeros(n+1, 1));
for i = 0:m
A=factorial(m)/(factorial(i)*factorial(mi))*(x/L)^i*(1x/L)^(mi);
B11(i+1)=A;
end
for i = 0:n
B=factorial(n)/(factorial(i)*factorial(ni))*(y/W)^i*(1y/W)^(ni);
B12(i+1)=B;
end
%Construct displacement vectors using the basis polynomials in x/y
zero_vector=zeros((m+1)*(n+1),1);
counter=0;
Vx1 = sym(zeros((m+1)*(n+1),1));
Vy1 = sym(zeros((m+1)*(n+1),1));
Vz1 = sym(zeros((m+1)*(n+1),1));
Vrx1 = sym(zeros((m+1)*(n+1),1));
Vry1 = sym(zeros((m+1)*(n+1),1));
for i = 0:m
for j=0:n
counter=counter+1;
C=B11(i+1)*B12(j+1);
Vx1(counter,1)=C; %Bernstein bases for xdisplacement
Vy1(counter,1)=C; %Bernstein bases for ydisplacement
Vz1(counter,1)=C; %Bernstein bases for zdisplacement
Vrx1(counter,1)=C; %Bernstein bases for xrotation
Vry1(counter,1)=C; %Bernstein bases for yrotation
end
end
Vx1=[Vx1;zero_vector;zero_vector;zero_vector;zero_vector]; %Displacement in x
Vy1=[zero_vector;Vy1;zero_vector;zero_vector;zero_vector]; %Displacement in y
Vz1=[zero_vector;zero_vector;Vz1;zero_vector;zero_vector]; %Displacement in z
Vrx1=[zero_vector;zero_vector;zero_vector;Vrx1;zero_vector]; %Rotation in x
Vry1=[zero_vector;zero_vector;zero_vector;zero_vector;Vry1]; %Rotation in y
%Define the strain tensor
Vxx1=diff(Vx1,x);
Vyy1=diff(Vy1,y);
Vxy1=diff(Vx1,y);
Vyx1=diff(Vy1,x);
Vzx1=diff(Vz1,x);
Vzy1=diff(Vz1,y);
Vrxx1=diff(Vrx1,x);
Vryy1=diff(Vry1,y);
Vrxy1=diff(Vrx1,y);
Vryx1=diff(Vry1,x);
B_epsilon_1=[Vxx1.';Vyy1.';Vxy1.'+Vyx1.']; %inplane strains
B_kappa_1=[Vrxx1.';Vryy1.';Vrxy1.'+Vryx1.']; %curvatures
B_gamma_1=[Vzy1.'Vry1.';Vzx1.'Vrx1.']; %shear strains
B_1=[B_epsilon_1;B_kappa_1;B_gamma_1]; %overall strain tensor
%Calculate stiffness matrix
%**This is where the longest part of the code is**
tic
K=int(int(B_1.'*C1*B_1,x,0,L),y,0,W); %Perform area integral over geometry
K=double(K); %Convert the matrix obtained to doubles
toc
%**This is where the longest part of the code ends**
%Faster stiffness matrix calculation
if matlabpool('size') == 0
matlabpool open 4 %Assign 4 cores to task
end
tic
K1_unintegrated=triu(B_1.'*C1*B_1);
K1_unintegrated = reshape(K1_unintegrated,1,(5*(m+1)*(n+1))^2);
Batch_length=(5*(m+1)*(n+1))^2;
addendum=0;
while mod(Batch_length,4)>0
addendum=addendum+1;
Batch_length=Batch_length+1;
end
Batch_in=sym(zeros(4,Batch_length/4));
Batch_in(1,:)=K1_unintegrated(1:Batch_length/4);
Batch_in(2,:)=K1_unintegrated(Batch_length/4+1:2*Batch_length/4);
Batch_in(3,:)=K1_unintegrated(2*Batch_length/4+1:3*Batch_length/4);
Batch_in(4,1:Batch_length/4addendum)=K1_unintegrated(3*Batch_length/4+1:4*Batch_length/4addendum);
Batch_out=sym(zeros(4,Batch_length/4));
parfor i=1:4
Batch_out(i,:)=int(int(Batch_in(i,:),x,0,L),y,0,W);
end
K1_fast =[Batch_out(1,:),Batch_out(2,:),Batch_out(3,:),Batch_out(4,1:Batch_length/4addendum)];
K1_fast = reshape(K1_fast,5*(m+1)*(n+1),5*(m+1)*(n+1));
K1_fast = K1_fast+triu(K1_fast,1).';
K1_fast=double(K1_fast);
toc
