function E=energy1nl(vdoft)
global L B pp nel
global npe nodg ndof conecg conecdof coordg rr
vdof=vdoft';
xh=-1:0.1:1;#integration points
wh=[1/20 1/10*ones(1,19) 1/20];#integr. weights
nhat=size(xh,2);
##[g dg d2g] = master(1,xh);
[g dg d2g]=masterlag(1,xh);
[sf dsf ddsf]=masterlag(pp,xh);
E=0;
for k=1:nel
 XX(:,1)=coordg(:,conecg(k,1));
 XX(:,2)=coordg(:,conecg(k,2));
 x=zeros(1,nhat);
 dx=zeros(1,nhat);
 for i=1:2
  x(1,:)=x(1,:)+XX(1,i)*g(i,:);
  dx=dx+XX(1,i)*dg(i,:);
 end
 uh=vdof(conecdof(k,:))*sf;
 duh=vdof(conecdof(k,:))*dsf./dx;
####data
 beta=ones(1,nhat)+0*(1+sign(x-2*L/3)).*ones(1,nhat);
 gamma=0*ones(1,nhat);
 ff=ones(1,nhat);
####
 E=E+0.5*((abs(duh).^(rr-2)).*duh)*(duh.*beta.*wh.*dx)'+0.5*uh*(uh.*gamma.*wh.*dx)'-ff*(uh.*wh.*dx)';
end
E=E-B*vdoft(ndof);
end