From 0625ca788d811afd5d195595ed50cf72f74c1ace Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Mon, 6 May 2019 17:48:19 -0400 Subject: [PATCH] Modernize step-62 --- examples/step-62/step-62.cc | 338 +++++++++++++++++------------------- 1 file changed, 155 insertions(+), 183 deletions(-) diff --git a/examples/step-62/step-62.cc b/examples/step-62/step-62.cc index 07930b0a7f..c87d257c6f 100644 --- a/examples/step-62/step-62.cc +++ b/examples/step-62/step-62.cc @@ -477,9 +477,7 @@ namespace step62 } } - return std::complex(1, - std::max(calculated_pml_x_coeff, - calculated_pml_y_coeff)); + return {1, std::max(calculated_pml_x_coeff, calculated_pml_y_coeff)}; } @@ -665,14 +663,12 @@ namespace step62 , dimension_y(data.get_attribute("dimension_y")) , nb_probe_points(data.get_attribute("nb_probe_points")) , grid_level(data.get_attribute("grid_level")) - , probe_start_point( - Point(data.get_attribute("probe_pos_x"), - data.get_attribute("probe_pos_y") - - data.get_attribute("probe_width_y") / 2)) - , probe_stop_point( - Point(data.get_attribute("probe_pos_x"), - data.get_attribute("probe_pos_y") + - data.get_attribute("probe_width_y") / 2)) + , probe_start_point(data.get_attribute("probe_pos_x"), + data.get_attribute("probe_pos_y") - + data.get_attribute("probe_width_y") / 2) + , probe_stop_point(data.get_attribute("probe_pos_x"), + data.get_attribute("probe_pos_y") + + data.get_attribute("probe_width_y") / 2) , right_hand_side(data) , pml(data) , rho(data) @@ -822,179 +818,158 @@ namespace step62 const FEValuesExtractors::Vector displacement(0); for (const auto &cell : dof_handler.active_cell_iterators()) - { - if (cell->is_locally_owned()) - { - cell_matrix = 0; - cell_rhs = 0; - - // We have to calculate the values of the right hand side, rho and - // the PML only if we are going to calculate the mass and the - // stiffness matrices. Otherwise we can skip this calculation which - // considerably reduces the total calculation time. - if (calculate_quadrature_data) - { - fe_values.reinit(cell); - - parameters.right_hand_side.vector_value_list( - fe_values.get_quadrature_points(), rhs_values); - parameters.rho.value_list(fe_values.get_quadrature_points(), - rho_values); - parameters.pml.vector_value_list( - fe_values.get_quadrature_points(), pml_values); - } - - // We have done this in step-18. Get a pointer to the quadrature - // cache data local to the present cell, and, as a defensive - // measure, make sure that this pointer is within the bounds of the - // global array: - QuadratureCache *local_quadrature_points_data = - reinterpret_cast *>(cell->user_pointer()); - Assert(local_quadrature_points_data >= &quadrature_cache.front(), - ExcInternalError()); - Assert(local_quadrature_points_data <= &quadrature_cache.back(), - ExcInternalError()); - for (unsigned int q = 0; q < n_q_points; ++q) - { - // The quadrature_data variable is used to store the mass and - // stiffness matrices, the right hand side vector and the value - // of `JxW`. - QuadratureCache &quadrature_data = - local_quadrature_points_data[q]; - - // Below we declare the force vector and the parameters of the - // PML $s$ and $\xi$. - Tensor<1, dim> force; - Tensor<1, dim, std::complex> s; - std::complex xi(1, 0); - - // The following block is calculated only in the first frequency - // step. - if (calculate_quadrature_data) - { - // Store the value of `JxW`. - quadrature_data.JxW = fe_values.JxW(q); - - for (unsigned int component = 0; component < dim; - ++component) - { - // Convert vectors to tensors and calculate xi - force[component] = rhs_values[q][component]; - s[component] = pml_values[q][component]; - xi *= s[component]; - } - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - const Tensor<1, dim> phi_i = - fe_values[displacement].value(i, q); - const Tensor<2, dim> grad_phi_i = - fe_values[displacement].gradient(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 Tensor<2, dim> grad_phi_j = - fe_values[displacement].gradient(j, q); - - // calculate the values of the mass matrix. - quadrature_data.mass_coefficient[i][j] = - rho_values[q] * xi * phi_i * phi_j; - - // Loop over the $mnkl$ indices of the stiffness - // tensor. - std::complex stiffness_coefficient = 0; - for (unsigned int m = 0; m < dim; ++m) - { - for (unsigned int n = 0; n < dim; ++n) + if (cell->is_locally_owned()) + { + cell_matrix = 0; + cell_rhs = 0; + + // We have to calculate the values of the right hand side, rho and + // the PML only if we are going to calculate the mass and the + // stiffness matrices. Otherwise we can skip this calculation which + // considerably reduces the total calculation time. + if (calculate_quadrature_data) + { + fe_values.reinit(cell); + + parameters.right_hand_side.vector_value_list( + fe_values.get_quadrature_points(), rhs_values); + parameters.rho.value_list(fe_values.get_quadrature_points(), + rho_values); + parameters.pml.vector_value_list( + fe_values.get_quadrature_points(), pml_values); + } + + // We have done this in step-18. Get a pointer to the quadrature + // cache data local to the present cell, and, as a defensive + // measure, make sure that this pointer is within the bounds of the + // global array: + QuadratureCache *local_quadrature_points_data = + reinterpret_cast *>(cell->user_pointer()); + Assert(local_quadrature_points_data >= &quadrature_cache.front(), + ExcInternalError()); + Assert(local_quadrature_points_data <= &quadrature_cache.back(), + ExcInternalError()); + for (unsigned int q = 0; q < n_q_points; ++q) + { + // The quadrature_data variable is used to store the mass and + // stiffness matrices, the right hand side vector and the value + // of `JxW`. + QuadratureCache &quadrature_data = + local_quadrature_points_data[q]; + + // Below we declare the force vector and the parameters of the + // PML $s$ and $\xi$. + Tensor<1, dim> force; + Tensor<1, dim, std::complex> s; + std::complex xi(1, 0); + + // The following block is calculated only in the first frequency + // step. + if (calculate_quadrature_data) + { + // Store the value of `JxW`. + quadrature_data.JxW = fe_values.JxW(q); + + for (unsigned int component = 0; component < dim; ++component) + { + // Convert vectors to tensors and calculate xi + force[component] = rhs_values[q][component]; + s[component] = pml_values[q][component]; + xi *= s[component]; + } + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const Tensor<1, dim> phi_i = + fe_values[displacement].value(i, q); + const Tensor<2, dim> grad_phi_i = + fe_values[displacement].gradient(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 Tensor<2, dim> grad_phi_j = + fe_values[displacement].gradient(j, q); + + // calculate the values of the mass matrix. + quadrature_data.mass_coefficient[i][j] = + rho_values[q] * xi * phi_i * phi_j; + + // Loop over the $mnkl$ indices of the stiffness + // tensor. + std::complex stiffness_coefficient = 0; + for (unsigned int m = 0; m < dim; ++m) + for (unsigned int n = 0; n < dim; ++n) + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = 0; l < dim; ++l) { - for (unsigned int k = 0; k < dim; ++k) - { - for (unsigned int l = 0; l < dim; ++l) - { - // Here we calculate the tensors - // $\alpha_{mnkl}$ and - // $\beta_{mnkl}$ - const std::complex alpha = - xi * - stiffness_tensor[m][n][k][l] / - (2.0 * s[n] * s[k]); - const std::complex beta = - xi * - stiffness_tensor[m][n][k][l] / - (2.0 * s[n] * s[l]); - - // Here we calculate the stiffness - // matrix. Note that the stiffness - // matrix is not symmetric because - // of the PMLs. We use the gradient - // function (see the - // [documentation](https://www.dealii.org/current/doxygen/deal.II/group__vector__valued.html)) - // which is a - // Tensor@<2,dim@>. - // The matrix $G_{ij}$ - // consists of entries - // @f[ - // G_{ij}= - // \frac{\partial\phi_i}{\partial - // x_j} - // =\partial_j \phi_i - // @f] - // Note the position of the indices - // $i$ and $j$ and the notation that - // we use in this tutorial: - // $\partial_j\phi_i$. As the - // stiffness tensor is not - // symmetric, it is very easy to - // make a mistake. - stiffness_coefficient += - grad_phi_i[m][n] * - (alpha * grad_phi_j[l][k] + - beta * grad_phi_j[k][l]); - } - } + // Here we calculate the tensors + // $\alpha_{mnkl}$ and $\beta_{mnkl}$ + const std::complex alpha = + xi * stiffness_tensor[m][n][k][l] / + (2.0 * s[n] * s[k]); + const std::complex beta = + xi * stiffness_tensor[m][n][k][l] / + (2.0 * s[n] * s[l]); + + // Here we calculate the stiffness matrix. + // Note that the stiffness matrix is not + // symmetric because of the PMLs. We use the + // gradient function (see the + // [documentation](https://www.dealii.org/current/doxygen/deal.II/group__vector__valued.html)) + // which is a Tensor@<2,dim@>. + // The matrix $G_{ij}$ consists of entries + // @f[ + // G_{ij}= + // \frac{\partial\phi_i}{\partial x_j} + // =\partial_j \phi_i + // @f] + // Note the position of the indices $i$ and + // $j$ and the notation that we use in this + // tutorial: $\partial_j\phi_i$. As the + // stiffness tensor is not symmetric, it is + // very easy to make a mistake. + stiffness_coefficient += + grad_phi_i[m][n] * + (alpha * grad_phi_j[l][k] + + beta * grad_phi_j[k][l]); } - } - - // We save the value of the stiffness matrix in - // quadrature_data - quadrature_data.stiffness_coefficient[i][j] = - stiffness_coefficient; - } - - // and the value of the right hand side in - // quadrature_data. - quadrature_data.right_hand_side[i] = - phi_i * force * fe_values.JxW(q); - } - } - // We loop again over the degrees of freedom of the cells to - // calculate the system matrix. These loops are really quick - // because we have already calculated the stiffness and mass - // matrices, only the value of $\omega$ changes. - for (unsigned int i = 0; i < dofs_per_cell; ++i) + // We save the value of the stiffness matrix in + // quadrature_data + quadrature_data.stiffness_coefficient[i][j] = + stiffness_coefficient; + } + + // and the value of the right hand side in + // quadrature_data. + quadrature_data.right_hand_side[i] = + phi_i * force * fe_values.JxW(q); + } + } + + // We loop again over the degrees of freedom of the cells to + // calculate the system matrix. These loops are really quick + // because we have already calculated the stiffness and mass + // matrices, only the value of $\omega$ changes. + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) { - for (unsigned int j = 0; j < dofs_per_cell; ++j) - { - std::complex matrix_sum = 0; - matrix_sum += -std::pow(omega, 2) * - quadrature_data.mass_coefficient[i][j]; - matrix_sum += - quadrature_data.stiffness_coefficient[i][j]; - cell_matrix(i, j) += matrix_sum * quadrature_data.JxW; - } - cell_rhs(i) += quadrature_data.right_hand_side[i]; + std::complex matrix_sum = 0; + matrix_sum += -std::pow(omega, 2) * + quadrature_data.mass_coefficient[i][j]; + matrix_sum += quadrature_data.stiffness_coefficient[i][j]; + cell_matrix(i, j) += matrix_sum * quadrature_data.JxW; } - } - cell->get_dof_indices(local_dof_indices); - constraints.distribute_local_to_global(cell_matrix, - cell_rhs, - local_dof_indices, - system_matrix, - system_rhs); - } - } + cell_rhs(i) += quadrature_data.right_hand_side[i]; + } + 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); @@ -1251,10 +1226,7 @@ namespace step62 quadrature_formula.size(), QuadratureCache(fe.dofs_per_cell)); unsigned int cache_index = 0; - for (typename Triangulation::active_cell_iterator cell = - triangulation.begin_active(); - cell != triangulation.end(); - ++cell) + for (const auto &cell : triangulation.active_cell_iterators()) if (cell->is_locally_owned()) { cell->set_user_pointer(&quadrature_cache[cache_index]); -- 2.39.5