diff --git a/Parallel_Vibro-Acoustic_Solver/CMakeLists.txt b/Parallel_Vibro-Acoustic_Solver/CMakeLists.txt new file mode 100644 index 00000000..3d1a9143 --- /dev/null +++ b/Parallel_Vibro-Acoustic_Solver/CMakeLists.txt @@ -0,0 +1,53 @@ +# Set the name of the project and target: +SET(TARGET "parallel_vibroacoustic_solver") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 3.13.4) + +FIND_PACKAGE(deal.II 9.6.0 + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +# +# Are all dependencies fulfilled? +# +if(NOT (DEAL_II_WITH_PETSC AND DEAL_II_PETSC_WITH_COMPLEX AND DEAL_II_WITH_P4EST)) + message(FATAL_ERROR " +Error! This programm requires a deal.II library that was configured with the following options: + DEAL_II_WITH_PETSC = ON + DEAL_II_PETSC_WITH_COMPLEX = ON + DEAL_II_WITH_P4EST = ON +However, the deal.II library found at ${DEAL_II_PATH} was configured with these options: + DEAL_II_WITH_PETSC = ${DEAL_II_WITH_PETSC} + DEAL_II_PETSC_WITH_COMPLEX = ${DEAL_II_PETSC_WITH_COMPLEX} + DEAL_II_WITH_P4EST = ${DEAL_II_WITH_P4EST} +This conflicts with the requirements. +One or both of the aforementioned combinations of prerequisites are not met by your installation, but at least one is required for this programm." + ) +endif() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/Parallel_Vibro-Acoustic_Solver/Readme.md b/Parallel_Vibro-Acoustic_Solver/Readme.md new file mode 100644 index 00000000..d8eb48f4 --- /dev/null +++ b/Parallel_Vibro-Acoustic_Solver/Readme.md @@ -0,0 +1,181 @@ +Overview +-------- + +The Sound Transmission Loss (STL) is a measure of the sound insulation of components and is therefore an important acoustic quantity. Here, the STL of a concrete wall is calculated by means of vibroacoustic FEM simulations. In the frequency domain, the elastodynamic equations are solved in the structural domain, while the Helmholtz equation is solved in the acoustic domain, with full coupling at the interface. The deal.II implementation is parallelized using MPI and MUMPS, employs the hp finite element method, and is mainly based on step-62 step-46, and step-27. + +Elastodynamic equations +----------------------- + +The behavior of the structure $ \Omega_S $ is modelled by means of the linear elastodynamic equations, which are given in strong form as +@f[ +\rho_S \ddot{u}_i-\left( c_{ijkl} \epsilon_{kl} \right)_{,j} = f_i \qquad \text{in } \Omega_S \\ +\qquad \qquad \qquad u_i = \bar{u}_i \qquad \text{on } \partial \Omega_{SD} \\ +\qquad \qquad \left( c_{ijkl} \epsilon_{kl} \right) n_j= \bar{t}_i \qquad \text{on } \partial \Omega_{SN} \\ +@f] + +Here, $u_i $ denotes the displacement, $ \rho $ is the mass density, $ f_i $ is the body force density, $ c_{ijkl} $ is the stiffness tensor, $ n_i $ is the outwards pointing normal unit vector, $ \bar{t}_i $ is the prescribed traction and $ \bar{u}_i $ is the prescribed displacement and the double dot represents the second derivative with respect to time. The linear strain tensor $ \epsilon_{ij} $ is defined as +@f[ \epsilon_{ij} = \frac{1}{2}\left( u_{i,j} + u_{j,i} \right) @f] +Multiplying by the test function $ \phi_i $, integration over the domain, applying integration by parts and the divergence theorem and introducing the boundary conditions results in the weak form +@f{align*}{ +\int_{\Omega_S} \phi_j \rho_S \ddot{u}_j dV + \int_{\Omega_S} \phi_{j,i} c_{ijkl} \epsilon_{kl} dV - \int_{\partial \Omega_{SN} } \phi_j \bar{t}_j = 0 +@f} + +The time domain equation above is transferred into the frequency domain by performing a Fourier transform with regard to the time variable. Then, the weak form becomes +@f{align*}{ +-\int_{\Omega_S} \phi_j \omega^2 \rho_S {u}_j dV + \left[ 1 + i\mathcal{I} \right] \int_{\Omega_S} \phi_{j,i} c_{ijkl} \epsilon_{kl} dV - \int_{\partial \Omega_{SN} } \phi_j \bar{t}_j = 0 +@f} +whereby $ \omega $ is the angular frequency. Here, the isotropic loss factor $ \mathcal{I} $ is introduced to account for structural damping and $ i $ is the imaginary unit defined as $ i^2=-1 $. + +Acoustic equations +------------------ + +The acoustic equation considered is the wave equation for pressure, given in strong form as +@f[ +\frac{1}{c^2} \ddot{p} - p_{,ii} = 0 \qquad \text{in } \Omega_A \\ +\qquad \qquad p = \bar{p} \qquad \text{on } \partial \Omega_{AD} \\ +\qquad \qquad p_{,i} n_i = -\rho_A \dot{\bar{v}}_n \qquad \text{on } \partial \Omega_{AN} \\ +@f] +Here, $ c $ is the sound speed and $ p $ is the acoustic pressure. As boundary conditions the prescribed normal velocity $ \bar{v}_n $ and the prescribed acoustic pressure $\bar p $ are considered. +Similar as for the elastodynamic equations, the weak form results as +@f{align*}{ +\int_{\Omega_A} \eta \frac{1}{c^2} \ddot{p} dV + \int_{\Omega_A} \eta_{,i} p_{,i} dV - \int_{\partial \Omega_{AN}} \eta \rho_A \dot{\bar{v}}_n dA = 0 +@f} +whereby $ \eta $ denotes the test function. +Performing a Fourier transform with regard to the time variable results in the frequency domain weak form, given with the wave number $ k = \omega / c $ as +@f{align*}{ +-\int_{\Omega_A} \eta k^2 {p} dV + \int_{\Omega_A} \eta_{,i} p_{,i} dV - \int_{\partial \Omega_{AN}} \eta \rho_A i \omega {\bar{v}}_n dA = 0 +@f} + +A perfectly matched layer (PML) is used to simulate an open boundary. A PML is a complex coordinate stretch in the form of +@f{align*}{ +x \to x + \frac{i}{\omega} \int^{x} \sigma(x') dx' +@f} +where $ \sigma $ is the PML absorption function. It vanishes outside of the PML, where no absorption is desired but is turned on inside the PML resulting in an absorption of outgoing waves while minimizing reflections back into the domain. Considering the PML, the weak form becomes +@f[ +-\int_{\Omega_A} \eta k^2 {p} J dV + \int_{\Omega_A} \eta_{,i} p_{,i} \Lambda_i^2 J dV - \int_{\partial \Omega_{AN}} \eta \rho_A i \omega {\bar{v}}_n \frac{J}{\Lambda_n} dA = 0 \\ +J = \prod_{i=1}^{\mathrm{dim}} s_i \\ +\Lambda_i = \frac{1}{s_i} \\ +s_i = 1+\frac{i}{\omega} \sigma_i +@f] +In the equations above the [transformation of both volume and surface elements](https://en.wikipedia.org/wiki/Finite_strain_theory#:~:text=value%20decomposition.-,Transformation%20of%20a%20surface%20and%20volume%20element,-%5Bedit%5D) is taken into account. + +Coupled Mechanic-Acoustic equations +----------------------------------- + +For the structure-acoustic coupling, the normal vector $ n_{i} $ is chosen to conform with the normal vector of the mechanical domain given as +@f[ +n_{i} = n_{i_S} = -n_{i_A} +@f] +Here, the normal vector $ n_{i_S} $ referes to the mechanical domain, while $ n_{i_A} $ denotes the normal vector of the acoustic domain. + +The coupling between the mechanic and the acoustic domain takes place at the interface $ \partial \Omega_{SA} $ and is modelled via boundary conditions, specified as +@f[ +-p n_{i} = \bar{t}_i \qquad \text{on } \partial \Omega_{SA} \\ +p_{,i} n_i = \omega^2 \rho_A u_i n_i \qquad \text{on } \partial \Omega_{SA} +@f] + +Then, the coupled system becomes +@f[ +-\int_{\Omega_S} \phi_j \omega^2 \rho_S {u}_j dV + \left[ 1 + i\mathcal{I} \right] \int_{\Omega_S} \phi_{j,i} c_{ijkl} \epsilon_{kl} dV - \int_{\partial \Omega_{SN} } \phi_j \bar{t}_j \\ +- \int_{\Omega_A} \eta k^2 {p} J dV + \int_{\Omega_A} \eta_{,i} p_{,i} \Lambda_i^2 J dV - \int_{\partial \Omega_{AN}} \eta \rho_A i \omega {\bar{v}}_n \frac{J}{\Lambda_n} dA \\ ++ \int_{\partial \Omega_{SA}} \phi_i p n_{i} dA + \int_{\partial \Omega_{SA}} \eta \rho_A \omega^2 u_i n_{i} dA = 0 +@f] +The coupled system above is implemented in the code. + +Problem description and modeling +-------------------------------- + +The STL of a concrete wall is calculated. Conceptually, several methods exist for calculating the STL, including modeling the source and receiver rooms. Here, however, in the interest of computational efficiency the source room is not modeled and the wall is directly excited by a prescribed diffuse sound field with sound power $ P_s$ introducing structural vibrations. These vibrations radiate airborne sound with sound power $ P_r$ into the receiver side, which is considered anechoic and modeled using PMLs, as illustrated in the figure below.
+The same problem and modeling setup as in [1] and [2] are considered here. + +![image](./doc/images/problem.png) + +The geometric dimension of the wall are: height $H = 4.37\,{m}$, width $W = 2.84\,\text{m}$ and thickness $T = 0.203\,\text{m}$. The density of the concrete is $ \rho_S = 2275\,\frac{kg}{m^3} $, its Young’s modulus is $ E = 31.6\,GPa $ and the Poisson’s ratio is $\nu = 0.2$. As damping, an isotropic loss factor of $\mathcal{I} = 0.01$ is assumed. The outer boundary of the wall is fixed, and the surrounding wall is considered sound-hard and therefore does not influence the STL. The density of the air is $ \rho_A = 1.204\,\frac{kg}{m^3} $ and the speed of sound in air is $ c = 343\,\frac{m}{s} $. + + +The diffuse sound pressure field $ p_{diffuse,wall} $ acting on the source side of the wall surface must be admissible, that is, it must satisfy the wave equation. Assuming $ N $ uncorrelated plane waves traveling in $ x$-direction, an admissible diffuse sound pressure field $ p_{diffuse,wall} $ on the wall surface can be obtained as +@f[ +p_{diffuse,wall} = \underbrace{ \frac{A}{\sqrt{2N}} \sum^{N}_{n=1} e^{-i\left( k_{n,x}x + k_{n,y}y + k_{n,z}z \right)} e^{i \Phi_n}}_{p_{incident}} + \underbrace{\frac{A}{\sqrt{2N}} \sum^{N}_{n=1} e^{-i\left( -k_{n,x}x + k_{n,y}y + k_{n,z}z \right)} e^{i \Phi_n}}_{p_{reflection}} \\ +k_{n,x} = \text{cos}(\Theta_n) \\ +k_{n,y} = \text{sin}(\Theta_n) \text{cos}(\phi_n) \\ +k_{n,z} = \text{sin}(\Theta_n) \text{sin}(\phi_n) +@f] +Here, $ A $ is the amplitude of the plane waves, the polar angles $\Theta_n$, $\phi_n$ and the phase angle $\Phi_n$ are uniformly distributed random variables in the interval $ \left[0,2\pi \right] $, the polar angle $ \Theta_n $ follows from $ \Theta_n = \text{acos}(q_n)$, whereby $ q_n$ is a uniformly distributed variable in the interval $ \left[0,1 \right] $. +On the wall surface the diffuse sound pressure field $ p_{diffuse,wall} $ is a superposition of the incoming waves $ p_{incident}$ and their reflections $ p_{reflection}$. In the equations above, a perfectly reflecting wall is assumed, so no reduction of amplitude or phase shift occurs. To overcome this simplification, the source room would need to be modeled, which is not done here. + +The STL is defined as +@f[ +STL = 10 \text{log}_{10} \left(\frac{P_s}{P_r} \right) +@f] +and therefore requires the calculation of the sound powers on the source $ P_s$ and receiver $ P_r$ sides. + +The sound power on the source side $ P_s$ can be calculated as +@f[ +P_s = \frac{1}{2}\int_{\partial A_{wall,source} } \text{Re}\left[ p_{incident} v_{i}^{*} n_i \right] dA +@f] +Here, the asterix ($ {}^*$) indicates complex conjugation, $ A_{wall,source} $ is the surface of the wall on the source side and $ v_i $ is the sound particle velocity, which can be calculated as +@f[ +v_{i} = \frac{-1}{i\omega \rho_a}p_{incident,i} +@f] +The sound power on the receiver side $ P_r$ can be calculated as +@f[ +P_r = \frac{1}{2}\int_{\partial A_{wall,receiver} } \text{Re}\left[ p \left(i\omega u_i\right)^{*} n_i \right] dA +@f] +Here, $ A_{wall,receiver} $ is the surface of the wall on the receiver side. + +Results +------- + +In the figure below, the STL obtained in the present work is compared with measurement data from [3], with results from a commercial FEM software reported in [1] for a single random seed in the diffuse sound field, and with the corresponding average STL over multiple random seeds as presented in [2]. Overall, the STL calculations are sensitive to the choice of random seeds for the diffuse sound field. The STL ranges between 30 and 60 dB, which corresponds to the incident sound power being between $10^3$ and $10^6$ times greater than the transmitted sound power, depending on frequency. + +![image](./doc/images/plot.png) + +In addition, the figure below shows the magnitude of the solution for the first resonance frequency at around 112 Hz. On the left, the magnitude of the displacement $u_i$ is presented, while the magnitude of the sound pressure $p$ directly behind the concrete wall is illustrated on the right side. As enforced by the PML, the sound pressure $p$ vanishes at the boundaries of the domain. + +![image](./doc/images/solution_magnitude_112_Hz.png) + + +Ideas for extensions +-------------------- + +* Implementation of a QuadratureCache class to store calculated data for single cells to be reused for different frequencies, similar to that in step-62. +* Application of an iterative solver. +* Consideration of non-matching grids between structural and acoustic domains. + +References +---------- + +[1] COMSOL Multiphysics 6.4, "Sound Transmission Loss Through a Concrete Wall", available at [https://www.comsol.com/model/sound-transmission-loss-through-a-concrete-wall-73371](https://www.comsol.com/model/sound-transmission-loss-through-a-concrete-wall-73371). + +[2] M.H. Jensen, "Modeling Sound Transmission Loss Through a Concrete Wall", COMSOL Blog, 2020, available at [https://www.comsol.com/blogs/modeling-sound-transmission-loss-through-a-concrete-wall](https://www.comsol.com/blogs/modeling-sound-transmission-loss-through-a-concrete-wall). + +[3] A. Litvin and H.W. Belliston, “Sound Transmission Loss Through Concrete and +Concrete Masonry Walls,” American Concrete Institute, Journal Proceedings, vol. 45, +pp. 641–646, 1978. + +Compilation and running the code +--------------- + +To run the code, deal.II must be installed together with PETSc (configured with complex number support) and the p4est library. + +Similar to the deal.ii examples, run +``` +cmake -DDEAL_II_DIR=/path/to/deal.II . +``` +in this directory to configure the project. +You can switch between debug and release mode by calling either +``` +make debug +``` +or +``` +make release +``` +To execute the program in serial run +``` +make run +``` +and for parallel execution (in this case, on `j` processors) run +``` +mpirun -np j ./parallel_vibroacoustic_solver +``` \ No newline at end of file diff --git a/Parallel_Vibro-Acoustic_Solver/doc/author b/Parallel_Vibro-Acoustic_Solver/doc/author new file mode 100644 index 00000000..416ba11f --- /dev/null +++ b/Parallel_Vibro-Acoustic_Solver/doc/author @@ -0,0 +1 @@ +Andreas Hegendörfer diff --git a/Parallel_Vibro-Acoustic_Solver/doc/builds-on b/Parallel_Vibro-Acoustic_Solver/doc/builds-on new file mode 100644 index 00000000..f2c2c064 --- /dev/null +++ b/Parallel_Vibro-Acoustic_Solver/doc/builds-on @@ -0,0 +1 @@ +step-62 step-46 step-27 diff --git a/Parallel_Vibro-Acoustic_Solver/doc/dependencies b/Parallel_Vibro-Acoustic_Solver/doc/dependencies new file mode 100644 index 00000000..40c2860b --- /dev/null +++ b/Parallel_Vibro-Acoustic_Solver/doc/dependencies @@ -0,0 +1 @@ +DEAL_II_WITH_PETSC DEAL_II_WITH_MPI DEAL_II_WITH_P4EST DEAL_II_WITH_COMPLEX_VALUES DEAL_II_PETSC_WITH_COMPLEX \ No newline at end of file diff --git a/Parallel_Vibro-Acoustic_Solver/doc/entry-name b/Parallel_Vibro-Acoustic_Solver/doc/entry-name new file mode 100644 index 00000000..dae29cb1 --- /dev/null +++ b/Parallel_Vibro-Acoustic_Solver/doc/entry-name @@ -0,0 +1 @@ +Parallel Vibroacoustic Solver diff --git a/Parallel_Vibro-Acoustic_Solver/doc/images/plot.png b/Parallel_Vibro-Acoustic_Solver/doc/images/plot.png new file mode 100644 index 00000000..1245bc6d Binary files /dev/null and b/Parallel_Vibro-Acoustic_Solver/doc/images/plot.png differ diff --git a/Parallel_Vibro-Acoustic_Solver/doc/images/problem.png b/Parallel_Vibro-Acoustic_Solver/doc/images/problem.png new file mode 100644 index 00000000..28034b03 Binary files /dev/null and b/Parallel_Vibro-Acoustic_Solver/doc/images/problem.png differ diff --git a/Parallel_Vibro-Acoustic_Solver/doc/images/solution_magnitude_112_Hz.png b/Parallel_Vibro-Acoustic_Solver/doc/images/solution_magnitude_112_Hz.png new file mode 100644 index 00000000..c0523ed0 Binary files /dev/null and b/Parallel_Vibro-Acoustic_Solver/doc/images/solution_magnitude_112_Hz.png differ diff --git a/Parallel_Vibro-Acoustic_Solver/doc/tooltip b/Parallel_Vibro-Acoustic_Solver/doc/tooltip new file mode 100644 index 00000000..32d66c2e --- /dev/null +++ b/Parallel_Vibro-Acoustic_Solver/doc/tooltip @@ -0,0 +1 @@ +Implementation of a parallel frequency-domain vibroacoustic solver calculating the sound transmission loss of a concrete wall. diff --git a/Parallel_Vibro-Acoustic_Solver/parallel_vibroacoustic_solver.cc b/Parallel_Vibro-Acoustic_Solver/parallel_vibroacoustic_solver.cc new file mode 100644 index 00000000..5bf7d4fe --- /dev/null +++ b/Parallel_Vibro-Acoustic_Solver/parallel_vibroacoustic_solver.cc @@ -0,0 +1,1357 @@ +/* ----------------------------------------------------------------------------- + * + * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later + * Copyright (C) 2026 by Andreas Hegendörfer + * + * This file is part of the deal.II code gallery. + * + * ----------------------------------------------------------------------------- + */ + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +namespace VibroAcousticProblem +{ + using namespace dealii; + // geometric tolerance + constexpr double geom_tol = 1e-10; + + enum class SurfaceID : dealii::types::boundary_id + { + // By default, boundary indicators are 0 + Default, + SourceSide, + ReceiverSide, + FixedBoundary, + ZeroPressure + }; + + enum class MaterialID : dealii::types::material_id + { + Concrete, + Air + }; + + // return the isotropic loss factor + std::complex + get_iso_loss() + { + return std::complex{1., 0.01}; + } + + // calculation of the linear strain tensor + template + inline SymmetricTensor<2, dim> + get_strain(const FEValues &fe_values, + const unsigned int shape_func, + const unsigned int q_point) + { + SymmetricTensor<2, dim> tmp; + + for (unsigned int i = 0; i < dim; ++i) + tmp[i][i] = fe_values.shape_grad_component(shape_func, q_point, i)[i]; + + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = i + 1; j < dim; ++j) + tmp[i][j] = + (fe_values.shape_grad_component(shape_func, q_point, i)[j] + + fe_values.shape_grad_component(shape_func, q_point, j)[i]) / + 2; + + return tmp; + } + + // Returning the stiffness tensor of the wall + template + SymmetricTensor<4, dim> + get_stiffness_tensor() + { + const double E = 31600. * 1e6; + const double v = 0.2; + double lambda = v / (1 - 2 * v) * 1 / (1 + v) * E; + double mu = 0.5 * 1 / (1 + v) * E; + + SymmetricTensor<4, dim> stiffness_tensor; + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = 0; l < dim; ++l) + stiffness_tensor[i][j][k][l] = + (((i == k) && (j == l) ? mu : 0.0) + + ((i == l) && (j == k) ? mu : 0.0) + + ((i == j) && (k == l) ? lambda : 0.0)); + return stiffness_tensor; + } + + // Returning the density of the wall + double + get_density_structure() + { + return 2.275 * 1e-3; + } + + // Returning the density of air + double + get_density_air() + { + return 1.204 * 1e-6; + } + + // Returning the speed of sound + double + get_sound_speed() + { + return 343. * 1e3; + } + + // Implementation of a perfectly matched layer + template + class PML : public Function> + { + public: + explicit PML(double omega) + : omega(omega) + , b_pos{} + , b_neg{} + , t_pos{} + , t_neg{} + , pml_coeff_degree(2.0) + , pml_coeff(1.e4 / omega) + { + Assert(dim == 3, ExcNotImplemented()); + + // x-direction + b_pos[0] = 812.0; + b_neg[0] = std::numeric_limits::lowest(); + t_pos[0] = 500.0; + t_neg[0] = 500.0; + + // y-direction + b_pos[1] = 1826.0; + b_neg[1] = -1826.0; + t_pos[1] = 500.0; + t_neg[1] = 500.0; + + // z-direction + b_pos[2] = 2591.0; + b_neg[2] = -2591.0; + t_pos[2] = 500.0; + t_neg[2] = 500.0; + } + + // calculates the value of the complex coordinate stretch + void + vector_value(const Point &p, + Vector> &value) const override; + + private: + double omega; + + std::array b_pos, b_neg; + std::array t_pos, t_neg; + + double pml_coeff_degree; + double pml_coeff; + }; + + template + void + PML::vector_value(const Point &p, + Vector> &value) const + { + value.reinit(dim); + + for (unsigned int d = 0; d < dim; ++d) + { + double coeff = 0.0; + + if (p[d] > b_pos[d]) + { + const double x_prime = p[d] - b_pos[d]; + const double a_coeff = + pml_coeff / std::pow(t_pos[d], pml_coeff_degree); + coeff = a_coeff * std::pow(x_prime, pml_coeff_degree); + } + else if (p[d] < b_neg[d]) + { + const double x_prime = b_neg[d] - p[d]; + const double a_coeff = + pml_coeff / std::pow(t_neg[d], pml_coeff_degree); + coeff = a_coeff * std::pow(x_prime, pml_coeff_degree); + } + + // complex coordinate stretching: s = 1 + i * sigma(x) + value[d] = std::complex(1.0, coeff); + } + } + + // Creation of a diffuse sound field for excitation of the wall on the source + // side + template + class DiffuseSoundField : public Function> + { + public: + DiffuseSoundField(unsigned int N, double omega, MPI_Comm mpi_communicator) + : Function>() + , N(N) + , dist_0_2PI(0.0, 2. * numbers::PI) + , dist_0_1(0.0, 1.) + , generator(static_cast(omega)) + , omega(omega) + { + // Only let one rank 0 create the random variables of the diffuse sound + // field... + const unsigned int rank = + Utilities::MPI::this_mpi_process(mpi_communicator); + if (rank == 0) + { + phi.resize(N); + Phi.resize(N); + phase.resize(N); + for (unsigned int n = 0; n < N; ++n) + { + phi[n] = dist_0_2PI(generator); + Phi[n] = std::acos(dist_0_1(generator)); + phase[n] = dist_0_2PI(generator); + } + } + // ... and broadcast from rank 0 to all ranks. + phi = Utilities::MPI::broadcast(mpi_communicator, phi, 0); + Phi = Utilities::MPI::broadcast(mpi_communicator, Phi, 0); + phase = Utilities::MPI::broadcast(mpi_communicator, phase, 0); + + // Generation of kn on all ranks + for (unsigned int n = 0; n < N; n++) + { + double scale = (omega / get_sound_speed()); + Tensor<1, dim> k; + k[0] = std::cos(Phi[n]) * scale; + k[1] = std::sin(Phi[n]) * std::cos(phi[n]) * scale; + k[2] = std::sin(Phi[n]) * std::sin(phi[n]) * scale; + k_vector.push_back(k); + } + } + // Calculation of the total pressure on the source side + virtual void + value_list(const std::vector> &points, + std::vector> &values, + const unsigned int component = 0) const override; + + // Calculation of the gradient of the total pressure on the source side + virtual void + gradient_list(const std::vector> &points, + std::vector>> &gradients, + const unsigned int component = 0) const override; + + // Calculation of the incident sound pressure on the source side + void + value_list_incidence(const std::vector> &points, + std::vector> &values, + const unsigned int component = 0) const; + + private: + // Number of plane waves + unsigned int N; + // For randmoness + std::uniform_real_distribution dist_0_2PI; + std::uniform_real_distribution dist_0_1; + std::mt19937 generator; + + // variables of the diffuse field + std::vector phi, Phi, phase, kn_x, kn_y, kn_z; + std::vector> k_vector; + + // angular velocity + double omega; + }; + + template + void + DiffuseSoundField::value_list(const std::vector> &points, + std::vector> &values, + const unsigned int component) const + { + AssertThrow(component == 0, ExcMessage("only component 0 is implemented")); + values.resize(points.size()); + std::complex j{0., 1.}; + + for (unsigned int q = 0; q < points.size(); ++q) + { + const auto &q_p = points[q]; + values[q] = {0., 0.}; + for (unsigned int n = 0; n < N; n++) + { + // sound waves traveling towards the wall... + double dot_towards = 0.; + double dot_reflection = 0.; + + for (unsigned int d = 0; d < dim; d++) + { + // sound waves traveling towards the wall... + dot_towards += k_vector[n][d] * q_p[d]; + // ... and reflections. + if (d == 0) + { + dot_reflection -= k_vector[n][d] * q_p[d]; + } + else + { + dot_reflection += k_vector[n][d] * q_p[d]; + } + } + // sound waves traveling towards the wall... + values[q] += std::exp(-j * (dot_towards) + j * phase[n]); + // ... and reflections. + values[q] += std::exp(-j * (dot_reflection) + j * phase[n]); + } + values[q] *= 1. / (std::sqrt(2. * static_cast(N))) * 1.e6; + } + } + + template + void + DiffuseSoundField::value_list_incidence( + const std::vector> &points, + std::vector> &values, + const unsigned int component) const + { + AssertThrow(component == 0, ExcMessage("only component 0 is implemented")); + + values.resize(points.size()); + std::complex j{0., 1.}; + for (unsigned int q = 0; q < points.size(); ++q) + { + const auto &q_p = points[q]; + values[q] = {0., 0.}; + for (unsigned int n = 0; n < N; n++) + { + // Only sound waves traveling towards the wall. + double dot_towards = 0.; + for (unsigned int d = 0; d < dim; d++) + { + // Sound waves traveling towards the wall. + dot_towards += k_vector[n][d] * q_p[d]; + } + // Only sound waves traveling towards the wall. + values[q] += std::exp(-j * (dot_towards) + j * phase[n]); + } + values[q] *= 1. / (std::sqrt(2. * static_cast(N))) * 1.e6; + } + } + + template + void + DiffuseSoundField::gradient_list( + const std::vector> &points, + std::vector>> &gradients, + const unsigned int component) const + { + AssertThrow(component == 0, ExcMessage("only component 0 is implemented")); + + gradients.resize(points.size()); + const std::complex j{0.0, 1.0}; + for (unsigned int q = 0; q < points.size(); ++q) + { + const auto &q_p = points[q]; + for (unsigned int d = 0; d < dim; ++d) + gradients[q][d] = 0.; + + for (unsigned int n = 0; n < N; ++n) + { + double dot = 0.; + for (unsigned int d = 0; d < dim; ++d) + dot += k_vector[n][d] * q_p[d]; + + const std::complex exp_term = + std::exp(-j * dot + j * phase[n]); + + for (unsigned int d = 0; d < dim; ++d) + gradients[q][d] += (-j * k_vector[n][d]) * exp_term * + (1. / std::sqrt(2. * static_cast(N))) * + 1.e6; + } + } + } + + template + class HarmonicResponse + { + public: + HarmonicResponse(double omega); + void + run(bool write_output = false); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + + // calculate the magnitude of u and p + void + calculate_magnitude(); + + // calculation of sound power at receiver side for one cell + double + cell_receiver_sound_power( + const FEFaceValuesBase &elasticity_fe_face_values, + const FEFaceValuesBase &air_fe_face_values, + const double &omega); + + // calculation of sound power at sender side + double + incident_sound_power(); + // calculation of sound power at receiver side + double + receiver_sound_power(); + // Assemble the air-structure coupling terms + // This function implements the assembly of the air-structure interface. + void + assemble_air_structure_interface_term( + const FEFaceValuesBase &elasticity_fe_face_values, + const FEFaceValuesBase &air_fe_face_values, + FullMatrix> &local_interface_matrix, + const double &omega); + MPI_Comm mpi_communicator; + parallel::distributed::Triangulation triangulation; + + const FESystem fe_structure, fe_air; + const QGauss quadrature_formula_structure, quadrature_formula_air; + const QGauss face_quadrature_formula_structure, + face_quadrature_formula_air; + + hp::FECollection fe_collection; + hp::QCollection q_collection; + hp::QCollection q_face_collection; + + DoFHandler dof_handler; + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + AffineConstraints> constraints; + LinearAlgebraPETSc::MPI::SparseMatrix system_matrix; + LinearAlgebraPETSc::MPI::Vector locally_relevant_solution; + LinearAlgebraPETSc::MPI::Vector locally_relevant_magnitude; + LinearAlgebraPETSc::MPI::Vector system_rhs; + + ConditionalOStream pcout; + const double omega; + DiffuseSoundField field; + PML pml; + }; + + template + HarmonicResponse::HarmonicResponse(double omega) + : mpi_communicator(MPI_COMM_WORLD) + , triangulation(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , fe_structure(FE_Q(2), dim, FE_Nothing(), 1) + , fe_air(FE_Nothing(), dim, FE_Q(1), 1) + , quadrature_formula_structure(fe_structure.degree + 1) + , quadrature_formula_air(fe_air.degree + 1) + , face_quadrature_formula_structure(fe_structure.degree + 1) + , face_quadrature_formula_air(fe_air.degree + 1) + , dof_handler(triangulation) + , pcout(std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + , omega(omega) + , field(1.e3, omega, mpi_communicator) + , pml(omega) + { + static_assert(dim == 3, + "HarmonicResponse is only implemented for dim == 3"); + fe_collection.push_back(fe_structure); + fe_collection.push_back(fe_air); + q_collection.push_back(quadrature_formula_structure); + q_collection.push_back(quadrature_formula_air); + q_face_collection.push_back(face_quadrature_formula_structure); + q_face_collection.push_back(face_quadrature_formula_air); + } + + template + void + HarmonicResponse::setup_system() + { + // Set material id and active FE indices. + for (const auto &cell : dof_handler.cell_iterators()) + { + cell->set_material_id(static_cast(MaterialID::Concrete)); + + if ((cell->center()[0] - 203.) > 0.) + { + cell->set_material_id(static_cast(MaterialID::Air)); + } + } + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + cell->set_active_fe_index( + static_cast(MaterialID::Concrete)); + } + if ((cell->center()[0] - 203.) > 0.) + { + if (cell->is_locally_owned()) + { + cell->set_active_fe_index( + static_cast(MaterialID::Air)); + } + } + } + // Definition of FE space + dof_handler.distribute_dofs(fe_collection); + pcout << " Number of degrees of freedom = " << dof_handler.n_dofs() + << std::endl; + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + locally_relevant_dofs.clear(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + locally_relevant_solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + locally_relevant_magnitude.reinit(locally_owned_dofs, mpi_communicator); + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + // set up contraints + constraints.clear(); + constraints.reinit(locally_owned_dofs, locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + // Set boundary ids + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + for (unsigned int f = 0; f < cell->n_faces(); ++f) + { + const auto face = cell->face(f); + const auto p = face->center(); + if (std::abs(p[0] - 203.) < geom_tol && + (cell->material_id() == + static_cast(MaterialID::Concrete))) + { + cell->face(f)->set_user_index( + static_cast(SurfaceID::ReceiverSide)); + } + if (std::abs(p[0]) < geom_tol) + { + cell->face(f)->set_user_index( + static_cast(SurfaceID::SourceSide) && + (cell->material_id() == + static_cast(MaterialID::Concrete))); + } + if (cell->face(f)->at_boundary()) + { + if ((p[0] + geom_tol) < 203. && p[0] > geom_tol) + { + cell->face(f)->set_boundary_id( + static_cast(SurfaceID::FixedBoundary)); + } + if (std::abs(p[0] - 1312.) < geom_tol && + cell->material_id() == + static_cast(MaterialID::Air)) + { + cell->face(f)->set_boundary_id( + static_cast(SurfaceID::ZeroPressure)); + } + } + } + } + } + // The wall is fixed at its outer boundary + const FEValuesExtractors::Vector displacement(0); + ComponentMask component_mask_displacement = + fe_collection.component_mask(displacement); + VectorTools::interpolate_boundary_values( + dof_handler, + static_cast(SurfaceID::FixedBoundary), + Functions::ZeroFunction>(4), + constraints, + component_mask_displacement); + constraints.close(); + + DynamicSparsityPattern dsp(locally_relevant_dofs); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, true); + DoFTools::make_flux_sparsity_pattern(dof_handler, dsp, constraints); + SparsityTools::distribute_sparsity_pattern(dsp, + locally_owned_dofs, + mpi_communicator, + locally_relevant_dofs); + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); + } + + template + void + HarmonicResponse::assemble_system() + { + // FE values for volume integration + hp::FEValues hp_fe_values(fe_collection, + q_collection, + update_values | update_gradients | + update_quadrature_points | + update_JxW_values); + hp::FEFaceValues hp_fe_face_values(fe_collection, + q_face_collection, + update_values | update_JxW_values | + update_normal_vectors | + update_quadrature_points); + + // Common face quadrature + const QGauss common_face_quadrature(fe_collection.max_degree() + + 1); + // FE face values for surface integration + FEFaceValues air_fe_face_values(fe_air, + common_face_quadrature, + update_values | update_JxW_values | + update_normal_vectors | + update_quadrature_points); + FEFaceValues elasticity_fe_face_values(fe_structure, + common_face_quadrature, + update_values | + update_JxW_values | + update_normal_vectors | + update_quadrature_points); + FESubfaceValues air_fe_sub_face_values(fe_air, + common_face_quadrature, + update_values | + update_JxW_values | + update_normal_vectors | + update_quadrature_points); + FESubfaceValues elasticity_fe_sub_face_values( + fe_structure, + common_face_quadrature, + update_values | update_JxW_values | update_normal_vectors | + update_quadrature_points); + + // Local element matrix for a cell + FullMatrix> cell_matrix; + // Local interface matrix between air and structure DoFs + FullMatrix> local_interface_matrix( + fe_air.n_dofs_per_cell(), fe_structure.n_dofs_per_cell()); + + // Right-hand side + Vector> cell_rhs; + std::vector local_dof_indices; + std::vector neighbor_dof_indices; + const FEValuesExtractors::Vector displacement(0); + const FEValuesExtractors::Scalar pressure(dim); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + // Assemble air cells contributions + if (cell->material_id() == + static_cast(MaterialID::Air)) + { + cell_matrix = 0; + cell_rhs = 0; + hp_fe_values.reinit(cell); + const FEValues &fe_values = + hp_fe_values.get_present_fe_values(); + + const unsigned int dofs_per_cell = fe_air.n_dofs_per_cell(); + const unsigned int n_q_points = fe_values.n_quadrature_points; + cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_rhs.reinit(dofs_per_cell); + local_dof_indices.resize(dofs_per_cell); + neighbor_dof_indices.resize(fe_structure.n_dofs_per_cell()); + const double sound_speed = get_sound_speed(); + + // Wave number k + const std::complex k = (omega / sound_speed); + + for (unsigned int q = 0; q < n_q_points; ++q) + { + const auto JxW = fe_values.JxW(q); + const Point &q_point = fe_values.quadrature_point(q); + // calucalte lambda and J for PML + Vector> lambda(dim); + pml.vector_value(q_point, lambda); + const std::complex J = + lambda[0] * lambda[1] * lambda[2]; + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const double phi_i = fe_values[pressure].value(i, q); + + // Gradient of phi_i, which is multiplied by 1/lambda + Tensor<1, dim, std::complex> phi_i_i = + fe_values[pressure].gradient(i, q); + for (unsigned int d = 0; d < dim; ++d) + phi_i_i[d] *= 1.0 / lambda[d]; + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + const double phi_j = + fe_values[pressure].value(j, q); + + // Gradient of phi_j, which is multiplied by + // 1/lambda + Tensor<1, dim, std::complex> phi_j_j = + fe_values[pressure].gradient(j, q); + for (unsigned int d = 0; d < dim; ++d) + phi_j_j[d] *= 1.0 / lambda[d]; + + cell_matrix[i][j] += phi_i_i * phi_j_j * JxW * J; + cell_matrix[i][j] -= Utilities::fixed_power<2>(k) * + phi_i * phi_j * JxW * J; + } + } + } + // assemble local system to global system. + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + + // Here, the air-structure interface is considered. Similar to + // step 46, 3 possibilities exist: The neighbor is + // at the same refinement level and has no children, the + // neighbor has children and the neighbor is coarser. + for (const auto f : cell->face_indices()) + { + if (!cell->at_boundary(f)) + { + const auto neighbor = cell->neighbor(f); + if (neighbor->material_id() == + static_cast(MaterialID::Concrete)) + + // if (neighbor->center()[0]<203.) + { + // The neighbor is at the same refinement level and + // has no children. + if ((cell->neighbor(f)->level() == cell->level()) && + (cell->neighbor(f)->has_children() == false)) + { + air_fe_face_values.reinit(cell, f); + elasticity_fe_face_values.reinit( + cell->neighbor(f), + cell->neighbor_of_neighbor(f)); + assemble_air_structure_interface_term( + elasticity_fe_face_values, + air_fe_face_values, + local_interface_matrix, + omega); + cell->neighbor(f)->get_dof_indices( + neighbor_dof_indices); + constraints.distribute_local_to_global( + local_interface_matrix, + local_dof_indices, + neighbor_dof_indices, + system_matrix); + } + // The neighbor has children. + else if ((cell->neighbor(f)->level() == + cell->level()) && + (cell->neighbor(f)->has_children() == + true)) + { + for (unsigned int subface = 0; + subface < cell->face(f)->n_children(); + ++subface) + { + air_fe_sub_face_values.reinit(cell, + f, + subface); + elasticity_fe_face_values.reinit( + cell->neighbor_child_on_subface(f, + subface), + cell->neighbor_of_neighbor(f)); + assemble_air_structure_interface_term( + elasticity_fe_face_values, + air_fe_sub_face_values, + local_interface_matrix, + omega); + cell->neighbor_child_on_subface(f, subface) + ->get_dof_indices(neighbor_dof_indices); + constraints.distribute_local_to_global( + local_interface_matrix, + local_dof_indices, + neighbor_dof_indices, + system_matrix); + } + } + // The neighbor is coarser. + else if (cell->neighbor_is_coarser(f)) + { + air_fe_face_values.reinit(cell, f); + elasticity_fe_sub_face_values.reinit( + cell->neighbor(f), + cell->neighbor_of_coarser_neighbor(f).first, + cell->neighbor_of_coarser_neighbor(f).second); + assemble_air_structure_interface_term( + elasticity_fe_sub_face_values, + air_fe_face_values, + local_interface_matrix, + omega); + cell->neighbor(f)->get_dof_indices( + neighbor_dof_indices); + constraints.distribute_local_to_global( + local_interface_matrix, + local_dof_indices, + neighbor_dof_indices, + system_matrix); + } + } + } + } + } + // Assemble structure cells contributions + else if (cell->is_locally_owned() && + cell->material_id() == + static_cast(MaterialID::Concrete)) + { + cell_matrix = 0; + cell_rhs = 0; + hp_fe_values.reinit(cell); + const FEValues &fe_values = + hp_fe_values.get_present_fe_values(); + + const unsigned int dofs_per_cell = + fe_structure.n_dofs_per_cell(); + const unsigned int n_q_points = fe_values.n_quadrature_points; + cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_rhs.reinit(dofs_per_cell); + local_dof_indices.resize(dofs_per_cell); + neighbor_dof_indices.resize(fe_air.n_dofs_per_cell()); + + const SymmetricTensor<4, dim> stiffness_tensor = + get_stiffness_tensor(); + const double density = get_density_structure(); + const std::complex iso_loss = get_iso_loss(); + + for (unsigned int q = 0; q < n_q_points; ++q) + { + const auto JxW = fe_values.JxW(q); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const Tensor<1, dim> phi_i = + fe_values[displacement].value(i, q); + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + const Tensor<1, dim> phi_j = + fe_values[displacement].value(j, q); + const SymmetricTensor<2, dim> eps_phi_i = + get_strain(fe_values, i, q); + const SymmetricTensor<2, dim> eps_phi_j = + get_strain(fe_values, j, q); + cell_matrix[i][j] += eps_phi_i * stiffness_tensor * + iso_loss * eps_phi_j * JxW; + cell_matrix[i][j] -= + Utilities::fixed_power<2>(omega) * phi_j * phi_i * + density * JxW; + } + } + } + // Here the diffuse sound field is considered on the source side + // of the wall. + for (const auto face_no : cell->face_indices()) + { + if (cell->face(face_no)->user_index() == + static_cast(SurfaceID::SourceSide)) + { + hp_fe_face_values.reinit(cell, face_no); + const FEFaceValues &fe_face_values = + hp_fe_face_values.get_present_fe_values(); + const unsigned int n_face_q_points = + fe_face_values.n_quadrature_points; + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + if (fe_structure.has_support_on_face(i, face_no)) + { + std::vector> pressure( + n_face_q_points); + const auto face_quadrature_points = + fe_face_values.get_quadrature_points(); + field.value_list(face_quadrature_points, + pressure); + + for (unsigned int q = 0; q < n_face_q_points; + ++q) + { + const auto JxW = fe_face_values.JxW(q); + const Tensor<1, dim> phi_i = + fe_face_values[displacement].value(i, q); + const Tensor<1, dim, double> NormalVector = + -fe_face_values.normal_vector(q); + cell_rhs[i] += pressure[q] * + (NormalVector * phi_i) * JxW; + } + } + } + } + } + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } + } + } + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + template + void + HarmonicResponse::assemble_air_structure_interface_term( + const FEFaceValuesBase &elasticity_fe_face_values, + const FEFaceValuesBase &air_fe_face_values, + FullMatrix> &local_interface_matrix, + const double &omega) + { + const FEValuesExtractors::Vector displacement(0); + const FEValuesExtractors::Scalar pressure(dim); + local_interface_matrix = 0; + const auto density_air = get_density_air(); + const unsigned int n_face_quadrature_points = + air_fe_face_values.n_quadrature_points; + + for (unsigned int q = 0; q < n_face_quadrature_points; ++q) + { + const Tensor<1, dim, double> normalVectorStructure = + -air_fe_face_values.normal_vector(q); + for (unsigned int i = 0; i < air_fe_face_values.dofs_per_cell; ++i) + { + const double phi_i = air_fe_face_values[pressure].value(i, q); + for (unsigned int j = 0; + j < elasticity_fe_face_values.dofs_per_cell; + ++j) + { + const Tensor<1, dim> phi_j = + elasticity_fe_face_values[displacement].value(j, q); + local_interface_matrix[i][j] += + density_air * Utilities::fixed_power<2>(omega) * phi_i * + normalVectorStructure * phi_j * air_fe_face_values.JxW(q); + local_interface_matrix[i][j] += + phi_j * normalVectorStructure * phi_i * + elasticity_fe_face_values.JxW(q); + } + } + } + } + + // calculation of incident sound power on the source side + template + double + HarmonicResponse::incident_sound_power() + { + FEFaceValues structure_fe_face_values( + fe_structure, + face_quadrature_formula_structure, + update_values | update_JxW_values | update_normal_vectors | + update_quadrature_points); + std::vector local_dof_indices; + const FEValuesExtractors::Vector displacement(0); + double sound_power = 0.; + const double density_air = get_density_air(); + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + for (const auto face_no : cell->face_indices()) + { + if (cell->face(face_no)->at_boundary() && + (cell->face(face_no)->user_index() == + static_cast(SurfaceID::SourceSide)) && + (cell->material_id() == + static_cast(MaterialID::Concrete))) + { + structure_fe_face_values.reinit(cell, face_no); + const unsigned int n_face_q_points = + structure_fe_face_values.n_quadrature_points; + std::vector> pressure(n_face_q_points); + const auto face_quadrature_points = + structure_fe_face_values.get_quadrature_points(); + field.value_list_incidence(face_quadrature_points, + pressure); + std::vector>> + pressure_gradients(n_face_q_points); + field.gradient_list(face_quadrature_points, + pressure_gradients); + for (unsigned int q = 0; q < n_face_q_points; ++q) + { + const Tensor<1, dim, double> NormalVector = + -structure_fe_face_values.normal_vector(q); + const auto JxW = structure_fe_face_values.JxW(q); + std::complex v_n = + (-1.0 / (omega * std::complex(0., 1.) * + density_air)) * + (pressure_gradients[q] * NormalVector); + sound_power += + 0.5 * std::real(pressure[q] * std::conj(v_n) * JxW); + } + } + } + } + } + return sound_power; + } + + // Calculating the sound power on the receiver side + template + double + HarmonicResponse::receiver_sound_power() + { + // Common face quadrature + const QGauss common_face_quadrature(fe_collection.max_degree() + + 1); + // FE face values for surface integration + FEFaceValues air_fe_face_values(fe_air, + common_face_quadrature, + update_values | update_JxW_values | + update_normal_vectors | + update_quadrature_points | + update_gradients); + FEFaceValues elasticity_fe_face_values(fe_structure, + common_face_quadrature, + update_values | + update_JxW_values | + update_normal_vectors | + update_quadrature_points | + update_gradients); + FESubfaceValues air_fe_sub_face_values(fe_air, + common_face_quadrature, + update_values | + update_JxW_values | + update_normal_vectors | + update_quadrature_points | + update_gradients); + FESubfaceValues elasticity_fe_sub_face_values( + fe_structure, + common_face_quadrature, + update_values | update_JxW_values | update_normal_vectors | + update_quadrature_points | update_gradients); + + double sound_power = 0.; + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + for (const auto f : cell->face_indices()) + { + if ((cell->face(f)->user_index() == + static_cast(SurfaceID::ReceiverSide)) && + !cell->at_boundary(f) && + (cell->material_id() == + static_cast(MaterialID::Concrete))) + { + const auto neighbor = cell->neighbor(f); + // The neighbor is at the same refinement level and has no + // children. + if ((neighbor->level() == cell->level()) && + (neighbor->has_children() == false) && + (neighbor->material_id() == + static_cast(MaterialID::Air))) + { + elasticity_fe_face_values.reinit(cell, f); + air_fe_face_values.reinit( + neighbor, cell->neighbor_of_neighbor(f)); + sound_power += + cell_receiver_sound_power(elasticity_fe_face_values, + air_fe_face_values, + omega); + } + // The neighbor has children. + else if ((neighbor->level() == cell->level()) && + (neighbor->has_children() == true)) + { + for (unsigned int subface = 0; + subface < cell->face(f)->n_children(); + ++subface) + { + elasticity_fe_sub_face_values.reinit(cell, + f, + subface); + air_fe_face_values.reinit( + cell->neighbor_child_on_subface(f, subface), + cell->neighbor_of_neighbor(f)); + sound_power += cell_receiver_sound_power( + elasticity_fe_sub_face_values, + air_fe_face_values, + omega); + } + } + // The neighbor is coarser. + else if (cell->neighbor_is_coarser(f)) + { + elasticity_fe_face_values.reinit(cell, f); + air_fe_sub_face_values.reinit( + neighbor, + cell->neighbor_of_coarser_neighbor(f).first, + cell->neighbor_of_coarser_neighbor(f).second); + sound_power += + cell_receiver_sound_power(elasticity_fe_face_values, + air_fe_sub_face_values, + omega); + } + } + } + } + } + return sound_power; + } + + template + double + HarmonicResponse::cell_receiver_sound_power( + const FEFaceValuesBase &elasticity_fe_face_values, + const FEFaceValuesBase &air_fe_face_values, + const double &omega) + { + std::vector local_dof_indices; + const FEValuesExtractors::Vector displacement(0); + const FEValuesExtractors::Scalar pressure(dim); + const unsigned int n_q_points = air_fe_face_values.n_quadrature_points; + double sound_power = 0.; + ; + + std::vector> local_dof_values_pressure( + air_fe_face_values.n_quadrature_points); + air_fe_face_values[pressure].get_function_values(locally_relevant_solution, + local_dof_values_pressure); + + std::vector>> + local_dof_values_displacement(n_q_points); + elasticity_fe_face_values[displacement].get_function_values( + locally_relevant_solution, local_dof_values_displacement); + for (unsigned int q = 0; q < n_q_points; ++q) + { + const Tensor<1, dim, double> NormalVector = + air_fe_face_values.normal_vector(q); + const auto JxW = air_fe_face_values.JxW(q); + std::complex normal_velocity = + local_dof_values_displacement[q] * NormalVector * + std::complex(0., 1.) * omega; + sound_power += + 0.5 * + std::real(local_dof_values_pressure[q] * std::conj(normal_velocity)) * + JxW; + } + return sound_power; + } + + template + void + HarmonicResponse::solve() + { + LinearAlgebraPETSc::MPI::Vector completely_distributed_solution( + locally_owned_dofs, mpi_communicator); + SolverControl solver_control; + PETScWrappers::SparseDirectMUMPS solver(solver_control, mpi_communicator); + solver.solve(system_matrix, completely_distributed_solution, system_rhs); + constraints.distribute(completely_distributed_solution); + locally_relevant_solution = completely_distributed_solution; + } + + template + void + HarmonicResponse::calculate_magnitude() + { + for (const auto i : locally_owned_dofs) + { + const std::complex value = locally_relevant_solution[i]; + // For postprocessing, the real part represents the magnitude, while the + // imaginary part vanishes. + locally_relevant_magnitude[i] = std::abs(value); + } + } + + template + void + HarmonicResponse::run(bool write_output) + { + static bool mode_output = false; + if (!mode_output) + { +#ifdef DEBUG + pcout << "Debug mode" << std::endl; +#else + pcout << "Release mode" << std::endl; +#endif + mode_output = true; + } + + const double frequency = omega / (numbers::PI * 2.); // frequency in Hz + pcout << std::endl << "Frequency = " << frequency << " Hz" << std::endl; + + GridIn grid_in; + grid_in.attach_triangulation(triangulation); + std::ifstream input_file("./tria.inp"); + grid_in.read_abaqus(input_file); + triangulation.refine_global(1); + for (const auto &cell : triangulation.active_cell_iterators()) + { + if (cell->center()[0] > 203.) + { + cell->set_refine_flag(); + } + } + triangulation.execute_coarsening_and_refinement(); + + pcout << " setup_system" << std::endl; + setup_system(); + pcout << " assemble_system" << std::endl; + assemble_system(); + pcout << " solve" << std::endl; + solve(); + // Calculation of the magnitude of the solution. For postprocessing, the + // magnitude of the solutions corresponds to the real part and the imaginary + // part vaishes. + calculate_magnitude(); + const auto sound_power_source_side_local = incident_sound_power(); + const auto sound_power_receiver_side_local = receiver_sound_power(); + const auto sound_power_source_side = + dealii::Utilities::MPI::sum(sound_power_source_side_local, + mpi_communicator); + const auto sound_power_receiver_side = + dealii::Utilities::MPI::sum(sound_power_receiver_side_local, + mpi_communicator); + + const auto stl = + 10. * std::log10(sound_power_source_side / sound_power_receiver_side); + pcout << " STL = " << stl << " dB" << std::endl; + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + std::ofstream outfile("STL_result.txt", std::ios::app); + outfile << frequency << " " << stl << std::endl; + } + + if (write_output) + { + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + + std::vector solution_names_magnitude(dim, "magnitude_u"); + solution_names_magnitude.push_back("magnitude_p"); + std::vector solution_names(dim, "u"); + solution_names.push_back("p"); + + std::vector + interpretation = { + DataComponentInterpretation::component_is_part_of_vector, + DataComponentInterpretation::component_is_part_of_vector, + DataComponentInterpretation::component_is_part_of_vector, + DataComponentInterpretation::component_is_scalar}; + data_out.add_data_vector(locally_relevant_solution, + solution_names, + DataOut::type_dof_data, + interpretation); + + + data_out.add_data_vector(locally_relevant_magnitude, + solution_names_magnitude, + DataOut::type_dof_data, + interpretation); + + // subdomain visualization (unchanged) + Vector subdomain(triangulation.n_active_cells()); + for (unsigned int i = 0; i < subdomain.size(); ++i) + subdomain(i) = triangulation.locally_owned_subdomain(); + + data_out.add_data_vector(subdomain, "subdomain"); + + data_out.build_patches(); + + std::ostringstream ss; + ss << std::fixed << std::setprecision(2) << frequency; + std::filesystem::create_directories("./output/"); + data_out.write_vtu_with_pvtu_record( + "./output/", "solution_" + ss.str(), 0, mpi_communicator, 2, 8); + } + } +} // namespace VibroAcousticProblem + +int +main(int argc, char *argv[]) +{ + try + { + using namespace dealii; + const unsigned int dim = 3; + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + double f0 = 50.; // starting frequency in Hz + unsigned int n = 0; // number of frequency steps + double f_c = f0; // evaluation frequency + + while (f_c < 1000.) + { + f_c = f0 * std::pow(10., static_cast(n) / 60.); + double omega = f_c * 2. * numbers::PI; + VibroAcousticProblem::HarmonicResponse elastic_problem(omega); + elastic_problem.run(true); + n++; + } + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} \ No newline at end of file diff --git a/Parallel_Vibro-Acoustic_Solver/tria.inp b/Parallel_Vibro-Acoustic_Solver/tria.inp new file mode 100644 index 00000000..9c2997f7 --- /dev/null +++ b/Parallel_Vibro-Acoustic_Solver/tria.inp @@ -0,0 +1,850 @@ +********************************** N O D E S ********************************** +*NODE, NSET=ALLNODES + 22307, 2.030000e+02, -1.826000e+03, -2.591000e+03 + 22308, 2.030000e+02, 1.826000e+03, -2.591000e+03 + 22313, 2.030000e+02, -1.826000e+03, 2.591000e+03 + 22320, 2.030000e+02, 1.826000e+03, 2.591000e+03 + 22331, 2.030000e+02, -1.420000e+03, -2.185000e+03 + 22332, 2.030000e+02, 1.420000e+03, -2.185000e+03 + 22336, 2.030000e+02, 1.420000e+03, 2.185000e+03 + 22342, 2.030000e+02, -1.420000e+03, 2.185000e+03 + 22351, 2.030000e+02, -1.369500e+03, -2.591000e+03 + 22352, 2.030000e+02, -9.130000e+02, -2.591000e+03 + 22353, 2.030000e+02, -4.565000e+02, -2.591000e+03 + 22354, 2.030000e+02, 0.000000e+00, -2.591000e+03 + 22355, 2.030000e+02, 4.565000e+02, -2.591000e+03 + 22356, 2.030000e+02, 9.130000e+02, -2.591000e+03 + 22357, 2.030000e+02, 1.369500e+03, -2.591000e+03 + 22358, 2.030000e+02, -1.826000e+03, 2.119909e+03 + 22359, 2.030000e+02, -1.826000e+03, 1.648818e+03 + 22360, 2.030000e+02, -1.826000e+03, 1.177727e+03 + 22361, 2.030000e+02, -1.826000e+03, 7.066364e+02 + 22362, 2.030000e+02, -1.826000e+03, 2.355455e+02 + 22363, 2.030000e+02, -1.826000e+03, -2.355455e+02 + 22364, 2.030000e+02, -1.826000e+03, -7.066364e+02 + 22365, 2.030000e+02, -1.826000e+03, -1.177727e+03 + 22366, 2.030000e+02, -1.826000e+03, -1.648818e+03 + 22367, 2.030000e+02, -1.826000e+03, -2.119909e+03 + 22368, 2.030000e+02, 1.369500e+03, 2.591000e+03 + 22369, 2.030000e+02, 9.130000e+02, 2.591000e+03 + 22370, 2.030000e+02, 4.565000e+02, 2.591000e+03 + 22371, 2.030000e+02, 0.000000e+00, 2.591000e+03 + 22372, 2.030000e+02, -4.565000e+02, 2.591000e+03 + 22373, 2.030000e+02, -9.130000e+02, 2.591000e+03 + 22374, 2.030000e+02, -1.369500e+03, 2.591000e+03 + 22375, 2.030000e+02, 1.826000e+03, -2.119909e+03 + 22376, 2.030000e+02, 1.826000e+03, -1.648818e+03 + 22377, 2.030000e+02, 1.826000e+03, -1.177727e+03 + 22378, 2.030000e+02, 1.826000e+03, -7.066364e+02 + 22379, 2.030000e+02, 1.826000e+03, -2.355455e+02 + 22380, 2.030000e+02, 1.826000e+03, 2.355455e+02 + 22381, 2.030000e+02, 1.826000e+03, 7.066364e+02 + 22382, 2.030000e+02, 1.826000e+03, 1.177727e+03 + 22383, 2.030000e+02, 1.826000e+03, 1.648818e+03 + 22384, 2.030000e+02, 1.826000e+03, 2.119909e+03 + 22385, 2.030000e+02, -9.466667e+02, -2.185000e+03 + 22386, 2.030000e+02, -4.733333e+02, -2.185000e+03 + 22387, 2.030000e+02, 0.000000e+00, -2.185000e+03 + 22388, 2.030000e+02, 4.733333e+02, -2.185000e+03 + 22389, 2.030000e+02, 9.466667e+02, -2.185000e+03 + 22390, 2.030000e+02, 1.420000e+03, -1.699444e+03 + 22391, 2.030000e+02, 1.420000e+03, -1.213889e+03 + 22392, 2.030000e+02, 1.420000e+03, -7.283333e+02 + 22393, 2.030000e+02, 1.420000e+03, -2.427778e+02 + 22394, 2.030000e+02, 1.420000e+03, 2.427778e+02 + 22395, 2.030000e+02, 1.420000e+03, 7.283333e+02 + 22396, 2.030000e+02, 1.420000e+03, 1.213889e+03 + 22397, 2.030000e+02, 1.420000e+03, 1.699444e+03 + 22398, 2.030000e+02, -9.466667e+02, 2.185000e+03 + 22399, 2.030000e+02, -4.733333e+02, 2.185000e+03 + 22400, 2.030000e+02, 0.000000e+00, 2.185000e+03 + 22401, 2.030000e+02, 4.733333e+02, 2.185000e+03 + 22402, 2.030000e+02, 9.466667e+02, 2.185000e+03 + 22403, 2.030000e+02, -1.420000e+03, -1.699444e+03 + 22404, 2.030000e+02, -1.420000e+03, -1.213889e+03 + 22405, 2.030000e+02, -1.420000e+03, -7.283333e+02 + 22406, 2.030000e+02, -1.420000e+03, -2.427778e+02 + 22407, 2.030000e+02, -1.420000e+03, 2.427778e+02 + 22408, 2.030000e+02, -1.420000e+03, 7.283333e+02 + 22409, 2.030000e+02, -1.420000e+03, 1.213889e+03 + 22410, 2.030000e+02, -1.420000e+03, 1.699444e+03 + 22411, 0.000000e+00, -1.420000e+03, 2.185000e+03 + 22412, 0.000000e+00, 1.420000e+03, 2.185000e+03 + 22413, 0.000000e+00, -9.466667e+02, 2.185000e+03 + 22414, 0.000000e+00, -4.733333e+02, 2.185000e+03 + 22415, 0.000000e+00, 0.000000e+00, 2.185000e+03 + 22416, 0.000000e+00, 4.733333e+02, 2.185000e+03 + 22417, 0.000000e+00, 9.466667e+02, 2.185000e+03 + 22418, 0.000000e+00, 1.420000e+03, -2.185000e+03 + 22419, 0.000000e+00, 1.420000e+03, -1.699444e+03 + 22420, 0.000000e+00, 1.420000e+03, -1.213889e+03 + 22421, 0.000000e+00, 1.420000e+03, -7.283333e+02 + 22422, 0.000000e+00, 1.420000e+03, -2.427778e+02 + 22423, 0.000000e+00, 1.420000e+03, 2.427778e+02 + 22424, 0.000000e+00, 1.420000e+03, 7.283333e+02 + 22425, 0.000000e+00, 1.420000e+03, 1.213889e+03 + 22426, 0.000000e+00, 1.420000e+03, 1.699444e+03 + 22427, 0.000000e+00, -1.420000e+03, -2.185000e+03 + 22428, 0.000000e+00, -9.466667e+02, -2.185000e+03 + 22429, 0.000000e+00, -4.733333e+02, -2.185000e+03 + 22430, 0.000000e+00, 0.000000e+00, -2.185000e+03 + 22431, 0.000000e+00, 4.733333e+02, -2.185000e+03 + 22432, 0.000000e+00, 9.466667e+02, -2.185000e+03 + 22433, 0.000000e+00, -1.420000e+03, -1.699444e+03 + 22434, 0.000000e+00, -1.420000e+03, -1.213889e+03 + 22435, 0.000000e+00, -1.420000e+03, -7.283333e+02 + 22436, 0.000000e+00, -1.420000e+03, -2.427778e+02 + 22437, 0.000000e+00, -1.420000e+03, 2.427778e+02 + 22438, 0.000000e+00, -1.420000e+03, 7.283333e+02 + 22439, 0.000000e+00, -1.420000e+03, 1.213889e+03 + 22440, 0.000000e+00, -1.420000e+03, 1.699444e+03 + 22441, 0.000000e+00, -9.466667e+02, 1.699444e+03 + 22442, 0.000000e+00, -4.733333e+02, 1.699444e+03 + 22443, 0.000000e+00, 0.000000e+00, 1.699444e+03 + 22444, 0.000000e+00, 4.733333e+02, 1.699444e+03 + 22445, 0.000000e+00, 9.466667e+02, 1.699444e+03 + 22446, 0.000000e+00, -9.466667e+02, 1.213889e+03 + 22447, 0.000000e+00, -4.733333e+02, 1.213889e+03 + 22448, 0.000000e+00, 0.000000e+00, 1.213889e+03 + 22449, 0.000000e+00, 4.733333e+02, 1.213889e+03 + 22450, 0.000000e+00, 9.466667e+02, 1.213889e+03 + 22451, 0.000000e+00, -9.466667e+02, 7.283333e+02 + 22452, 0.000000e+00, -4.733333e+02, 7.283333e+02 + 22453, 0.000000e+00, 0.000000e+00, 7.283333e+02 + 22454, 0.000000e+00, 4.733333e+02, 7.283333e+02 + 22455, 0.000000e+00, 9.466667e+02, 7.283333e+02 + 22456, 0.000000e+00, -9.466667e+02, 2.427778e+02 + 22457, 0.000000e+00, -4.733333e+02, 2.427778e+02 + 22458, 0.000000e+00, 2.273737e-13, 2.427778e+02 + 22459, 0.000000e+00, 4.733333e+02, 2.427778e+02 + 22460, 0.000000e+00, 9.466667e+02, 2.427778e+02 + 22461, 0.000000e+00, -9.466667e+02, -2.427778e+02 + 22462, 0.000000e+00, -4.733333e+02, -2.427778e+02 + 22463, 0.000000e+00, -2.273737e-13, -2.427778e+02 + 22464, 0.000000e+00, 4.733333e+02, -2.427778e+02 + 22465, 0.000000e+00, 9.466667e+02, -2.427778e+02 + 22466, 0.000000e+00, -9.466667e+02, -7.283333e+02 + 22467, 0.000000e+00, -4.733333e+02, -7.283333e+02 + 22468, 0.000000e+00, 2.273737e-13, -7.283333e+02 + 22469, 0.000000e+00, 4.733333e+02, -7.283333e+02 + 22470, 0.000000e+00, 9.466667e+02, -7.283333e+02 + 22471, 0.000000e+00, -9.466667e+02, -1.213889e+03 + 22472, 0.000000e+00, -4.733333e+02, -1.213889e+03 + 22473, 0.000000e+00, 0.000000e+00, -1.213889e+03 + 22474, 0.000000e+00, 4.733333e+02, -1.213889e+03 + 22475, 0.000000e+00, 9.466667e+02, -1.213889e+03 + 22476, 0.000000e+00, -9.466667e+02, -1.699444e+03 + 22477, 0.000000e+00, -4.733333e+02, -1.699444e+03 + 22478, 0.000000e+00, 0.000000e+00, -1.699444e+03 + 22479, 0.000000e+00, 4.733333e+02, -1.699444e+03 + 22480, 0.000000e+00, 9.466667e+02, -1.699444e+03 + 22481, 2.030000e+02, 2.326000e+03, -3.091000e+03 + 22482, 2.030000e+02, -2.326000e+03, -3.091000e+03 + 22483, 2.030000e+02, 1.860800e+03, -3.091000e+03 + 22484, 2.030000e+02, 1.395600e+03, -3.091000e+03 + 22485, 2.030000e+02, 9.304000e+02, -3.091000e+03 + 22486, 2.030000e+02, 4.652000e+02, -3.091000e+03 + 22487, 2.030000e+02, 0.000000e+00, -3.091000e+03 + 22488, 2.030000e+02, -4.652000e+02, -3.091000e+03 + 22489, 2.030000e+02, -9.304000e+02, -3.091000e+03 + 22490, 2.030000e+02, -1.395600e+03, -3.091000e+03 + 22491, 2.030000e+02, -1.860800e+03, -3.091000e+03 + 22492, 2.030000e+02, -2.326000e+03, 3.091000e+03 + 22493, 2.030000e+02, -2.326000e+03, -2.615462e+03 + 22494, 2.030000e+02, -2.326000e+03, -2.139923e+03 + 22495, 2.030000e+02, -2.326000e+03, -1.664385e+03 + 22496, 2.030000e+02, -2.326000e+03, -1.188846e+03 + 22497, 2.030000e+02, -2.326000e+03, -7.133077e+02 + 22498, 2.030000e+02, -2.326000e+03, -2.377692e+02 + 22499, 2.030000e+02, -2.326000e+03, 2.377692e+02 + 22500, 2.030000e+02, -2.326000e+03, 7.133077e+02 + 22501, 2.030000e+02, -2.326000e+03, 1.188846e+03 + 22502, 2.030000e+02, -2.326000e+03, 1.664385e+03 + 22503, 2.030000e+02, -2.326000e+03, 2.139923e+03 + 22504, 2.030000e+02, -2.326000e+03, 2.615462e+03 + 22505, 2.030000e+02, 2.326000e+03, 3.091000e+03 + 22506, 2.030000e+02, -1.860800e+03, 3.091000e+03 + 22507, 2.030000e+02, -1.395600e+03, 3.091000e+03 + 22508, 2.030000e+02, -9.304000e+02, 3.091000e+03 + 22509, 2.030000e+02, -4.652000e+02, 3.091000e+03 + 22510, 2.030000e+02, 0.000000e+00, 3.091000e+03 + 22511, 2.030000e+02, 4.652000e+02, 3.091000e+03 + 22512, 2.030000e+02, 9.304000e+02, 3.091000e+03 + 22513, 2.030000e+02, 1.395600e+03, 3.091000e+03 + 22514, 2.030000e+02, 1.860800e+03, 3.091000e+03 + 22515, 2.030000e+02, 2.326000e+03, 2.615462e+03 + 22516, 2.030000e+02, 2.326000e+03, 2.139923e+03 + 22517, 2.030000e+02, 2.326000e+03, 1.664385e+03 + 22518, 2.030000e+02, 2.326000e+03, 1.188846e+03 + 22519, 2.030000e+02, 2.326000e+03, 7.133077e+02 + 22520, 2.030000e+02, 2.326000e+03, 2.377692e+02 + 22521, 2.030000e+02, 2.326000e+03, -2.377692e+02 + 22522, 2.030000e+02, 2.326000e+03, -7.133077e+02 + 22523, 2.030000e+02, 2.326000e+03, -1.188846e+03 + 22524, 2.030000e+02, 2.326000e+03, -1.664385e+03 + 22525, 2.030000e+02, 2.326000e+03, -2.139923e+03 + 22526, 2.030000e+02, 2.326000e+03, -2.615462e+03 + 22527, 8.120000e+02, 2.326000e+03, -3.091000e+03 + 22528, 8.120000e+02, -2.326000e+03, -3.091000e+03 + 22529, 8.120000e+02, 1.860800e+03, -3.091000e+03 + 22530, 8.120000e+02, 1.395600e+03, -3.091000e+03 + 22531, 8.120000e+02, 9.304000e+02, -3.091000e+03 + 22532, 8.120000e+02, 4.652000e+02, -3.091000e+03 + 22533, 8.120000e+02, 0.000000e+00, -3.091000e+03 + 22534, 8.120000e+02, -4.652000e+02, -3.091000e+03 + 22535, 8.120000e+02, -9.304000e+02, -3.091000e+03 + 22536, 8.120000e+02, -1.395600e+03, -3.091000e+03 + 22537, 8.120000e+02, -1.860800e+03, -3.091000e+03 + 22538, 8.120000e+02, 2.326000e+03, 3.091000e+03 + 22539, 8.120000e+02, 2.326000e+03, 2.615462e+03 + 22540, 8.120000e+02, 2.326000e+03, 2.139923e+03 + 22541, 8.120000e+02, 2.326000e+03, 1.664385e+03 + 22542, 8.120000e+02, 2.326000e+03, 1.188846e+03 + 22543, 8.120000e+02, 2.326000e+03, 7.133077e+02 + 22544, 8.120000e+02, 2.326000e+03, 2.377692e+02 + 22545, 8.120000e+02, 2.326000e+03, -2.377692e+02 + 22546, 8.120000e+02, 2.326000e+03, -7.133077e+02 + 22547, 8.120000e+02, 2.326000e+03, -1.188846e+03 + 22548, 8.120000e+02, 2.326000e+03, -1.664385e+03 + 22549, 8.120000e+02, 2.326000e+03, -2.139923e+03 + 22550, 8.120000e+02, 2.326000e+03, -2.615462e+03 + 22551, 8.120000e+02, -2.326000e+03, 3.091000e+03 + 22552, 8.120000e+02, -2.326000e+03, -2.615462e+03 + 22553, 8.120000e+02, -2.326000e+03, -2.139923e+03 + 22554, 8.120000e+02, -2.326000e+03, -1.664385e+03 + 22555, 8.120000e+02, -2.326000e+03, -1.188846e+03 + 22556, 8.120000e+02, -2.326000e+03, -7.133077e+02 + 22557, 8.120000e+02, -2.326000e+03, -2.377692e+02 + 22558, 8.120000e+02, -2.326000e+03, 2.377692e+02 + 22559, 8.120000e+02, -2.326000e+03, 7.133077e+02 + 22560, 8.120000e+02, -2.326000e+03, 1.188846e+03 + 22561, 8.120000e+02, -2.326000e+03, 1.664385e+03 + 22562, 8.120000e+02, -2.326000e+03, 2.139923e+03 + 22563, 8.120000e+02, -2.326000e+03, 2.615462e+03 + 22564, 1.312000e+03, 2.326000e+03, -3.091000e+03 + 22565, 1.312000e+03, -2.326000e+03, -3.091000e+03 + 22566, 1.312000e+03, -1.860800e+03, -3.091000e+03 + 22567, 1.312000e+03, -1.395600e+03, -3.091000e+03 + 22568, 1.312000e+03, -9.304000e+02, -3.091000e+03 + 22569, 1.312000e+03, -4.652000e+02, -3.091000e+03 + 22570, 1.312000e+03, 0.000000e+00, -3.091000e+03 + 22571, 1.312000e+03, 4.652000e+02, -3.091000e+03 + 22572, 1.312000e+03, 9.304000e+02, -3.091000e+03 + 22573, 1.312000e+03, 1.395600e+03, -3.091000e+03 + 22574, 1.312000e+03, 1.860800e+03, -3.091000e+03 + 22575, 8.120000e+02, -1.860800e+03, 3.091000e+03 + 22576, 8.120000e+02, -1.395600e+03, 3.091000e+03 + 22577, 8.120000e+02, -9.304000e+02, 3.091000e+03 + 22578, 8.120000e+02, -4.652000e+02, 3.091000e+03 + 22579, 8.120000e+02, 0.000000e+00, 3.091000e+03 + 22580, 8.120000e+02, 4.652000e+02, 3.091000e+03 + 22581, 8.120000e+02, 9.304000e+02, 3.091000e+03 + 22582, 8.120000e+02, 1.395600e+03, 3.091000e+03 + 22583, 8.120000e+02, 1.860800e+03, 3.091000e+03 + 22584, 1.312000e+03, 2.326000e+03, 3.091000e+03 + 22585, 1.312000e+03, 2.326000e+03, -2.615462e+03 + 22586, 1.312000e+03, 2.326000e+03, -2.139923e+03 + 22587, 1.312000e+03, 2.326000e+03, -1.664385e+03 + 22588, 1.312000e+03, 2.326000e+03, -1.188846e+03 + 22589, 1.312000e+03, 2.326000e+03, -7.133077e+02 + 22590, 1.312000e+03, 2.326000e+03, -2.377692e+02 + 22591, 1.312000e+03, 2.326000e+03, 2.377692e+02 + 22592, 1.312000e+03, 2.326000e+03, 7.133077e+02 + 22593, 1.312000e+03, 2.326000e+03, 1.188846e+03 + 22594, 1.312000e+03, 2.326000e+03, 1.664385e+03 + 22595, 1.312000e+03, 2.326000e+03, 2.139923e+03 + 22596, 1.312000e+03, 2.326000e+03, 2.615462e+03 + 22597, 1.312000e+03, -2.326000e+03, 3.091000e+03 + 22598, 1.312000e+03, -2.326000e+03, 2.615462e+03 + 22599, 1.312000e+03, -2.326000e+03, 2.139923e+03 + 22600, 1.312000e+03, -2.326000e+03, 1.664385e+03 + 22601, 1.312000e+03, -2.326000e+03, 1.188846e+03 + 22602, 1.312000e+03, -2.326000e+03, 7.133077e+02 + 22603, 1.312000e+03, -2.326000e+03, 2.377692e+02 + 22604, 1.312000e+03, -2.326000e+03, -2.377692e+02 + 22605, 1.312000e+03, -2.326000e+03, -7.133077e+02 + 22606, 1.312000e+03, -2.326000e+03, -1.188846e+03 + 22607, 1.312000e+03, -2.326000e+03, -1.664385e+03 + 22608, 1.312000e+03, -2.326000e+03, -2.139923e+03 + 22609, 1.312000e+03, -2.326000e+03, -2.615462e+03 + 22610, 1.312000e+03, 1.860800e+03, 3.091000e+03 + 22611, 1.312000e+03, 1.395600e+03, 3.091000e+03 + 22612, 1.312000e+03, 9.304000e+02, 3.091000e+03 + 22613, 1.312000e+03, 4.652000e+02, 3.091000e+03 + 22614, 1.312000e+03, 0.000000e+00, 3.091000e+03 + 22615, 1.312000e+03, -4.652000e+02, 3.091000e+03 + 22616, 1.312000e+03, -9.304000e+02, 3.091000e+03 + 22617, 1.312000e+03, -1.395600e+03, 3.091000e+03 + 22618, 1.312000e+03, -1.860800e+03, 3.091000e+03 + 22619, 8.120000e+02, 9.466667e+02, -2.185000e+03 + 22620, 8.120000e+02, 1.420000e+03, -2.185000e+03 + 22621, 8.120000e+02, 4.733333e+02, -2.185000e+03 + 22622, 8.120000e+02, -2.042929e-30, -2.185000e+03 + 22623, 8.120000e+02, -4.733333e+02, -2.185000e+03 + 22624, 8.120000e+02, -9.466667e+02, -2.185000e+03 + 22625, 8.120000e+02, -1.420000e+03, -2.185000e+03 + 22626, 8.120000e+02, 1.420000e+03, 1.699444e+03 + 22627, 8.120000e+02, 1.420000e+03, 2.185000e+03 + 22628, 8.120000e+02, 1.420000e+03, 1.213889e+03 + 22629, 8.120000e+02, 1.420000e+03, 7.283333e+02 + 22630, 8.120000e+02, 1.420000e+03, 2.427778e+02 + 22631, 8.120000e+02, 1.420000e+03, -2.427778e+02 + 22632, 8.120000e+02, 1.420000e+03, -7.283333e+02 + 22633, 8.120000e+02, 1.420000e+03, -1.213889e+03 + 22634, 8.120000e+02, 1.420000e+03, -1.699444e+03 + 22635, 8.120000e+02, -1.420000e+03, -1.699444e+03 + 22636, 8.120000e+02, -1.420000e+03, -1.213889e+03 + 22637, 8.120000e+02, -1.420000e+03, -7.283333e+02 + 22638, 8.120000e+02, -1.420000e+03, -2.427778e+02 + 22639, 8.120000e+02, -1.420000e+03, 2.427778e+02 + 22640, 8.120000e+02, -1.420000e+03, 7.283333e+02 + 22641, 8.120000e+02, -1.420000e+03, 1.213889e+03 + 22642, 8.120000e+02, -1.420000e+03, 1.699444e+03 + 22643, 8.120000e+02, -1.420000e+03, 2.185000e+03 + 22644, 8.120000e+02, -9.466667e+02, 2.185000e+03 + 22645, 8.120000e+02, -4.733333e+02, 2.185000e+03 + 22646, 8.120000e+02, -2.042929e-30, 2.185000e+03 + 22647, 8.120000e+02, 4.733333e+02, 2.185000e+03 + 22648, 8.120000e+02, 9.466667e+02, 2.185000e+03 + 22649, 8.120000e+02, 1.826000e+03, -2.119909e+03 + 22650, 8.120000e+02, 1.826000e+03, -2.591000e+03 + 22651, 8.120000e+02, 1.369500e+03, -2.591000e+03 + 22652, 8.120000e+02, 9.130000e+02, -2.591000e+03 + 22653, 8.120000e+02, 4.565000e+02, -2.591000e+03 + 22654, 8.120000e+02, -2.816835e-14, -2.591000e+03 + 22655, 8.120000e+02, -4.565000e+02, -2.591000e+03 + 22656, 8.120000e+02, -9.130000e+02, -2.591000e+03 + 22657, 8.120000e+02, -1.369500e+03, -2.591000e+03 + 22658, 8.120000e+02, 1.826000e+03, 2.119909e+03 + 22659, 8.120000e+02, 1.826000e+03, 1.648818e+03 + 22660, 8.120000e+02, 1.826000e+03, 1.177727e+03 + 22661, 8.120000e+02, 1.826000e+03, 7.066364e+02 + 22662, 8.120000e+02, 1.826000e+03, 2.355455e+02 + 22663, 8.120000e+02, 1.826000e+03, -2.355455e+02 + 22664, 8.120000e+02, 1.826000e+03, -7.066364e+02 + 22665, 8.120000e+02, 1.826000e+03, -1.177727e+03 + 22666, 8.120000e+02, 1.826000e+03, -1.648818e+03 + 22667, 8.120000e+02, -1.826000e+03, -2.591000e+03 + 22668, 8.120000e+02, -1.826000e+03, -2.119909e+03 + 22669, 8.120000e+02, -1.826000e+03, -1.648818e+03 + 22670, 8.120000e+02, -1.826000e+03, -1.177727e+03 + 22671, 8.120000e+02, -1.826000e+03, -7.066364e+02 + 22672, 8.120000e+02, -1.826000e+03, -2.355455e+02 + 22673, 8.120000e+02, -1.826000e+03, 2.355455e+02 + 22674, 8.120000e+02, -1.826000e+03, 7.066364e+02 + 22675, 8.120000e+02, -1.826000e+03, 1.177727e+03 + 22676, 8.120000e+02, -1.826000e+03, 1.648818e+03 + 22677, 8.120000e+02, -1.826000e+03, 2.119909e+03 + 22678, 8.120000e+02, -1.826000e+03, 2.591000e+03 + 22679, 8.120000e+02, -1.369500e+03, 2.591000e+03 + 22680, 8.120000e+02, -9.130000e+02, 2.591000e+03 + 22681, 8.120000e+02, -4.565000e+02, 2.591000e+03 + 22682, 8.120000e+02, 2.816835e-14, 2.591000e+03 + 22683, 8.120000e+02, 4.565000e+02, 2.591000e+03 + 22684, 8.120000e+02, 9.130000e+02, 2.591000e+03 + 22685, 8.120000e+02, 1.369500e+03, 2.591000e+03 + 22686, 8.120000e+02, 1.826000e+03, 2.591000e+03 + 22687, 2.030000e+02, -9.466667e+02, 1.699444e+03 + 22688, 2.030000e+02, -4.733333e+02, 1.699444e+03 + 22689, 2.030000e+02, -2.826403e-15, 1.699444e+03 + 22690, 2.030000e+02, 4.733333e+02, 1.699444e+03 + 22691, 2.030000e+02, 9.466667e+02, 1.699444e+03 + 22692, 2.030000e+02, -9.466667e+02, 1.213889e+03 + 22693, 2.030000e+02, -4.733333e+02, 1.213889e+03 + 22694, 2.030000e+02, -2.018859e-15, 1.213889e+03 + 22695, 2.030000e+02, 4.733333e+02, 1.213889e+03 + 22696, 2.030000e+02, 9.466667e+02, 1.213889e+03 + 22697, 2.030000e+02, -9.466667e+02, 7.283333e+02 + 22698, 2.030000e+02, -4.733333e+02, 7.283333e+02 + 22699, 2.030000e+02, -1.211316e-15, 7.283333e+02 + 22700, 2.030000e+02, 4.733333e+02, 7.283333e+02 + 22701, 2.030000e+02, 9.466667e+02, 7.283333e+02 + 22702, 2.030000e+02, -9.466667e+02, 2.427778e+02 + 22703, 2.030000e+02, -4.733333e+02, 2.427778e+02 + 22704, 2.030000e+02, 2.269699e-13, 2.427778e+02 + 22705, 2.030000e+02, 4.733333e+02, 2.427778e+02 + 22706, 2.030000e+02, 9.466667e+02, 2.427778e+02 + 22707, 2.030000e+02, -9.466667e+02, -2.427778e+02 + 22708, 2.030000e+02, -4.733333e+02, -2.427778e+02 + 22709, 2.030000e+02, -2.269699e-13, -2.427778e+02 + 22710, 2.030000e+02, 4.733333e+02, -2.427778e+02 + 22711, 2.030000e+02, 9.466667e+02, -2.427778e+02 + 22712, 2.030000e+02, -9.466667e+02, -7.283333e+02 + 22713, 2.030000e+02, -4.733333e+02, -7.283333e+02 + 22714, 2.030000e+02, 2.285850e-13, -7.283333e+02 + 22715, 2.030000e+02, 4.733333e+02, -7.283333e+02 + 22716, 2.030000e+02, 9.466667e+02, -7.283333e+02 + 22717, 2.030000e+02, -9.466667e+02, -1.213889e+03 + 22718, 2.030000e+02, -4.733333e+02, -1.213889e+03 + 22719, 2.030000e+02, 2.018859e-15, -1.213889e+03 + 22720, 2.030000e+02, 4.733333e+02, -1.213889e+03 + 22721, 2.030000e+02, 9.466667e+02, -1.213889e+03 + 22722, 2.030000e+02, -9.466667e+02, -1.699444e+03 + 22723, 2.030000e+02, -4.733333e+02, -1.699444e+03 + 22724, 2.030000e+02, 2.826403e-15, -1.699444e+03 + 22725, 2.030000e+02, 4.733333e+02, -1.699444e+03 + 22726, 2.030000e+02, 9.466667e+02, -1.699444e+03 + 22727, 8.120000e+02, -9.466667e+02, 1.699444e+03 + 22728, 8.120000e+02, -4.733333e+02, 1.699444e+03 + 22729, 8.120000e+02, -1.130561e-14, 1.699444e+03 + 22730, 8.120000e+02, 4.733333e+02, 1.699444e+03 + 22731, 8.120000e+02, 9.466667e+02, 1.699444e+03 + 22732, 8.120000e+02, -9.466667e+02, 1.213889e+03 + 22733, 8.120000e+02, -4.733333e+02, 1.213889e+03 + 22734, 8.120000e+02, -8.075438e-15, 1.213889e+03 + 22735, 8.120000e+02, 4.733333e+02, 1.213889e+03 + 22736, 8.120000e+02, 9.466667e+02, 1.213889e+03 + 22737, 8.120000e+02, -9.466667e+02, 7.283333e+02 + 22738, 8.120000e+02, -4.733333e+02, 7.283333e+02 + 22739, 8.120000e+02, -4.845263e-15, 7.283333e+02 + 22740, 8.120000e+02, 4.733333e+02, 7.283333e+02 + 22741, 8.120000e+02, 9.466667e+02, 7.283333e+02 + 22742, 8.120000e+02, -9.466667e+02, 2.427778e+02 + 22743, 8.120000e+02, -4.733333e+02, 2.427778e+02 + 22744, 8.120000e+02, 2.257586e-13, 2.427778e+02 + 22745, 8.120000e+02, 4.733333e+02, 2.427778e+02 + 22746, 8.120000e+02, 9.466667e+02, 2.427778e+02 + 22747, 8.120000e+02, -9.466667e+02, -2.427778e+02 + 22748, 8.120000e+02, -4.733333e+02, -2.427778e+02 + 22749, 8.120000e+02, -2.257586e-13, -2.427778e+02 + 22750, 8.120000e+02, 4.733333e+02, -2.427778e+02 + 22751, 8.120000e+02, 9.466667e+02, -2.427778e+02 + 22752, 8.120000e+02, -9.466667e+02, -7.283333e+02 + 22753, 8.120000e+02, -4.733333e+02, -7.283333e+02 + 22754, 8.120000e+02, 2.322189e-13, -7.283333e+02 + 22755, 8.120000e+02, 4.733333e+02, -7.283333e+02 + 22756, 8.120000e+02, 9.466667e+02, -7.283333e+02 + 22757, 8.120000e+02, -9.466667e+02, -1.213889e+03 + 22758, 8.120000e+02, -4.733333e+02, -1.213889e+03 + 22759, 8.120000e+02, 8.075438e-15, -1.213889e+03 + 22760, 8.120000e+02, 4.733333e+02, -1.213889e+03 + 22761, 8.120000e+02, 9.466667e+02, -1.213889e+03 + 22762, 8.120000e+02, -9.466667e+02, -1.699444e+03 + 22763, 8.120000e+02, -4.733333e+02, -1.699444e+03 + 22764, 8.120000e+02, 1.130561e-14, -1.699444e+03 + 22765, 8.120000e+02, 4.733333e+02, -1.699444e+03 + 22766, 8.120000e+02, 9.466667e+02, -1.699444e+03 + 22767, 1.312000e+03, 1.826000e+03, -2.119909e+03 + 22768, 1.312000e+03, 1.826000e+03, -2.591000e+03 + 22769, 1.312000e+03, 1.369500e+03, -2.591000e+03 + 22770, 1.312000e+03, 9.130000e+02, -2.591000e+03 + 22771, 1.312000e+03, 4.565000e+02, -2.591000e+03 + 22772, 1.312000e+03, 8.778924e-15, -2.591000e+03 + 22773, 1.312000e+03, -4.565000e+02, -2.591000e+03 + 22774, 1.312000e+03, -9.130000e+02, -2.591000e+03 + 22775, 1.312000e+03, -1.369500e+03, -2.591000e+03 + 22776, 1.312000e+03, 1.826000e+03, 2.119909e+03 + 22777, 1.312000e+03, 1.826000e+03, 1.648818e+03 + 22778, 1.312000e+03, 1.826000e+03, 1.177727e+03 + 22779, 1.312000e+03, 1.826000e+03, 7.066364e+02 + 22780, 1.312000e+03, 1.826000e+03, 2.355455e+02 + 22781, 1.312000e+03, 1.826000e+03, -2.355455e+02 + 22782, 1.312000e+03, 1.826000e+03, -7.066364e+02 + 22783, 1.312000e+03, 1.826000e+03, -1.177727e+03 + 22784, 1.312000e+03, 1.826000e+03, -1.648818e+03 + 22785, 1.312000e+03, -1.826000e+03, -2.591000e+03 + 22786, 1.312000e+03, -1.826000e+03, -2.119909e+03 + 22787, 1.312000e+03, -1.826000e+03, -1.648818e+03 + 22788, 1.312000e+03, -1.826000e+03, -1.177727e+03 + 22789, 1.312000e+03, -1.826000e+03, -7.066364e+02 + 22790, 1.312000e+03, -1.826000e+03, -2.355455e+02 + 22791, 1.312000e+03, -1.826000e+03, 2.355455e+02 + 22792, 1.312000e+03, -1.826000e+03, 7.066364e+02 + 22793, 1.312000e+03, -1.826000e+03, 1.177727e+03 + 22794, 1.312000e+03, -1.826000e+03, 1.648818e+03 + 22795, 1.312000e+03, -1.826000e+03, 2.119909e+03 + 22796, 1.312000e+03, -1.826000e+03, 2.591000e+03 + 22797, 1.312000e+03, -1.369500e+03, 2.591000e+03 + 22798, 1.312000e+03, -9.130000e+02, 2.591000e+03 + 22799, 1.312000e+03, -4.565000e+02, 2.591000e+03 + 22800, 1.312000e+03, -8.778924e-15, 2.591000e+03 + 22801, 1.312000e+03, 4.565000e+02, 2.591000e+03 + 22802, 1.312000e+03, 9.130000e+02, 2.591000e+03 + 22803, 1.312000e+03, 1.369500e+03, 2.591000e+03 + 22804, 1.312000e+03, 1.826000e+03, 2.591000e+03 + 22805, 1.312000e+03, 9.466667e+02, -2.185000e+03 + 22806, 1.312000e+03, 1.420000e+03, -2.185000e+03 + 22807, 1.312000e+03, 4.733333e+02, -2.185000e+03 + 22808, 1.312000e+03, 7.403299e-15, -2.185000e+03 + 22809, 1.312000e+03, -4.733333e+02, -2.185000e+03 + 22810, 1.312000e+03, -9.466667e+02, -2.185000e+03 + 22811, 1.312000e+03, -1.420000e+03, -2.185000e+03 + 22812, 1.312000e+03, 1.420000e+03, 1.699444e+03 + 22813, 1.312000e+03, 1.420000e+03, 2.185000e+03 + 22814, 1.312000e+03, 1.420000e+03, 1.213889e+03 + 22815, 1.312000e+03, 1.420000e+03, 7.283333e+02 + 22816, 1.312000e+03, 1.420000e+03, 2.427778e+02 + 22817, 1.312000e+03, 1.420000e+03, -2.427778e+02 + 22818, 1.312000e+03, 1.420000e+03, -7.283333e+02 + 22819, 1.312000e+03, 1.420000e+03, -1.213889e+03 + 22820, 1.312000e+03, 1.420000e+03, -1.699444e+03 + 22821, 1.312000e+03, -1.420000e+03, -1.699444e+03 + 22822, 1.312000e+03, -1.420000e+03, -1.213889e+03 + 22823, 1.312000e+03, -1.420000e+03, -7.283333e+02 + 22824, 1.312000e+03, -1.420000e+03, -2.427778e+02 + 22825, 1.312000e+03, -1.420000e+03, 2.427778e+02 + 22826, 1.312000e+03, -1.420000e+03, 7.283333e+02 + 22827, 1.312000e+03, -1.420000e+03, 1.213889e+03 + 22828, 1.312000e+03, -1.420000e+03, 1.699444e+03 + 22829, 1.312000e+03, -1.420000e+03, 2.185000e+03 + 22830, 1.312000e+03, -9.466667e+02, 2.185000e+03 + 22831, 1.312000e+03, -4.733333e+02, 2.185000e+03 + 22832, 1.312000e+03, -7.403299e-15, 2.185000e+03 + 22833, 1.312000e+03, 4.733333e+02, 2.185000e+03 + 22834, 1.312000e+03, 9.466667e+02, 2.185000e+03 + 22835, 1.312000e+03, -9.466667e+02, 1.699444e+03 + 22836, 1.312000e+03, -4.733333e+02, 1.699444e+03 + 22837, 1.312000e+03, 1.554476e-15, 1.699444e+03 + 22838, 1.312000e+03, 4.733333e+02, 1.699444e+03 + 22839, 1.312000e+03, 9.466667e+02, 1.699444e+03 + 22840, 1.312000e+03, -9.466667e+02, 1.213889e+03 + 22841, 1.312000e+03, -4.733333e+02, 1.213889e+03 + 22842, 1.312000e+03, -1.055124e-15, 1.213889e+03 + 22843, 1.312000e+03, 4.733333e+02, 1.213889e+03 + 22844, 1.312000e+03, 9.466667e+02, 1.213889e+03 + 22845, 1.312000e+03, -9.466667e+02, 7.283333e+02 + 22846, 1.312000e+03, -4.733333e+02, 7.283333e+02 + 22847, 1.312000e+03, -3.664723e-15, 7.283333e+02 + 22848, 1.312000e+03, 4.733333e+02, 7.283333e+02 + 22849, 1.312000e+03, 9.466667e+02, 7.283333e+02 + 22850, 1.312000e+03, -9.466667e+02, 2.427778e+02 + 22851, 1.312000e+03, -4.733333e+02, 2.427778e+02 + 22852, 1.312000e+03, 2.210994e-13, 2.427778e+02 + 22853, 1.312000e+03, 4.733333e+02, 2.427778e+02 + 22854, 1.312000e+03, 9.466667e+02, 2.427778e+02 + 22855, 1.312000e+03, -9.466667e+02, -2.427778e+02 + 22856, 1.312000e+03, -4.733333e+02, -2.427778e+02 + 22857, 1.312000e+03, -2.362576e-13, -2.427778e+02 + 22858, 1.312000e+03, 4.733333e+02, -2.427778e+02 + 22859, 1.312000e+03, 9.466667e+02, -2.427778e+02 + 22860, 1.312000e+03, -9.466667e+02, -7.283333e+02 + 22861, 1.312000e+03, -4.733333e+02, -7.283333e+02 + 22862, 1.312000e+03, 2.158802e-13, -7.283333e+02 + 22863, 1.312000e+03, 4.733333e+02, -7.283333e+02 + 22864, 1.312000e+03, 9.466667e+02, -7.283333e+02 + 22865, 1.312000e+03, -9.466667e+02, -1.213889e+03 + 22866, 1.312000e+03, -4.733333e+02, -1.213889e+03 + 22867, 1.312000e+03, -1.410312e-14, -1.213889e+03 + 22868, 1.312000e+03, 4.733333e+02, -1.213889e+03 + 22869, 1.312000e+03, 9.466667e+02, -1.213889e+03 + 22870, 1.312000e+03, -9.466667e+02, -1.699444e+03 + 22871, 1.312000e+03, -4.733333e+02, -1.699444e+03 + 22872, 1.312000e+03, -1.671272e-14, -1.699444e+03 + 22873, 1.312000e+03, 4.733333e+02, -1.699444e+03 + 22874, 1.312000e+03, 9.466667e+02, -1.699444e+03 +********************************** E L E M E N T S **************************** +*ELEMENT, TYPE=C3D8R, ELSET=EB1 + 1, 22726, 22390, 22332, 22389, 22480, 22419, 22418, 22432 + 2, 22725, 22726, 22389, 22388, 22479, 22480, 22432, 22431 + 3, 22724, 22725, 22388, 22387, 22478, 22479, 22431, 22430 + 4, 22723, 22724, 22387, 22386, 22477, 22478, 22430, 22429 + 5, 22722, 22723, 22386, 22385, 22476, 22477, 22429, 22428 + 6, 22403, 22722, 22385, 22331, 22433, 22476, 22428, 22427 + 7, 22721, 22391, 22390, 22726, 22475, 22420, 22419, 22480 + 8, 22720, 22721, 22726, 22725, 22474, 22475, 22480, 22479 + 9, 22719, 22720, 22725, 22724, 22473, 22474, 22479, 22478 + 10, 22718, 22719, 22724, 22723, 22472, 22473, 22478, 22477 + 11, 22717, 22718, 22723, 22722, 22471, 22472, 22477, 22476 + 12, 22404, 22717, 22722, 22403, 22434, 22471, 22476, 22433 + 13, 22716, 22392, 22391, 22721, 22470, 22421, 22420, 22475 + 14, 22715, 22716, 22721, 22720, 22469, 22470, 22475, 22474 + 15, 22714, 22715, 22720, 22719, 22468, 22469, 22474, 22473 + 16, 22713, 22714, 22719, 22718, 22467, 22468, 22473, 22472 + 17, 22712, 22713, 22718, 22717, 22466, 22467, 22472, 22471 + 18, 22405, 22712, 22717, 22404, 22435, 22466, 22471, 22434 + 19, 22711, 22393, 22392, 22716, 22465, 22422, 22421, 22470 + 20, 22710, 22711, 22716, 22715, 22464, 22465, 22470, 22469 + 21, 22709, 22710, 22715, 22714, 22463, 22464, 22469, 22468 + 22, 22708, 22709, 22714, 22713, 22462, 22463, 22468, 22467 + 23, 22707, 22708, 22713, 22712, 22461, 22462, 22467, 22466 + 24, 22406, 22707, 22712, 22405, 22436, 22461, 22466, 22435 + 25, 22706, 22394, 22393, 22711, 22460, 22423, 22422, 22465 + 26, 22705, 22706, 22711, 22710, 22459, 22460, 22465, 22464 + 27, 22704, 22705, 22710, 22709, 22458, 22459, 22464, 22463 + 28, 22703, 22704, 22709, 22708, 22457, 22458, 22463, 22462 + 29, 22702, 22703, 22708, 22707, 22456, 22457, 22462, 22461 + 30, 22407, 22702, 22707, 22406, 22437, 22456, 22461, 22436 + 31, 22701, 22395, 22394, 22706, 22455, 22424, 22423, 22460 + 32, 22700, 22701, 22706, 22705, 22454, 22455, 22460, 22459 + 33, 22699, 22700, 22705, 22704, 22453, 22454, 22459, 22458 + 34, 22698, 22699, 22704, 22703, 22452, 22453, 22458, 22457 + 35, 22697, 22698, 22703, 22702, 22451, 22452, 22457, 22456 + 36, 22408, 22697, 22702, 22407, 22438, 22451, 22456, 22437 + 37, 22696, 22396, 22395, 22701, 22450, 22425, 22424, 22455 + 38, 22695, 22696, 22701, 22700, 22449, 22450, 22455, 22454 + 39, 22694, 22695, 22700, 22699, 22448, 22449, 22454, 22453 + 40, 22693, 22694, 22699, 22698, 22447, 22448, 22453, 22452 + 41, 22692, 22693, 22698, 22697, 22446, 22447, 22452, 22451 + 42, 22409, 22692, 22697, 22408, 22439, 22446, 22451, 22438 + 43, 22691, 22397, 22396, 22696, 22445, 22426, 22425, 22450 + 44, 22690, 22691, 22696, 22695, 22444, 22445, 22450, 22449 + 45, 22689, 22690, 22695, 22694, 22443, 22444, 22449, 22448 + 46, 22688, 22689, 22694, 22693, 22442, 22443, 22448, 22447 + 47, 22687, 22688, 22693, 22692, 22441, 22442, 22447, 22446 + 48, 22410, 22687, 22692, 22409, 22440, 22441, 22446, 22439 + 49, 22402, 22336, 22397, 22691, 22417, 22412, 22426, 22445 + 50, 22401, 22402, 22691, 22690, 22416, 22417, 22445, 22444 + 51, 22400, 22401, 22690, 22689, 22415, 22416, 22444, 22443 + 52, 22399, 22400, 22689, 22688, 22414, 22415, 22443, 22442 + 53, 22398, 22399, 22688, 22687, 22413, 22414, 22442, 22441 + 54, 22342, 22398, 22687, 22410, 22411, 22413, 22441, 22440 + 55, 22766, 22634, 22620, 22619, 22726, 22390, 22332, 22389 + 56, 22765, 22766, 22619, 22621, 22725, 22726, 22389, 22388 + 57, 22764, 22765, 22621, 22622, 22724, 22725, 22388, 22387 + 58, 22763, 22764, 22622, 22623, 22723, 22724, 22387, 22386 + 59, 22762, 22763, 22623, 22624, 22722, 22723, 22386, 22385 + 60, 22635, 22762, 22624, 22625, 22403, 22722, 22385, 22331 + 61, 22761, 22633, 22634, 22766, 22721, 22391, 22390, 22726 + 62, 22760, 22761, 22766, 22765, 22720, 22721, 22726, 22725 + 63, 22759, 22760, 22765, 22764, 22719, 22720, 22725, 22724 + 64, 22758, 22759, 22764, 22763, 22718, 22719, 22724, 22723 + 65, 22757, 22758, 22763, 22762, 22717, 22718, 22723, 22722 + 66, 22636, 22757, 22762, 22635, 22404, 22717, 22722, 22403 + 67, 22756, 22632, 22633, 22761, 22716, 22392, 22391, 22721 + 68, 22755, 22756, 22761, 22760, 22715, 22716, 22721, 22720 + 69, 22754, 22755, 22760, 22759, 22714, 22715, 22720, 22719 + 70, 22753, 22754, 22759, 22758, 22713, 22714, 22719, 22718 + 71, 22752, 22753, 22758, 22757, 22712, 22713, 22718, 22717 + 72, 22637, 22752, 22757, 22636, 22405, 22712, 22717, 22404 + 73, 22751, 22631, 22632, 22756, 22711, 22393, 22392, 22716 + 74, 22750, 22751, 22756, 22755, 22710, 22711, 22716, 22715 + 75, 22749, 22750, 22755, 22754, 22709, 22710, 22715, 22714 + 76, 22748, 22749, 22754, 22753, 22708, 22709, 22714, 22713 + 77, 22747, 22748, 22753, 22752, 22707, 22708, 22713, 22712 + 78, 22638, 22747, 22752, 22637, 22406, 22707, 22712, 22405 + 79, 22746, 22630, 22631, 22751, 22706, 22394, 22393, 22711 + 80, 22745, 22746, 22751, 22750, 22705, 22706, 22711, 22710 + 81, 22744, 22745, 22750, 22749, 22704, 22705, 22710, 22709 + 82, 22743, 22744, 22749, 22748, 22703, 22704, 22709, 22708 + 83, 22742, 22743, 22748, 22747, 22702, 22703, 22708, 22707 + 84, 22639, 22742, 22747, 22638, 22407, 22702, 22707, 22406 + 85, 22741, 22629, 22630, 22746, 22701, 22395, 22394, 22706 + 86, 22740, 22741, 22746, 22745, 22700, 22701, 22706, 22705 + 87, 22739, 22740, 22745, 22744, 22699, 22700, 22705, 22704 + 88, 22738, 22739, 22744, 22743, 22698, 22699, 22704, 22703 + 89, 22737, 22738, 22743, 22742, 22697, 22698, 22703, 22702 + 90, 22640, 22737, 22742, 22639, 22408, 22697, 22702, 22407 + 91, 22736, 22628, 22629, 22741, 22696, 22396, 22395, 22701 + 92, 22735, 22736, 22741, 22740, 22695, 22696, 22701, 22700 + 93, 22734, 22735, 22740, 22739, 22694, 22695, 22700, 22699 + 94, 22733, 22734, 22739, 22738, 22693, 22694, 22699, 22698 + 95, 22732, 22733, 22738, 22737, 22692, 22693, 22698, 22697 + 96, 22641, 22732, 22737, 22640, 22409, 22692, 22697, 22408 + 97, 22731, 22626, 22628, 22736, 22691, 22397, 22396, 22696 + 98, 22730, 22731, 22736, 22735, 22690, 22691, 22696, 22695 + 99, 22729, 22730, 22735, 22734, 22689, 22690, 22695, 22694 + 100, 22728, 22729, 22734, 22733, 22688, 22689, 22694, 22693 + 101, 22727, 22728, 22733, 22732, 22687, 22688, 22693, 22692 + 102, 22642, 22727, 22732, 22641, 22410, 22687, 22692, 22409 + 103, 22648, 22627, 22626, 22731, 22402, 22336, 22397, 22691 + 104, 22647, 22648, 22731, 22730, 22401, 22402, 22691, 22690 + 105, 22646, 22647, 22730, 22729, 22400, 22401, 22690, 22689 + 106, 22645, 22646, 22729, 22728, 22399, 22400, 22689, 22688 + 107, 22644, 22645, 22728, 22727, 22398, 22399, 22688, 22687 + 108, 22643, 22644, 22727, 22642, 22342, 22398, 22687, 22410 + 109, 22537, 22528, 22552, 22667, 22491, 22482, 22493, 22307 + 110, 22536, 22537, 22667, 22657, 22490, 22491, 22307, 22351 + 111, 22535, 22536, 22657, 22656, 22489, 22490, 22351, 22352 + 112, 22534, 22535, 22656, 22655, 22488, 22489, 22352, 22353 + 113, 22533, 22534, 22655, 22654, 22487, 22488, 22353, 22354 + 114, 22532, 22533, 22654, 22653, 22486, 22487, 22354, 22355 + 115, 22531, 22532, 22653, 22652, 22485, 22486, 22355, 22356 + 116, 22530, 22531, 22652, 22651, 22484, 22485, 22356, 22357 + 117, 22529, 22530, 22651, 22650, 22483, 22484, 22357, 22308 + 118, 22650, 22550, 22527, 22529, 22308, 22526, 22481, 22483 + 119, 22649, 22549, 22550, 22650, 22375, 22525, 22526, 22308 + 120, 22666, 22548, 22549, 22649, 22376, 22524, 22525, 22375 + 121, 22665, 22547, 22548, 22666, 22377, 22523, 22524, 22376 + 122, 22664, 22546, 22547, 22665, 22378, 22522, 22523, 22377 + 123, 22663, 22545, 22546, 22664, 22379, 22521, 22522, 22378 + 124, 22662, 22544, 22545, 22663, 22380, 22520, 22521, 22379 + 125, 22661, 22543, 22544, 22662, 22381, 22519, 22520, 22380 + 126, 22660, 22542, 22543, 22661, 22382, 22518, 22519, 22381 + 127, 22659, 22541, 22542, 22660, 22383, 22517, 22518, 22382 + 128, 22658, 22540, 22541, 22659, 22384, 22516, 22517, 22383 + 129, 22686, 22539, 22540, 22658, 22320, 22515, 22516, 22384 + 130, 22563, 22551, 22575, 22678, 22504, 22492, 22506, 22313 + 131, 22562, 22563, 22678, 22677, 22503, 22504, 22313, 22358 + 132, 22561, 22562, 22677, 22676, 22502, 22503, 22358, 22359 + 133, 22560, 22561, 22676, 22675, 22501, 22502, 22359, 22360 + 134, 22559, 22560, 22675, 22674, 22500, 22501, 22360, 22361 + 135, 22558, 22559, 22674, 22673, 22499, 22500, 22361, 22362 + 136, 22557, 22558, 22673, 22672, 22498, 22499, 22362, 22363 + 137, 22556, 22557, 22672, 22671, 22497, 22498, 22363, 22364 + 138, 22555, 22556, 22671, 22670, 22496, 22497, 22364, 22365 + 139, 22554, 22555, 22670, 22669, 22495, 22496, 22365, 22366 + 140, 22553, 22554, 22669, 22668, 22494, 22495, 22366, 22367 + 141, 22552, 22553, 22668, 22667, 22493, 22494, 22367, 22307 + 142, 22583, 22538, 22539, 22686, 22514, 22505, 22515, 22320 + 143, 22582, 22583, 22686, 22685, 22513, 22514, 22320, 22368 + 144, 22581, 22582, 22685, 22684, 22512, 22513, 22368, 22369 + 145, 22580, 22581, 22684, 22683, 22511, 22512, 22369, 22370 + 146, 22579, 22580, 22683, 22682, 22510, 22511, 22370, 22371 + 147, 22578, 22579, 22682, 22681, 22509, 22510, 22371, 22372 + 148, 22577, 22578, 22681, 22680, 22508, 22509, 22372, 22373 + 149, 22576, 22577, 22680, 22679, 22507, 22508, 22373, 22374 + 150, 22575, 22576, 22679, 22678, 22506, 22507, 22374, 22313 + 151, 22685, 22686, 22658, 22627, 22368, 22320, 22384, 22336 + 152, 22684, 22685, 22627, 22648, 22369, 22368, 22336, 22402 + 153, 22683, 22684, 22648, 22647, 22370, 22369, 22402, 22401 + 154, 22682, 22683, 22647, 22646, 22371, 22370, 22401, 22400 + 155, 22681, 22682, 22646, 22645, 22372, 22371, 22400, 22399 + 156, 22680, 22681, 22645, 22644, 22373, 22372, 22399, 22398 + 157, 22679, 22680, 22644, 22643, 22374, 22373, 22398, 22342 + 158, 22643, 22677, 22678, 22679, 22342, 22358, 22313, 22374 + 159, 22642, 22676, 22677, 22643, 22410, 22359, 22358, 22342 + 160, 22641, 22675, 22676, 22642, 22409, 22360, 22359, 22410 + 161, 22640, 22674, 22675, 22641, 22408, 22361, 22360, 22409 + 162, 22639, 22673, 22674, 22640, 22407, 22362, 22361, 22408 + 163, 22638, 22672, 22673, 22639, 22406, 22363, 22362, 22407 + 164, 22637, 22671, 22672, 22638, 22405, 22364, 22363, 22406 + 165, 22636, 22670, 22671, 22637, 22404, 22365, 22364, 22405 + 166, 22635, 22669, 22670, 22636, 22403, 22366, 22365, 22404 + 167, 22625, 22668, 22669, 22635, 22331, 22367, 22366, 22403 + 168, 22657, 22667, 22668, 22625, 22351, 22307, 22367, 22331 + 169, 22666, 22649, 22620, 22634, 22376, 22375, 22332, 22390 + 170, 22665, 22666, 22634, 22633, 22377, 22376, 22390, 22391 + 171, 22664, 22665, 22633, 22632, 22378, 22377, 22391, 22392 + 172, 22663, 22664, 22632, 22631, 22379, 22378, 22392, 22393 + 173, 22662, 22663, 22631, 22630, 22380, 22379, 22393, 22394 + 174, 22661, 22662, 22630, 22629, 22381, 22380, 22394, 22395 + 175, 22660, 22661, 22629, 22628, 22382, 22381, 22395, 22396 + 176, 22659, 22660, 22628, 22626, 22383, 22382, 22396, 22397 + 177, 22658, 22659, 22626, 22627, 22384, 22383, 22397, 22336 + 178, 22624, 22656, 22657, 22625, 22385, 22352, 22351, 22331 + 179, 22623, 22655, 22656, 22624, 22386, 22353, 22352, 22385 + 180, 22622, 22654, 22655, 22623, 22387, 22354, 22353, 22386 + 181, 22621, 22653, 22654, 22622, 22388, 22355, 22354, 22387 + 182, 22619, 22652, 22653, 22621, 22389, 22356, 22355, 22388 + 183, 22620, 22651, 22652, 22619, 22332, 22357, 22356, 22389 + 184, 22649, 22650, 22651, 22620, 22375, 22308, 22357, 22332 + 185, 22874, 22820, 22806, 22805, 22766, 22634, 22620, 22619 + 186, 22873, 22874, 22805, 22807, 22765, 22766, 22619, 22621 + 187, 22872, 22873, 22807, 22808, 22764, 22765, 22621, 22622 + 188, 22871, 22872, 22808, 22809, 22763, 22764, 22622, 22623 + 189, 22870, 22871, 22809, 22810, 22762, 22763, 22623, 22624 + 190, 22821, 22870, 22810, 22811, 22635, 22762, 22624, 22625 + 191, 22869, 22819, 22820, 22874, 22761, 22633, 22634, 22766 + 192, 22868, 22869, 22874, 22873, 22760, 22761, 22766, 22765 + 193, 22867, 22868, 22873, 22872, 22759, 22760, 22765, 22764 + 194, 22866, 22867, 22872, 22871, 22758, 22759, 22764, 22763 + 195, 22865, 22866, 22871, 22870, 22757, 22758, 22763, 22762 + 196, 22822, 22865, 22870, 22821, 22636, 22757, 22762, 22635 + 197, 22864, 22818, 22819, 22869, 22756, 22632, 22633, 22761 + 198, 22863, 22864, 22869, 22868, 22755, 22756, 22761, 22760 + 199, 22862, 22863, 22868, 22867, 22754, 22755, 22760, 22759 + 200, 22861, 22862, 22867, 22866, 22753, 22754, 22759, 22758 + 201, 22860, 22861, 22866, 22865, 22752, 22753, 22758, 22757 + 202, 22823, 22860, 22865, 22822, 22637, 22752, 22757, 22636 + 203, 22859, 22817, 22818, 22864, 22751, 22631, 22632, 22756 + 204, 22858, 22859, 22864, 22863, 22750, 22751, 22756, 22755 + 205, 22857, 22858, 22863, 22862, 22749, 22750, 22755, 22754 + 206, 22856, 22857, 22862, 22861, 22748, 22749, 22754, 22753 + 207, 22855, 22856, 22861, 22860, 22747, 22748, 22753, 22752 + 208, 22824, 22855, 22860, 22823, 22638, 22747, 22752, 22637 + 209, 22854, 22816, 22817, 22859, 22746, 22630, 22631, 22751 + 210, 22853, 22854, 22859, 22858, 22745, 22746, 22751, 22750 + 211, 22852, 22853, 22858, 22857, 22744, 22745, 22750, 22749 + 212, 22851, 22852, 22857, 22856, 22743, 22744, 22749, 22748 + 213, 22850, 22851, 22856, 22855, 22742, 22743, 22748, 22747 + 214, 22825, 22850, 22855, 22824, 22639, 22742, 22747, 22638 + 215, 22849, 22815, 22816, 22854, 22741, 22629, 22630, 22746 + 216, 22848, 22849, 22854, 22853, 22740, 22741, 22746, 22745 + 217, 22847, 22848, 22853, 22852, 22739, 22740, 22745, 22744 + 218, 22846, 22847, 22852, 22851, 22738, 22739, 22744, 22743 + 219, 22845, 22846, 22851, 22850, 22737, 22738, 22743, 22742 + 220, 22826, 22845, 22850, 22825, 22640, 22737, 22742, 22639 + 221, 22844, 22814, 22815, 22849, 22736, 22628, 22629, 22741 + 222, 22843, 22844, 22849, 22848, 22735, 22736, 22741, 22740 + 223, 22842, 22843, 22848, 22847, 22734, 22735, 22740, 22739 + 224, 22841, 22842, 22847, 22846, 22733, 22734, 22739, 22738 + 225, 22840, 22841, 22846, 22845, 22732, 22733, 22738, 22737 + 226, 22827, 22840, 22845, 22826, 22641, 22732, 22737, 22640 + 227, 22839, 22812, 22814, 22844, 22731, 22626, 22628, 22736 + 228, 22838, 22839, 22844, 22843, 22730, 22731, 22736, 22735 + 229, 22837, 22838, 22843, 22842, 22729, 22730, 22735, 22734 + 230, 22836, 22837, 22842, 22841, 22728, 22729, 22734, 22733 + 231, 22835, 22836, 22841, 22840, 22727, 22728, 22733, 22732 + 232, 22828, 22835, 22840, 22827, 22642, 22727, 22732, 22641 + 233, 22834, 22813, 22812, 22839, 22648, 22627, 22626, 22731 + 234, 22833, 22834, 22839, 22838, 22647, 22648, 22731, 22730 + 235, 22832, 22833, 22838, 22837, 22646, 22647, 22730, 22729 + 236, 22831, 22832, 22837, 22836, 22645, 22646, 22729, 22728 + 237, 22830, 22831, 22836, 22835, 22644, 22645, 22728, 22727 + 238, 22829, 22830, 22835, 22828, 22643, 22644, 22727, 22642 + 239, 22566, 22565, 22609, 22785, 22537, 22528, 22552, 22667 + 240, 22567, 22566, 22785, 22775, 22536, 22537, 22667, 22657 + 241, 22568, 22567, 22775, 22774, 22535, 22536, 22657, 22656 + 242, 22569, 22568, 22774, 22773, 22534, 22535, 22656, 22655 + 243, 22570, 22569, 22773, 22772, 22533, 22534, 22655, 22654 + 244, 22571, 22570, 22772, 22771, 22532, 22533, 22654, 22653 + 245, 22572, 22571, 22771, 22770, 22531, 22532, 22653, 22652 + 246, 22573, 22572, 22770, 22769, 22530, 22531, 22652, 22651 + 247, 22574, 22573, 22769, 22768, 22529, 22530, 22651, 22650 + 248, 22768, 22585, 22564, 22574, 22650, 22550, 22527, 22529 + 249, 22767, 22586, 22585, 22768, 22649, 22549, 22550, 22650 + 250, 22784, 22587, 22586, 22767, 22666, 22548, 22549, 22649 + 251, 22783, 22588, 22587, 22784, 22665, 22547, 22548, 22666 + 252, 22782, 22589, 22588, 22783, 22664, 22546, 22547, 22665 + 253, 22781, 22590, 22589, 22782, 22663, 22545, 22546, 22664 + 254, 22780, 22591, 22590, 22781, 22662, 22544, 22545, 22663 + 255, 22779, 22592, 22591, 22780, 22661, 22543, 22544, 22662 + 256, 22778, 22593, 22592, 22779, 22660, 22542, 22543, 22661 + 257, 22777, 22594, 22593, 22778, 22659, 22541, 22542, 22660 + 258, 22776, 22595, 22594, 22777, 22658, 22540, 22541, 22659 + 259, 22804, 22596, 22595, 22776, 22686, 22539, 22540, 22658 + 260, 22598, 22597, 22618, 22796, 22563, 22551, 22575, 22678 + 261, 22599, 22598, 22796, 22795, 22562, 22563, 22678, 22677 + 262, 22600, 22599, 22795, 22794, 22561, 22562, 22677, 22676 + 263, 22601, 22600, 22794, 22793, 22560, 22561, 22676, 22675 + 264, 22602, 22601, 22793, 22792, 22559, 22560, 22675, 22674 + 265, 22603, 22602, 22792, 22791, 22558, 22559, 22674, 22673 + 266, 22604, 22603, 22791, 22790, 22557, 22558, 22673, 22672 + 267, 22605, 22604, 22790, 22789, 22556, 22557, 22672, 22671 + 268, 22606, 22605, 22789, 22788, 22555, 22556, 22671, 22670 + 269, 22607, 22606, 22788, 22787, 22554, 22555, 22670, 22669 + 270, 22608, 22607, 22787, 22786, 22553, 22554, 22669, 22668 + 271, 22609, 22608, 22786, 22785, 22552, 22553, 22668, 22667 + 272, 22610, 22584, 22596, 22804, 22583, 22538, 22539, 22686 + 273, 22611, 22610, 22804, 22803, 22582, 22583, 22686, 22685 + 274, 22612, 22611, 22803, 22802, 22581, 22582, 22685, 22684 + 275, 22613, 22612, 22802, 22801, 22580, 22581, 22684, 22683 + 276, 22614, 22613, 22801, 22800, 22579, 22580, 22683, 22682 + 277, 22615, 22614, 22800, 22799, 22578, 22579, 22682, 22681 + 278, 22616, 22615, 22799, 22798, 22577, 22578, 22681, 22680 + 279, 22617, 22616, 22798, 22797, 22576, 22577, 22680, 22679 + 280, 22618, 22617, 22797, 22796, 22575, 22576, 22679, 22678 + 281, 22803, 22804, 22776, 22813, 22685, 22686, 22658, 22627 + 282, 22802, 22803, 22813, 22834, 22684, 22685, 22627, 22648 + 283, 22801, 22802, 22834, 22833, 22683, 22684, 22648, 22647 + 284, 22800, 22801, 22833, 22832, 22682, 22683, 22647, 22646 + 285, 22799, 22800, 22832, 22831, 22681, 22682, 22646, 22645 + 286, 22798, 22799, 22831, 22830, 22680, 22681, 22645, 22644 + 287, 22797, 22798, 22830, 22829, 22679, 22680, 22644, 22643 + 288, 22829, 22795, 22796, 22797, 22643, 22677, 22678, 22679 + 289, 22828, 22794, 22795, 22829, 22642, 22676, 22677, 22643 + 290, 22827, 22793, 22794, 22828, 22641, 22675, 22676, 22642 + 291, 22826, 22792, 22793, 22827, 22640, 22674, 22675, 22641 + 292, 22825, 22791, 22792, 22826, 22639, 22673, 22674, 22640 + 293, 22824, 22790, 22791, 22825, 22638, 22672, 22673, 22639 + 294, 22823, 22789, 22790, 22824, 22637, 22671, 22672, 22638 + 295, 22822, 22788, 22789, 22823, 22636, 22670, 22671, 22637 + 296, 22821, 22787, 22788, 22822, 22635, 22669, 22670, 22636 + 297, 22811, 22786, 22787, 22821, 22625, 22668, 22669, 22635 + 298, 22775, 22785, 22786, 22811, 22657, 22667, 22668, 22625 + 299, 22784, 22767, 22806, 22820, 22666, 22649, 22620, 22634 + 300, 22783, 22784, 22820, 22819, 22665, 22666, 22634, 22633 + 301, 22782, 22783, 22819, 22818, 22664, 22665, 22633, 22632 + 302, 22781, 22782, 22818, 22817, 22663, 22664, 22632, 22631 + 303, 22780, 22781, 22817, 22816, 22662, 22663, 22631, 22630 + 304, 22779, 22780, 22816, 22815, 22661, 22662, 22630, 22629 + 305, 22778, 22779, 22815, 22814, 22660, 22661, 22629, 22628 + 306, 22777, 22778, 22814, 22812, 22659, 22660, 22628, 22626 + 307, 22776, 22777, 22812, 22813, 22658, 22659, 22626, 22627 + 308, 22810, 22774, 22775, 22811, 22624, 22656, 22657, 22625 + 309, 22809, 22773, 22774, 22810, 22623, 22655, 22656, 22624 + 310, 22808, 22772, 22773, 22809, 22622, 22654, 22655, 22623 + 311, 22807, 22771, 22772, 22808, 22621, 22653, 22654, 22622 + 312, 22805, 22770, 22771, 22807, 22619, 22652, 22653, 22621 + 313, 22806, 22769, 22770, 22805, 22620, 22651, 22652, 22619 + 314, 22767, 22768, 22769, 22806, 22649, 22650, 22651, 22620