bld_master_full.sci

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

function [Ef]=bld_master_full(vtx,%simp,mat,elec,zc)
Ef=[];
//function [Ef] = bld_master_full(vtx,simp,mat,elec,zc);
// 
//System matrix assembling based on the complete electrode model. 
//This function is called within fem_master_full.
// 
// 
// 
//Ef   = The UNreferenced system matrix.
//vtx  = The vertices matrix. The coordinates of the nodes in 3D.
//simp = The simplices matrix. Unstructured tetrahedral.
//mat  = As for MATerial information. The conductivity vector.(isotropic)
//elec = The matrix that holds the boundary faces assigned as electrodes. Its typical
//       dimensions are [number of electrodes x 3*number of faces per electrode].
//zc   = The array of electrode contact impedances. 
 
 
[vr,vc] = size(vtx);
[sr,sc] = size(%simp);
[er,ec] = size(elec);
 
 
 
//! unknown arg type, using mtlb_length 
if mtlb_length(mat)~=sr then
  error('Invalid conductivity information for this mesh');
end
 
 
 
//! unknown arg type, using mtlb_length 
if mtlb_length(zc)==er then
   
   
  //The column vector zc with the contact impedances in [Ohms] is required
   
   
  Ef = bld_master(vtx,%simp,mat);
   
  vertos = sparse(zeros(vr,er));
  horos = sparse(zeros(er,vr+er));
   
  Ef = [[Ef,vertos];horos];
   
   
  //Up to this point we have calculated the master matrix without the influence of contact impedance.
   
  while %t then
   
  //! unknown arg type, using mtlb_length 
  if mtlb_length(zc)~=er then break;end
    disp(mtlb_sprintf('Please enter the correct zc column vector with length: %d',er));
  end
   
   
  for q = 1:er
     
    tang_area = 0;
     
    q_th_ele = elec(q,:);
    // Select the row of nodes corresponding to the current electrode
     
     
    //! unknown arg type, using mtlb_length 
    for w = 1:3:mtlb_length(q_th_ele)
       
      m = q_th_ele(w);
      n = q_th_ele(w+1);
      l = q_th_ele(w+2);
       
       
      // This way m & n nodes belong to the edge tangential to the electrode and also at the same simplex.
       
      // We will now evaluate the distance """"tangential contact area"""" between m & n 
       
      xm = vtx(m,1);
      ym = vtx(m,2);
      // m node coords
      zm = vtx(m,3);
       
      xn = vtx(n,1);
      yn = vtx(n,2);
      // n node coords
      zn = vtx(n,3);
       
      xl = vtx(l,1);
      yl = vtx(l,2);
      // l node coords
      zl = vtx(l,3);
       
       
      p1 = [ym,zm,1;yn,zn,1;yl,zl,1];
       
      p2 = [zm,xm,1;zn,xn,1;zl,xl,1];
       
      p3 = [xm,ym,1;xn,yn,1;xl,yl,1];
       
      Are = 0.5*sqrt(det(p1)^2+det(p2)^2+det(p3)^2);
      // area mnl
       
      cali_area = 2*Are ./ zc(q);
      // coefficient for the area mnl
       
       
      tang_area = tang_area+cali_area;
       
      // Start modifying """"expanding"""" the E master matrix
       
      Ef(m,vr+q) = Ef(m,vr+q)-cali_area/6;
      // Kv -> Ec  -> Vertical bar
      Ef(n,vr+q) = Ef(n,vr+q)-cali_area/6;
      Ef(l,vr+q) = Ef(l,vr+q)-cali_area/6;
       
      Ef(vr+q,m) = Ef(vr+q,m)-cali_area/6;
      // Kv' -> Ec' -> Horizontal bar
      Ef(vr+q,n) = Ef(vr+q,n)-cali_area/6;
      Ef(vr+q,l) = Ef(vr+q,l)-cali_area/6;
       
      Ef(m,m) = Ef(m,m)+cali_area/12;
      // Kz -> E -> Main bar
      Ef(m,n) = Ef(m,n)+cali_area/24;
      Ef(m,l) = Ef(m,l)+cali_area/24;
       
      Ef(n,m) = Ef(n,m)+cali_area/24;
      Ef(n,n) = Ef(n,n)+cali_area/12;
      Ef(n,l) = Ef(n,l)+cali_area/24;
       
      Ef(l,m) = Ef(l,m)+cali_area/24;
      Ef(l,n) = Ef(l,n)+cali_area/24;
      Ef(l,l) = Ef(l,l)+cali_area/12;
       
       
    end
    // dealing with this electrode
     
    Ef(vr+q,vr+q) = Ef(vr+q,vr+q)+tang_area;
     
  end
  //for the whole set of electrodes
   
end
 
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// 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
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%