include "cube.idp"
load "medit"
int[int, int] LL = [[2,2], [1,1], [1,1]];
real[int, int] BB = [[0,1], [0,1], [0,1]];
int[int] CC = [4, 4, 4];
load "splitmesh12"
mesh3 Th3=Cube(CC, BB, LL);
mesh3 Th12=splitmesh12(Th3);
mesh3 Th12WF=splitmesh12WorseyFarin(Th3);
medit("Th12WF",Th12WF);
assert(Th12.nt == Th3.nt * 12);
cout << Th12.nt << " == 12 x " << Th3.nt << endl;
assert(Th12.nbe == Th3.nbe * 3);
cout << Th12.nbe << " == 3 x " << Th3.nbe << endl;
plot(Th3, cmm="Th3", wait=1);
plot(Th12, cmm="Th12", wait=1);
fespace Vh(Th12, P1);
Vh xyz = x + y + z;
medit("splitmesh12",Th12,xyz);
