load "scotch"
load "metis"
border C(t=0, 2*pi){ x=cos(t); y=sin(t); }
mesh Th = buildmesh(C(1000)); // first mesh: a simple circle
int lpart = 10; // how many clusters must be created ?
fespace Ph(Th, P0); // using P0 FE to store the partition
int[int] nupart(Th.nt); // the actual partition
Ph part = nuTriangle; // number of each triangles in the mesh
int[int] weight(Th.nt); // array of integers to store weights
for(int i = 0; i < Th.nt;++i)
    weight[i] = Th[i][1].x > 0. && Th[i][1].y > 0. ? 10 : 1; // if we want to partition has much as possible the region x > 0 && y > 0
plot(Th, part, value = 1, fill = 1, wait = 1, cmm = "original");
scotch(nupart, Th, lpart, weight = weight); // with weight
Ph balance;
for(int i = 0; i < Th.nt; ++i)
    balance[][i] = nupart[i];
plot(Th, balance, value = 1, fill = 1, wait = 1, cmm = "SCOTCH with weights");
scotch(nupart, Th, lpart); // without weight
for(int i = 0; i < Th.nt; ++i)
    balance[][i] = nupart[i];
plot(Th, balance, value = 1, fill = 1, wait = 1, cmm = "SCOTCH");
metisdual(nupart, Th, lpart);
for(int i = 0; i < Th.nt; ++i)
    balance[][i] = nupart[i];
plot(Th, balance, value = 1, fill = 1, wait = 1, cmm = "Metis");
Th = square(400, 400); // second mesh: a simple square
balance = 0;
part = nuTriangle;
weight.resize(Th.nt);
nupart.resize(Th.nt);
for(int i = 0; i < Th.nt;++i)
    weight[i] = Th[i][1].x > 0.5 && Th[i][1].y > 0.5 ? 100 : 1; // if we want to partition has much as possible the region x > 0.5 && y > 0.5
plot(Th, part, value = 1, fill = 1, wait = 1, cmm = "original");
scotch(nupart, Th, lpart, weight = weight);
for(int i = 0; i < Th.nt; ++i)
    balance[][i] = nupart[i];
plot(Th, balance, value = 1, fill = 1, wait = 1, cmm = "SCOTCH with weights");
scotch(nupart, Th, lpart);
for(int i = 0; i < Th.nt; ++i)
    balance[][i] = nupart[i];
plot(Th, balance, value = 1, fill = 1, wait = 1, cmm = "SCOTCH");
metisdual(nupart, Th, lpart);
for(int i = 0; i < Th.nt; ++i)
    balance[][i] = nupart[i];
plot(Th, balance, value = 1, fill = 1, wait = 1, cmm = "Metis");
