function [J,g] = costFunc(p)

% J = costFunc(p)
%        calculates the following cost functional:
%        J = (1/meas(Omega_p)) * \int_{Omega_p} |curl(u)-Breq|^2 dx

sprintf('%s','costFunc')

global xinit e ep ey ec ID IN Ishape Ip Iy Ic IDD IDD1 IDD2 shapeDisplacement gridDisplacement A u B J g p_last h_last costSolved;


xl = -0.2;
xt = 0.1;
xpl = -0.02;
xpt = 0.02;
xyl = -0.16;
xyr = -0.04;
xyt = 0.04;
xcl = -0.12;
xcr = -0.08;
xct = 0.08;

h = 0.01;


if max(size(p_last)~=size(p))==1 | size(h_last)~=size(h)
    costSolved = [];
elseif max(p_last~=p) | h_last~=h
    costSolved = [];
end


if ~in(costSolved,'e')
    [xinit,e,ep,ey,ec,ID,IN,Ishape,Ip,Iy,Ic,IDD,IDD1,IDD2] = discret(xl,xt,xpl,xpt,xyl,xyr,xyt,xcl,xcr,xct,h);
    costSolved = [costSolved 'e'];
end

if ~in(costSolved,'s')
    shapeDisplacement = designMap(Ishape,p);
    costSolved = [costSolved 's'];
end

if ~in(costSolved,'g')
    gridDisplacement = shapeMap(xinit,e,IDD,IDD1,IDD2,Ishape,shapeDisplacement);
    costSolved = [costSolved 'g'];
end

if ~in(costSolved,'A')
    [A,f] = assemble(ey,ec,xinit+gridDisplacement,e,ID,IN);
    costSolved = [costSolved 'A'];
end

if ~in(costSolved,'u')
    u = A\f;
    costSolved = [costSolved 'u'];
end

if ~in(costSolved,'B')
    B = curl(xinit+gridDisplacement,e,u);
    costSolved = [costSolved 'B'];
end

if ~in(costSolved,'J')
    Breq = [0.025;0];
    J = objectiveFunc(xpl,xpt,xinit+gridDisplacement,e,ep,B,Breq);
    costSolved = [costSolved 'J'];
end

g = -1;

p'
J

p_last = p;
h_last = h;
