bld_master.sci~

用来实现三维阻抗及光学断层成像重建的matlab程序

function [Ef]=bld_master(vtx,%simp,mat_ref)
Ef=[];
//function [Ef] = bld_master(vtx,simp,mat_ref);
// 
//Builds up the main compartment (GAP-SHUNT) of the system matrix 
//for the complete electrode model. It is called within the function 
//fem_master_full.
// 
// 
// 
//Ef      = The UNreferenced GAP based system matrix
//vtx     = The vertices matrix.
//simp    = The simplices matrix.
//mat_ref = The reference CONDUCTIVITY at each element.
 
 
[nv,dimen] = size(vtx);
//Number of vertices and dimension
[ns,dimen_p1] = size(%simp);
//Number of simplices
a = mat_ref;
dimen2 = 2*dimen;
 
ils = 1:dimen*ns;
ilst(1,1:2:dimen2*ns) = ils;
 
//! ilst = mtlb_i(ilst,2:2:dimen2*ns,ils) may be replaced by ilst(2:2:dimen2*ns) = ils
//!   if ilst is not a  vector
//!   or if ilst and ils are both row or column vectors
ilst = mtlb_i(ilst,2:2:dimen2*ns,ils)
 
 
patt = 1:dimen2:ns*dimen2;
 
jlst(1,patt) = %simp(:,1);
 
//! jlst = mtlb_i(jlst,patt+1,%simp(:,2)) may be replaced by jlst(patt+1) = %simp(:,2)
//!   if jlst is not a  vector
//!   or if jlst and %simp(:,2) are both row or column vectors
jlst = mtlb_i(jlst,patt+1,%simp(:,2))
 
//! jlst = mtlb_i(jlst,patt+2,%simp(:,1)) may be replaced by jlst(patt+2) = %simp(:,1)
//!   if jlst is not a  vector
//!   or if jlst and %simp(:,1) are both row or column vectors
jlst = mtlb_i(jlst,patt+2,%simp(:,1))
 
//! jlst = mtlb_i(jlst,patt+3,%simp(:,3)) may be replaced by jlst(patt+3) = %simp(:,3)
//!   if jlst is not a  vector
//!   or if jlst and %simp(:,3) are both row or column vectors
jlst = mtlb_i(jlst,patt+3,%simp(:,3))
 
//! jlst = mtlb_i(jlst,patt+4,%simp(:,1)) may be replaced by jlst(patt+4) = %simp(:,1)
//!   if jlst is not a  vector
//!   or if jlst and %simp(:,1) are both row or column vectors
jlst = mtlb_i(jlst,patt+4,%simp(:,1))
 
//! jlst = mtlb_i(jlst,patt+5,%simp(:,4)) may be replaced by jlst(patt+5) = %simp(:,4)
//!   if jlst is not a  vector
//!   or if jlst and %simp(:,4) are both row or column vectors
jlst = mtlb_i(jlst,patt+5,%simp(:,4))
 
 
sls = ones(1,dimen*ns);
slst(1,1:2:dimen2*ns) = -sls;
 
//! slst = mtlb_i(slst,2:2:dimen2*ns,sls) may be replaced by slst(2:2:dimen2*ns) = sls
//!   if slst is not a  vector
//!   or if slst and sls are both row or column vectors
slst = mtlb_i(slst,2:2:dimen2*ns,sls)
 
 
D0 = sparse([ilst(:),jlst(:)],slst,[dimen*ns,nv]);
//Correct up to here
 
 
J1 = D0*vtx(:,1);
J2 = D0*vtx(:,2);
J3 = D0*vtx(:,3);
 
for w = 1:3
  r = 1;
  JJ(r,w) = J1(w:dimen:ns*dimen);
  JJ(r+1,w) = J2(w:dimen:ns*dimen);
  JJ(r+2,w) = J3(w:dimen:ns*dimen);
end
 
 
//!! Unknown function squeeze ,the original calling sequence is used
dj = 1;
 
 
ilst = (1:dimen*ns).*.ones(1,dimen);
jlst = zeros(1,ns*(dimen^2));
for d = 1:dimen
  jlst(d:dimen:ns*(dimen^2)) = (d:dimen:dimen*ns).*.ones(1,dimen)
end
slst = zeros(1,ns*(dimen^2));
 
pat = 1:dimen^2:ns*(dimen^2);
 
 
//!! Unknown function squeeze ,the original calling sequence is used
%v = squeeze(sum([-JJ(3,:);prod([3;JJ(2,:)])])) ./ dj
slst(1,2) = %v(:).';
slst = pat;
 
//!! Unknown function squeeze ,the original calling sequence is used
%v = squeeze(sum([-JJ(1,:);prod([3;JJ(3,:)])])) ./ dj
slst(1,2) = %v(:).';
slst = pat+1;
 
//!! Unknown function squeeze ,the original calling sequence is used
%v = squeeze(sum([-JJ(1,:);prod([2;JJ(2,:)])])) ./ dj
slst(1,3) = %v(:).';
slst = pat+2;
 
//!! Unknown function squeeze ,the original calling sequence is used
%v = squeeze(sum([-JJ(2,:);prod([3;JJ(3,:)])])) ./ dj
slst(1,1) = %v(:).';
slst = pat+3;
 
//!! Unknown function squeeze ,the original calling sequence is used
%v = squeeze(sum([-JJ(3,:);prod([3;JJ(1,:)])])) ./ dj
slst(1,1) = %v(:).';
slst = pat+4;
 
//!! Unknown function squeeze ,the original calling sequence is used
%v = squeeze(sum([-JJ(1,:);prod([3;JJ(2,:)])])) ./ dj
slst(1,1) = %v(:).';
slst = pat+5;
 
//!! Unknown function squeeze ,the original calling sequence is used
%v = squeeze(sum([-JJ(2,:);prod([1;JJ(3,:)])])) ./ dj
slst(1,2) = %v(:).';
slst = pat+6;
 
 
//!! Unknown function squeeze ,the original calling sequence is used
%v = squeeze(sum([-JJ(1,:);prod([2;JJ(3,:)])])) ./ dj
slst(1,1) = %v(:).';
slst = pat+7;
 
//!! Unknown function squeeze ,the original calling sequence is used
%v = squeeze(sum([-JJ(1,:);prod([1;JJ(2,:)])])) ./ dj
slst(1,2) = %v(:).';
slst = pat+8;
 
 
LocJac = sparse([ilst(:),jlst(:)],slst,[dimen*ns,dimen*ns]);
 
A = LocJac*D0;
 
 
Vols = abs(dj)/prod(1:dimen);
 
 
//! unknown arg type, using mtlb_length 
materials = mtlb_length(a);
volumes = size(Vols);
 
if materials~=volumes then
  error('Some elements have not been assigned');
end
 
%v$1 = 1:dimen*ns;%v$1=%v$1(:)
%v$2 = 1:dimen*ns;%v$2=%v$2(:)
Ela = sparse([%v$1,%v$2],(a .* Vols)'.*.ones(1,dimen));
 
Ef = A'*Ela*A;
 
 
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// This is part of the EIDORS suite.
// Copyright (c) N. Polydorides 2001
// Copying permitted under terms of GNU GPL
// See enclosed file gpl.html for details.
// EIDORS 3D version 1.0
// MATLAB version 5.3 R11
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%