function [error,ncp] = IGA_laplace1d(p,nel,k)

% A (very simple) 1D IGA code
%
% Solve -u'' = f with homogeneous Dirichlet b.c.'s and with a manufactured
%                body load s.t. the analytical solution is available
% 
% Input:  p    = approximation degree
%         nel  = number of (nonzero) elements
%         k    = parameter s.t. regularity is C^(p-k) [0<k<p+1]
%
% Output: error = (approximation of the) L2-norm error
%         ncp   = number of control points
%                 [ncp = nel*k + p + 1 - k]
%
%
%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         % By:                                            %
%         %          Alessandro Reali                      %
%         %          University of Pavia                   %
%         % email:   alessandro.reali@unipv.it             %
%         % webpage: http://www.unipv.it/alereali          %
%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set tolerance for checks
tol = 1e-15;

% set number of Gauss pts for integration
ngauss = p + 1;
if p > 9
    ngauss = 15;
    display('ngauss = 15... check Gauss rule for very high degrees!')
end

% set mesh [linear parameterization, continuity C^(p-k)]

%knot      = augknt([0 1],p+1);
knot      = [zeros(1,p+1),ones(1,p+1)];
h         = 1/nel;
ins       = sort(reshape((h:h:1-h)'*ones(1,k),1,k*(nel-1)));
[cp,knot] = bspkntins(p,0:1/p:1,knot,ins);

ncp = length(cp);
nspan = ncp - p;

% intializations
f_gl = zeros(ncp,1);
row  = zeros(1,ncp^2);
col  = zeros(1,ncp^2);
val  = zeros(1,ncp^2);

% precompute values of shape functions and derivatives at Gauss pts
[gp,gw] = gauss(ngauss);

gausspt = zeros(1,ngauss*nel);
n = 0;
for n0 = p+1:p+nspan;
    if (knot(n0+1)-knot(n0))>tol
        gausspt(n*ngauss+1:(n+1)*ngauss) =...
            ((knot(n0+1) - knot(n0))*gp + knot(n0+1) + knot(n0))/2;
        n = n + 1;
    end
end
N  = spcol(knot,p+1,sort([gausspt,gausspt]));
dN = N(2:2:end,:);
N  = N(1:2:end,:);
clear gausspt

% loop over elements for local array construction
icount = 0;
for iel = 1:nel
    dl = h/2;
    f_el = zeros(p+1,1);
    k_el = zeros(p+1,p+1);
    % loop over Gauss pts
    for igauss = 1:ngauss
        % compute shape functions and their derivatives with respect to
        % parametric variable, evaluated at current gauss point
        sh  =  N((iel-1)*ngauss+igauss,...
                    (k-1)*(iel-1)+iel+p:-1:(k-1)*(iel-1)+iel);
        dsh = dN((iel-1)*ngauss+igauss,...
                    (k-1)*(iel-1)+iel+p:-1:(k-1)*(iel-1)+iel);
        cpi = cp((k-1)*(iel-1)+iel+p:-1:(k-1)*(iel-1)+iel);
        
        % compute Gauss pt coordinate in phys space and jacobian @ Gauss pt       
        x      = sh*cpi';        
        J      = dsh*cpi';
        % compute shape fcn derivatives with respect to physical variable
        dshape = dsh/J;
        
        % compute weigth for integration
        gwt = gw(igauss)*J*dl;
        
        % manufactured body load
        f = 4*pi^2*sin(2*pi*x);
            
        % compute local r.h.s. vector and stiffness
        f_el = f_el + sh'*f*gwt;
        k_el = k_el + dsh'*dsh*gwt;
    end
    
    % assemble r.h.s. and prepare vectors for sparse stiffness assembly
    for a = 1:p+1
        i = (k-1)*(iel-1) + iel + p - (a-1);
        if (i>1 && i<ncp)
            f_gl(i) = f_gl(i) + f_el(a);
            for b = 1:p+1
                j = (k-1)*(iel-1) + iel + p - (b-1);
                icount = icount + 1;                
                row(icount) = i;
                col(icount) = j;
                val(icount) = k_el(a,b);
            end
        else % apply homogeneous Dirichlet b.c.'s
            icount = icount + 1;
            row(icount) = i;
            col(icount) = i;
            val(icount) = 1;
        end
    end
end

row = row(1:icount);
col = col(1:icount);
val = val(1:icount);

% assemble stiffness
k_gl = sparse(row,col,val,ncp,ncp);
clear row col val

% solve linear system
sol = k_gl\f_gl;
clear k_gl f_gl


% approx relative L^2-norm error vs the manufactured analytical solution
xx = linspace(0,1,2*ncp);
approx = bspeval(p,sol',knot,xx);
exact = sin(2*pi*xx);

error = sqrt(sum((exact-approx).^2)/sum(exact.^2));

return 


function [GP,GW] = gauss(ng)
% Gauss pts' coordinates and weights

if (ng==1)
   GP(1) = 0;
   GW(1) = 2;
elseif (ng==2)
   GP(1) = -.5773502691896259;
   GP(2) = .5773502691896259;
   GW(1) = 1;
   GW(2) = 1;
elseif (ng==3)
   GP(1) = -.7745966692414834;
   GP(2) = 0;
   GP(3) = .7745966692414834;
   GW(1) = .5555555555555556;
   GW(2) = .8888888888888889;
   GW(3) = .5555555555555556;
elseif (ng==4)
   GP(1) = -.86113631159405257524;
   GP(2) = -.33998104358485626481;
   GP(3) = .33998104358485626481;
   GP(4) = .86113631159405257524;
   GW(1) = .34785484513745385736;
   GW(2) = .65214515486254614264;
   GW(3) = .65214515486254614264;
   GW(4) = .34785484513745385736;
elseif (ng==5)
   GP(1) = -.90617984593866399282;
   GP(2) = -.53846931010568309105;
   GP(3) = 0;
   GP(4) = .53846931010568309105;
   GP(5) = .90617984593866399282;
   GW(1) = .23692688505618908749;
   GW(2) = .47862867049936646808;
   GW(3) = .56888888888888888888;
   GW(4) = .47862867049936646808;
   GW(5) = .23692688505618908749;
elseif (ng==6)
   GP(1) = -.9324695142031520;
   GP(2) = -.6612093864662645;
   GP(3) = -.2386191860831969;
   GP(4) = .2386191860831969;
   GP(5) = .6612093864662645;
   GP(6) = .9324695142031520;
   GW(1) = .1713244923791703;
   GW(2) = .3607615730481386;
   GW(3) = .4679139345726911;
   GW(4) = .4679139345726911;
   GW(5) = .3607615730481386;
   GW(6) = .1713244923791703;
elseif (ng==7)
   GP(1) = -.9491079123427585;
   GP(2) = -.7415311855993944;
   GP(3) = -.4058451513773972;
   GP(4) = 0;
   GP(5) = .4058451513773972;
   GP(6) = .7415311855993944;
   GP(7) = .9491079123427585;
   GW(1) = .1294849661688697;
   GW(2) = .2797053914892767;
   GW(3) = .3818300505051189;
   GW(4) = .4179591836734694;
   GW(5) = .3818300505051189;
   GW(6) = .2797053914892767;
   GW(7) = .1294849661688697;
elseif (ng==8)
   GP(1) = -.9602898564975362;
   GP(2) = -.7966664774136267;
   GP(3) = -.5255324099163290;
   GP(4) = -.1834346424956498;
   GP(5) = .1834346424956498;
   GP(6) = .5255324099163290;
   GP(7) = .7966664774136267;
   GP(8) = .9602898564975362;
   GW(1) = .1012285362903763;
   GW(2) = .2223810344533745;
   GW(3) = .3137066458778873;
   GW(4) = .3626837833783620;
   GW(5) = .3626837833783620;
   GW(6) = .3137066458778873;
   GW(7) = .2223810344533745;
   GW(8) = .1012285362903763;
elseif (ng==9)
   GP(1) = -.9681602395076261;
   GP(2) = -.8360311073266358;
   GP(3) = -.6133714327005904;
   GP(4) = -.3242534234038089;
   GP(5) = 0;
   GP(6) = .3242534234038089;
   GP(7) = .6133714327005904;
   GP(8) = .8360311073266358;
   GP(9) = .9681602395076261;
   GW(1) = .0812743883615744;
   GW(2) = .1806481606948574;
   GW(3) = .2606106964029354;
   GW(4) = .3123470770400029;
   GW(5) = .3302393550012598;
   GW(6) = .3123470770400028;
   GW(7) = .2606106964029355;
   GW(8) = .1806481606948574;
   GW(9) = .0812743883615744;
elseif (ng==10)
   GP(1) = -.973906528517172;
   GP(2) = -.865063366688985;
   GP(3) = -.679409568299024;
   GP(4) = -.433395394129247;
   GP(5) = -.148874338981631;
   GP(6) = .148874338981631;
   GP(7) = .433395394129247;
   GP(8) = .679409568299024;
   GP(9) = .865063366688985;
   GP(10) = .973906528517172;
   GW(1) = .066671344308688;
   GW(2) = .149451349150581;
   GW(3) = .219086362515982;
   GW(4) = .269266719309996;
   GW(5) = .295524224714753;
   GW(6) = .295524224714753;
   GW(7) = .269266719309996;
   GW(8) = .219086362515982;
   GW(9) = .149451349150581;
   GW(10) = .066671344308688;
elseif (ng==15)
   GP(1) = -.9879925180204854;
   GP(2) = -.9372733924007059;
   GP(3) = -.8482065834104272;
   GP(4) = -.7244177313601700;
   GP(5) = -.5709721726085388;
   GP(6) = -.3941513470775634;
   GP(7) = -.2011940939974345;
   GP(8) = 0;
   GP(9) = .2011940939974345;
   GP(10) = .3941513470775634;
   GP(11) = .5709721726085388;
   GP(12) = .7244177313601700;
   GP(13) = .8482065834104272;
   GP(14) = .9372733924007059;
   GP(15) = .9879925180204854;
   GW(1) = .03075324199611807;
   GW(2) = .07036604748811134;
   GW(3) = .1071592204671351;
   GW(4) = .1395706779261761;
   GW(5) = .1662692058169852;
   GW(6) = .1861610000155741;
   GW(7) = .1984314853271374;
   GW(8) = .2025782419255562;
   GW(9) = .1984314853271374;
   GW(10) = .1861610000155741;
   GW(11) = .1662692058169852;
   GW(12) = .1395706779261761;
   GW(13) = .1071592204671351;
   GW(14) = .07036604748811134;
   GW(15) = .03075324199611807;
end

return