1021 * The speed of sound is defined by
1023 * c = \frac{K_e}{\rho}
1025 * where @f$K_e@f$ is the effective elastic
constant and @f$\rho@f$ the density.
1026 * Here we consider the
case in which the waveguide width is much smaller
1027 * than the wavelength. In
this case it can be shown that
for the two
1032 * and
for the three dimensional
case @f$K_e@f$ is
equal to the Young
's modulus.
1034 * K_e = \mu\frac{3\lambda +2\mu}{\lambda+\mu}
1038 * double elastic_constant;
1041 * elastic_constant = 4 * mu * (lambda + mu) / (lambda + 2 * mu);
1043 * else if (dim == 3)
1045 * elastic_constant = mu * (3 * lambda + 2 * mu) / (lambda + mu);
1049 * Assert(false, ExcInternalError());
1051 * const double material_a_speed_of_sound =
1052 * std::sqrt(elastic_constant / material_a_rho);
1053 * const double material_a_wavelength =
1054 * material_a_speed_of_sound / cavity_resonance_frequency;
1055 * const double material_b_speed_of_sound =
1056 * std::sqrt(elastic_constant / material_b_rho);
1057 * const double material_b_wavelength =
1058 * material_b_speed_of_sound / cavity_resonance_frequency;
1062 * The density @f$\rho@f$ takes the following form
1063 * <img alt="Phononic superlattice cavity"
1064 * src="https://www.dealii.org/images/steps/developer/step-62.04.svg"
1066 * where the brown color represents material_a and the green color
1067 * represents material_b.
1070 * for (unsigned int idx = 0; idx < nb_mirror_pairs; ++idx)
1072 * const double layer_transition_center =
1073 * material_a_wavelength / 2 +
1074 * idx * (material_b_wavelength / 4 + material_a_wavelength / 4);
1075 * if (std::abs(p[0]) >=
1076 * (layer_transition_center - average_rho_width / 2) &&
1077 * std::abs(p[0]) <= (layer_transition_center + average_rho_width / 2))
1079 * const double coefficient =
1081 * (layer_transition_center - average_rho_width / 2)) /
1082 * average_rho_width;
1083 * return (1 - coefficient) * material_a_rho +
1084 * coefficient * material_b_rho;
1090 * Here we define the
1092 * smoothing](https://meep.readthedocs.io/en/latest/Subpixel_Smoothing/)
1093 * which improves the precision of the simulation.
1096 * for (unsigned int idx = 0; idx < nb_mirror_pairs; ++idx)
1098 * const double layer_transition_center =
1099 * material_a_wavelength / 2 +
1100 * idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1101 * material_b_wavelength / 4;
1102 * if (std::abs(p[0]) >=
1103 * (layer_transition_center - average_rho_width / 2) &&
1104 * std::abs(p[0]) <= (layer_transition_center + average_rho_width / 2))
1106 * const double coefficient =
1108 * (layer_transition_center - average_rho_width / 2)) /
1109 * average_rho_width;
1110 * return (1 - coefficient) * material_b_rho +
1111 * coefficient * material_a_rho;
1120 * if (std::abs(p[0]) <= material_a_wavelength / 2)
1122 * return material_a_rho;
1127 * the material_a layers
1130 * for (unsigned int idx = 0; idx < nb_mirror_pairs; ++idx)
1132 * const double layer_center =
1133 * material_a_wavelength / 2 +
1134 * idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1135 * material_b_wavelength / 4 + material_a_wavelength / 8;
1136 * const double layer_width = material_a_wavelength / 4;
1137 * if (std::abs(p[0]) >= (layer_center - layer_width / 2) &&
1138 * std::abs(p[0]) <= (layer_center + layer_width / 2))
1140 * return material_a_rho;
1146 * the material_b layers
1149 * for (unsigned int idx = 0; idx < nb_mirror_pairs; ++idx)
1151 * const double layer_center =
1152 * material_a_wavelength / 2 +
1153 * idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1154 * material_b_wavelength / 8;
1155 * const double layer_width = material_b_wavelength / 4;
1156 * if (std::abs(p[0]) >= (layer_center - layer_width / 2) &&
1157 * std::abs(p[0]) <= (layer_center + layer_width / 2))
1159 * return material_b_rho;
1165 * and finally the default is material_a.
1168 * return material_a_rho;
1176 * <a name="TheParametersclassimplementation"></a>
1177 * <h4>The `Parameters` class implementation</h4>
1181 * The constructor reads all the parameters from the HDF5::Group `data` using
1182 * the HDF5::Group::get_attribute() function.
1185 * template <int dim>
1186 * Parameters<dim>::Parameters(HDF5::Group &data)
1188 * , simulation_name(data.get_attribute<std::string>("simulation_name"))
1189 * , save_vtu_files(data.get_attribute<bool>("save_vtu_files"))
1190 * , start_frequency(data.get_attribute<double>("start_frequency"))
1191 * , stop_frequency(data.get_attribute<double>("stop_frequency"))
1192 * , nb_frequency_points(data.get_attribute<int>("nb_frequency_points"))
1193 * , lambda(data.get_attribute<double>("lambda"))
1194 * , mu(data.get_attribute<double>("mu"))
1195 * , dimension_x(data.get_attribute<double>("dimension_x"))
1196 * , dimension_y(data.get_attribute<double>("dimension_y"))
1197 * , nb_probe_points(data.get_attribute<int>("nb_probe_points"))
1198 * , grid_level(data.get_attribute<int>("grid_level"))
1199 * , probe_start_point(data.get_attribute<double>("probe_pos_x"),
1200 * data.get_attribute<double>("probe_pos_y") -
1201 * data.get_attribute<double>("probe_width_y") / 2)
1202 * , probe_stop_point(data.get_attribute<double>("probe_pos_x"),
1203 * data.get_attribute<double>("probe_pos_y") +
1204 * data.get_attribute<double>("probe_width_y") / 2)
1205 * , right_hand_side(data)
1215 * <a name="TheQuadratureCacheclassimplementation"></a>
1216 * <h4>The `QuadratureCache` class implementation</h4>
1220 * We need to reserve enough space for the mass and stiffness matrices and the
1221 * right hand side vector.
1224 * template <int dim>
1225 * QuadratureCache<dim>::QuadratureCache(const unsigned int dofs_per_cell)
1226 * : dofs_per_cell(dofs_per_cell)
1227 * , mass_coefficient(dofs_per_cell, dofs_per_cell)
1228 * , stiffness_coefficient(dofs_per_cell, dofs_per_cell)
1229 * , right_hand_side(dofs_per_cell)
1237 * <a name="ImplementationoftheElasticWaveclass"></a>
1238 * <h3>Implementation of the `ElasticWave` class</h3>
1243 * <a name="Constructor"></a>
1244 * <h4>Constructor</h4>
1248 * This is very similar to the constructor of @ref step_40 "step-40". In addition we create
1249 * the HDF5 datasets `frequency_dataset`, `position_dataset` and
1250 * `displacement`. Note the use of the `template` keyword for the creation of
1251 * the HDF5 datasets. It is a C++ requirement to use the `template` keyword in
1252 * order to treat `create_dataset` as a dependent template name.
1255 * template <int dim>
1256 * ElasticWave<dim>::ElasticWave(const Parameters<dim> ¶meters)
1257 * : parameters(parameters)
1258 * , mpi_communicator(MPI_COMM_WORLD)
1259 * , triangulation(mpi_communicator,
1260 * typename Triangulation<dim>::MeshSmoothing(
1261 * Triangulation<dim>::smoothing_on_refinement |
1262 * Triangulation<dim>::smoothing_on_coarsening))
1263 * , quadrature_formula(2)
1264 * , fe(FE_Q<dim>(1), dim)
1265 * , dof_handler(triangulation)
1266 * , frequency(parameters.nb_frequency_points)
1267 * , probe_positions(parameters.nb_probe_points, dim)
1268 * , frequency_dataset(parameters.data.template create_dataset<double>(
1270 * std::vector<hsize_t>{parameters.nb_frequency_points}))
1271 * , probe_positions_dataset(parameters.data.template create_dataset<double>(
1273 * std::vector<hsize_t>{parameters.nb_probe_points, dim}))
1275 * parameters.data.template create_dataset<std::complex<double>>(
1277 * std::vector<hsize_t>{parameters.nb_probe_points,
1278 * parameters.nb_frequency_points}))
1279 * , pcout(std::cout,
1280 * (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
1281 * , computing_timer(mpi_communicator,
1283 * TimerOutput::never,
1284 * TimerOutput::wall_times)
1292 * <a name="ElasticWavesetup_system"></a>
1293 * <h4>ElasticWave::setup_system</h4>
1297 * There is nothing new in this function, the only difference with @ref step_40 "step-40" is
1298 * that we don't have to
apply boundary conditions because we use the PMLs to
1299 * truncate the domain.
1302 *
template <
int dim>
1303 *
void ElasticWave<dim>::setup_system()
1307 * dof_handler.distribute_dofs(fe);
1309 * locally_owned_dofs = dof_handler.locally_owned_dofs();
1310 * locally_relevant_dofs =
1313 * locally_relevant_solution.reinit(locally_owned_dofs,
1314 * locally_relevant_dofs,
1315 * mpi_communicator);
1317 * system_rhs.reinit(locally_owned_dofs, mpi_communicator);
1319 * constraints.clear();
1320 * constraints.reinit(locally_relevant_dofs);
1323 * constraints.close();
1329 * locally_owned_dofs,
1331 * locally_relevant_dofs);
1333 * system_matrix.reinit(locally_owned_dofs,
1334 * locally_owned_dofs,
1336 * mpi_communicator);
1344 * <a name=
"ElasticWaveassemble_system"></a>
1345 * <h4>ElasticWave::assemble_system</h4>
1349 * This function is also very similar to @ref step_40
"step-40", though there are notable
1350 * differences. We
assemble the system
for each frequency/omega step. In the
1351 *
first step we
set `calculate_quadrature_data = True` and we calculate the
1352 * mass and stiffness matrices and the right hand side vector. In the
1353 * subsequent steps we will use that data to accelerate the calculation.
1356 *
template <
int dim>
1357 *
void ElasticWave<dim>::assemble_system(
const double omega,
1358 *
const bool calculate_quadrature_data)
1363 * quadrature_formula,
1366 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1367 *
const unsigned int n_q_points = quadrature_formula.size();
1372 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1376 * Here we store the
value of the right hand side, rho and the PML.
1379 * std::vector<Vector<double>> rhs_values(n_q_points,
Vector<double>(dim));
1380 * std::vector<double> rho_values(n_q_points);
1381 * std::vector<Vector<std::complex<double>>> pml_values(
1382 * n_q_points,
Vector<std::complex<double>>(dim));
1386 * We calculate the stiffness tensor
for the @f$\lambda@f$ and @f$\mu@f$ that have
1387 * been defined in the jupyter notebook. Note that contrary to @f$\rho@f$ the
1388 * stiffness is
constant among
for the whole domain.
1392 * get_stiffness_tensor<dim>(parameters.lambda, parameters.mu);
1396 * We use the same method of @ref step_20
"step-20" for vector-valued problems.
1401 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1402 * if (cell->is_locally_owned())
1409 * We have to calculate the
values of the right hand side, rho and
1410 * the PML only
if we are going to calculate the mass and the
1411 * stiffness matrices. Otherwise we can skip
this calculation which
1412 * considerably reduces the total calculation time.
1415 *
if (calculate_quadrature_data)
1417 * fe_values.reinit(cell);
1419 * parameters.right_hand_side.vector_value_list(
1420 * fe_values.get_quadrature_points(), rhs_values);
1421 * parameters.rho.value_list(fe_values.get_quadrature_points(),
1423 * parameters.pml.vector_value_list(
1424 * fe_values.get_quadrature_points(), pml_values);
1429 * We have done
this in @ref step_18
"step-18". Get a pointer to the quadrature
1430 * cache data local to the present cell, and, as a defensive
1431 * measure, make sure that
this pointer is within the bounds of the
1435 * QuadratureCache<dim> *local_quadrature_points_data =
1436 *
reinterpret_cast<QuadratureCache<dim> *
>(cell->user_pointer());
1437 *
Assert(local_quadrature_points_data >= &quadrature_cache.front(),
1438 * ExcInternalError());
1439 *
Assert(local_quadrature_points_data <= &quadrature_cache.back(),
1440 * ExcInternalError());
1441 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1445 * The quadrature_data variable is used to store the mass and
1446 * stiffness matrices, the right hand side vector and the
value
1450 * QuadratureCache<dim> &quadrature_data =
1451 * local_quadrature_points_data[q];
1455 * Below we declare the force vector and the parameters of the
1456 * PML @f$s@f$ and @f$\xi@f$.
1461 * std::complex<double> xi(1, 0);
1465 * The following block is calculated only in the
first frequency
1469 *
if (calculate_quadrature_data)
1473 * Store the
value of `JxW`.
1476 * quadrature_data.JxW = fe_values.JxW(q);
1478 *
for (
unsigned int component = 0; component < dim; ++component)
1482 * Convert vectors to tensors and calculate xi
1485 * force[component] = rhs_values[q][component];
1486 * s[component] = pml_values[q][component];
1487 * xi *= s[component];
1492 * Here we calculate the @f$\alpha_{mnkl}@f$ and @f$\beta_{mnkl}@f$
1498 *
for (
unsigned int m = 0; m < dim; ++m)
1499 *
for (
unsigned int n = 0; n < dim; ++n)
1500 *
for (
unsigned int k = 0; k < dim; ++k)
1501 *
for (
unsigned int l = 0;
l < dim; ++
l)
1503 * alpha[m][n][k][
l] = xi *
1504 * stiffness_tensor[m][n][k][
l] /
1505 * (2.0 * s[n] * s[k]);
1506 * beta[m][n][k][
l] = xi *
1507 * stiffness_tensor[m][n][k][
l] /
1508 * (2.0 * s[n] * s[l]);
1511 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1514 * fe_values[displacement].value(i, q);
1516 * fe_values[displacement].gradient(i, q);
1518 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1521 * fe_values[displacement].value(j, q);
1523 * fe_values[displacement].gradient(j, q);
1527 * calculate the
values of the @ref GlossMassMatrix
"mass matrix".
1530 * quadrature_data.mass_coefficient[i][j] =
1531 * rho_values[q] * xi * phi_i * phi_j;
1535 * Loop over the @f$mnkl@f$ indices of the stiffness
1539 * std::complex<double> stiffness_coefficient = 0;
1540 *
for (
unsigned int m = 0; m < dim; ++m)
1541 *
for (
unsigned int n = 0; n < dim; ++n)
1542 *
for (
unsigned int k = 0; k < dim; ++k)
1543 *
for (
unsigned int l = 0;
l < dim; ++
l)
1547 * Here we calculate the stiffness
matrix.
1548 * Note that the stiffness
matrix is not
1549 *
symmetric because of the PMLs. We use the
1551 * [documentation](https:
1552 * which is a <code>
Tensor@<2,dim@></code>.
1553 * The
matrix @f$G_{ij}@f$ consists of entries
1556 * \frac{\partial\phi_i}{\partial x_j}
1557 * =\partial_j \phi_i
1559 * Note the position of the indices @f$i@f$ and
1560 * @f$j@f$ and the notation that we use in
this
1561 * tutorial: @f$\partial_j\phi_i@f$. As the
1562 * stiffness tensor is not
symmetric, it is
1563 * very easy to make a mistake.
1566 * stiffness_coefficient +=
1567 * grad_phi_i[m][n] *
1568 * (alpha[m][n][k][l] * grad_phi_j[l][k] +
1569 * beta[m][n][k][l] * grad_phi_j[k][l]);
1578 * quadrature_data.stiffness_coefficient[i][j] =
1579 * stiffness_coefficient;
1584 * and the
value of the right hand side in
1588 * quadrature_data.right_hand_side[i] =
1589 * phi_i * force * fe_values.JxW(q);
1595 * We
loop again over the degrees of freedom of the cells to
1596 * calculate the system
matrix. These loops are really quick
1597 * because we have already calculated the stiffness and mass
1598 * matrices, only the
value of @f$\omega@f$ changes.
1601 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1603 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1605 * std::complex<double> matrix_sum = 0;
1606 * matrix_sum += -
std::pow(omega, 2) *
1607 * quadrature_data.mass_coefficient[i][j];
1608 * matrix_sum += quadrature_data.stiffness_coefficient[i][j];
1609 *
cell_matrix(i, j) += matrix_sum * quadrature_data.JxW;
1611 * cell_rhs(i) += quadrature_data.right_hand_side[i];
1614 * cell->get_dof_indices(local_dof_indices);
1615 * constraints.distribute_local_to_global(cell_matrix,
1617 * local_dof_indices,
1629 * <a name=
"ElasticWavesolve"></a>
1630 * <h4>ElasticWave::solve</h4>
1634 * This is even more simple than in @ref step_40
"step-40". We use the
parallel direct solver
1635 * MUMPS which
requires less options than an iterative solver. The drawback is
1636 * that it does not
scale very well. It is not straightforward to solve the
1637 * Helmholtz equation with an iterative solver. The shifted Laplacian
1638 * multigrid method is a well known approach to precondition
this system, but
1639 *
this is beyond the scope of
this tutorial.
1642 *
template <
int dim>
1643 *
void ElasticWave<dim>::solve()
1647 * locally_owned_dofs, mpi_communicator);
1651 * solver.solve(system_matrix, completely_distributed_solution, system_rhs);
1653 * pcout <<
" Solved in " << solver_control.last_step() <<
" iterations."
1655 * constraints.distribute(completely_distributed_solution);
1656 * locally_relevant_solution = completely_distributed_solution;
1662 * <a name=
"ElasticWaveinitialize_position_vector"></a>
1663 * <h4>ElasticWave::initialize_position_vector</h4>
1667 * We use
this function to calculate the
values of the position vector.
1670 *
template <
int dim>
1671 *
void ElasticWave<dim>::initialize_probe_positions_vector()
1673 *
for (
unsigned int position_idx = 0;
1674 * position_idx < parameters.nb_probe_points;
1679 * Because of the way the
operator + and - are overloaded to subtract
1680 * two points, the following has to be done:
1681 * `Point_b<dim> + (-Point_a<dim>)`
1685 * (position_idx / ((
double)(parameters.nb_probe_points - 1))) *
1686 * (parameters.probe_stop_point + (-parameters.probe_start_point)) +
1687 * parameters.probe_start_point;
1688 * probe_positions[position_idx][0] = p[0];
1689 * probe_positions[position_idx][1] = p[1];
1692 * probe_positions[position_idx][2] = p[2];
1700 * <a name=
"ElasticWavestore_frequency_step_data"></a>
1701 * <h4>ElasticWave::store_frequency_step_data</h4>
1705 * This function stores in the
HDF5 file the measured energy by the probe.
1708 *
template <
int dim>
1710 * ElasticWave<dim>::store_frequency_step_data(
const unsigned int frequency_idx)
1716 * We store the displacement in the @f$x@f$ direction; the displacement in the
1717 * @f$y@f$ direction is negligible.
1720 *
const unsigned int probe_displacement_component = 0;
1724 * The vector coordinates contains the coordinates in the
HDF5 file of the
1725 * points of the probe that are located in locally owned cells. The vector
1726 * displacement_data contains the
value of the displacement at these points.
1729 * std::vector<hsize_t> coordinates;
1730 * std::vector<std::complex<double>> displacement_data;
1735 * std::vector<bool> marked_vertices = {};
1736 *
const double tolerance = 1.e-10;
1738 *
for (
unsigned int position_idx = 0;
1739 * position_idx < parameters.nb_probe_points;
1743 *
for (
unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1745 *
point[dim_idx] = probe_positions[position_idx][dim_idx];
1747 *
bool point_in_locally_owned_cell =
false;
1750 * cache, point, cell_hint, marked_vertices, tolerance);
1753 * cell_hint = cell_and_ref_point.first;
1754 * point_in_locally_owned_cell =
1755 * cell_and_ref_point.first->is_locally_owned();
1758 *
if (point_in_locally_owned_cell)
1762 * Then we can store the
values of the displacement in the points of
1763 * the probe in `displacement_data`.
1768 * locally_relevant_solution,
1771 * coordinates.emplace_back(position_idx);
1772 * coordinates.emplace_back(frequency_idx);
1773 * displacement_data.emplace_back(
1774 * tmp_vector(probe_displacement_component));
1780 * We write the displacement data in the
HDF5 file. The
call
1782 * the processes have to participate.
1785 * if (coordinates.size() > 0)
1787 * displacement.write_selection(displacement_data, coordinates);
1791 * Therefore even
if the process has no data to write it has to participate
1793 * Note that we have to specify the data type, in
this case
1794 * `std::complex<double>`.
1799 * displacement.write_none<std::complex<double>>();
1804 * If the variable `save_vtu_files` in the input file equals `True` then all
1805 * the data will be saved as
vtu. The procedure to write `
vtu` files has
1806 * been described in @ref step_40
"step-40".
1809 *
if (parameters.save_vtu_files)
1811 * std::vector<std::string> solution_names(dim,
"displacement");
1812 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
1818 * locally_relevant_solution,
1822 *
for (
unsigned int i = 0; i < subdomain.size(); ++i)
1824 * data_out.add_data_vector(subdomain,
"subdomain");
1826 * std::vector<Vector<double>> force(
1828 * std::vector<Vector<double>> pml(
1834 * if (cell->is_locally_owned())
1836 * for (unsigned
int dim_idx = 0; dim_idx < dim; ++dim_idx)
1838 * force[dim_idx](cell->active_cell_index()) =
1839 * parameters.right_hand_side.value(cell->center(), dim_idx);
1840 * pml[dim_idx](cell->active_cell_index()) =
1841 * parameters.pml.value(cell->center(), dim_idx).imag();
1843 * rho(cell->active_cell_index()) =
1844 * parameters.rho.value(cell->center());
1848 * And on the cells that we are not interested in,
set the
1849 * respective
value to a bogus
value in order to make sure that
if
1850 * we were somehow wrong about our assumption we would find out by
1851 * looking at the graphical output:
1856 *
for (
unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1858 * force[dim_idx](cell->active_cell_index()) = -1e+20;
1859 * pml[dim_idx](cell->active_cell_index()) = -1e+20;
1861 * rho(cell->active_cell_index()) = -1
e+20;
1865 *
for (
unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1867 * data_out.add_data_vector(force[dim_idx],
1868 *
"force_" + std::to_string(dim_idx));
1869 * data_out.add_data_vector(pml[dim_idx],
1870 *
"pml_" + std::to_string(dim_idx));
1872 * data_out.add_data_vector(rho,
"rho");
1874 * data_out.build_patches();
1876 * std::stringstream frequency_idx_stream;
1877 *
const unsigned int nb_number_positions =
1878 * ((
unsigned int)std::log10(parameters.nb_frequency_points)) + 1;
1879 * frequency_idx_stream << std::setw(nb_number_positions)
1880 * << std::setfill(
'0') << frequency_idx;
1881 * std::string filename = (parameters.simulation_name +
"_" +
1882 * frequency_idx_stream.str() +
".vtu");
1883 * data_out.write_vtu_in_parallel(filename.c_str(), mpi_communicator);
1892 * <a name=
"ElasticWaveoutput_results"></a>
1893 * <h4>ElasticWave::output_results</h4>
1897 * This function writes the datasets that have not already been written.
1900 *
template <
int dim>
1901 *
void ElasticWave<dim>::output_results()
1905 * The vectors `frequency` and `position` are the same
for all the
1906 * processes. Therefore any of the processes can write the corresponding
1913 * frequency_dataset.write(frequency);
1914 * probe_positions_dataset.write(probe_positions);
1918 * frequency_dataset.write_none<
double>();
1919 * probe_positions_dataset.write_none<
double>();
1928 * <a name=
"ElasticWavesetup_quadrature_cache"></a>
1929 * <h4>ElasticWave::setup_quadrature_cache</h4>
1933 * We use
this function at the beginning of our computations to
set up
initial
1934 *
values of the cache variables. This function has been described in @ref step_18
"step-18".
1935 * There are no differences with the function of @ref step_18
"step-18".
1938 *
template <
int dim>
1939 *
void ElasticWave<dim>::setup_quadrature_cache()
1944 * std::vector<QuadratureCache<dim>> tmp;
1945 * quadrature_cache.swap(tmp);
1948 * quadrature_cache.resize(
triangulation.n_locally_owned_active_cells() *
1949 * quadrature_formula.size(),
1950 * QuadratureCache<dim>(fe.n_dofs_per_cell()));
1951 *
unsigned int cache_index = 0;
1952 *
for (
const auto &cell :
triangulation.active_cell_iterators())
1953 * if (cell->is_locally_owned())
1955 * cell->set_user_pointer(&quadrature_cache[cache_index]);
1956 * cache_index += quadrature_formula.size();
1958 *
Assert(cache_index == quadrature_cache.size(), ExcInternalError());
1966 * <a name=
"ElasticWavefrequency_sweep"></a>
1967 * <h4>ElasticWave::frequency_sweep</h4>
1971 * For clarity we divide the function `
run` of @ref step_40
"step-40" into the
functions
1972 * `
run` and `frequency_sweep`. In the function `frequency_sweep` we place the
1973 * iteration over the frequency vector.
1976 *
template <
int dim>
1977 *
void ElasticWave<dim>::frequency_sweep()
1979 *
for (
unsigned int frequency_idx = 0;
1980 * frequency_idx < parameters.nb_frequency_points;
1983 * pcout << parameters.simulation_name +
" frequency idx: "
1984 * << frequency_idx <<
'/' << parameters.nb_frequency_points - 1
1990 *
if (frequency_idx == 0)
1992 * pcout <<
" Number of active cells : "
1994 * pcout <<
" Number of degrees of freedom : "
1995 * << dof_handler.n_dofs() << std::endl;
1998 *
if (frequency_idx == 0)
2002 * Write the simulation parameters only once
2005 * parameters.data.set_attribute(
"active_cells",
2007 * parameters.data.set_attribute(
"degrees_of_freedom",
2008 * dof_handler.n_dofs());
2013 * We calculate the frequency and omega
values for this particular step.
2016 *
const double current_loop_frequency =
2017 * (parameters.start_frequency +
2019 * (parameters.stop_frequency - parameters.start_frequency) /
2020 * (parameters.nb_frequency_points - 1));
2021 *
const double current_loop_omega =
2026 * In the
first frequency step we calculate the mass and stiffness
2027 * matrices and the right hand side. In the subsequent frequency steps
2028 * we will use those
values. This improves considerably the calculation
2032 * assemble_system(current_loop_omega,
2033 * (frequency_idx == 0) ?
true : false);
2036 * frequency[frequency_idx] = current_loop_frequency;
2037 * store_frequency_step_data(frequency_idx);
2039 * computing_timer.print_summary();
2040 * computing_timer.reset();
2041 * pcout << std::endl;
2050 * <a name=
"ElasticWaverun"></a>
2051 * <h4>ElasticWave::run</h4>
2055 * This function is very similar to the one in @ref step_40
"step-40".
2058 *
template <
int dim>
2059 *
void ElasticWave<dim>::run()
2062 * pcout <<
"Debug mode" << std::endl;
2064 * pcout <<
"Release mode" << std::endl;
2069 * p1(0) = -parameters.dimension_x / 2;
2070 * p1(1) = -parameters.dimension_y / 2;
2073 * p1(2) = -parameters.dimension_y / 2;
2076 * p2(0) = parameters.dimension_x / 2;
2077 * p2(1) = parameters.dimension_y / 2;
2080 * p2(2) = parameters.dimension_y / 2;
2082 * std::vector<unsigned int> divisions(dim);
2083 * divisions[0] =
int(parameters.dimension_x / parameters.dimension_y);
2097 * setup_quadrature_cache();
2099 * initialize_probe_positions_vector();
2101 * frequency_sweep();
2112 * <a name=
"Themainfunction"></a>
2113 * <h4>The main function</h4>
2117 * The main function is very similar to the one in @ref step_40
"step-40".
2120 *
int main(
int argc,
char *argv[])
2124 *
using namespace dealii;
2125 *
const unsigned int dim = 2;
2132 *
auto data = data_file.create_group(
"data");
2136 * Each of the simulations (displacement and calibration) is stored in a
2140 *
const std::array<std::string, 2> group_names{
2141 * {
"displacement",
"calibration"}};
2142 *
for (
const std::string &group_name : group_names)
2146 * For each of these two
group names, we now create the
group and put
2147 * attributes into these groups.
2148 * Specifically, these are:
2149 * - The dimensions of the waveguide (in @f$x@f$ and @f$y@f$ directions)
2150 * - The position of the probe (in @f$x@f$ and @f$y@f$ directions)
2151 * - The number of points in the probe
2152 * - The global refinement
level
2153 * - The cavity resonance frequency
2154 * - The number of mirror pairs
2155 * - The material properties
2156 * - The force parameters
2157 * - The PML parameters
2158 * - The frequency parameters
2166 *
group.set_attribute<
double>(
"dimension_x", 2
e-5);
2167 *
group.set_attribute<
double>(
"dimension_y", 2
e-8);
2168 *
group.set_attribute<
double>(
"probe_pos_x", 8
e-6);
2169 *
group.set_attribute<
double>(
"probe_pos_y", 0);
2170 *
group.set_attribute<
double>(
"probe_width_y", 2
e-08);
2171 *
group.set_attribute<
unsigned int>(
"nb_probe_points", 5);
2172 *
group.set_attribute<
unsigned int>(
"grid_level", 1);
2173 *
group.set_attribute<
double>(
"cavity_resonance_frequency", 20e9);
2174 *
group.set_attribute<
unsigned int>(
"nb_mirror_pairs", 15);
2176 *
group.set_attribute<
double>(
"poissons_ratio", 0.27);
2177 *
group.set_attribute<
double>(
"youngs_modulus", 270000000000.0);
2178 *
group.set_attribute<
double>(
"material_a_rho", 3200);
2180 *
if (group_name ==
"displacement")
2181 *
group.set_attribute<
double>(
"material_b_rho", 2000);
2183 *
group.set_attribute<
double>(
"material_b_rho", 3200);
2185 *
group.set_attribute(
2187 *
group.get_attribute<
double>(
"youngs_modulus") *
2188 *
group.get_attribute<
double>(
"poissons_ratio") /
2189 * ((1 +
group.get_attribute<
double>(
"poissons_ratio")) *
2190 * (1 - 2 *
group.get_attribute<
double>(
"poissons_ratio"))));
2191 *
group.set_attribute(
"mu",
2192 *
group.get_attribute<
double>(
"youngs_modulus") /
2193 * (2 * (1 +
group.get_attribute<
double>(
2194 *
"poissons_ratio"))));
2196 *
group.set_attribute<
double>(
"max_force_amplitude", 1e26);
2197 *
group.set_attribute<
double>(
"force_sigma_x", 1
e-7);
2198 *
group.set_attribute<
double>(
"force_sigma_y", 1);
2199 *
group.set_attribute<
double>(
"max_force_width_x", 3
e-7);
2200 *
group.set_attribute<
double>(
"max_force_width_y", 2
e-8);
2201 *
group.set_attribute<
double>(
"force_x_pos", -8
e-6);
2202 *
group.set_attribute<
double>(
"force_y_pos", 0);
2204 *
group.set_attribute<
bool>(
"pml_x",
true);
2205 *
group.set_attribute<
bool>(
"pml_y",
false);
2206 *
group.set_attribute<
double>(
"pml_width_x", 1.8e-6);
2207 *
group.set_attribute<
double>(
"pml_width_y", 5
e-7);
2208 *
group.set_attribute<
double>(
"pml_coeff", 1.6);
2209 *
group.set_attribute<
unsigned int>(
"pml_coeff_degree", 2);
2211 *
group.set_attribute<
double>(
"center_frequency", 20e9);
2212 *
group.set_attribute<
double>(
"frequency_range", 0.5e9);
2213 *
group.set_attribute<
double>(
2214 *
"start_frequency",
2215 *
group.get_attribute<
double>(
"center_frequency") -
2216 *
group.get_attribute<
double>(
"frequency_range") / 2);
2217 *
group.set_attribute<
double>(
2219 *
group.get_attribute<
double>(
"center_frequency") +
2220 *
group.get_attribute<
double>(
"frequency_range") / 2);
2221 *
group.set_attribute<
unsigned int>(
"nb_frequency_points", 400);
2223 *
if (group_name == std::string(
"displacement"))
2224 *
group.set_attribute<std::string>(
2225 *
"simulation_name", std::string(
"phononic_cavity_displacement"));
2227 *
group.set_attribute<std::string>(
2228 *
"simulation_name", std::string(
"phononic_cavity_calibration"));
2230 *
group.set_attribute<
bool>(
"save_vtu_files",
false);
2236 * Displacement simulation. The parameters are read from the
2237 * displacement
HDF5 group and the results are saved in the same
HDF5
2241 *
auto displacement = data.open_group(
"displacement");
2242 * step62::Parameters<dim> parameters(displacement);
2244 * step62::ElasticWave<dim> elastic_problem(parameters);
2245 * elastic_problem.run();
2251 * Calibration simulation. The parameters are read from the calibration
2255 *
auto calibration = data.open_group(
"calibration");
2256 * step62::Parameters<dim> parameters(calibration);
2258 * step62::ElasticWave<dim> elastic_problem(parameters);
2259 * elastic_problem.run();
2262 *
catch (std::exception &exc)
2264 * std::cerr << std::endl
2266 * <<
"----------------------------------------------------"
2268 * std::cerr <<
"Exception on processing: " << std::endl
2269 * << exc.what() << std::endl
2270 * <<
"Aborting!" << std::endl
2271 * <<
"----------------------------------------------------"
2278 * std::cerr << std::endl
2280 * <<
"----------------------------------------------------"
2282 * std::cerr <<
"Unknown exception!" << std::endl
2283 * <<
"Aborting!" << std::endl
2284 * <<
"----------------------------------------------------"
2292<a name=
"Results"></a><h1>Results</h1>
2295<a name=
"Resonancefrequencyandbandgap"></a><h3>Resonance frequency and bandgap</h3>
2298The results are analyzed in the
2299[jupyter notebook](https:
2300with the following code
2302h5_file = h5py.File(
'results.h5',
'r')
2303data = h5_file[
'data']
2305# Gaussian function that we use to fit the resonance
2306def resonance_f(freq, freq_m, quality_factor, max_amplitude):
2307 omega = 2 * constants.pi * freq
2308 omega_m = 2 * constants.pi * freq_m
2309 gamma = omega_m / quality_factor
2310 return max_amplitude * omega_m**2 *
gamma**2 / (((omega_m**2 - omega**2)**2 +
gamma**2 * omega**2))
2312frequency = data[
'displacement'][
'frequency'][...]
2313# Average the probe points
2314displacement = np.
mean(data[
'displacement'][
'displacement'], axis=0)
2315calibration_displacement = np.
mean(data[
'calibration'][
'displacement'], axis=0)
2316reflection_coefficient = displacement / calibration_displacement
2317reflectivity = (np.
abs(np.
mean(data[
'displacement'][
'displacement'][...]**2, axis=0))/
2318 np.
abs(np.
mean(data[
'calibration'][
'displacement'][...]**2, axis=0)))
2322 y_data = reflectivity
2323 quality_factor_guess = 1e3
2324 freq_guess = x_data[np.argmax(y_data)]
2325 amplitude_guess = np.
max(y_data)
2326 fit_result, covariance = scipy.optimize.curve_fit(resonance_f, x_data, y_data,
2327 [freq_guess, quality_factor_guess, amplitude_guess])
2328 freq_m = fit_result[0]
2329 quality_factor = np.
abs(fit_result[1])
2330 max_amplitude = fit_result[2]
2331 y_data_fit = resonance_f(x_data, freq_m, quality_factor, max_amplitude)
2334 plt.plot(frequency / 1e9, reflectivity, frequency / 1e9, y_data_fit)
2335 plt.xlabel(
'frequency (GHz)')
2336 plt.ylabel(
'amplitude^2 (a.u.)')
2337 plt.title(
'Transmission\n' +
'freq = ' +
"%.7g" % (freq_guess / 1e9) +
'GHz Q = ' +
"%.6g" % quality_factor)
2340 plt.plot(frequency / 1e9, reflectivity)
2341 plt.xlabel(
'frequency (GHz)')
2342 plt.ylabel(
'amplitude^2 (a.u.)')
2343 plt.title(
'Transmission')
2346plt.plot(frequency / 1e9, np.
angle(reflection_coefficient))
2347plt.xlabel(
'frequency (GHz)')
2348plt.ylabel(
'phase (rad)')
2349plt.title(
'Phase (transmission coefficient)\n')
2355A phononic cavity is characterized by the
2356[resonance frequency](https:
2357[the quality factor](https:
2358The quality factor is
equal to the ratio between the stored energy in the resonator and the energy
2359dissipated energy per cycle, which is approximately equivalent to the ratio between the
2360resonance frequency and the
2361[full width at half maximum (FWHM)](https:
2362The FWHM is
equal to the bandwidth over which the power of vibration is greater than half the
2363power at the resonant frequency.
2365Q = \frac{f_r}{\Delta f} = \frac{\omega_r}{\Delta \omega} =
23662 \pi \times \frac{\text{energy stored}}{\text{energy dissipated per cycle}}
2369The square of the amplitude of the mechanical resonance @f$a^2@f$ as a function of the frequency
2372a^2 = a_\textrm{
max}^2\frac{\omega^2\Gamma^2}{(\omega_r^2-\omega^2)^2+\Gamma^2\omega^2}
2374where @f$f_r = \frac{\omega_r}{2\pi}@f$ is the resonance frequency and @f$\Gamma=\frac{\omega_r}{Q}@f$ is the dissipation rate.
2375We used the previous equation in the jupyter notebook to fit the mechanical resonance.
2377Given the
values we have chosen
for the parameters, one could estimate the resonance frequency
2378analytically. Indeed,
this is then confirmed by what we get in
this program:
2379the phononic superlattice cavity exhibits a mechanical resonance at 20GHz and a quality factor of 5046.
2380The following images show the transmission amplitude and phase as a function of frequency in the
2381vicinity of the resonance frequency:
2383<img alt=
"Phononic superlattice cavity" src=
"https://www.dealii.org/images/steps/developer/step-62.05.png" height=
"400" />
2384<img alt=
"Phononic superlattice cavity" src=
"https://www.dealii.org/images/steps/developer/step-62.06.png" height=
"400" />
2386The images above suggest that the periodic structure has its intended effect: It really only lets waves of a very
2387specific frequency pass through, whereas all other waves are reflected. This is of course precisely what one builds
2388these sorts of devices
for.
2389But it is not quite
this easy. In practice, there is really only a
"band gap", i.e., the device blocks waves other than
2390the desired one at 20GHz only within a certain frequency range. Indeed, to find out how large
this "gap" is within
2391which waves are blocked, we can extend the frequency range to 16 GHz through the appropriate parameters in the
2392input file. We then obtain the following image:
2394<img alt=
"Phononic superlattice cavity" src=
"https://www.dealii.org/images/steps/developer/step-62.07.png" height=
"400" />
2396What
this image suggests is that in the range of around 18 to around 22 GHz, really only the waves with a frequency
2397of 20 GHz are allowed to pass through, but beyond
this range, there are plenty of other frequencies that can pass
2400<a name=
"Modeprofile"></a><h3>Mode profile</h3>
2403We can inspect the mode profile with Paraview or VisIt.
2404As we have discussed, at resonance all the mechanical
2405energy is transmitted and the amplitude of motion is amplified
inside the cavity.
2406It can be observed that the PMLs are quite effective to truncate the solution.
2407The following image shows the mode profile at resonance:
2409<img alt=
"Phononic superlattice cavity" src=
"https://www.dealii.org/images/steps/developer/step-62.08.png" height=
"400" />
2411On the other hand, out of resonance all the mechanical energy is
2412reflected. The following image shows the profile at 19.75 GHz.
2413Note the interference between the force pulse and the reflected wave
2414at the position @f$x=-8\mu\textrm{m}@f$.
2416<img alt=
"Phononic superlattice cavity" src=
"https://www.dealii.org/images/steps/developer/step-62.09.png" height=
"400" />
2418<a name=
"Experimentalapplications"></a><h3>Experimental applications</h3>
2421Phononic superlattice cavities find application in
2422[quantum optomechanics](https:
2423Here we have presented the simulation of a 2D superlattice cavity,
2424but
this code can be used as well to simulate
"real world" 3D devices such as
2425[micropillar superlattice cavities](https:
2426which are promising candidates to study macroscopic quantum phenomena.
2427The 20GHz mode of a micropillar superlattice cavity is essentially a mechanical harmonic oscillator that is very well isolated
2428from the environment. If the device is cooled down to 20mK in a dilution fridge, the mode would then become a
2429macroscopic quantum harmonic oscillator.
2432<a name=
"Possibilitiesforextensions"></a><h3>Possibilities
for extensions</h3>
2435Instead of setting the parameters in the
C++ file we could
set the parameters
2436using a python script and save them in the
HDF5 file that we will use
for
2437the simulations. Then the deal.II program will read the parameters from the
2443import matplotlib.pyplot as plt
2445import scipy.constants as constants
2446import scipy.optimize
2448# This considerably reduces the size of the svg data
2449plt.rcParams[
'svg.fonttype'] =
'none'
2451h5_file = h5py.
File(
'results.h5',
'w')
2457for group in [displacement, calibration]:
2458 # Dimensions of the domain
2459 # The waveguide length is equal to dimension_x
2460 group.attrs[
'dimension_x'] = 2
e-5
2461 # The waveguide width is equal to dimension_y
2462 group.attrs[
'dimension_y'] = 2
e-8
2464 # Position of the probe that we use to measure the flux
2465 group.attrs[
'probe_pos_x'] = 8
e-6
2466 group.attrs[
'probe_pos_y'] = 0
2467 group.attrs[
'probe_width_y'] = 2
e-08
2469 # Number of points in the probe
2470 group.attrs[
'nb_probe_points'] = 5
2473 group.attrs[
'grid_level'] = 1
2476 group.attrs[
'cavity_resonance_frequency'] = 20e9
2477 group.attrs[
'nb_mirror_pairs'] = 15
2480 group.attrs[
'poissons_ratio'] = 0.27
2481 group.attrs[
'youngs_modulus'] = 270000000000.0
2482 group.attrs[
'material_a_rho'] = 3200
2483 if group == displacement:
2484 group.attrs[
'material_b_rho'] = 2000
2486 group.attrs[
'material_b_rho'] = 3200
2487 group.attrs[
'lambda'] = (
group.attrs[
'youngs_modulus'] *
group.attrs[
'poissons_ratio'] /
2488 ((1 +
group.attrs[
'poissons_ratio']) *
2489 (1 - 2 *
group.attrs[
'poissons_ratio'])))
2490 group.attrs[
'mu']= (
group.attrs[
'youngs_modulus'] / (2 * (1 +
group.attrs[
'poissons_ratio'])))
2493 group.attrs[
'max_force_amplitude'] = 1e26
2494 group.attrs[
'force_sigma_x'] = 1e-7
2495 group.attrs[
'force_sigma_y'] = 1
2496 group.attrs[
'max_force_width_x'] = 3e-7
2497 group.attrs[
'max_force_width_y'] = 2e-8
2498 group.attrs[
'force_x_pos'] = -8e-6
2499 group.attrs[
'force_y_pos'] = 0
2502 group.attrs[
'pml_x'] = True
2503 group.attrs[
'pml_y'] = False
2504 group.attrs[
'pml_width_x'] = 1.8e-6
2505 group.attrs[
'pml_width_y'] = 5e-7
2506 group.attrs[
'pml_coeff'] = 1.6
2507 group.attrs[
'pml_coeff_degree'] = 2
2510 group.attrs[
'center_frequency'] = 20e9
2511 group.attrs[
'frequency_range'] = 0.5e9
2512 group.attrs[
'start_frequency'] =
group.attrs[
'center_frequency'] -
group.attrs[
'frequency_range'] / 2
2513 group.attrs[
'stop_frequency'] =
group.attrs[
'center_frequency'] +
group.attrs[
'frequency_range'] / 2
2514 group.attrs[
'nb_frequency_points'] = 400
2517 if group == displacement:
2518 group.attrs[
'simulation_name'] =
'phononic_cavity_displacement'
2520 group.attrs[
'simulation_name'] =
'phononic_cavity_calibration'
2521 group.attrs[
'save_vtu_files'] = False
2526In order to read the
HDF5 parameters we have to use the
2527HDF5::File::FileAccessMode::open flag.
2529 HDF5::File data_file(
"results.h5",
2530 HDF5::File::FileAccessMode::open,
2532 auto data = data_file.open_group(
"data");
2536<a name=
"PlainProg"></a>
2537<h1> The plain program</h1>
2538@include
"step-62.cc"
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
void write(const Container &data)
void write_selection(const Container &data, const std::vector< hsize_t > &coordinates)
Group create_group(const std::string &name) const
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
void loop(ITERATOR begin, std_cxx20::type_identity_t< ITERATOR > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
@ component_is_part_of_vector
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
int create_group(const MPI_Comm comm, const MPI_Group &group, const int tag, MPI_Comm *new_comm)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
long double gamma(const unsigned int n)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
static constexpr double PI
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation