From 661410cd668a8a966aa55e78e3b1e37704e8545d Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Mon, 20 Dec 2010 10:40:33 +0100 Subject: [PATCH] ProcessFields now writes TD data into HDF5 group: /FieldData/TD + adapted matlab HDF5 Signed-off-by: Thorsten Liebig --- Common/processfields.cpp | 4 +++- matlab/GetField_Interpolation.m | 16 +++++++++---- matlab/GetField_TD2FD.m | 21 +++++++++++------ matlab/PlotHDF5FieldData.m | 16 ++++++------- matlab/ReadHDF5FieldData.m | 42 ++++++++++++++++++++++++++++----- 5 files changed, 73 insertions(+), 26 deletions(-) diff --git a/Common/processfields.cpp b/Common/processfields.cpp index f9f25ff..761a5da 100644 --- a/Common/processfields.cpp +++ b/Common/processfields.cpp @@ -87,6 +87,8 @@ void ProcessFields::InitProcess() delete group; group = new H5::Group( file->createGroup( "/FieldData/FD" )); delete group; + group = new H5::Group( file->createGroup( "/FieldData/TD" )); + delete group; delete file; } } @@ -423,7 +425,7 @@ bool ProcessFields::DumpVectorArray2HDF5(string filename, string name, FDTD_FLOA H5::H5File file( FILE_NAME, H5F_ACC_RDWR ); - H5::Group group( file.openGroup( "/FieldData" )); + H5::Group group( file.openGroup( "/FieldData/TD" )); hsize_t dimsf[4]; // dataset dimensions diff --git a/matlab/GetField_Interpolation.m b/matlab/GetField_Interpolation.m index abc8490..13511a2 100644 --- a/matlab/GetField_Interpolation.m +++ b/matlab/GetField_Interpolation.m @@ -49,11 +49,19 @@ mesh_i.lines{1} = x_i; mesh_i.lines{2} = y_i; mesh_i.lines{3} = z_i; -NULL = zeros(numel(x_i),numel(y_i),numel(z_i),3); -for n=1:numel(field.values) - field_i.values{n} = NULL; +if (isfield(field,'TD')) + field_i.TD = interpolate_fields(field.TD,x,y,z, x_i, y_i, z_i); + field_i.TD.time = field.TD.time; end -clear NULL; + +if (isfield(field,'FD')) + field_i.FD = interpolate_fields(field.FD,x,y,z, x_i, y_i, z_i); + field_i.FD.freq = field.FD.freq; +end + +return + +function field_i = interpolate_fields(field, x,y,z, x_i, y_i, z_i) % matlab cannot handle 3D data to be 2D data, workaround for these cases if (numel(x)==1) diff --git a/matlab/GetField_TD2FD.m b/matlab/GetField_TD2FD.m index 4cf8835..4bf6c91 100644 --- a/matlab/GetField_TD2FD.m +++ b/matlab/GetField_TD2FD.m @@ -1,5 +1,5 @@ -function field_FD = GetField_TD2FD(field_TD, freq) -% function field_FD = GetField_TD2FD(field_TD, freq) +function field = GetField_TD2FD(field, freq) +% function field = GetField_TD2FD(field, freq) % % Transforms time-domain field data into the frequency domain % Autocorrects the half-timestep offset of the H-field @@ -15,20 +15,27 @@ function field_FD = GetField_TD2FD(field_TD, freq) % % See also ReadHDF5FieldData -t = field_TD.time; +if (~isfield(field,'TD')) + warning('openEMS:GetField_TD2FD','field has no time domain data... skipping FD transformation...'); + return +end -field_FD.freq = freq; +t = field.TD.time; + +clear field.FD + +field.FD.freq = freq; for nf = 1:numel(freq) - field_FD.values{nf} = 0; + field.FD.values{nf} = 0; end -numTS = numel(field_TD.values); +numTS = numel(field.TD.values); for n=1:numTS for nf = 1:numel(freq) f = freq(nf); - field_FD.values{nf} = field_FD.values{nf} + 2/numTS * field_TD.values{n}.*exp(-1i*2*pi*f*t(n)); + field.FD.values{nf} = field.FD.values{nf} + 2/numTS * field.TD.values{n}.*exp(-1i*2*pi*f*t(n)); % t(n) is absolute time and therefore the half-timestep offset of % the H-field is automatically compensated % openEMS output: E-fields start at t=0 diff --git a/matlab/PlotHDF5FieldData.m b/matlab/PlotHDF5FieldData.m index 624d30d..553bee4 100644 --- a/matlab/PlotHDF5FieldData.m +++ b/matlab/PlotHDF5FieldData.m @@ -27,9 +27,9 @@ fields = ReadHDF5FieldData(file); if (mesh.type==0) % cartesian mesh [X Y Z] = meshgrid(mesh.lines{1},mesh.lines{2},mesh.lines{3}); - for n=1:numel(fields.values) + for n=1:numel(fields.TD.values) % since Matlab 7.1SP3 the field needs to be reordered - fields.values{n} = permute(fields.values{n},[2 1 3 4]); % reorder: y,x,z (or y,x) + fields.TD.values{n} = permute(fields.TD.values{n},[2 1 3 4]); % reorder: y,x,z (or y,x) end else disp(['PlotHDF5FieldData:: Error: unknown mesh type ' num2str(mesh.type)]); @@ -38,14 +38,14 @@ end max_amp = 0; if (component>0) - for n=1:numel(fields.values) - Field{n} = fields.values{n}(:,:,:,component); + for n=1:numel(fields.TD.values) + Field{n} = fields.TD.values{n}(:,:,:,component); end else - for n=1:numel(fields.values) - fx = fields.values{n}(:,:,:,1); - fy = fields.values{n}(:,:,:,2); - fz = fields.values{n}(:,:,:,3); + for n=1:numel(fields.TD.values) + fx = fields.TD.values{n}(:,:,:,1); + fy = fields.TD.values{n}(:,:,:,2); + fz = fields.TD.values{n}(:,:,:,3); Field{n} = sqrt(fx.^2 + fy.^2 + fz.^2); end end diff --git a/matlab/ReadHDF5FieldData.m b/matlab/ReadHDF5FieldData.m index 18bb038..08879c6 100644 --- a/matlab/ReadHDF5FieldData.m +++ b/matlab/ReadHDF5FieldData.m @@ -25,21 +25,51 @@ if isOctave end info = hdf5info(file); +TD.names = {}; +FD.names = {}; +hdf_fielddata = []; for n=1:numel(info.GroupHierarchy.Groups) if strcmp(info.GroupHierarchy.Groups(n).Name,'/FieldData') - for m=1:numel(info.GroupHierarchy.Groups(n).Datasets) - names{m} = info.GroupHierarchy.Groups(n).Datasets(m).Name; - hdf_fielddata.time(m) = double(info.GroupHierarchy.Groups(n).Datasets(m).Attributes.Value); + %found /FieldData, look for either TD or FD data + for nGroup=1:numel(info.GroupHierarchy.Groups(n).Groups) + %search and read TD data + if strcmp(info.GroupHierarchy.Groups(n).Groups(nGroup).Name,'/FieldData/TD') + for m=1:numel(info.GroupHierarchy.Groups(n).Groups(nGroup).Datasets) + TD.names{m} = info.GroupHierarchy.Groups(n).Groups(nGroup).Datasets(m).Name; + hdf_fielddata.TD.time(m) = double(info.GroupHierarchy.Groups(n).Groups(nGroup).Datasets(m).Attributes.Value); + end + end + %search and read FD data + if strcmp(info.GroupHierarchy.Groups(n).Groups(nGroup).Name,'/FieldData/FD') + for m=1:numel(info.GroupHierarchy.Groups(n).Groups(nGroup).Datasets) + FD.names{m} = info.GroupHierarchy.Groups(n).Groups(nGroup).Datasets(m).Name; + FD.freq(m) = double(info.GroupHierarchy.Groups(n).Groups(nGroup).Datasets(m).Attributes.Value); + end + end end + end end -hdf_fielddata.names = names; -for n=1:numel(names) - hdf_fielddata.values{n} = double(hdf5read(file,names{n})); +if (numel(TD.names)>0) + for n=1:numel(TD.names) + hdf_fielddata.TD.values{n} = double(hdf5read(file,TD.names{n})); + end end +if (numel(FD.names)>0) + Nr_freq = numel(FD.names)/2; + for n=1:Nr_freq + name = ['/FieldData/FD/f' int2str(n-1) '_real']; + ind = find(strcmp(FD.names,name)); + hdf_fielddata.FD.values{n} = double(hdf5read(file,FD.names{ind})); + hdf_fielddata.FD.freq(n) = FD.freq(ind); + name = ['/FieldData/FD/f' int2str(n-1) '_imag']; + ind = find(strcmp(FD.names,name)); + hdf_fielddata.FD.values{n} = hdf_fielddata.FD.values{n} + 1j*double(hdf5read(file,FD.names{ind})); + end +end function hdf_fielddata = ReadHDF5FieldData_octave(file)