function [C,C1,C2,C3] = interpolate(Invariant,x)
global ini_Val 



p = ini_Val.p;
Nshape = Invariant.Nshape;
dNsh = Invariant.dNsh;
ddNsh = Invariant.ddNsh;
spv = Invariant.spv;
ngp_total=length(spv);

C = zeros(ngp_total,2);
C1 = zeros(ngp_total,2);
C2 = zeros(ngp_total,2);
if nargout==4
    C3 = zeros(ngp_total,2);
    dddNsh = Invariant.dddNsh;
else
    C3 = [];
end
for k=1:(p+1)
    C(:,1) = C(:,1)+ Nshape(:,k).*x(2*(spv(:)+k-1)-1);
    C(:,2) = C(:,2)+ Nshape(:,k).*x(2*(spv(:)+k-1));
    C1(:,1) = C1(:,1)+ dNsh(:,k).*x(2*(spv(:)+k-1)-1);
    C1(:,2) = C1(:,2)+ dNsh(:,k).*x(2*(spv(:)+k-1));
    C2(:,1) = C2(:,1)+ ddNsh(:,k).*x(2*(spv(:)+k-1)-1);
    C2(:,2) = C2(:,2)+ ddNsh(:,k).*x(2*(spv(:)+k-1));
    if nargout==4
        C3(:,1) = C3(:,1)+ dddNsh(:,k).*x(2*(spv(:)+k-1)-1);
        C3(:,2) = C3(:,2)+ dddNsh(:,k).*x(2*(spv(:)+k-1));
    end

end
