femwind/ndt_boundary_conditions_fortran.m

   1 function Kb=ndt_boundary_conditions_fortran(K,params);
   2 % Kb=ndt_boundary_conditions_fortran(K,params);
   3 % call fortran version and compare results
   4 % in:
   5 %     K global stiffness matrix ndt14 format
   6
   7 [n1,n2,n3,st] = size(K);
   8 if st~=14,
   9     error('K must be in ndt 14 storage scheme')
  10 end
  11
  12 %Writing to file for use by fortran tester
  13 write_array_nd(swap23(K),'K');
  14
  15 system('./fortran/ndt_boundary_conditions_test.exe');
  16
  17 Kb_f = swap23(read_array_nd('Kb')); % boundary condition from fortran
  18
  19 % the same in matlab
  20 Ks = ndt_convert(K,'sparse');
  21 [Kb_ms,~]=apply_boundary_conditions(Ks,[],{K(:,:,:,1)});
  22
  23 % convert matlab output to nd14 and compare
  24 Kb_m =  ndt_convert(Kb_ms,n1,n2,n3,st);
  25 err = big(Kb_m - Kb_f)
  26
  27 % convert fortran output to sparse and compare
  28 Kb_fs = ndt_convert(Kb_f,'sparse');
  29 err_s = big(Kb_ms - Kb_fs)
  30 if err_s > eps(single(1.0))*big(Kb_fs)*10
  31     error('error too large')
  32 end
  33 end
  34