load "exactpartition"
mesh Th=square(10,10);
mesh3 Th3=buildlayers(Th,10); 
fespace Vh(Th,P1);
fespace Vh3(Th3,P1);
int n=10;
Vh[int] p(n); 
Vh3[int] p3(n); 
real[int][int] pp(n);

for(int i=0; i<n; ++i)
p[i]= lrint(x*y*n/3)== lrint(i*1./3.);
for(int i=0; i<n; ++i)
  p3[i]= lrint(x*y*(1-z)*n/3)== lrint(i*1./3.);
for(int i=0;i<n;++i)
{
    pp[i].resize(Vh.ndof);
    pp[i]=p[i][]; 
}
plot(p, wait=1);
int[int] pj=[3,4,1,8,9,100,8,7,10,11,12]; 
exactpartition(p,pj); 
exactpartition(p3,pj); 
exactpartition(pp,pj); 
Vh s=1;
for(int i=0; i<n;++i)
 s[] -= p[i][];
 
 assert(s[].linfty ==0); 
