vector n = (20, 20, 20);
vertex a = (-1,-1,-1);
vertex b = ( 1, 1, 1);
double pi = 4*atan(1);
 
scene Sc = pov("example1.pov"); // the pov-ray file for the geometry

mesh M = structured(n,a,b);

domain O = domain(Sc,inside(<1,0,0>));

function KiO = one(<1,0,0>);
 
function uexact = (x*(x-1) + y*(y-1) + z*(z-1));
function vexact = sin(pi*x)*sin(pi*y)*sin(pi*z);

solve(u,v) in O by M method(type=penalty)
{
   pde(u)
     - dx(KiO*dx(u)) - dy(KiO*dy(u)) - dz(KiO*dz(u)) + KiO*v
       = -6*KiO + KiO*vexact;
     u = uexact on <1,0,0>;
   
   pde(v)
     KiO*v - div(KiO*grad(v)) + KiO*u
        = KiO*(3*pi^2 + 1)*vexact + KiO*uexact;
     v = vexact on <1,0,0>;
}

double I=int(M)(KiO*(uexact-u)^2);
double J=int(M)(KiO*uexact^2);
cout << sqrt(I/J) << "\n";
I = int(M)(KiO*(vexact-v)^2);
J = int(M)(KiO*vexact^2);
cout << sqrt(I/J) << "\n";

save(dx,"v.dat",KiO*v,M);
save(dx,"u.dat",KiO*u,M);