Splines basis of order p=2

Contents

Knot vector on the reference element and "bspline"-command

knots_in = [0 1 2 3];
knots_le = [0 0 1 2];
knots_ri = [1 2 3 3];

ref_spl_in = bspline(knots_in);
ref_spl_le = bspline(knots_le);
ref_spl_ri = bspline(knots_ri);

The "bspline"-command yields the piece-wise polynomials of the spline function defined by the knot vector, with coefficients such that each piece is defined on [0,1] (shape functions)!!

phi1 = ref_spl_in.coefs(1,:);
phi0 = ref_spl_in.coefs(2,:);
phi2 = ref_spl_in.coefs(3,:);

psi0 = ref_spl_le.coefs(1,:);
psi2 = ref_spl_le.coefs(2,:);

chi1 = ref_spl_ri.coefs(1,:);
chi0 = ref_spl_ri.coefs(2,:);

Dphi1 = polyder(phi1);
Dphi0 = polyder(phi0);
Dphi2 = polyder(phi2);

Dpsi0 = polyder(psi0);
Dpsi2 = polyder(psi2);

Dchi1 = polyder(chi1);
Dchi0 = polyder(chi0);

bsp_in = @(t) ppval(ref_spl_in,t);
bsp_le = @(t) ppval(ref_spl_le,t);
bsp_ri = @(t) ppval(ref_spl_ri,t);

absz_in = linspace(0,3,100);
absz_le = linspace(0,2,100);
absz_ri = linspace(1,3,100);
plot( absz_in, bsp_in(absz_in),'b', absz_le, bsp_le(absz_le),'g', absz_ri, bsp_ri(absz_ri),'r' )
title([swi_meth])
legend('interior','left','right')

Integrals on the reference element

% Mass interior:
pr = polyint(conv(phi0,phi0));
M00 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(phi1,phi1));
M11 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(phi2,phi2));
M22 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(phi0,phi1));
M01 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(phi0,phi2));
M02 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(phi1,phi2));
M12 = polyval(pr,1) - polyval(pr,0);

% Mass left:
pr = polyint(conv(psi0,psi0));
le_M00 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(psi2,psi2));
le_M22 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(psi0,phi1));
le_M01 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(psi2,phi0));
le_M02 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(psi2,phi1));
le_M12 = polyval(pr,1) - polyval(pr,0);

% Mass right:
pr = polyint(conv(chi0,chi0));
ri_M00 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(chi1,chi1));
ri_M11 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(chi1,phi0));
ri_M01 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(chi0,phi2));
ri_M02 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(chi1,phi2));
ri_M12 = polyval(pr,1) - polyval(pr,0);

% Stiffness interior:
pr = polyint(conv(Dphi0,Dphi0));
S00 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dphi1,Dphi1));
S11 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dphi2,Dphi2));
S22 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dphi0,Dphi1));
S01 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dphi0,Dphi2));
S02 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dphi1,Dphi2));
S12 = polyval(pr,1) - polyval(pr,0);

% Stiffness left:
pr = polyint(conv(Dpsi0,Dpsi0));
le_S00 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dpsi2,Dpsi2));
le_S22 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dpsi0,Dphi1));
le_S01 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dpsi2,Dphi0));
le_S02 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dpsi2,Dphi1));
le_S12 = polyval(pr,1) - polyval(pr,0);

% Stiffness right:
pr = polyint(conv(Dchi0,Dchi0));
ri_S00 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dchi1,Dchi1));
ri_S11 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dchi1,Dphi0));
ri_S01 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dchi0,Dphi2));
ri_S02 = polyval(pr,1) - polyval(pr,0);

pr = polyint(conv(Dchi1,Dphi2));
ri_S12 = polyval(pr,1) - polyval(pr,0);

N = M; % dimension of the FEM space (= M+p-2, because of p-1 regularity).
% A basis function can be associated to each element.

C = x(2:end) - x(1:end-1); % vector of cell volumes

nodes = zeros(1,M-1); % node vector
nodes = x(2:end-1);

Sparse matrix assembly

nz = 1; % counter of non-zero elements (+1)

% boundary conditions on the left:
i = 1; % all integrales named "le_..."
ivec(nz) = i; % vector of row indices
jvec(nz) = i; % vector of column indices
Svec(nz) = 1/C(i)*le_S00 + 1/C(i+1)*le_S22 ; % vector of stiffness matrix entries
Mvec(nz) = C(i)*le_M00 + C(i+1)*le_M22; % vector of mass matrix entries
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+1;
Svec(nz) = 1/C(i)*le_S01 + 1/C(i+1)*le_S02;;
Mvec(nz) = C(i)*le_M01 + C(i+1)*le_M02;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+2;
Svec(nz) = 1/C(i+1)*le_S12;
Mvec(nz) = C(i+1)*le_M12;
nz = nz+1;

i=2; % boundary
ivec(nz) = i;
jvec(nz) = i-1;
Svec(nz) = 1/C(i-1)*le_S01 + 1/C(i)*le_S02; % here are integrales named "le_..."
Mvec(nz) = C(i-1)*le_M01 + C(i)*le_M02; % here are integrales named "le_..."
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i;
Svec(nz) = 1/C(i-1)*S11 + 1/C(i)*S00 + 1/C(i+1)*S22;
Mvec(nz) = C(i-1)*M11 + C(i)*M00 + C(i+1)*M22;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+1;
Svec(nz) = 1/C(i)*S01 + 1/C(i+1)*S02;
Mvec(nz) = C(i)*M01 + C(i+1)*M02;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+2;
Svec(nz) = 1/C(i+1)*S12;
Mvec(nz) = C(i+1)*M12;
nz = nz+1;

i=3; % boundary
ivec(nz) = i;
jvec(nz) = i-2;
Svec(nz) = 1/C(i-1)*le_S12; % here are integrales named "le_..."
Mvec(nz) = C(i-1)*le_M12; % here are integrales named "le_..."
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i-1;
Svec(nz) = 1/C(i-1)*S01 + 1/C(i)*S02;
Mvec(nz) = C(i-1)*M01 + C(i)*M02;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i;
Svec(nz) = 1/C(i-1)*S11 + 1/C(i)*S00 + 1/C(i+1)*S22;
Mvec(nz) = C(i-1)*M11 + C(i)*M00 + C(i+1)*M22;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+1;
Svec(nz) = 1/C(i)*S01 + 1/C(i+1)*S02;
Mvec(nz) = C(i)*M01 + C(i+1)*M02;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+2;
Svec(nz) = 1/C(i+1)*S12;
Mvec(nz) = C(i+1)*M12;
nz = nz+1;

for i=4:N-3

    ivec(nz) = i;
    jvec(nz) = i-2;
    Svec(nz) = 1/C(i-1)*S12;
    Mvec(nz) = C(i-1)*M12;
    nz = nz+1;
    ivec(nz) = i;
    jvec(nz) = i-1;
    Svec(nz) = 1/C(i-1)*S01 + 1/C(i)*S02;
    Mvec(nz) = C(i-1)*M01 + C(i)*M02;
    nz = nz+1;
    ivec(nz) = i;
    jvec(nz) = i;
    Svec(nz) = 1/C(i-1)*S11 + 1/C(i)*S00 + 1/C(i+1)*S22;
    Mvec(nz) = C(i-1)*M11 + C(i)*M00 + C(i+1)*M22;
    nz = nz+1;
    ivec(nz) = i;
    jvec(nz) = i+1;
    Svec(nz) = 1/C(i)*S01 + 1/C(i+1)*S02;
    Mvec(nz) = C(i)*M01 + C(i+1)*M02;
    nz = nz+1;
    ivec(nz) = i;
    jvec(nz) = i+2;
    Svec(nz) = 1/C(i+1)*S12;
    Mvec(nz) = C(i+1)*M12;
    nz = nz+1;

end

% boundary conditions on the right:
i = N-2; % boundary
ivec(nz) = i;
jvec(nz) = i-2;
Svec(nz) = 1/C(i-1)*S12;
Mvec(nz) = C(i-1)*M12;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i-1;
Svec(nz) = 1/C(i-1)*S01 + 1/C(i)*S02;
Mvec(nz) = C(i-1)*M01 + C(i)*M02;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i;
Svec(nz) = 1/C(i-1)*S11 + 1/C(i)*S00 + 1/C(i+1)*S22;
Mvec(nz) = C(i-1)*M11 + C(i)*M00 + C(i+1)*M22;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+1;
Svec(nz) = 1/C(i)*S01 + 1/C(i+1)*S02;
Mvec(nz) = C(i)*M01 + C(i+1)*M02;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+2;
Svec(nz) = 1/C(i+1)*ri_S12; % here are integrales named "ri_..."
Mvec(nz) = C(i+1)*ri_M12; % here are integrales named "ri_..."
nz = nz+1;

i = N-1; % boundary
ivec(nz) = i;
jvec(nz) = i-2;
Svec(nz) = 1/C(i-1)*S12;
Mvec(nz) = C(i-1)*M12;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i-1;
Svec(nz) = 1/C(i-1)*S01 + 1/C(i)*S02;
Mvec(nz) = C(i-1)*M01 + C(i)*M02;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i;
Svec(nz) = 1/C(i-1)*S11 + 1/C(i)*S00 + 1/C(i+1)*S22;
Mvec(nz) = C(i-1)*M11 + C(i)*M00 + C(i+1)*M22;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+1;
Svec(nz) = 1/C(i)*ri_S01 + 1/C(i+1)*ri_S02; % here are integrales named "ri_..."
Mvec(nz) = C(i)*ri_M01 + C(i+1)*ri_M02; % here are integrales named "ri_..."
nz = nz+1;

i = N; % all integrales named "ri_..."
ivec(nz) = i;
jvec(nz) = i-2;
Svec(nz) = 1/C(i-1)*ri_S12;
Mvec(nz) = C(i-1)*ri_M12;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i-1;
Svec(nz) = 1/C(i-1)*ri_S01 + 1/C(i)*ri_S02;
Mvec(nz) = C(i-1)*ri_M01 + C(i)*ri_M02;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i;
Svec(nz) = 1/C(i-1)*ri_S11 + 1/C(i)*ri_S00;
Mvec(nz) = C(i-1)*ri_M11 + C(i)*ri_M00;
nz = nz+1;

A = sparse(ivec,jvec,Svec + Mvec); % system matrix
figure
spy(A)
title('A')
Mass = sparse(ivec,jvec,Mvec); % mass matrix
figure
spy(Mass)
title('Mass matrix')

Right-hand side vector

syms xi
if swi_rhs=='exact'

    i = 1;
    clear z igrand han int
    z = x(i) + xi*C(i); % is the variable transformation R_{i}
    igrand(xi) = eval(rho)*poly2sym(psi0,xi); % the integrand in xi-space
    han = matlabFunction(igrand(xi));
    int = C(i)*integral(han,0,1); % numerical integration

    rhovec(i) = int;

    clear z igrand han int
    z = x(i+1) + xi*C(i+1); % is the variable transformation R_{i+1}
    igrand(xi) = eval(rho)*poly2sym(psi2,xi); % the integrand in xi-space
    han = matlabFunction(igrand(xi));
    int = C(i+1)*integral(han,0,1); % numerical integration

    rhovec(i) = rhovec(i) + int;

    for i = 2:N-1

        clear z igrand han int
        z = x(i-1) + xi*C(i-1); % is the variable transformation R_{i-1}
        igrand(xi) = eval(rho)*poly2sym(phi1,xi); % the integrand in xi-space
        han = matlabFunction(igrand(xi));

        int = C(i-1)*integral(han,0,1); % numerical integration

        rhovec(i) = int;

        clear z igrand han int
        z = x(i) + xi*C(i); % is the variable transformation R_{i}
        igrand(xi) = eval(rho)*poly2sym(phi0,xi); % the integrand in xi-space
        han = matlabFunction(igrand(xi));
        int = C(i)*integral(han,0,1); % numerical integration

        rhovec(i) = rhovec(i) + int;

        clear z igrand han int
        z = x(i+1) + xi*C(i+1); % is the variable transformation R_{i+1}
        igrand(xi) = eval(rho)*poly2sym(phi2,xi); % the integrand in xi-space
        han = matlabFunction(igrand(xi));
        int = C(i+1)*integral(han,0,1); % numerical integration

        rhovec(i) = rhovec(i) + int;

    end

    i = N;
    clear z igrand han int
    z = x(i-1) + xi*C(i-1); % is the variable transformation R_{i-1}
    igrand(xi) = eval(rho)*poly2sym(chi1,xi); % the integrand in xi-space
    han = matlabFunction(igrand(xi));

    int = C(i-1)*integral(han,0,1); % numerical integration

    rhovec(i) = int;

    clear z igrand han int
    z = x(i) + xi*C(i); % is the variable transformation R_{i}
    igrand(xi) = eval(rho)*poly2sym(chi0,xi); % the integrand in xi-space
    han = matlabFunction(igrand(xi));
    int = C(i)*integral(han,0,1); % numerical integration

    rhovec(i) = rhovec(i) + int;

elseif swi_rhs=='basis'

    z = nodes; % set the symbolic variable
    rhovec = (Mass * eval(rho)')';

end

Solve with direct solver

u = A \ rhovec';

Get the values of $u_h$ and $u_h'$ at the vertices and at the midpoints:

% boundary:
i = 1;
u_000(i) = u(i)*polyval(psi0,0.00) + u(i+1)*polyval(phi1,0.00);
u_025(i) = u(i)*polyval(psi0,0.25) + u(i+1)*polyval(phi1,0.25);
u_050(i) = u(i)*polyval(psi0,0.50) + u(i+1)*polyval(phi1,0.50);
u_075(i) = u(i)*polyval(psi0,0.75) + u(i+1)*polyval(phi1,0.75);
u_100(i) = u(i)*polyval(psi0,1.00) + u(i+1)*polyval(phi1,1.00);

Du_000(i) = 1/C(i)*(u(i)*polyval(Dpsi0,0.00) + u(i+1)*polyval(Dphi1,0.00));
Du_025(i) = 1/C(i)*(u(i)*polyval(Dpsi0,0.25) + u(i+1)*polyval(Dphi1,0.25));
Du_050(i) = 1/C(i)*(u(i)*polyval(Dpsi0,0.50) + u(i+1)*polyval(Dphi1,0.50));
Du_075(i) = 1/C(i)*(u(i)*polyval(Dpsi0,0.75) + u(i+1)*polyval(Dphi1,0.75));
Du_100(i) = 1/C(i)*(u(i)*polyval(Dpsi0,1.00) + u(i+1)*polyval(Dphi1,1.00));

% boundary:
i = 2;
u_000(i) = u(i-1)*polyval(psi2,0.00) + u(i)*polyval(phi0,0.00) + u(i+1)*polyval(phi1,0.00);
u_025(i) = u(i-1)*polyval(psi2,0.25) + u(i)*polyval(phi0,0.25) + u(i+1)*polyval(phi1,0.25);
u_050(i) = u(i-1)*polyval(psi2,0.50) + u(i)*polyval(phi0,0.50) + u(i+1)*polyval(phi1,0.50);
u_075(i) = u(i-1)*polyval(psi2,0.75) + u(i)*polyval(phi0,0.75) + u(i+1)*polyval(phi1,0.75);
u_100(i) = u(i-1)*polyval(psi2,1.00) + u(i)*polyval(phi0,1.00) + u(i+1)*polyval(phi1,1.00);

Du_000(i) = 1/C(i)*(u(i-1)*polyval(Dpsi2,0.00) + u(i)*polyval(Dphi0,0.00) + u(i+1)*polyval(Dphi1,0.00));
Du_025(i) = 1/C(i)*(u(i-1)*polyval(Dpsi2,0.25) + u(i)*polyval(Dphi0,0.25) + u(i+1)*polyval(Dphi1,0.25));
Du_050(i) = 1/C(i)*(u(i-1)*polyval(Dpsi2,0.50) + u(i)*polyval(Dphi0,0.50) + u(i+1)*polyval(Dphi1,0.50));
Du_075(i) = 1/C(i)*(u(i-1)*polyval(Dpsi2,0.75) + u(i)*polyval(Dphi0,0.75) + u(i+1)*polyval(Dphi1,0.75));
Du_100(i) = 1/C(i)*(u(i-1)*polyval(Dpsi2,1.00) + u(i)*polyval(Dphi0,1.00) + u(i+1)*polyval(Dphi1,1.00));

for i=3:N-2

    u_000(i) = u(i-1)*polyval(phi2,0.00) + u(i)*polyval(phi0,0.00) + u(i+1)*polyval(phi1,0.00);
    u_025(i) = u(i-1)*polyval(phi2,0.25) + u(i)*polyval(phi0,0.25) + u(i+1)*polyval(phi1,0.25);
    u_050(i) = u(i-1)*polyval(phi2,0.50) + u(i)*polyval(phi0,0.50) + u(i+1)*polyval(phi1,0.50);
    u_075(i) = u(i-1)*polyval(phi2,0.75) + u(i)*polyval(phi0,0.75) + u(i+1)*polyval(phi1,0.75);
    u_100(i) = u(i-1)*polyval(phi2,1.00) + u(i)*polyval(phi0,1.00) + u(i+1)*polyval(phi1,1.00);

    Du_000(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,0.00) + u(i)*polyval(Dphi0,0.00) + u(i+1)*polyval(Dphi1,0.00));
    Du_025(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,0.25) + u(i)*polyval(Dphi0,0.25) + u(i+1)*polyval(Dphi1,0.25));
    Du_050(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,0.50) + u(i)*polyval(Dphi0,0.50) + u(i+1)*polyval(Dphi1,0.50));
    Du_075(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,0.75) + u(i)*polyval(Dphi0,0.75) + u(i+1)*polyval(Dphi1,0.75));
    Du_100(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,1.00) + u(i)*polyval(Dphi0,1.00) + u(i+1)*polyval(Dphi1,1.00));

end

% boundary:
i = N-1;
u_000(i) = u(i-1)*polyval(phi2,0.00) + u(i)*polyval(phi0,0.00) + u(i+1)*polyval(chi1,0.00);
u_025(i) = u(i-1)*polyval(phi2,0.25) + u(i)*polyval(phi0,0.25) + u(i+1)*polyval(chi1,0.25);
u_050(i) = u(i-1)*polyval(phi2,0.50) + u(i)*polyval(phi0,0.50) + u(i+1)*polyval(chi1,0.50);
u_075(i) = u(i-1)*polyval(phi2,0.75) + u(i)*polyval(phi0,0.75) + u(i+1)*polyval(chi1,0.75);
u_100(i) = u(i-1)*polyval(phi2,1.00) + u(i)*polyval(phi0,1.00) + u(i+1)*polyval(chi1,1.00);

Du_000(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,0.00) + u(i)*polyval(Dphi0,0.00) + u(i+1)*polyval(Dchi1,0.00));
Du_025(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,0.25) + u(i)*polyval(Dphi0,0.25) + u(i+1)*polyval(Dchi1,0.25));
Du_050(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,0.50) + u(i)*polyval(Dphi0,0.50) + u(i+1)*polyval(Dchi1,0.50));
Du_075(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,0.75) + u(i)*polyval(Dphi0,0.75) + u(i+1)*polyval(Dchi1,0.75));
Du_100(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,1.00) + u(i)*polyval(Dphi0,1.00) + u(i+1)*polyval(Dchi1,1.00));

% boundary:
i = N;
u_000(i) = u(i-1)*polyval(phi2,0.00) + u(i)*polyval(chi0,0.00);
u_025(i) = u(i-1)*polyval(phi2,0.25) + u(i)*polyval(chi0,0.25);
u_050(i) = u(i-1)*polyval(phi2,0.50) + u(i)*polyval(chi0,0.50);
u_075(i) = u(i-1)*polyval(phi2,0.75) + u(i)*polyval(chi0,0.75);
u_100(i) = u(i-1)*polyval(phi2,1.00) + u(i)*polyval(chi0,1.00);

Du_000(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,0.00) + u(i)*polyval(Dchi0,0.00));
Du_025(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,0.25) + u(i)*polyval(Dchi0,0.25));
Du_050(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,0.50) + u(i)*polyval(Dchi0,0.50));
Du_075(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,0.75) + u(i)*polyval(Dchi0,0.75));
Du_100(i) = 1/C(i)*(u(i-1)*polyval(Dphi2,1.00) + u(i)*polyval(Dchi0,1.00));

u = [u_000 u_100(end)];
Du_plot = [Du_000 Du_100(end)];