test_3d_resistor

 Create 3D model of a Rectangular resistor
 $Id: test_3d_resistor.m 7100 2024-12-26 12:05:50Z aadler $

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