test_3d_resistor

 Create 3D model of a Rectangular resistor
 $Id: test_3d_resistor.m 6246 2022-03-23 13:45:26Z aadler $

0001 % Create 3D model of a Rectangular resistor
0002 % $Id: test_3d_resistor.m 6246 2022-03-23 13:45:26Z aadler $
0003 function test_3d_resistor
0004 
0005 ll=5*1; % length
0006 ww=1*2; % width
0007 hh=1*3; % height
0008 conduc= .13;  % conductivity in Ohm-meters
0009 current= 4;  % Amps
0010 z_contact= 1e-1;
0011 scale = .46;
0012 mdl= eidors_obj('fwd_model','3D rectangle');
0013 nn=0;
0014 for z=0:ll; for x=0:ww; for y=0:hh
0015    nn=nn+1;
0016    mdl.nodes(nn,:) = [x,y,z];
0017 end; end; end
0018 
0019 if 0
0020 % Matlab's delaunayn triangularization is screwed up
0021    mdl.elems = delaunayn(mdl.nodes);
0022 elseif 0
0023 % Matlab's delaunay3 triangularization is screwed up
0024    mdl.elems = delaunay3(mdl.nodes(:,1), mdl.nodes(:,2), mdl.nodes(:,3));
0025 elseif 0
0026    elem1= [ 1 2 3 5; 5 6 2 3; 5 6 7 3; ...
0027             4 2 3 8; 8 6 2 3; 8 6 7 3;];
0028    mdl.elems=[];
0029    for i=0:ll-1;
0030        mdl.elems= [mdl.elems; elem1+4*i];
0031    end
0032 else
0033    mdl= mk_grid_model([],0:ww,0:hh,0:ll);
0034    mdl.nodes= mdl.nodes*scale;
0035    mdl= rmfield(mdl,'coarse2fine');
0036 end
0037 mdl.boundary= find_boundary(mdl.elems);
0038 mdl.gnd_node = 1;
0039 elec_nodes= [1:(ww+1)*(hh+1)];
0040 elec(1).nodes= elec_nodes;      elec(1).z_contact= z_contact;
0041 elec(2).nodes= nn-elec_nodes+1; elec(2).z_contact= z_contact;
0042 stim.stim_pattern= [-1;1]*current;
0043 stim.meas_pattern= [-1,1];
0044 mdl.stimulation= stim;
0045 mdl.electrode= elec;
0046 show_fem(mdl);
0047 mdl = mdl_normalize(mdl,0);
0048 
0049 
0050 
0051 mdl.solve = @fwd_solve_1st_order;
0052 mdl.system_mat = @system_mat_1st_order;
0053 
0054 n_el = size(mdl.elems,1);
0055 img= eidors_obj('image','3D rectangle', ...
0056       'elem_data', ones(n_el,1) * conduc, ...
0057       'fwd_model', mdl); 
0058 
0059 fsol= fwd_solve(img);
0060 fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0061 
0062 % analytical solution
0063 Block_R =  ll / ww / hh / scale/ conduc;
0064 Contact_R = z_contact/(ww*hh)/scale^2;
0065 % Contact R reflects z_contact / (width/scale)^2. Here we need to use
0066 %  the scale, since this is not reflected in the size of the
0067 %  FEM as created by the grid. This is different to the test_2d_resistor,
0068 %  where the FEM is created scaled, so that the ww
0069 %  don't need to be scaled by the scale parameter.
0070 R = Block_R + 2*Contact_R;
0071 
0072 V= current*R;
0073 fprintf('Solver %s: %f\n', 'analytic', V);
0074 fprintf('Solver %s: %f\n', 'analytic (no z_contact)', V - 2*Contact_R*current);
0075 
0076 if 0 % OLD CODE
0077    mdl.solve = @np_fwd_solve;
0078    mdl.jacobian = @np_calc_jacobian;
0079 else
0080    mdl.solve = @fwd_solve_1st_order;
0081    mdl.jacobian= @jacobian_adjoint;
0082 end
0083 mdl.misc.perm_sym= '{n}';
0084 mdl.normalize_measurements = 0;
0085 img= eidors_obj('image','3D rectangle', ...
0086       'elem_data', ones(size(mdl.elems,1),1) * conduc, ...
0087       'fwd_model', mdl); 
0088 
0089 fsol= fwd_solve(img);
0090 fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0091 
0092 
0093 % NOW CALCULATE THE ANALYTICAL JACOBIAN
0094 if 0 % OLD CODE
0095    mdl.solve = @np_fwd_solve;
0096    mdl.jacobian = @np_calc_jacobian;
0097 else
0098    mdl.solve = @fwd_solve_1st_order;
0099    mdl.jacobian= @jacobian_adjoint;
0100 end
0101 img.fwd_model = mdl;
0102 img.elem_data= ones(n_el,1) * conduc ;
0103 Jnp= calc_jacobian(img);
0104 
0105 mdl.jacobian = @jacobian_perturb;
0106 Jp1= calc_jacobian(img);
0107 
0108 img.elem_data= ones(n_el,1) * conduc ;
0109 Jp2= zeros(size(Jnp));
0110 delta= 1e-8;
0111 for i=1:n_el
0112    fsol_h= fwd_solve(img);
0113    img.elem_data(i) = img.elem_data(i) + delta;
0114    fsol_i= fwd_solve(img);
0115    Jp2(:,i) = (fsol_i.meas - fsol_h.meas)/delta; 
0116 end
0117    
0118 mdl.solve = @fwd_solve_1st_order;
0119 mdl.jacobian = @jacobian_adjoint;
0120 Jaa= calc_jacobian(img);
0121 
0122 fprintf('Jacobians: Cols by Jaa, Jnp, Jp1, Jp2\n');
0123 JJ = [Jaa;Jnp;Jp1;Jp2];
0124 disp(JJ(:,1:6))
0125 
0126 % UNIT_TEST
0127 unit_test_cmp('Jaa vs Jnp', Jaa, Jnp, 1e-6);
0128 unit_test_cmp('Jaa vs Jp1', Jaa, Jp1, 1e-6);
0129 unit_test_cmp('Jaa vs Jp2', Jaa, Jp2, 1e-3);
0130 
0131 
0132 function cleaner_code
0133 % Cleaner code using newer EIDORS helper functions (AA - dec 2018)
0134 ll=5; ww=1; hh=1; z_contact= 1e-1; conduc= .13; current= 4; scale = .46;
0135 mdl= rmfield( mk_grid_model([],0:ww,0:hh,0:ll), 'coarse2fine');
0136 mdl.boundary= find_boundary(mdl.elems); mdl.gnd_node = 1;
0137 mdl.electrode(1)=struct('z_contact', z_contact, 'nodes', ...
0138      find(mdl.nodes(:,3)==0));
0139 mdl.electrode(2)=struct('z_contact', z_contact, 'nodes', ...
0140      find(mdl.nodes(:,3)==ll));
0141 mdl.stimulation= stim_meas_list([1,2,1,2],2,current);
0142 mdl.nodes = mdl.nodes * scale;
0143 show_fem(mdl);
0144 Block_R =  ll / ww / hh / scale/ conduc;
0145 Contact_R = z_contact/(ww*hh)/scale^2;
0146 Vs = current * ( Block_R + 2*Contact_R );
0147 vh = fwd_solve(mk_image(mdl, conduc));
0148 unit_test_cmp('3D Resistor',Vs, vh.meas, 1e-10);