function AssembleCorrelationOperator(iscale,nmod,dflag)
if nargin<3, dflag=0;end
load(fullfile('TMP','params'));
onflight=0;
if isfield(param,'onflight')
onflight=param.onflight;
end
if onflight
global phidf wdetJ
end
param0=param;
load(fullfile('TMP',sprintf('%d_params',nmod)),'param');
tic();
if iscell(param0.reference_image)
ncams=length(param0.reference_image);
else
ncams=1;
end
load(fullfile('TMP',sprintf('sample0_%d',iscale-1)),'sizeim');
lc=(min(sizeim)/10);
lm=lc;
if isfield(param0,'regularization_parameter')
lm=param0.regularization_parameter/2^(iscale-1);
if isfield(param0,'detect')&&iscale>1
if param0.detect
lm=0.1*lm;
end
end
end
reg_type='none';
if isfield(param0,'regularization_type')
if strcmp(param0.regularization_type,'none')
reg_type=param0.regularization_type;
else
if (iscale==1)
reg_type=param0.regularization_type;
else
if strcmp(param0.regularization_type,'median')
reg_type='median';
else
reg_type='tiko';
end
end
end
end
if (strcmp(param.basis,'nurbs')||strcmp(param.basis,'btri'))&&(iscale>1)
param.basis='fem';
end
disp(sprintf(' Correlation operator...'));
for ncam=1:ncams
load(fullfile('TMP',sprintf('%d_mask_%d',nmod,iscale-1)),'mask');
if ~onflight
load(fullfile('TMP',sprintf('%d_phidf_%d',nmod*10^(ncam-1),iscale-1)),'phidf');
load(fullfile('TMP',sprintf('%d_phix_%d',nmod*10^(ncam-1),iscale-1)),'wdetJ');
end
M=phidf'*mask*wdetJ*phidf;
nM=norm(diag(M));
R=sparse(size(M,1),size(M,2));
P=sparse(size(M,1),size(M,2));
if strcmp(param.basis,'fem')&&~isfield(param,'gp_file')
load(fullfile('TMP',sprintf('%d_mesh_%d',nmod*10^(ncam-1),iscale-1)),'xo','yo','selected','Nnodes');
load(fullfile('TMP',sprintf('sample0_%d',iscale-1)),'sizeim');
ki=1/lc;
V=[repmat(cos(2*pi*(xo)*ki).*cos(2*pi*(yo)*ki),2,1);zeros(size(M,1)-2*prod(Nnodes),1)];
if mean(Nnodes(1:2))<2
reg_type='none';
end
switch reg_type
case 'tiko'
% if iscale==1
% keyboard
% end
%if iscale>1
if ~((iscale==1)&&isfield(param,'coupling_parameter'))
Rall=AssembleRegularizationOperator([],iscale,nmod,2);
a=(V'*M*V)/(V'*Rall*V);
a=a*(2*lm/lc)^2;
R=R+a*Rall;
disp(sprintf(' Tikhonov regularization...'));
else
if dflag
Nddl_tot=size(M,1);
save(fullfile('TMP',sprintf('%d_3d_epsxx_%d',nmod,iscale-1)),'Nddl_tot');
Ri=AssembleRegularizationOperator3D(nmod,1);
else
selected_nodes=find(~selected(:));
Ri=AssembleRegularizationOperator(selected_nodes,iscale,nmod,2);
disp(sprintf(' Tikhonov boundary regularization...'));
end
a=(V'*M*V)/(V'*Ri*V);
a=0.1*a*(2*lm/lc)^2;
R=R+a*Ri;
end
load(fullfile('TMP',sprintf('%d_mask_%d',nmod,iscale-1)),'cut_nodes');
% if numel(cut_nodes)
% Ri=AssembleRegularizationForMaskedElements(cut_nodes,iscale,nmod);
% R=R+10*nM*Ri/norm(diag(Ri));
% end
if (iscale==1)&&isfield(param,'enrichment')&&(~isfield(param,'coupling_parameter'))
load(fullfile('TMP',sprintf('%d_mesh_%d',nmod*10^(ncam-1),iscale-1)),'Nnodes');
CreateGradBasisFunction(iscale,nmod*10^(ncam-1));
load(fullfile('TMP',sprintf('%d_xepsxx_%d',nmod,iscale-1)),'epsxx','wdetJ');
load(fullfile('TMP',sprintf('%d_xepsyy_%d',nmod,iscale-1)),'epsyy');
load(fullfile('TMP',sprintf('%d_xepsxy_%d',nmod,iscale-1)),'Uxy','Uyx');
load(fullfile('TMP',sprintf('%d_mask_%d',nmod,iscale-1)),'masks');
Nddl_tot=size(epsxx,2);
load(fullfile('TMP',sprintf('%d_mask_%d',nmod,iscale-1)),'unmasked_nodes');
if isempty(unmasked_nodes)
Nddls=2*prod(Nnodes);
else
Nddls=2*numel(unmasked_nodes);
end
Ri=epsxx'*masks*wdetJ*epsxx+epsyy'*masks*wdetJ*epsyy+Uxy'*masks*wdetJ*Uxy+Uyx'*masks*wdetJ*Uyx;
[indi,indj,val]=find(Ri);
keep=~((indi<=Nddls)&(indj<=Nddls));
Ri=sparse(indi,indj,val.*keep,Nddl_tot,Nddl_tot);
% R=R+Po*Ri*Po;
R=R+a*Ri;
% if ~iscell(param0.levelset_file)
% nbfis=1;
% else
% nbfis=numel(param0.levelset_file);
% end
% for ic=1:nbfis
% load(fullfile('TMP',sprintf('%d_enrichment_%d',nmod,ic)),'face_nodes');
% Ri=AssembleRegularizationForMaskedElements(face_nodes,iscale,nmod);
% size(Ri)
% display(sprintf('Assembling penlaty terms for enriched nodes...'));
% %ind=[face_nodes,face_nodes+prod(Nnodes)];
% %Ri=sparse(ind,ind,1,size(M,1),size(M,2));
%
%
%
%
%
% % R=R+1.e1*nM*Ri/norm(diag(Ri));
% end
end
case {'equilibrium_gap','constitutive_gap'}
if (param0.detect==0)
selected_nodes=find(~selected(:));
Ri=AssembleRegularizationOperator(selected_nodes,iscale,nmod,2);
disp(sprintf(' Tikhonov boundary regularization...'));
% if (norm(diag(Ri)))>0
% R=R+1*nM*Ri/norm(diag(Ri));
% end
S=diag(sparse([repmat(selected(:),2,1);ones(size(M,1)-2*prod(Nnodes),1)]));
a=((S*V)'*M*(S*V))/(V'*Ri*V);
a=a*(0.1*2*lm/lc)^2;
R=R+a*Ri;
end
if ~isfield(param,'coupling_parameter')
wdetJo=wdetJ;
CreateGradBasisFunction(iscale,nmod*10^(ncam-1));
load(fullfile('TMP',sprintf('%d_mask_%d',nmod,iscale-1)),'unmasked_nodes');
load(fullfile('TMP',sprintf('%d_mesh_%d',nmod*10^(ncam-1),iscale-1)),'Nnodes','selected','conn');
load(fullfile('TMP',sprintf('%d_epsxx_%d',nmod,10*(iscale-1))),'epsxx','wdetJ');
load(fullfile('TMP',sprintf('%d_epsyy_%d',nmod,10*(iscale-1))),'epsyy');
load(fullfile('TMP',sprintf('%d_epsxy_%d',nmod,10*(iscale-1))),'epsxy');
load(fullfile('TMP',sprintf('%d_mask_%d',nmod,iscale-1)),'maskg');
% mesh_file=fullfile('TMP',sprintf('%d_mesh_%d',nmod,iscale-1));
% [phig]=CreateFiniteElementBasis(mesh_file,1,1,[],'Gauss_points');
% outer_nodes=~selected(:);
% fouter_nodes=AddOneNodeLayer(conn,find(outer_nodes));
% outer_nodes(fouter_nodes)=1;
% maskg=diag(sparse((phig*(1-outer_nodes))>0.1));
if isfield(param,'material_model')
[sxx,syy,sxy]=ComputeTangentStress(epsxx,epsyy,epsxy,nmod);
else
nuo=0.3;
mu=0.5/(1+nuo);
lambda = nuo/((1.-nuo^2));%CP
ltr=lambda*(epsxx+epsyy);
sxx=2*mu*epsxx+ltr;
syy=2*mu*epsyy+ltr;
sxy=mu*epsxy;
end
Ri= epsxx'*maskg*(wdetJ*sxx)+epsyy'*maskg*(wdetJ*syy)+epsxy'*maskg*(wdetJ*sxy);
if isfield(param,'enrichment')&&(iscale==1)
load(fullfile('TMP',sprintf('%d_xepsxx_%d',nmod,iscale-1)),'epsxx','wdetJ');
load(fullfile('TMP',sprintf('%d_xepsyy_%d',nmod,iscale-1)),'epsyy');
load(fullfile('TMP',sprintf('%d_xepsxy_%d',nmod,iscale-1)),'epsxy');
load(fullfile('TMP',sprintf('%d_mask_%d',nmod,iscale-1)),'masks');
Nddl_tot=size(epsxx,2);
% phis=CreateFiniteElementBasis(mesh_file,1,1,enriched_nodes,'sub_cells',true);
% masks=diag(sparse((phis*(1-outer_nodes))>0.0001));
if isfield(param,'material_model')
[sxx,syy,sxy]=ComputeTangentStress(epsxx,epsyy,epsxy,nmod);
else
ltr=lambda*(epsxx+epsyy);
sxx=2*mu*epsxx+ltr;
syy=2*mu*epsyy+ltr;
sxy=mu*epsxy;
end
Rx= epsxx'*masks*wdetJ*sxx+epsyy'*masks*wdetJ*syy+epsxy'*masks*wdetJ*sxy;
[indi,indj,val]=find(Rx);
Nddls=size(Ri,1);
keep=~((indi<=Nddls)&(indj<=Nddls));
Rx=sparse(indi,indj,val.*keep,Nddl_tot,Nddl_tot);
if Nddl_tot>size(Ri,2)
Ro=sparse(Nddl_tot-size(Ri,1),Nddl_tot-size(Ri,2));
Ri=blkdiag(Ri,Ro);
end
Ri=Ri+Rx;
end
select=ones(size(Ri,1),1);
outer_nodes=~selected;
Nnods=prod(Nnodes);
if ~isempty(unmasked_nodes)
outer_nodes=outer_nodes(unmasked_nodes);
Nnods=length(unmasked_nodes);
end
reject=find(outer_nodes(:));
select(reject)=0;
select(prod(Nnodes)+reject)=0;
if isfield(param,'enrichment')&&(iscale==1)
if ~iscell(param0.levelset_file)
nbfis=1;
else
nbfis=numel(param0.levelset_file);
end
for ic=1:nbfis
load(fullfile('TMP',sprintf('%d_enrichment_%d',nmod,ic)),'face_nodes');
select(face_nodes)=1;
select(Nnods+face_nodes)=1;
end
end
select=diag(sparse(select));
switch reg_type
case 'equilibrium_gap'
disp(sprintf(' Equilibrium Gap regularization...'));
Ri=select*Ri;
Ri=Ri'*Ri;
Pi=sparse(size(M,1),size(M,2));
a=(V'*M*V)/(V'*Ri*V);
a=a*(2*lm/lc)^4;
case 'constitutive_gap'
disp(sprintf(' Constitutive Equation Gap regularization...'));
Pi=1e9*diag(diag(Ri).*(~diag(select)));
a=(V'*M*V)/(V'*Ri*V);
a=a*(2*lm/lc)^2;
end
R=R+a*Ri;
P=P+a*Pi;
% R=R+Po*Ri*Po;
wdetJ=wdetJo;
end
case 'none'
disp(sprintf(' No regularization...'));
case 'median'
disp(sprintf(' Median regularization...'));
otherwise
error('REGULARIZATION TYPE NOT CODED YET');
end
end
if strcmp(param.basis,'vic-nurbs')
if param.Nelems>1
load(fullfile('TMP',sprintf('sample0_%d',1-1)),'ls1');
load(fullfile('TMP',sprintf('sample0_%d',iscale-1)),'sizeim','nband','on');
ls1=ls1(nband(on));
load(fullfile('TMP',sprintf('%d_phi_%d',nmod,iscale-1)),'dphii','phii');
load(fullfile('TMP',sprintf('%d_vicmesh_%d',nmod,iscale-1)),'mesh_size');
if lm==lc
lm=mean(mesh_size);
end
Ri=dphii'*dphii;
ki=1/lc;
vi=cos(2*pi*ls1*ki);
L=phii'*wdetJ*phii;
b=phii'*wdetJ*vi;
V=L\b;
a=(V'*M*V)/(V'*Ri*V);
a=a*(2*lm/lc)^2;
R=R+a*Ri;
% val=(1*max(diag(M))-mean(diag(M)))/mean(diag(Ri));
% % Po=diag(sparse(ones(size(M,1),1)*sqrt(val)));
% Po=sqrt(val);
% R=R+penalty*Po*Ri*Po;
disp(sprintf(' Tikhonov regularization...'));
end
end
if strcmp(param.basis,'nurbs')||strcmp(param.basis,'btri')
switch reg_type
case 'tiko'
disp(sprintf(' Tikhonov regularization...'));
load(fullfile('TMP',sprintf('%d_mesh_%d',nmod*10^(ncam-1),iscale-1)),'xo','yo','Nnodes');
load(fullfile('TMP',sprintf('sample0_%d',iscale-1)),'sizeim');
load(fullfile('TMP',sprintf('%d_phio_%d',nmod,iscale-1)),'phio');
ki=1/lc;
Vn=cos(2*pi*(xo)*ki).*cos(2*pi*(yo)*ki);
L=phio'*phio;
b=phio'*Vn;
V=L\b;
V=repmat(V,2,1);
CreateGradBasisFunction(iscale,nmod);
load(fullfile('TMP',sprintf('%d_epsxx_%d',nmod,10*(iscale-1))),'epsxx','wdetJ');
load(fullfile('TMP',sprintf('%d_epsyy_%d',nmod,10*(iscale-1))),'epsyy');
load(fullfile('TMP',sprintf('%d_epsxy_%d',nmod,10*(iscale-1))),'Uxy','Uyx');
Ri= epsxx'*wdetJ*epsxx+Uxy'*wdetJ*Uxy...
+epsyy'*wdetJ*epsyy+Uyx'*wdetJ*Uyx;
a=(V'*M*V)/(V'*Ri*V);
a=a*(2*lm/lc)^2;
R=R+a*Ri;
case 'none'
disp(sprintf(' No regularization...'));
otherwise
error('REGULARIZATION TYPE NOT CODED YET');
end
end
save(fullfile('TMP',sprintf('%d_operator_%d',nmod*10^(ncam-1),iscale-1)),'M','R','P');
end
disp(sprintf('Computing correlation operator for model %d...%6.2f s',nmod,toc()));
end
