非均匀材料的多尺度均匀化matlab程序

function CH = homogenize（lx ly lambda mu phi x）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lx        = Unit cell length in x-direction.
% ly        = Unit cell length in y-direction.
% lambda    = Lame‘s first parameter for both materials. Two entries.
% mu        = Lame‘s second parameter for both materials. Two entries.
% phi       = Angle between horizontal and vertical cell wall. Degrees
% x         = Material indicator matrix. Size used to determine nelx/nely
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INITIALIZE
% Deduce discretization
[nely nelx] = size（x）;
% Stiffness matrix consists of two parts one belonging to lambda and
% one belonging to mu. Same goes for load vector
dx = lx/nelx; dy = ly/nely;
nel = nelx*nely;
[keLambda keMu feLambda feMu] = elementMatVec（dx/2 dy/2 phi）;
% Node numbers and element degrees of freedom for full （not periodic） mesh
nodenrs = reshape（1:（1+nelx）*（1+nely）1+nely1+nelx）;
edofVec = reshape（2*nodenrs（1:end-11:end-1）+1nel1）;
edofMat = repmat（edofVec18）+repmat（[0 1 2*nely+[2 3 0 1] -2 -1]nel1）;
%% IMPOSE PERIODIC BOUNDARY CONDITIONS
% Use original edofMat to index into list with the periodic dofs
nn = （nelx+1）*（nely+1）; % Total number of nodes
nnP = （nelx）*（nely）;    % Total number of unique nodes
nnPArray = reshape（1:nnP nely nelx）;
% Extend with a mirror of the top border
nnPArray（end+1:） = nnPArray（1:）;
% Extend with a mirror of the left border
nnPArray（:end+1） = nnPArray（:1）;
% Make a vector into which we can index using edofMat:
dofVector = zeros（2*nn 1）;
dofVector（1:2:end） = 2*nnPArray（:）-1;
dofVector（2:2:end） = 2*nnPArray（:）;
edofMat = dofVector（edofMat）;
ndof = 2*nnP; % Number of dofs
%% ASSEMBLE STIFFNESS MATRIX
% Indexing vectors
iK = kron（edofMatones（81））‘;
jK = kron（edofMatones（18））‘;
% Material properties in the different elements
lambda = lambda（1）*（x==1） + lambda（2）*（x==2）;
mu     = mu（1）*（x==1） + mu（2）*（x==2）;
% The corresponding stiffness matrix entries
sK = keLambda（:）*lambda（:）.‘ + keMu（:）*mu（:）.‘;
K  = sparse（iK（:） jK（:） sK（:） ndof ndof）; 
%% LOAD VECTORS AND SOLUTION
% Assembly three load cases corresponding to the three strain cases
sF = feLambda（:）*lambda（:）.‘+feMu（:）*mu（:）.‘;
iF = repmat（edofMat‘31）;
jF = [ones（8nel）; 2*ones（8nel）; 3*ones（8nel）];
F  = sparse（iF（:） jF（:） sF（:） ndof 3）;
% Solve （remember to constrain one node）
chi（3:ndof:） = K（3:ndof3:ndof）\F（3:ndof:）;
%% HOMOGENIZATION
% The displacement vectors corresponding to