% First building block
n = 4;
p1 = [0 0
    3 0
    0 1
    1.5 1
    0 2
    1.5 2
    0 3
    3 3];

f1 = [1 2 3 4
    3 4 5 6
    5 6 7 8
    2 8 4 6];

% second building block (square)
p2 = [0 0
    1 0
    0 1
    1 1];

f2 = [1 2 3 4];

% Build lauyers using the building blocks
% first "layer"
glist = {};
for i=1:3
    glist{end+1} = {p1+[10 9-3*i] f1};
end
for k=1:3
    for j=1:n
        for i=1:3
            glist{end+1} = {p2+[10-j 9-i-3*(k-1)] f2};
        end
    end
end
[p,f] = combineg(glist);

% second "layer"
pb = (p-[6 0])/3; % building block
for i=1:3
    glist{end+1} = {pb+[10-n-7/3 9-3*i] f};
end

% third layer
for i=1:3*3
    glist{end+1} = {pb/3+[10-n-7/3-7/9 9-i] f};
end
[p,f] = combineg(glist);

[pu,fu] = uniquevertices(p,f);
import com.comsol.model.util.*
ModelUtil.remove('MeshModel');
model = ModelUtil.create('MeshModel');
model.hist.disable
model.geom.create('geom1', 2);
model.mesh.create('mesh1', 'geom1');
model.mesh('mesh1').data().setVertex(pu');
model.mesh('mesh1').data().setElem('quad', fu'-1);
model.mesh('mesh1').data().createMesh;
mphmesh(model,'mesh1')

function [p,f]=combineg(glist)

p = [];
f = [];
for i=1:length(glist)
    idx = size(p,1);
    item = glist{i};
    p0 = item{1};
    f0 = item{2};
    p = [p; p0];
    f = [f; f0+idx];
end
end

function [pout,fout]=uniquevertices(pin,fin)

f = 50; % precision factor (repair setting)
[pout,~,idx2] = unique(round(pin*f), 'rows', 'stable');
pout = pout/f;
fout = reshape(idx2(fin(:)), [], 4);

end