function material_cell = mean_material_node(grid3d, gt, material_node)

chkarg(istypesizeof(grid3d, 'Grid3d'), '"grid3d" should be instance of Grid.');
chkarg(istypesizeof(gt, 'GT'), '"gt" should be instance of GT.');
chkarg(istypesizeof(material_node, 'complexcell', [1 Axis.count], grid3d.N), ...
	'"material_node" should be length-%d cell array, whose each element is %d-by-%d-by-%d array with complex elements.', ...
	Axis.count, grid3d.Ncell{:});

for w = Axis.elems
	material_node{w} = expand_node_array(grid3d, material_node{w});  % (Nx+2) x (Ny+2) x (Nz+2)
end

if gt == GT.prim
	% material parameters for fields on primary grid
	material_cell = arithmetic_mean_material_node(material_node);
else  % gt == GT.dual
	% material parameters for fields on dual grid
	material_cell = harmonic_mean_material_node(material_node);
end


function material_edge_cell = arithmetic_mean_material_node(material_node)
material_edge_cell = cell(1, Axis.count);

material_edge_cell{Axis.x} = (...
	material_node{Axis.x}(2:end-1, 2:end-1, 2:end-1) ...
	+ material_node{Axis.x}(2:end-1, 1:end-2, 2:end-1) ...
	+ material_node{Axis.x}(2:end-1, 2:end-1, 1:end-2) ...
	+ material_node{Axis.x}(2:end-1, 1:end-2, 1:end-2)...
) / 4;

material_edge_cell{Axis.y} = (...
	material_node{Axis.y}(2:end-1, 2:end-1, 2:end-1) ...
	+ material_node{Axis.y}(2:end-1, 2:end-1, 1:end-2) ...
	+ material_node{Axis.y}(1:end-2, 2:end-1, 2:end-1) ...
	+ material_node{Axis.y}(1:end-2, 2:end-1, 1:end-2)...
) / 4;

material_edge_cell{Axis.z} = (...
	material_node{Axis.z}(2:end-1, 2:end-1, 2:end-1) ...
	+ material_node{Axis.z}(1:end-2, 2:end-1, 2:end-1) ...
	+ material_node{Axis.z}(2:end-1, 1:end-2, 2:end-1) ...
	+ material_node{Axis.z}(1:end-2, 1:end-2, 2:end-1)...
) / 4;


function material_face_cell = harmonic_mean_material_node(material_node)
material_face_cell = cell(1, Axis.count);
material_face_cell{Axis.x} = 2./(1./material_node{Axis.x}(1:end-2, 2:end-1, 2:end-1) + 1./material_node{Axis.x}(2:end-1, 2:end-1, 2:end-1));
material_face_cell{Axis.y} = 2./(1./material_node{Axis.y}(2:end-1, 1:end-2, 2:end-1) + 1./material_node{Axis.y}(2:end-1, 2:end-1, 2:end-1));
material_face_cell{Axis.z} = 2./(1./material_node{Axis.z}(2:end-1, 2:end-1, 1:end-2) + 1./material_node{Axis.z}(2:end-1, 2:end-1, 2:end-1));
