Splines basis of order p=1

Contents

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

knots = [0 1 2];

ref_spl = bspline(knots);

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)!!

phiM1 = ref_spl.coefs(1,:);
phi0 = ref_spl.coefs(2,:);

Dphi0 = polyder(phi0);
DphiM1 = polyder(phiM1);

bsp = @(t) ppval(ref_spl,t);

absziss = linspace(0,2,100);
plot( absziss, bsp(absziss) )
title([swi_meth])

Integrals on the reference element

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

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

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

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

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

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

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

C = x(2:end) - x(1:end-1); % vector of cell volumes
% Attention: the first element of C is C_0 !! --> C(i) = C_{i-1}

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

Sparse matrix assembly

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

i = 1; % boundary conditions on the left:
ivec(nz) = i; % vector of row indices
jvec(nz) = i; % vector of column indices
Svec(nz) = 1/C(i)*S11 + 1/C(i+1)*S00 ; % vector of stiffness matrix entries
Mvec(nz) = C(i)*M11 + C(i+1)*M00; % vector of mass matrix entries
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i+1;
Svec(nz) = 1/C(i+1)*S01;
Mvec(nz) = C(i+1)*M01;
nz = nz+1;

for i=2:N-1

    ivec(nz) = i;
    jvec(nz) = i;
    Svec(nz) = 1/C(i)*S11 + 1/C(i+1)*S00;
    Mvec(nz) = C(i)*M11 + C(i+1)*M00;
    nz = nz+1;
    ivec(nz) = i;
    jvec(nz) = i+1;
    Svec(nz) = 1/C(i+1)*S01;
    Mvec(nz) = C(i+1)*M01;
    nz = nz+1;
    ivec(nz) = i;
    jvec(nz) = i-1;
    Svec(nz) = 1/C(i)*S01;
    Mvec(nz) = C(i)*M01;
    nz = nz+1;

end

% boundary conditions on the right:
i = N;
ivec(nz) = i;
jvec(nz) = i;
Svec(nz) = 1/C(i)*S11 + 1/C(i+1)*S00;
Mvec(nz) = C(i)*M11 + C(i+1)*M00;
nz = nz+1;
ivec(nz) = i;
jvec(nz) = i-1;
Svec(nz) = 1/C(i)*S01;
Mvec(nz) = C(i)*M01;
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('M')

Right-hand side vector

if swi_rhs=='exact'

    syms xi
    for i = 1:N

        clear z
        z = x(i) + xi*C(i); % set the symbolic variable z to new symbolic
        igrand1(xi) = eval(rho)*xi; % the integrand in xi-space
        han1 = matlabFunction(igrand1(xi));
        int1 = C(i)*integral(han1,0,1); % numerical integration
%             int1 = int(igrand,xi,0,1); % analytic integration (very slow)
%             int1 = C(k)*eval(int1);

        z = x(i+1) + xi*C(i+1);
        igrand2(xi) = eval(rho)*(1-xi);
        han2 = matlabFunction(igrand2(xi));
        int2 = C(i+1)*integral(han2,0,1);
%             int2 = int(igrand2,xi,0,1);
%             int2 = C(k+1)*eval(int2);

        rhovec(i) = int1 + int2;

    end

elseif swi_rhs=='basis'

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

end

Solve with direct solver

u = A \ rhovec';
u = [0 u' 0]; % adding boundary conditions (possible because the basis is in fact a Lagrange basis).

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

u_000 = u(1:end-1);
u_025 = 3/4*u(1:end-1) + 1/4*u(2:end);
u_050 = 1/2*u(1:end-1) + 1/2*u(2:end);
u_075 = 1/4*u(1:end-1) + 3/4*u(2:end);
u_100 = u(2:end);

Du = (u(2:end) - u(1:end-1))./C; % independent of xi
Du_000 = Du;
Du_025 = Du;
Du_050 = Du;
Du_075 = Du;
Du_100 = Du;

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