function precal_gauss
%Calculating the values at gauss points
%
global ini_Val Invariant
U = ini_Val.U;
p = ini_Val.p;
d = ini_Val.d;
gaussno = ini_Val.gaussno;
a = U';
a(1:p)=[];
a(end-p+1:end)=[];
b=a;
a(length(a)) = [];
b(1) = [];
[gpoint,gweight_local]=gaquad(gaussno);
h = ((b-a)/2);

Nshape = zeros(length(a)*gaussno,(p+1));
dNsh = zeros(length(a)*gaussno,(p+1));
ddNsh = zeros(length(a)*gaussno,(p+1));
dddNsh = zeros(length(a)*gaussno,(p+1));
gweight=zeros(length(a)*gaussno,1);
connect=zeros(length(a),p+1);
iglobal=0;
for i=1:length(a)
    for k=1:gaussno
        iglobal=iglobal+1;
        gweight(iglobal)=gweight_local(k)*h(i);
        gaussp(iglobal) =((b(i)-a(i))/2)*gpoint(k)+((a(i)+b(i))/2);
        [Nsh,sp]=evaluate_BSp(U,p,d,gaussp(iglobal));
        Nshape(iglobal,:)=Nsh(1,:);
        dNsh(iglobal,:)=Nsh(2,:);
        ddNsh(iglobal,:)=Nsh(3,:);
        if ini_Val.d == 3
            dddNsh(iglobal,:)=Nsh(4,:);
        end
        spv(iglobal)=sp-p;
    end
    connect(i,:)=[i:i+p];
end


ngp=length(spv);
ncp=ini_Val.ncontrol;
spvm=spv'*ones(1,p+1)+ones(ngp,1)*(0:p);
BNshape=zeros(ngp,ncp);
for g=1:ngp
    BNshape(g,spvm(g,:))=Nshape(g,:);
end
Invariant.BNshape = BNshape;


Invariant.Nshape = Nshape;
Invariant.dNsh = dNsh;
Invariant.ddNsh = ddNsh;
if ini_Val.d == 3
    Invariant.dddNsh = dddNsh;
end
Invariant.spv = spv;
Invariant.h = h;
Invariant.gweight = gweight;
Invariant.connect=connect;
Invariant.gaussp=gaussp;



