diff --git a/docker/Dockerfile b/docker/Dockerfile index ba97fd92..7ee5033b 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -63,7 +63,7 @@ ENV LIBROM_DIR=$LIB_DIR/libROM # install python RUN sudo apt-get update && sudo apt-get install -yq python3 python3-pip -RUN sudo pip3 install --upgrade pip && sudo pip3 install numpy scipy argparse tables PyYAML h5py +RUN sudo pip3 install --upgrade pip && sudo pip3 install numpy scipy argparse tables PyYAML h5py matplotlib # install h5dump and gmsh RUN sudo apt-get install -yq hdf5-tools gmsh diff --git a/examples/linelast/CMakeLists.txt b/examples/linelast/CMakeLists.txt index 573c89c5..b7bc8bb9 100644 --- a/examples/linelast/CMakeLists.txt +++ b/examples/linelast/CMakeLists.txt @@ -7,9 +7,11 @@ file(COPY config/linelast.simpleL.h5 DESTINATION ${CMAKE_BINARY_DIR}/examples/li file(COPY linelast.lattice.yml DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast) file(COPY config/linelast.lattice.h5 DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/config) file(COPY linelast.force_cantilever.yml DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast) +file(COPY linelast.comp_train.yml DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast) +file(COPY config/linelast.comp_train.h5 DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/config) file(COPY meshes/beam-tri.mesh DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/meshes) file(COPY meshes/beam-tet.mesh DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/meshes) file(COPY meshes/joint2D.mesh DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/meshes) file(COPY meshes/rod2D_H.mesh DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/meshes) -file(COPY meshes/rod2D_V.mesh DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/meshes) \ No newline at end of file +file(COPY meshes/rod2D_V.mesh DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/meshes) diff --git a/examples/linelast/config/linelast.comp_train.h5 b/examples/linelast/config/linelast.comp_train.h5 new file mode 100644 index 00000000..91c34837 Binary files /dev/null and b/examples/linelast/config/linelast.comp_train.h5 differ diff --git a/examples/linelast/config/linelast.lattice.h5 b/examples/linelast/config/linelast.lattice.h5 index 93da8e1e..e020e024 100644 Binary files a/examples/linelast/config/linelast.lattice.h5 and b/examples/linelast/config/linelast.lattice.h5 differ diff --git a/examples/linelast/linelast.comp_train.yml b/examples/linelast/linelast.comp_train.yml new file mode 100644 index 00000000..14b568db --- /dev/null +++ b/examples/linelast/linelast.comp_train.yml @@ -0,0 +1,227 @@ +main: + mode: single_run + use_rom: true + solver: linelast + +mesh: + type: component-wise + component-wise: + global_config: "config/linelast.comp_train.h5" + components: + - name: "joint2D" + file: "meshes/joint2D.mesh" + - name: "rod2D_H" + file: "meshes/rod2D_H.mesh" + - name: "rod2D_V" + file: "meshes/rod2D_V.mesh" + uniform_refinement: 0 + +domain-decomposition: + type: interior_penalty + +solver: + direct_solve: true + +discretization: + order: 1 + full-discrete-galerkin: true + +visualization: + enabled: false + unified_paraview: true + file_path: + prefix: paraview_output + +parameterized_problem: + name: linelast_cwtrain + +single_run: + linelast_cwtrain: + l_ux: 0.0 + l_uy: 0.0 + r_fx: 0.0 + r_fy: -2.0e-2 + u_fx: 0.0 + u_fy: 0.0 + d_fx: 0.0 + d_fy: 0.0 + lx: 0.5 + ly: 0.5 + rx: 0.0 + ry: 0.5 + dx: 0.0 + dy: 0.0 + ux: 0.0 + uy: 0.0 + xu_amp: 1.0 + xu_freq: 1.0 + xu_offset: 0.0 + xf_amp: 1.0 + xf_freq: 1.0 + xf_offset: 0.0 + yu_amp: 1.0 + yu_freq: 1.0 + yu_offset: 0.0 + yf_amp: 1.0 + yf_freq: 1.0 + yf_offset: 0.0 + +sample_generation: + maximum_number_of_snapshots: 70000 + component_sampling: true + type: random + random_sample_generator: + number_of_samples: 100 + file_path: + prefix: "linelast_cwtrain" + parameters: + - key: single_run/linelast_cwtrain/l_ux + type: double + minimum: 0.0 + maximum: 10.0 + - key: single_run/linelast_cwtrain/l_uy + type: double + minimum: 0.0 + maximum: 10.0 + - key: single_run/linelast_cwtrain/r_fx + type: double + minimum: -20.0 + maximum: 20.0 + - key: single_run/linelast_cwtrain/r_fy + type: double + minimum: -20.0 + maximum: 20.0 + - key: single_run/linelast_cwtrain/u_fx + type: double + minimum: -20.0 + maximum: 20.0 + - key: single_run/linelast_cwtrain/u_fy + type: double + minimum: -20.0 + maximum: 20.0 + - key: single_run/linelast_cwtrain/d_fx + type: double + minimum: -20.0 + maximum: 20.0 + - key: single_run/linelast_cwtrain/d_fy + type: double + minimum: -20.0 + maximum: 20.0 + - key: single_run/linelast_cwtrain/rx + type: double + minimum: 0.2 + maximum: 1.0 + - key: single_run/linelast_cwtrain/ry + type: double + minimum: 0.2 + maximum: 1.0 + - key: single_run/linelast_cwtrain/lx + type: double + minimum: 0.5 + maximum: 1.0 + - key: single_run/linelast_cwtrain/ly + type: double + minimum: 0.5 + maximum: 1.0 + - key: single_run/linelast_cwtrain/dx + type: double + minimum: 0.3 + maximum: 1.0 + - key: single_run/linelast_cwtrain/dy + type: double + minimum: 0.3 + maximum: 1.0 + - key: single_run/linelast_cwtrain/ux + type: double + minimum: 0.3 + maximum: 1.0 + - key: single_run/linelast_cwtrain/uy + type: double + minimum: 0.3 + maximum: 1.0 + - key: single_run/linelast_cwtrain/xu_offset + type: double + minimum: 0.0 + maximum: 1.0 + - key: single_run/linelast_cwtrain/yu_offset + type: double + minimum: 0.0 + maximum: 1.0 + - key: single_run/linelast_cwtrain/xf_offset + type: double + minimum: 0.0 + maximum: 1.0 + - key: single_run/linelast_cwtrain/yf_offset + type: double + minimum: 0.0 + maximum: 1.0 + - key: single_run/linelast_cwtrain/xu_amp + type: double + minimum: -10.0 + maximum: 10.0 + - key: single_run/linelast_cwtrain/yu_amp + type: double + minimum: -10.0 + maximum: 10.0 + - key: single_run/linelast_cwtrain/xf_amp + type: double + minimum: -10.0 + maximum: 10.0 + - key: single_run/linelast_cwtrain/yf_amp + type: double + minimum: -10.0 + maximum: 10.0 + - key: single_run/linelast_cwtrain/xu_freq + type: double + minimum: 0.1 + maximum: 1.0 + - key: single_run/linelast_cwtrain/yu_freq + type: double + minimum: 0.1 + maximum: 1.0 + - key: single_run/linelast_cwtrain/xf_freq + type: double + minimum: 0.1 + maximum: 1.0 + - key: single_run/linelast_cwtrain/yf_freq + type: double + minimum: 0.1 + maximum: 1.0 + +basis: + prefix: "linelast_cwtrain" + number_of_basis: 16 + tags: + - name: "joint2D" + number_of_basis: 16 + snapshot_format: + minimum: 0 + maximum: 3 + format: "linelast_cwtrain_sample_joint2D_snapshot%01d" + - name: "rod2D_H" + number_of_basis: 64 + snapshot_format: + minimum: 0 + maximum: 3 + format: "linelast_cwtrain_sample_rod2D_H_snapshot%01d" + - name: "rod2D_V" + number_of_basis: 64 + snapshot_format: + minimum: 0 + maximum: 3 + format: "linelast_cwtrain_sample_rod2D_V_snapshot%01d" + svd: + save_spectrum: true + update_right_sv: false + visualization: + enabled: false + +model_reduction: + rom_handler_type: mfem + subdomain_training: universal + linear_solver_type: direct + save_operator: + level: component + prefix: "linelast_cw.rom_elem" + compare_solution: + enabled: true diff --git a/examples/linelast/linelast.lattice.yml b/examples/linelast/linelast.lattice.yml index e0013c6c..26f386bf 100644 --- a/examples/linelast/linelast.lattice.yml +++ b/examples/linelast/linelast.lattice.yml @@ -26,7 +26,7 @@ visualization: enabled: false unified_paraview: true file_path: - prefix: lattice_output + prefix: paraview_output parameterized_problem: name: linelast_disp_lattice diff --git a/include/parameterized_problem.hpp b/include/parameterized_problem.hpp index dccca6f0..0800b497 100644 --- a/include/parameterized_problem.hpp +++ b/include/parameterized_problem.hpp @@ -92,6 +92,47 @@ extern double rforce_f; void tip_force(const Vector &x, Vector &u); } +namespace linelast_cwtrain +{ +extern double l_ux; +extern double l_uy; +extern double r_fx; +extern double r_fy; +extern double u_fx; +extern double u_fy; +extern double d_fx; +extern double d_fy; + +extern double lx; +extern double ly; +extern double rx; +extern double ry; +extern double dx; +extern double dy; +extern double ux; +extern double uy; + +extern double xu_amp; +extern double xu_freq; +extern double xu_offset; +extern double xf_amp; +extern double xf_freq; +extern double xf_offset; +extern double yu_amp; +extern double yu_freq; +extern double yu_offset; +extern double yf_amp; +extern double yf_freq; +extern double yf_offset; + +double perturb_func(const double x, const double amp, const double freq, const double offset); +void left_disp(const Vector &x, Vector &u); +void up_disp(const Vector &x, Vector &u); +void down_disp(const Vector &x, Vector &u); +void right_disp(const Vector &x, Vector &u); + +} + namespace advdiff_problem { @@ -305,6 +346,12 @@ class AdvDiffFlowPastArray : public StokesFlowPastArray } }; +class LinElastComponentWiseTrain : public LinElastProblem +{ +public: + LinElastComponentWiseTrain(); +}; + ParameterizedProblem* InitParameterizedProblem(); #endif diff --git a/src/parameterized_problem.cpp b/src/parameterized_problem.cpp index 59a7b73c..fbe4c804 100644 --- a/src/parameterized_problem.cpp +++ b/src/parameterized_problem.cpp @@ -226,6 +226,82 @@ double qbdr(const Vector &x) } // namespace advdiff_problem +namespace linelast_cwtrain +{ + +double l_ux; +double l_uy; +double r_fx; +double r_fy; +double u_fx; +double u_fy; +double d_fx; +double d_fy; + +double xu_amp; +double xu_freq; +double xu_offset; + +double xf_amp; +double xf_freq; +double xf_offset; + +double yu_amp; +double yu_freq; +double yu_offset; + +double yf_amp; +double yf_freq; +double yf_offset; + +double lx; +double ly; +double rx; +double ry; +double dx; +double dy; +double ux; +double uy; + +double perturb_func(const double x, const double amp, const double freq, const double offset) +{ + return amp * sin(pi * freq *( x / 3.0 + 2 * offset) ); +} + +void left_disp(const Vector &x, Vector &u) +{ + if (lx >= 0.5) + u(0) = l_ux + perturb_func(x(0), xu_amp, xu_freq, xu_offset); + if (ly >= 0.5) + u(1) = l_uy + perturb_func(x(1), yu_amp, yu_freq, yu_offset); +} + +void right_disp(const Vector &x, Vector &u) +{ + if (rx >= 0.5) + u(0) = r_fx + perturb_func(x(0), xf_amp, xf_freq, xf_offset); + if (ry >= 0.5) + u(1) = r_fy + perturb_func(x(1), yf_amp, yf_freq, yf_offset); +} + +void up_disp(const Vector &x, Vector &u) +{ + if (ux >= 0.5) + u(0) = u_fx + perturb_func(x(0), xf_amp, xf_freq, xf_offset); + if (uy >= 0.5) + u(1) = u_fy + perturb_func(x(1), yf_amp, yf_freq, yf_offset); +} + +void down_disp(const Vector &x, Vector &u) +{ + if (dx >= 0.5) + u(0) = d_fx + perturb_func(x(0), xf_amp, xf_freq, xf_offset); + if (dy >= 0.5) + u(1) = d_fy + perturb_func(x(1), yf_amp, yf_freq, yf_offset); +} + +} // namespace linelast_cwtrain + } // namespace function_factory ParameterizedProblem::ParameterizedProblem() @@ -329,6 +405,10 @@ ParameterizedProblem* InitParameterizedProblem() { problem = new LinElastForceCantilever(); } + else if (problem_name == "linelast_cwtrain") + { + problem = new LinElastComponentWiseTrain(); + } else if (problem_name == "advdiff_flow_past_array") { problem = new AdvDiffFlowPastArray(); @@ -746,7 +826,6 @@ LinElastDispLattice::LinElastDispLattice() } - LinElastForceCantilever::LinElastForceCantilever() : LinElastProblem() { @@ -840,3 +919,141 @@ AdvDiffFlowPastArray::~AdvDiffFlowPastArray() { delete flow_problem; } + +LinElastComponentWiseTrain::LinElastComponentWiseTrain() + : LinElastProblem() +{ + // pointer to static function. + bdr_type.SetSize(5); + battr.SetSize(5); + vector_bdr_ptr.SetSize(5); + + // Left + battr[0] = 1; + bdr_type[0] = BoundaryType::DIRICHLET; + vector_bdr_ptr[0] = &(function_factory::linelast_cwtrain::left_disp); + + // Right + battr[1] = 2; + bdr_type[1] = BoundaryType::NEUMANN; + vector_bdr_ptr[1] = &(function_factory::linelast_cwtrain::right_disp); + + // None + battr[2] = 5; + bdr_type[2] = BoundaryType::NEUMANN; + vector_bdr_ptr[2] = NULL; + + // Up + battr[3] = 4; + bdr_type[3] = BoundaryType::NEUMANN; + vector_bdr_ptr[3] = &(function_factory::linelast_cwtrain::up_disp); + + // Down + battr[4] = 3; + bdr_type[4] = BoundaryType::NEUMANN; + vector_bdr_ptr[4] = &(function_factory::linelast_cwtrain::down_disp); + + + // Set materials + general_scalar_ptr.SetSize(2); + general_scalar_ptr[0] = function_factory::linelast_problem::lambda; + general_scalar_ptr[1] = function_factory::linelast_problem::mu; + + // Default values. + function_factory::linelast_cwtrain::l_ux= 0.0; + function_factory::linelast_cwtrain::l_uy = 0.0; + function_factory::linelast_cwtrain::r_fx = 0.0; + function_factory::linelast_cwtrain::r_fy = 0.0; + function_factory::linelast_cwtrain::u_fx = 0.0; + function_factory::linelast_cwtrain::u_fy = 0.0; + function_factory::linelast_cwtrain::d_fx = 0.0; + function_factory::linelast_cwtrain::d_fy = 0.0; + function_factory::linelast_cwtrain::lx = 0.0; + function_factory::linelast_cwtrain::ly = 0.0; + function_factory::linelast_cwtrain::rx = 0.0; + function_factory::linelast_cwtrain::ry = 0.0; + function_factory::linelast_cwtrain::dx = 0.0; + function_factory::linelast_cwtrain::dy = 0.0; + function_factory::linelast_cwtrain::ux = 0.0; + function_factory::linelast_cwtrain::uy = 0.0; + function_factory::linelast_problem::_lambda = 1.0; + function_factory::linelast_problem::_mu = 1.0; + function_factory::linelast_cwtrain::xu_amp = 1.0; + function_factory::linelast_cwtrain::xu_freq = 1.0; + function_factory::linelast_cwtrain::xu_offset = 0.0; + function_factory::linelast_cwtrain::xf_amp = 1.0; + function_factory::linelast_cwtrain::xf_freq = 1.0; + function_factory::linelast_cwtrain::xf_offset = 0.0; + function_factory::linelast_cwtrain::yu_amp = 1.0; + function_factory::linelast_cwtrain::yu_freq = 1.0; + function_factory::linelast_cwtrain::yu_offset = 0.0; + function_factory::linelast_cwtrain::yf_amp = 1.0; + function_factory::linelast_cwtrain::yf_freq = 1.0; + function_factory::linelast_cwtrain::yf_offset = 0.0; + + param_map["l_ux"] = 0; + param_map["l_uy"] = 1; + param_map["r_fx"] = 2; + param_map["r_fy"] = 3; + param_map["u_fx"] = 4; + param_map["u_fy"] = 5; + param_map["d_fx"] = 6; + param_map["d_fy"] = 7; + param_map["lx"] = 8; + param_map["ly"] = 9; + param_map["rx"] = 10; + param_map["ry"] = 11; + param_map["dx"] = 12; + param_map["dy"] = 13; + param_map["ux"] = 14; + param_map["uy"] = 15; + param_map["lambda"] = 16; + param_map["mu"] = 17; + param_map["xu_amp"] = 18; + param_map["xu_freq"] = 19; + param_map["xu_offset"] = 20; + param_map["xf_amp"] = 21; + param_map["xf_freq"] = 22; + param_map["xf_offset"] = 23; + param_map["yu_amp"] = 24; + param_map["yu_freq"] = 25; + param_map["yu_offset"] = 26; + param_map["yf_amp"] = 27; + param_map["yf_freq"] = 28; + param_map["yf_offset"] = 29; + + param_ptr.SetSize(30); + param_ptr[0] = &(function_factory::linelast_cwtrain::l_ux); + param_ptr[1] = &(function_factory::linelast_cwtrain::l_uy); + param_ptr[2] = &(function_factory::linelast_cwtrain::r_fx); + param_ptr[3] = &(function_factory::linelast_cwtrain::r_fy); + param_ptr[4] = &(function_factory::linelast_cwtrain::u_fx); + param_ptr[5] = &(function_factory::linelast_cwtrain::u_fy); + param_ptr[6] = &(function_factory::linelast_cwtrain::d_fx); + param_ptr[7] = &(function_factory::linelast_cwtrain::d_fy); + param_ptr[8] = &(function_factory::linelast_cwtrain::lx); + param_ptr[9] = &(function_factory::linelast_cwtrain::ly); + param_ptr[10] = &(function_factory::linelast_cwtrain::rx); + param_ptr[11] = &(function_factory::linelast_cwtrain::ry); + param_ptr[12] = &(function_factory::linelast_cwtrain::dx); + param_ptr[13] = &(function_factory::linelast_cwtrain::dy); + param_ptr[14] = &(function_factory::linelast_cwtrain::ux); + param_ptr[15] = &(function_factory::linelast_cwtrain::uy); + param_ptr[16] = &(function_factory::linelast_problem::_lambda); + param_ptr[17] = &(function_factory::linelast_problem::_mu); + param_ptr[18] = &(function_factory::linelast_cwtrain::xu_amp); + param_ptr[19] = &(function_factory::linelast_cwtrain::xu_freq); + param_ptr[20] = &(function_factory::linelast_cwtrain::xu_offset); + param_ptr[21] = &(function_factory::linelast_cwtrain::xf_amp); + param_ptr[22] = &(function_factory::linelast_cwtrain::xf_freq); + param_ptr[23] = &(function_factory::linelast_cwtrain::xf_offset); + param_ptr[24] = &(function_factory::linelast_cwtrain::yu_amp); + param_ptr[25] = &(function_factory::linelast_cwtrain::yu_freq); + param_ptr[26] = &(function_factory::linelast_cwtrain::yu_offset); + param_ptr[27] = &(function_factory::linelast_cwtrain::yf_amp); + param_ptr[28] = &(function_factory::linelast_cwtrain::yf_freq); + param_ptr[29] = &(function_factory::linelast_cwtrain::yf_offset); + + general_vector_ptr.SetSize(1); + general_vector_ptr[0] = NULL; // for now, change if current params doesn't work well enough. +} \ No newline at end of file diff --git a/utils/python/linelast_comp.py b/utils/python/linelast_comp.py index 0a54b4f5..efbab729 100644 --- a/utils/python/linelast_comp.py +++ b/utils/python/linelast_comp.py @@ -148,7 +148,7 @@ def grid_mesh_configs(nx, ny, l, w): def LatticeCantilever(nx, ny): # nx and ny are the number of sections in x- and y-direction, respectively - l = 4.0 + l = 4.0 w = 1.0 mesh_configs, bdr_data, if_data, mesh_type, n_mesh = grid_mesh_configs(nx, ny, l, w) @@ -223,4 +223,4 @@ def LatticeCantilever(nx, ny): elif name == "lattice_cantilever": nx = int(sys.argv[2]) ny = int(sys.argv[3]) - LatticeCantilever(nx, ny) + LatticeCantilever(nx, ny) \ No newline at end of file diff --git a/utils/python/linelast_cwtrain.py b/utils/python/linelast_cwtrain.py new file mode 100644 index 00000000..02f0b3e8 --- /dev/null +++ b/utils/python/linelast_cwtrain.py @@ -0,0 +1,371 @@ +# Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +# +# SPDX-License-Identifier: MIT + +import numpy as np +import h5py +import os +import matplotlib.pyplot as plt + +def shift_mesh_battrs(num_attr, path): + # Read original file + f = open(path, "r") + all_lines = [x for x in f] + f.close() + f = open(path, "r") + bdr_start = next(x for x, val in enumerate(f) + if val == "boundary\n") + bdr_len = int(f.readline(bdr_start+1)) + bdr_txt = [f.readline(bdr_start+2 + i) for i in range(bdr_len)] + battrs = [int(txt.split(' ')[0]) for txt in bdr_txt] + first_lines = all_lines[0:bdr_start+2] + last_lines = all_lines[bdr_start+2+bdr_len:] + f.close() + prefix = path.split('.')[0] + for i in range(num_attr): + filename = prefix + "_"+ str(i) + ".mesh" + shift_battrs = [(b + i ) if (b + i) == num_attr else (b + i )%num_attr for b in battrs] + shift_battrs_txt = [' '.join([str(shift_battrs[i]), *txt.split(' ')[1:]]) for i, txt in enumerate(bdr_txt)] + f = open(filename, "w") + new_txt = ''.join([''.join(first_lines), ''.join(shift_battrs_txt), ''.join(last_lines)]) + f.write(new_txt) + f.close() + +def create_training_meshes(prefix): + shift_mesh_battrs(4, os.path.join(prefix, "joint2D.mesh")) + shift_mesh_battrs(4, os.path.join(prefix, "rod2D_H.mesh")) + shift_mesh_battrs(4, os.path.join(prefix, "rod2D_V.mesh")) + +def cwgrid_mesh_configs(nx, ny, l, w): + bdr_data = [] + if_data = [] + n_joint = (nx + 1) * (ny + 1) + n_rod_H = (ny + 1) * nx + n_rod_V = (nx + 1) * ny + n_rod = n_rod_H + n_rod_V + n_mesh = n_joint + n_rod + mesh_configs = np.zeros([n_mesh, 6]) + mesh_type = np.concatenate((np.full((int(n_joint),), 0.0),np.full((int(n_rod_H),), 1.0),np.full((int(n_rod_V),), 2.0))) + # Setup joints + for i in range(n_joint): + xi = (i + nx+1) % (nx+1) + yi = np.floor(i / (nx+1)) + mesh_configs[i,:] = [xi*(l+w), yi*(l+w), 0., 0., 0., 0.] + + # Boundary check + # global_battr / mesh_idx / comp_battr + if xi == 0.0: + bdr_data += [[1, i, 4]] # constrain one end + elif xi == nx: + bdr_data += [[2, i, 2]] # constrain one end + if yi==0.0: + bdr_data += [[4, i, 1]] # constrained ends + elif yi==ny: + bdr_data += [[5, i, 3]] # constrained ends + + # Interface check + # mesh1 / mesh2 / battr1 / battr2 / port_idx + # Case port idx = 0 + if xi < nx: + if_data += [[i, i + n_joint - yi, 2, 4, 0]] + # Case port idx = 1 + if xi > 0: + if_data += [[i, i + n_joint - yi - 1, 4, 2, 1]] + # Case port idx = 2 + if yi < ny: + if_data += [[i, n_joint + n_rod_H + i, 3, 1, 2]] + # Case port idx = 3 + if yi > 0: + if_data += [[i, n_joint + n_rod_H + i - (nx+1), 1, 3, 3]] + + # Setup horizontal rods + for i in range(n_rod_H): + xi = (i + nx) % (nx) + yi = np.floor(i / (nx)) + mesh_configs[n_joint + i,:] = [xi*(l+w) + w, yi*(l+w), 0., 0., 0., 0.] + # Boundary check + if yi==0.0: + bdr_data += [[3, n_joint + i, 3]] # free boundary + bdr_data += [[4, n_joint + i, 1]] # fixed boundary + + elif yi==ny: + bdr_data += [[5, n_joint + i, 3]] # fixed boundary + bdr_data += [[3, n_joint + i, 1]] # free boundary + else: + bdr_data += [[3, n_joint + i, 1]] # free boundary + bdr_data += [[3, n_joint + i, 3]] # free boundary + + # Setup vertical rods + for i in range(n_rod_V): + xi = (i + nx+1) % (nx+1) + yi = np.floor(i / (nx+1)) + mesh_configs[n_joint + n_rod_H + i,:] = [xi*(l+w), yi*(l+w) + w, 0., 0., 0., 0.] + + # Boundary check + if xi == 0.0: + bdr_data += [[1, n_joint + n_rod_H + i, 4]] # constrain one end + bdr_data += [[3, n_joint + n_rod_H + i, 2]] # free boundary + elif xi == nx: + bdr_data += [[2, n_joint + n_rod_H + i, 2]] # constrain one end + bdr_data += [[3, n_joint + n_rod_H + i, 4]] # free boundary + else: + bdr_data += [[3, n_joint + n_rod_H + i, 2]] # constrain one end + bdr_data += [[3, n_joint + n_rod_H + i, 4]] # free boundary + + return mesh_configs, bdr_data, if_data, mesh_type, n_mesh + +def CWTrain1x1(): + n_mesh = 8 + mesh_type = [0, 1, 0, 2, 2, 0, 1, 0] + mesh_configs = np.zeros([n_mesh, 6]) + mesh_configs[1,:] = [1., 0., 0., 0., 0., 0.] + mesh_configs[2,:] = [5., 0., 0., 0., 0., 0.] + mesh_configs[3,:] = [0., 1., 0., 0., 0., 0.] + mesh_configs[4,:] = [5., 1., 0., 0., 0., 0.] + mesh_configs[5,:] = [0., 5., 0., 0., 0., 0.] + mesh_configs[6,:] = [1., 5., 0., 0., 0., 0.] + mesh_configs[7,:] = [5., 5., 0., 0., 0., 0.] + + # interface data + # mesh1 / mesh2 / battr1 / battr2 / port_idx + if_data = np.zeros([8, 5]) + if_data[0, :] = [0, 1, 2, 4, 0] + if_data[1, :] = [2, 1, 4, 2, 1] + if_data[2, :] = [0, 3, 3, 1, 2] + if_data[3, :] = [2, 4, 3, 1, 2] + if_data[4, :] = [5, 3, 1, 3, 3] + if_data[5, :] = [7, 4, 1, 3, 3] + if_data[6, :] = [5, 6, 2, 4, 0] + if_data[7, :] = [7, 6, 4, 2, 1] + + # boundary attributes + # global_battr / mesh_idx / comp_battr + bdr_data = [] + + # applied loads and displacements + bdr_data += [[1, 0, 1]] + bdr_data += [[1, 1, 1]] + bdr_data += [[1, 2, 1]] + bdr_data += [[2, 2, 2]] + bdr_data += [[2, 4, 2]] + bdr_data += [[2, 7, 2]] + bdr_data += [[3, 5, 3]] + bdr_data += [[3, 6, 3]] + bdr_data += [[3, 7, 3]] + bdr_data += [[4, 0, 4]] + bdr_data += [[4, 3, 4]] + bdr_data += [[4, 5, 4]] + + # homogenous neumann + bdr_data += [[5, 1, 3]] + bdr_data += [[5, 4, 4]] + bdr_data += [[5, 6, 1]] + bdr_data += [[5, 3, 2]] + + # boundary attributes + bdr_data0 = np.array(bdr_data) + + for i in range(4): # Loop over all the variants + bdr_data = PermuteBdrData(bdr_data0, i) + filename = "linelast.comp_train" + str(i) + ".h5" + with h5py.File(filename, 'w') as f: + # c++ currently cannot read datasets of string. + # change to multiple attributes, only as a temporary implementation. + grp = f.create_group("components") + grp.attrs["number_of_components"] = 3 + grp.attrs["0"] = "joint2D" + grp.attrs["1"] = "rod2D_H" + grp.attrs["2"] = "rod2D_V" + # component index of each mesh + grp.create_dataset("meshes", (n_mesh,), data=mesh_type) + # 3-dimension vector for translation / rotation + grp.create_dataset("configuration", mesh_configs.shape, data=mesh_configs) + + grp = f.create_group("ports") + grp.attrs["number_of_references"] = 4 + grp.attrs["0"] = "port1" + grp.attrs["1"] = "port2" + grp.attrs["2"] = "port3" + grp.attrs["3"] = "port4" + grp.create_dataset("interface", if_data.shape, data=if_data) + + port = grp.create_group("port1") + port.attrs["comp1"] = "joint2D" + port.attrs["comp2"] = "rod2D_H" + port.attrs["attr1"] = 2 + port.attrs["attr2"] = 4 + port.create_dataset("comp2_configuration", (6,), data=[1., 0., 0., 0., 0., 0.]) + + port = grp.create_group("port2") + port.attrs["comp1"] = "joint2D" + port.attrs["comp2"] = "rod2D_H" + port.attrs["attr1"] = 4 + port.attrs["attr2"] = 2 + port.create_dataset("comp2_configuration", (6,), data=[-4., 0., 0., 0., 0., 0.]) + + port = grp.create_group("port3") + port.attrs["comp1"] = "joint2D" + port.attrs["comp2"] = "rod2D_V" + port.attrs["attr1"] = 3 + port.attrs["attr2"] = 1 + port.create_dataset("comp2_configuration", (6,), data=[0., 1., 0., 0., 0., 0.]) + + port = grp.create_group("port4") + port.attrs["comp1"] = "joint2D" + port.attrs["comp2"] = "rod2D_V" + port.attrs["attr1"] = 1 + port.attrs["attr2"] = 3 + port.create_dataset("comp2_configuration", (6,), data=[0., -4., 0., 0., 0., 0.]) + + # boundary attributes + f.create_dataset("boundary", bdr_data.shape, data=bdr_data) + return + +battr_map = np.array( +[ + [1, 2 ,3, 4], + [4, 1, 2 ,3], + [3, 4, 1, 2], + [2, 3, 4, 1] +] +) + +def PermuteBdrData(bdr_data0, j): + bdr_data = bdr_data0.copy() + for i in range(bdr_data.shape[0]): + battr = bdr_data[i, 0] + if battr < 5: + bdr_data[i, 0] = battr_map[j, battr-1] + return bdr_data + + +def CWTrain2x2(): + nx = 2 + ny = 2 + l = 4.0 + w = 1.0 + mesh_configs, bdr_data, if_data, mesh_type, n_mesh = cwgrid_mesh_configs(nx, ny, l, w) + + # interface data + if_data = np.array(if_data) + + # boundary attributes + bdr_data0 = np.array(bdr_data) + + for i in range(4): # Loop over all the variants + bdr_data = PermuteBdrData(bdr_data0, i) + filename = "linelast.comp_train" + str(i) + ".h5" + with h5py.File(filename, 'w') as f: + # c++ currently cannot read datasets of string. + # change to multiple attributes, only as a temporary implementation. + grp = f.create_group("components") + grp.attrs["number_of_components"] = 3 + grp.attrs["0"] = "joint2D" + grp.attrs["1"] = "rod2D_H" + grp.attrs["2"] = "rod2D_V" + # component index of each mesh + grp.create_dataset("meshes", (n_mesh,), data=mesh_type) + # 3-dimension vector for translation / rotation + grp.create_dataset("configuration", mesh_configs.shape, data=mesh_configs) + + grp = f.create_group("ports") + grp.attrs["number_of_references"] = 4 + grp.attrs["0"] = "port1" + grp.attrs["1"] = "port2" + grp.attrs["2"] = "port3" + grp.attrs["3"] = "port4" + grp.create_dataset("interface", if_data.shape, data=if_data) + + port = grp.create_group("port1") + port.attrs["comp1"] = "joint2D" + port.attrs["comp2"] = "rod2D_H" + port.attrs["attr1"] = 2 + port.attrs["attr2"] = 4 + port.create_dataset("comp2_configuration", (6,), data=[w, 0., 0., 0., 0., 0.]) + + port = grp.create_group("port2") + port.attrs["comp1"] = "joint2D" + port.attrs["comp2"] = "rod2D_H" + port.attrs["attr1"] = 4 + port.attrs["attr2"] = 2 + port.create_dataset("comp2_configuration", (6,), data=[-l, 0., 0., 0., 0., 0.]) + + port = grp.create_group("port3") + port.attrs["comp1"] = "joint2D" + port.attrs["comp2"] = "rod2D_V" + port.attrs["attr1"] = 3 + port.attrs["attr2"] = 1 + port.create_dataset("comp2_configuration", (6,), data=[0., w, 0., 0., 0., 0.]) + + port = grp.create_group("port4") + port.attrs["comp1"] = "joint2D" + port.attrs["comp2"] = "rod2D_V" + port.attrs["attr1"] = 1 + port.attrs["attr2"] = 3 + port.create_dataset("comp2_configuration", (6,), data=[0., -l, 0., 0., 0., 0.]) + + # boundary attributes + f.create_dataset("boundary", bdr_data.shape, data=bdr_data) + return + +def get_file_data(filename): + with h5py.File(filename, 'r') as f: + fom_t = f['fom_solve'][0] + rom_t = f['rom_solve'][0] + speedup = fom_t / rom_t + relerr = f['rel_error'][0] + return fom_t, rom_t, relerr, speedup + +def get_results(samples, prefix): + res=[] + + for name in samples: + filename = os.path.join(prefix, 'comparison' + str(name) + '.h5') + res.append(get_file_data(filename)) + + return np.array(res) + +def create_scaling_plot(samples, res, scale_prefix, plt_name = "plot.png"): + fig, axs = plt.subplots(1,2, figsize=(15,5)) + + axs[0].plot(samples, res[:,0], label='FOM solve time') + axs[0].plot(samples, res[:,1], label='ROM solve time') + axs[0].plot(samples, res[:,2], label='Relative error') + axs[0].set_xlabel(scale_prefix) + axs[0].set_yscale('log') + axs[0].legend() + + + axs[1].plot(samples, res[:,3], label='Speedup factor', ) + axs[1].set_xlabel(scale_prefix) + axs[1].legend() + + fig.suptitle(scale_prefix+'-scaling', fontsize = 24) + plt.tight_layout() + plt.savefig(plt_name) + +def basis_scaling_plot(prefix, plot_path): + scale_prefix = '$n_{basis}$' + samples = (2, 4, 6, 8, 10, 12, 16, 20, 24, 28, 32, 40, 48, 56, 64, 72, 80, 88, 96, 104, 112, 120, 128) + res = get_results(samples, prefix) + create_scaling_plot(samples, res, scale_prefix, plot_path) + +if __name__ == "__main__": + import sys + if len(sys.argv) == 1: + print("Please specify a component to generate") + else: + name = sys.argv[1] + if name == "cwtrain_1x1": + CWTrain1x1() + elif name == "cwtrain_2x2": + CWTrain2x2() + elif name == "bscale": + prefix = sys.argv[2] + plot_path = sys.argv[3] + basis_scaling_plot(prefix, plot_path) + elif name == "cwtrain_mesh": + prefix = sys.argv[2] + create_training_meshes(prefix) + else: + print("Option not found")