From 063f4a1780aab7dad5ec33f72863187669bf816e Mon Sep 17 00:00:00 2001 From: Sam Cunliffe Date: Tue, 13 Dec 2022 17:50:57 +0000 Subject: [PATCH] =?UTF-8?q?Rename=20phasorinc,=20add=20{C,D}Collection=20d?= =?UTF-8?q?ocumentation,=20and=20`main()`=20=E2=86=92=20`main.cpp`=20(#200?= =?UTF-8?q?)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Move main into its own file. * Typo and fix details tags. * Rename phasorinc. * Document the CCollection and DCollection classes. * Recursive glob everything in the source directory. Co-authored-by: Will Graham <32364977+willGraham01@users.noreply.github.com> --- .github/workflows/doxygen-gh-pages.yml | 2 +- README.md | 52 +++++++++++++++----- doc/developers.md | 6 ++- tdms/CMakeLists.txt | 35 +++---------- tdms/cmake/targets.cmake | 4 +- tdms/include/arrays.h | 32 ++++++++++-- tdms/include/mesh_base.h | 2 +- tdms/include/simulation_parameters.h | 68 +++++++++++++------------- tdms/src/argument_parser.cpp | 6 +-- tdms/src/iterator.cpp | 6 +-- tdms/src/main.cpp | 67 +++++++++++++++++++++++++ tdms/src/mesh_base.cpp | 4 +- tdms/src/openandorder.cpp | 63 +----------------------- tdms/src/simulation_parameters.cpp | 8 +-- 14 files changed, 199 insertions(+), 156 deletions(-) create mode 100644 tdms/src/main.cpp diff --git a/.github/workflows/doxygen-gh-pages.yml b/.github/workflows/doxygen-gh-pages.yml index 265181e7b..23a381e09 100644 --- a/.github/workflows/doxygen-gh-pages.yml +++ b/.github/workflows/doxygen-gh-pages.yml @@ -8,7 +8,7 @@ jobs: deploy: runs-on: ubuntu-latest steps: - - uses: DenverCoder1/doxygen-github-pages-action@v1.1.0 + - uses: DenverCoder1/doxygen-github-pages-action@v1.2.0 with: github_token: ${{ secrets.GITHUB_TOKEN }} folder: html diff --git a/README.md b/README.md index 6a6a6a365..657575d59 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ *** ## Introduction -TDMS (Time Domain Maxwell Solver) is a hybrid C++ and MATLAB for solving +TDMS, the Time Domain Maxwell Solver, is a hybrid C++ and MATLAB tool for solving Maxwell's equations to simulate light propagation through a medium. See the [pdf documentation](https://github.com/UCL/TDMS/blob/gh-doc/masterdoc.pdf) for further details. @@ -20,19 +20,19 @@ further details. *** ## Compilation -TDMS requires building against [FFTW](https://www.fftw.org/) and +TDMS needs to be built against [FFTW](https://www.fftw.org/) and [MATLAB](https://www.mathworks.com/products/matlab.html), thus both need to be downloaded and installed prior to compiling TDMS. Install with ```bash -cd tdms -mkdir build; cd build -cmake .. \ +$ cd tdms +$ mkdir build; cd build +$ cmake .. \ # -DMatlab_ROOT_DIR=/usr/local/MATLAB/R2019b/ \ # -DFFTW_ROOT=/usr/local/fftw3/ \ -# -DCMAKE_INSTALL_PREFIX=$HOME/.local/ +# -DCMAKE_INSTALL_PREFIX=$HOME/.local/ \ # -DBUILD_TESTING=ON -make install +$ make install ``` where lines need to be commented in and the paths modified if cmake cannot (1) find MATLAB, (2) find FFTW or (3) install to the default install prefix. @@ -62,16 +62,45 @@ where lines need to be commented in and the paths modified if cmake cannot *** ## Usage +Once the executable has been compiled and installed, `tdms` should be in the `PATH`. +Check that installation worked with + +```bash +$ tdms -h +``` + +You can invoke it directly or call it from a MATLAB script. +We recommend that beginners with MATLAB installed start with the demonstration MATLAB script. + #### To run the demonstration code -Once the executable has been compiled, move into directory _examples/arc_01_, -launch Matlab and run the Matlab script: +Move into directory [`examples/arc_01`](./examples/arc_01/), +launch MATLAB and run the MATLAB script: -_run_pstd_bscan.m_ +[`run_pstd_bscan.m`](./examples/arc_01/run_pstd_bscan.m) This script will generate the input to the executable, run the executable and display sample output. +#### To run standalone + +You can also run `tdms` from the command line... + +```bash +$ tdms --help +Usage: +tdms [options] infile outfile +tdms [options] infile gridfile outfile +Options: +-h: Display this help message +-fd, --finite-difference: Use the finite-difference solver, instead of the pseudo-spectral method. +-q: Quiet operation. Silence all logging +-m: Minimise output file size by not saving vertex and facet information +``` + +The basic workflow is with two arguments, an input file as specified by [`iterate_fdtd_matrix.m`](./tdms/matlab/iteratefdtd_matrix.m), and an output file name to be created. + +You can choose two possible solver methods: either pseudo-spectral time-domain (PSTD, the default) or finite-difference (FDTD, with option `--finite-difference`). #### Parallelism @@ -80,9 +109,10 @@ number of threads can be set with the `OMP_NUM_THREADS` environment variable. For example, to use 4 threads, in a bash shell, use: ```bash -export OMP_NUM_THREADS=4 +$ export OMP_NUM_THREADS=4 ``` +Before calling the `tdms` executable. ## Contributing diff --git a/doc/developers.md b/doc/developers.md index fbff434c8..8d5b201fd 100644 --- a/doc/developers.md +++ b/doc/developers.md @@ -99,7 +99,9 @@ spdlog::debug("Send help"); ``` ### Compiling on UCL's Myriad cluster +
+ > **Warning** > These instructions are a bit experimental. Please use with care (and report anything that's wrong here)! @@ -134,13 +136,13 @@ spdlog::debug("Send help"); CApath: none ``` it's because the MATLAB module is interfering with the SSL certificates (and we clone over https by default). This issue is known and reported. As a workaround, we've added the build option `-DGIT_SSH=ON` to switch to `git clone` over ssh instead. +
## Where's the main? -The C++ `main` function is in openandorder.cpp however this only really does file I/O and setup. -The main FDTD algorithm code is in iterator.cpp and classes included therein. +The C++ `main` function is in main.cpp however the main algorithm code is execute_simulation() in iterator.cpp ## Testing diff --git a/tdms/CMakeLists.txt b/tdms/CMakeLists.txt index e64b2ffe8..d5e8225a5 100644 --- a/tdms/CMakeLists.txt +++ b/tdms/CMakeLists.txt @@ -81,34 +81,15 @@ else() add_definitions(-DSPDLOG_ACTIVE_LEVEL=SPDLOG_LEVEL_INFO) endif() - # TDMS target ----------------------------------------------------------------- -set(SOURCES - src/fields/base.cpp - src/fields/electric.cpp - src/fields/magnetic.cpp - src/fields/split.cpp - src/fields/td_field_exporter_2d.cpp - src/argument_parser.cpp - src/array_init.cpp - src/arrays.cpp - src/dimensions.cpp - src/fdtd_grid_initialiser.cpp - src/grid_labels.cpp - src/interface.cpp - src/interpolate.cpp - src/iterator.cpp - src/interpolation_methods.cpp - src/matrix_collection.cpp - src/mesh_base.cpp - src/numerical_derivative.cpp - src/shapes.cpp - src/simulation_parameters.cpp - src/source.cpp - src/timer.cpp - src/utils.cpp - matlabio/matlabio.cpp - ) + +# wildcard to find all *.cpp files in the src directory except for main.cpp +# (which we don't want for the testing target) +file(GLOB_RECURSE SOURCES "src/*.cpp") +list(FILTER SOURCES EXCLUDE REGEX "src/main.cpp") + +# TODO: can delete this line when matlabio is removed. +set(SOURCES ${SOURCES} matlabio/matlabio.cpp) if (BUILD_TESTING) test_target() diff --git a/tdms/cmake/targets.cmake b/tdms/cmake/targets.cmake index e828cdac0..1f7caae94 100644 --- a/tdms/cmake/targets.cmake +++ b/tdms/cmake/targets.cmake @@ -1,7 +1,7 @@ function(release_target) add_executable(tdms) - target_sources(tdms PRIVATE ${SOURCES} src/openandorder.cpp) + target_sources(tdms PRIVATE ${SOURCES} src/main.cpp) target_link_libraries(tdms FFTW::Double @@ -36,7 +36,7 @@ function(test_target) add_library(tdms_lib SHARED) target_sources(tdms_lib PUBLIC ${SOURCES}) - add_executable(tdms "src/openandorder.cpp") + add_executable(tdms "src/main.cpp") target_link_libraries(tdms_lib LINK_PUBLIC FFTW::Double diff --git a/tdms/include/arrays.h b/tdms/include/arrays.h index d5ea06cad..c9e5353fd 100644 --- a/tdms/include/arrays.h +++ b/tdms/include/arrays.h @@ -91,6 +91,18 @@ class MaterialCollection{ static void init_xyz_vectors(const mxArray *ptr, XYZVectors &arrays, const std::string &prefix); }; +/** + * @brief A class to encapsulate collection of algebraic terms in the + * discretized forms of Maxwells equations for E fields. The symbol chosen in + * the original reference is \f$C\f$. + * + * @details Algebraic terms \f$C_{a,b,c}\f$ defined in Section 4.2 of Munro, P,. + * "Application of numerical methods to high numerical aperture imaging", 2006, + * PhD thesis, Imperial College London. + * + * The definitions are equations 4.13, 4.14 (pp 82-3). Part of Maxwell's E-field + * equations in equations 4.7-9. + */ class CCollectionBase { public: XYZVectors a; @@ -98,7 +110,7 @@ class CCollectionBase { XYZVectors c; }; -// TODO: docstring +/*! @copydoc CCollectionBase */ class CCollection : public CCollectionBase { private: void init_xyz_vectors(const mxArray *ptr, XYZVectors &arrays, const std::string &prefix); @@ -110,19 +122,31 @@ class CCollection : public CCollectionBase { explicit CCollection(const mxArray *ptr); }; -// TODO: docstring +/*! @copydoc CCollectionBase */ class CMaterial : public CCollectionBase, MaterialCollection { public: explicit CMaterial(const mxArray *ptr); }; +/** + * @brief A class to encapsulate collection of algebraic terms in the + * discretized forms of Maxwells equations for H fields. The symbol chosen in + * the original reference is \f$D\f$. + * + * @details Algebraic terms \f$D_{a,b}\f$ defined in Section 4.2 of Munro, P,. + * "Application of numerical methods to high numerical aperture imaging", 2006, + * PhD thesis, Imperial College London. + * + * The definitions are equations 4.15, 4.16 (pp 82-3). Part of Maxwell's H-field + * equations in equations 4.10-12. + */ class DCollectionBase { public: XYZVectors a; XYZVectors b; }; -// TODO: docstring +/*! @copydoc DCollectionBase */ class DCollection: public DCollectionBase{ private: static void init_xyz_vectors(const mxArray *ptr, XYZVectors &arrays, const std::string &prefix); @@ -131,7 +155,7 @@ class DCollection: public DCollectionBase{ explicit DCollection(const mxArray *ptr); }; -// TODO: docstring +/*! @copydoc DCollectionBase */ class DMaterial : public DCollectionBase, MaterialCollection { public: explicit DMaterial(const mxArray *ptr); diff --git a/tdms/include/mesh_base.h b/tdms/include/mesh_base.h index 5a507cee8..7e04eecf5 100644 --- a/tdms/include/mesh_base.h +++ b/tdms/include/mesh_base.h @@ -94,6 +94,6 @@ void triangulateCuboidSkip(int I0, int I1, int J0, int J1, int K0, int K1, mxArr void conciseTriangulateCuboid(int I0, int I1, int J0, int J1, int K0, int K1, mxArray **vertices, mxArray **facets); void conciseTriangulateCuboidSkip(int I0, int I1, int J0, int J1, int K0, int K1, - PhasorInc &phasorinc, mxArray **vertices, mxArray **facets); + SurfaceSpacingStride &spacing_stride, mxArray **vertices, mxArray **facets); void conciseCreateBoundary(int I0, int I1,int K0, int K1, mxArray **vertices, mxArray ** facets); diff --git a/tdms/include/simulation_parameters.h b/tdms/include/simulation_parameters.h index 230a48d1b..5c7d6a77d 100644 --- a/tdms/include/simulation_parameters.h +++ b/tdms/include/simulation_parameters.h @@ -34,9 +34,9 @@ enum RunMode{ }; enum Dimension{ - THREE, // Full dimensionality - compute all H and E components - TE, // Transverse electric - only compute Ex, Ey, and Hz components - TM // Transverse magnetic - only compute Hx, Hy, and Ez components + THREE, //< Full dimensionality - compute all H and E components + TE, //< Transverse electric - only compute Ex, Ey, and Hz components + TM //< Transverse magnetic - only compute Hx, Hy, and Ez components }; enum InterpolationMethod{ @@ -50,7 +50,7 @@ enum InterpolationMethod{ * to extract the phasors on. the surface i.e. x = 2 means extract from * every 2nd Yee cell. */ -struct PhasorInc{ +struct SurfaceSpacingStride{ int x = 1; int y = 1; int z = 1; @@ -74,36 +74,36 @@ class SimulationParameters{ public: SimulationParameters(); - double omega_an = 0.0; // Angular ω - unsigned int Nt = 0; // Number of simulation steps - unsigned int Np = 0; // Number of times to extract the phasors - unsigned int Npe = 0; // Extract the phasors every this number of iterations - unsigned int start_tind = 0; // Starting iteration number for the time steps - double dt = 0.0; // Time step - bool has_tdfdir = false; // Is the tdfdir (time domain field directory) defined? - bool is_multilayer = false; // Is this simulation of a multilayer? - bool is_disp_ml = false; // Is the material dispersive? - double to_l = 0.0; // Time delay of pulse - double hwhm = 0.0; // Half width at half max of pulse - PerfectlyMatchedLayer pml; // Perfectly matched layer struct with size attributes - bool exphasorsvolume = false; // Should phasors be extracted in the whole volume? - bool exphasorssurface = false; // Should phasors be extracted on a surface? - bool intphasorssurface = false; // Should phasors be extracted/interpolated? - RunMode run_mode = complete; // Run mode - SourceMode source_mode = pulsed; // Source mode - Dimension dimension = THREE; // Dimensions to calculate in - bool is_structure = false; // Has a grating structure been defined? - bool exdetintegral = false; // TODO: detector sensitivity evaluation ? - int k_det_obs = 0; // TODO: no idea what this is?! - double z_obs = 0.0; // Value of the input grid_labels_z at k_det_obs + double omega_an = 0.0; //< Angular ω + unsigned int Nt = 0; //< Number of simulation steps + unsigned int Np = 0; //< Number of times to extract the phasors + unsigned int Npe = 0; //< Extract the phasors every this number of iterations + unsigned int start_tind = 0; //< Starting iteration number for the time steps + double dt = 0.0; //< Time step + bool has_tdfdir = false; //< Is the tdfdir (time domain field directory) defined? + bool is_multilayer = false; //< Is this simulation of a multilayer? + bool is_disp_ml = false; //< Is the material dispersive? + double to_l = 0.0; //< Time delay of pulse + double hwhm = 0.0; //< Half width at half max of pulse + PerfectlyMatchedLayer pml; //< Perfectly matched layer struct with size attributes + bool exphasorsvolume = false; //< Should phasors be extracted in the whole volume? + bool exphasorssurface = false; //< Should phasors be extracted on a surface? + bool intphasorssurface = false; //< Should phasors be extracted/interpolated? + RunMode run_mode = complete; //< Run mode + SourceMode source_mode = pulsed; //< Source mode + Dimension dimension = THREE; //< Dimensions to calculate in + bool is_structure = false; //< Has a grating structure been defined? + bool exdetintegral = false; //< TODO: detector sensitivity evaluation ? + int k_det_obs = 0; //< TODO: no idea what this is?! + double z_obs = 0.0; //< Value of the input grid_labels_z at k_det_obs bool air_interface_present = false; - double air_interface = 0.0; // TODO: what is this?! - bool interp_mat_props = false; // Should the material properties be interpolated? - InterpolationMethod interp_method = cubic; // Type of surface field interpolation to do - bool exi_present = false; // Is the time dependent x incident field present? - bool eyi_present = false; // Is the time dependent y incident field present? - PhasorInc phasorinc; // Surface stride for extracting phasors - YeeCellDimensions delta; // Yee cell dimensions (dx, dy, dz) + double air_interface = 0.0; //< TODO: what is this?! + bool interp_mat_props = false; //< Should the material properties be interpolated? + InterpolationMethod interp_method = cubic; //< Type of surface field interpolation to do + bool exi_present = false; //< Is the time dependent x incident field present? + bool eyi_present = false; //< Is the time dependent y incident field present? + SurfaceSpacingStride spacing_stride; //< Surface stride for extracting phasors (in matlab this is called 'phasorinc') + YeeCellDimensions delta; //< Yee cell dimensions (dx, dy, dz) void set_run_mode(std::string mode_string); @@ -111,7 +111,7 @@ class SimulationParameters{ void set_dimension(std::string mode_string); - void set_phasorinc(const double* vector); + void set_spacing_stride(const double* vector); void set_Np(FrequencyExtractVector &f_ex_vec); }; diff --git a/tdms/src/argument_parser.cpp b/tdms/src/argument_parser.cpp index 4027b81fa..a95e8c724 100644 --- a/tdms/src/argument_parser.cpp +++ b/tdms/src/argument_parser.cpp @@ -47,11 +47,11 @@ ArgumentNamespace ArgumentParser::parse_args(int n_args, char *arg_ptrs[]) { void ArgumentParser::print_help_message(){ fprintf(stdout,"Usage:\n" - "openandorder [options] infile outfile\n" - "openandorder [options] infile gridfile outfile\n" + "tdms [options] infile outfile\n" + "tdms [options] infile gridfile outfile\n" "Options:\n" "-h:\tDisplay this help message\n" - "--finite-difference:\tUse the finite-difference solver, instead of the pseudo-spectral method.\n" + "-fd, --finite-difference:\tUse the finite-difference solver, instead of the pseudo-spectral method.\n" "-q:\tQuiet operation. Silence all logging\n" "-m:\tMinimise output file size by not saving vertex and facet information\n\n"); } diff --git a/tdms/src/iterator.cpp b/tdms/src/iterator.cpp index f767b13bb..d69ff8376 100644 --- a/tdms/src/iterator.cpp +++ b/tdms/src/iterator.cpp @@ -412,7 +412,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs input_counter++; //fprintf(stderr,"Got phasorsurface\n"); - params.set_phasorinc(mxGetPr(prhs[input_counter++])); + params.set_spacing_stride(mxGetPr(prhs[input_counter++])); params.set_dimension(string_in(prhs[input_counter++], "dimension")); /*Get conductive_aux */ @@ -582,7 +582,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs &mx_surface_facets); else conciseTriangulateCuboidSkip(cuboid[0], cuboid[1], cuboid[2], cuboid[3], cuboid[4], cuboid[5], - params.phasorinc, &mx_surface_vertices, + params.spacing_stride, &mx_surface_vertices, &mx_surface_facets); //fprintf(stderr,"Qos 00a:\n"); //we don't need the facets so destroy the matrix now to save memory @@ -4636,7 +4636,7 @@ void execute_simulation(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs &mx_surface_facets); else conciseTriangulateCuboidSkip(cuboid[0], cuboid[1], cuboid[2], cuboid[3], cuboid[4], cuboid[5], - params.phasorinc, &dummy_vertex_list, + params.spacing_stride, &dummy_vertex_list, &mx_surface_facets); mxDestroyArray(dummy_vertex_list); mxArray *vertex_list; diff --git a/tdms/src/main.cpp b/tdms/src/main.cpp new file mode 100644 index 000000000..3ef212339 --- /dev/null +++ b/tdms/src/main.cpp @@ -0,0 +1,67 @@ +/** + * @file main.cpp + * @brief The main function. Launches TDMS. + */ +#include + +#include "openandorder.h" +#include "iterator.h" +#include "mat_io.h" + +#define NMATRICES 49 //< number of input matrices +#define NOUTMATRICES_WRITE 23 //< number of output matrices to be written to output file +#define NOUTMATRICES_WRITE_ALL 25 //< number of output matrices to be written to output file +#define NOUTMATRICES_PASSED 31 //< number of output matrices passed by mexFunction + +int main(int nargs, char *argv[]){ + + // Set the logging level with a compile-time define for debugging + #if SPDLOG_ACTIVE_LEVEL == SPDLOG_LEVEL_DEBUG + spdlog::set_level(spdlog::level::debug); + #elif SPDLOG_ACTIVE_LEVEL == SPDLOG_LEVEL_INFO + spdlog::set_level(spdlog::level::info); + #endif + + /* + There are two cases to consider, when the fdtdgrid matrix is specified in a separate mat file + and when it is in the same file as the other matrices. + */ + + const char *matrixnames[] = {"fdtdgrid","Cmaterial","Dmaterial","C","D","freespace","disp_params","delta","interface","Isource","Jsource","Ksource","grid_labels","omega_an","to_l","hwhm","Dxl","Dxu","Dyl","Dyu","Dzl","Dzu","Nt","dt","tind","sourcemode","runmode","exphasorsvolume","exphasorssurface","intphasorssurface","phasorsurface","phasorinc","dimension","conductive_aux","dispersive_aux","structure","f_ex_vec","exdetintegral","f_vec","Pupil","D_tilde","k_det_obs_global","air_interface","intmatprops","intmethod","tdfield","tdfdir","fieldsample","campssample"}; + const char *matrixnames_infile[] = {"Cmaterial","Dmaterial","C","D","freespace","disp_params","delta","interface","Isource","Jsource","Ksource","grid_labels","omega_an","to_l","hwhm","Dxl","Dxu","Dyl","Dyu","Dzl","Dzu","Nt","dt","tind","sourcemode","runmode","exphasorsvolume","exphasorssurface","intphasorssurface","phasorsurface","phasorinc","dimension","conductive_aux","dispersive_aux","structure","f_ex_vec","exdetintegral","f_vec","Pupil","D_tilde","k_det_obs_global","air_interface","intmatprops","intmethod","tdfield","tdfdir","fieldsample","campssample"}; + const char *matrixnames_gridfile[] = {"fdtdgrid"}; + const char *outputmatrices_all[] = {"Ex_out","Ey_out","Ez_out","Hx_out","Hy_out","Hz_out","x_out","y_out","z_out","Ex_i","Ey_i","Ez_i","Hx_i","Hy_i","Hz_i","x_i","y_i","z_i","vertices","camplitudes","facets","maxresfield","Id","fieldsample","campssample"}; + const char *outputmatrices[] = {"Ex_out","Ey_out","Ez_out","Hx_out","Hy_out","Hz_out","x_out","y_out","z_out","Ex_i","Ey_i","Ez_i","Hx_i","Hy_i","Hz_i","x_i","y_i","z_i","camplitudes","maxresfield","Id","fieldsample","campssample"}; + int matricestosave_all[] = {0,1,2,3,4,5,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28}; + int matricestosave[] = {0,1,2,3,4,5,10,11,12,13,14,15,16,17,18,19,20,21,23,25,26,27,28}; + mxArray *plhs[NOUTMATRICES_PASSED]; + const mxArray *matrixptrs[NMATRICES]; + + auto args = ArgumentParser::parse_args(nargs, argv); + check_files_can_be_accessed(args); + + //now it is safe to use matlab routines to open the file and order the matrices + if(!args.has_grid_filename()) + openandorder(args.input_filename(), (char **)matrixnames, matrixptrs, NMATRICES); + else{ + openandorder(args.input_filename(), (char **)matrixnames_infile, matrixptrs+1, NMATRICES-1); + openandorder(args.grid_filename(), (char **)matrixnames_gridfile, matrixptrs, 1); + } + + // decide which derivative method to use (PSTD or FDTD) + SolverMethod method = PseudoSpectral; // default + if (args.finite_difference()) + method = SolverMethod::FiniteDifference; + + //now run the time propagation code + execute_simulation(NOUTMATRICES_PASSED, (mxArray **)plhs, NMATRICES, (const mxArray **)matrixptrs, method); + + if( !args.have_flag("-m") ){ //prints vertices and facets + saveoutput(plhs, matricestosave_all, (char **)outputmatrices_all, NOUTMATRICES_WRITE_ALL, args.output_filename()); + } + else{ // minimise the file size by not printing vertices and facets + saveoutput(plhs, matricestosave, (char **)outputmatrices, NOUTMATRICES_WRITE, args.output_filename()); + } + + return 0; +} diff --git a/tdms/src/mesh_base.cpp b/tdms/src/mesh_base.cpp index 09c36d6cb..5251bb1a9 100644 --- a/tdms/src/mesh_base.cpp +++ b/tdms/src/mesh_base.cpp @@ -537,9 +537,9 @@ void conciseCreateBoundary(int I0, int I1,int K0, int K1, } void conciseTriangulateCuboidSkip(int I0, int I1, int J0, int J1, int K0, int K1, - PhasorInc &phasorinc, mxArray **vertices, mxArray ** facets){ + SurfaceSpacingStride &spacing_stride, mxArray **vertices, mxArray ** facets){ - int dI = phasorinc.x, dJ = phasorinc.y, dK = phasorinc.z; + int dI = spacing_stride.x, dJ = spacing_stride.y, dK = spacing_stride.z; mxArray *triangles[6]; mxArray *index_map; int ndims, ***index_map_int, nindices, **vertices_int, i, **facets_int, j,k, **triangles_int,ii,jj,kk; diff --git a/tdms/src/openandorder.cpp b/tdms/src/openandorder.cpp index 5fc54c967..e3a565d5f 100644 --- a/tdms/src/openandorder.cpp +++ b/tdms/src/openandorder.cpp @@ -1,6 +1,6 @@ /** * @file openandorder.cpp - * @brief Launch and file IO + * @brief Passing arguments and file IO * * Code for processing command line arguments, opening input files, passing * matrices to the mexFunction and writing the output to the specified output @@ -11,73 +11,12 @@ #include #include -#include - #include "utils.h" #include "fdtd_grid_initialiser.h" #include "iterator.h" -#define NMATRICES 49 //< number of input matrices -#define NOUTMATRICES_WRITE 23 //< number of output matrices to be written to output file -#define NOUTMATRICES_WRITE_ALL 25 //< number of output matrices to be written to output file -#define NOUTMATRICES_PASSED 31 //< number of output matrices passed by mexFunction - using namespace std; - -int main(int nargs, char *argv[]){ - - // Set the logging level with a compile-time define for debugging - #if SPDLOG_ACTIVE_LEVEL == SPDLOG_LEVEL_DEBUG - spdlog::set_level(spdlog::level::debug); - #elif SPDLOG_ACTIVE_LEVEL == SPDLOG_LEVEL_INFO - spdlog::set_level(spdlog::level::info); - #endif - - /* - There are two cases to consider, when the fdtdgrid matrix is specified in a separate mat file - and when it is in the same file as the other matrices. - */ - - const char *matrixnames[] = {"fdtdgrid","Cmaterial","Dmaterial","C","D","freespace","disp_params","delta","interface","Isource","Jsource","Ksource","grid_labels","omega_an","to_l","hwhm","Dxl","Dxu","Dyl","Dyu","Dzl","Dzu","Nt","dt","tind","sourcemode","runmode","exphasorsvolume","exphasorssurface","intphasorssurface","phasorsurface","phasorinc","dimension","conductive_aux","dispersive_aux","structure","f_ex_vec","exdetintegral","f_vec","Pupil","D_tilde","k_det_obs_global","air_interface","intmatprops","intmethod","tdfield","tdfdir","fieldsample","campssample"}; - const char *matrixnames_infile[] = {"Cmaterial","Dmaterial","C","D","freespace","disp_params","delta","interface","Isource","Jsource","Ksource","grid_labels","omega_an","to_l","hwhm","Dxl","Dxu","Dyl","Dyu","Dzl","Dzu","Nt","dt","tind","sourcemode","runmode","exphasorsvolume","exphasorssurface","intphasorssurface","phasorsurface","phasorinc","dimension","conductive_aux","dispersive_aux","structure","f_ex_vec","exdetintegral","f_vec","Pupil","D_tilde","k_det_obs_global","air_interface","intmatprops","intmethod","tdfield","tdfdir","fieldsample","campssample"}; - const char *matrixnames_gridfile[] = {"fdtdgrid"}; - const char *outputmatrices_all[] = {"Ex_out","Ey_out","Ez_out","Hx_out","Hy_out","Hz_out","x_out","y_out","z_out","Ex_i","Ey_i","Ez_i","Hx_i","Hy_i","Hz_i","x_i","y_i","z_i","vertices","camplitudes","facets","maxresfield","Id","fieldsample","campssample"}; - const char *outputmatrices[] = {"Ex_out","Ey_out","Ez_out","Hx_out","Hy_out","Hz_out","x_out","y_out","z_out","Ex_i","Ey_i","Ez_i","Hx_i","Hy_i","Hz_i","x_i","y_i","z_i","camplitudes","maxresfield","Id","fieldsample","campssample"}; - int matricestosave_all[] = {0,1,2,3,4,5,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28}; - int matricestosave[] = {0,1,2,3,4,5,10,11,12,13,14,15,16,17,18,19,20,21,23,25,26,27,28}; - mxArray *plhs[NOUTMATRICES_PASSED]; - const mxArray *matrixptrs[NMATRICES]; - - auto args = ArgumentParser::parse_args(nargs, argv); - check_files_can_be_accessed(args); - - //now it is safe to use matlab routines to open the file and order the matrices - if(!args.has_grid_filename()) - openandorder(args.input_filename(), (char **)matrixnames, matrixptrs, NMATRICES); - else{ - openandorder(args.input_filename(), (char **)matrixnames_infile, matrixptrs+1, NMATRICES-1); - openandorder(args.grid_filename(), (char **)matrixnames_gridfile, matrixptrs, 1); - } - - // decide which derivative method to use (PSTD or FDTD) - SolverMethod method = PseudoSpectral; // default - if (args.finite_difference()) - method = SolverMethod::FiniteDifference; - - //now run the time propagation code - execute_simulation(NOUTMATRICES_PASSED, (mxArray **)plhs, NMATRICES, (const mxArray **)matrixptrs, method); - - if( !args.have_flag("-m") ){ //prints vertices and facets - saveoutput(plhs, matricestosave_all, (char **)outputmatrices_all, NOUTMATRICES_WRITE_ALL, args.output_filename()); - } - else{ // minimise the file size by not printing vertices and facets - saveoutput(plhs, matricestosave, (char **)outputmatrices, NOUTMATRICES_WRITE, args.output_filename()); - } - - return 0; -} - void openandorder(const char *mat_filename, char **matrix_names, const mxArray **matrix_ptrs, int n_matrices){ auto expected = MatrixCollection(matrix_names, n_matrices); diff --git a/tdms/src/simulation_parameters.cpp b/tdms/src/simulation_parameters.cpp index f3dfd2f04..381af78b6 100644 --- a/tdms/src/simulation_parameters.cpp +++ b/tdms/src/simulation_parameters.cpp @@ -43,10 +43,10 @@ void SimulationParameters::set_dimension(string mode_string) { } } -void SimulationParameters::set_phasorinc(const double* vector) { - phasorinc.x = (int) vector[0]; - phasorinc.y = (int) vector[1]; - phasorinc.z = (int) vector[2]; +void SimulationParameters::set_spacing_stride(const double* vector) { + spacing_stride.x = (int) vector[0]; + spacing_stride.y = (int) vector[1]; + spacing_stride.z = (int) vector[2]; } void SimulationParameters::set_Np(FrequencyExtractVector &f_ex_vec) {