Skip to content

Commit

Permalink
Use FEValues extractors, avoid explicit use of view.
Browse files Browse the repository at this point in the history
  • Loading branch information
@bangerth
bangerth committed Oct 24, 2023
1 parent 9050964 commit 59e7f8a
Showing 1 changed file with 36 additions and 29 deletions.
65 changes: 36 additions & 29 deletions examples/step-81/step-81.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -622,27 +622,28 @@ namespace Step81
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

// This is assembling the interior of the domain on the left hand side.
// So we are assembling
// Next, let us assemble on the interior of the domain on the left hand
// side. So we are computing
// \f{align*}{
// \int_\Omega (\mu_r^{-1}\nabla \times \varphi_i) \cdot
// (\nabla\times\bar{\varphi}_j)\text{d}x
// - \int_\Omega \varepsilon_r\varphi_i \cdot \bar{\varphi}_j\text{d}x
// \int_\Omega (\mu_r^{-1}\nabla \times \varphi_i) \cdot
// (\nabla\times\bar{\varphi}_j)\text{d}x
// -
// \int_\Omega \varepsilon_r\varphi_i \cdot \bar{\varphi}_j\text{d}x
// \f}
// and
// \f{align}{
// i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x
// - \int_\Omega \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i}) \text{d}x.
// i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x
// - \int_\Omega \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i}) \text{d}x.
// \f}
// In doing so, we need test functions $\varphi_i$ and $\varphi_j$, and the
// curl of these test variables. We must be careful with the signs of the
// imaginary parts of these complex test variables. Moreover, we have a
// conditional that changes the parameters if the cell is in the PML region.
const FEValuesExtractors::Vector real_part(0);
const FEValuesExtractors::Vector imag_part(dim);
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
FEValuesViews::Vector<dim> real_part(fe_values, 0);
FEValuesViews::Vector<dim> imag_part(fe_values, dim);

cell_matrix = 0.;
cell_rhs = 0.;
Expand Down Expand Up @@ -671,21 +672,25 @@ namespace Step81
{
constexpr std::complex<double> imag{0., 1.};

const auto phi_i = real_part.value(i, q_point) -
imag * imag_part.value(i, q_point);
const auto curl_phi_i = real_part.curl(i, q_point) -
imag * imag_part.curl(i, q_point);
const auto phi_i =
fe_values[real_part].value(i, q_point) -
imag * fe_values[imag_part].value(i, q_point);
const auto curl_phi_i =
fe_values[real_part].curl(i, q_point) -
imag * fe_values[imag_part].curl(i, q_point);

const auto rhs_value =
(imag * scalar_product(J_a, phi_i)) * fe_values.JxW(q_point);
cell_rhs(i) += rhs_value.real();

for (const auto j : fe_values.dof_indices())
{
const auto phi_j = real_part.value(j, q_point) +
imag * imag_part.value(j, q_point);
const auto curl_phi_j = real_part.curl(j, q_point) +
imag * imag_part.curl(j, q_point);
const auto phi_j =
fe_values[real_part].value(j, q_point) +
imag * fe_values[imag_part].value(j, q_point);
const auto curl_phi_j =
fe_values[real_part].curl(j, q_point) +
imag * fe_values[imag_part].curl(j, q_point);

const auto temp =
(scalar_product(mu_inv * curl_phi_j, curl_phi_i) -
Expand All @@ -706,6 +711,8 @@ namespace Step81
// \f}
// respectively. The test variables and the PML are implemented
// similarly as the domain.
const FEValuesExtractors::Vector real_part(0);
const FEValuesExtractors::Vector imag_part(dim);
for (const auto &face : cell->face_iterators())
{
if (face->at_boundary())
Expand All @@ -714,8 +721,6 @@ namespace Step81
if (id != 0)
{
fe_face_values.reinit(cell, face);
FEValuesViews::Vector<dim> real_part(fe_face_values, 0);
FEValuesViews::Vector<dim> imag_part(fe_face_values, dim);

for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
Expand All @@ -742,15 +747,17 @@ namespace Step81
constexpr std::complex<double> imag{0., 1.};

const auto phi_i =
real_part.value(i, q_point) -
imag * imag_part.value(i, q_point);
fe_face_values[real_part].value(i, q_point) -
imag *
fe_face_values[imag_part].value(i, q_point);
const auto phi_i_T = tangential_part(phi_i, normal);

for (const auto j : fe_face_values.dof_indices())
{
const auto phi_j =
real_part.value(j, q_point) +
imag * imag_part.value(j, q_point);
fe_face_values[real_part].value(j, q_point) +
imag *
fe_face_values[imag_part].value(j, q_point);
const auto phi_j_T =
tangential_part(phi_j, normal) *
fe_face_values.JxW(q_point);
Expand Down Expand Up @@ -779,8 +786,6 @@ namespace Step81
continue; /* skip this face */

fe_face_values.reinit(cell, face);
FEValuesViews::Vector<dim> real_part(fe_face_values, 0);
FEValuesViews::Vector<dim> imag_part(fe_face_values, dim);

for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
Expand All @@ -799,15 +804,17 @@ namespace Step81
{
constexpr std::complex<double> imag{0., 1.};

const auto phi_i = real_part.value(i, q_point) -
imag * imag_part.value(i, q_point);
const auto phi_i =
fe_face_values[real_part].value(i, q_point) -
imag * fe_face_values[imag_part].value(i, q_point);
const auto phi_i_T = tangential_part(phi_i, normal);

for (const auto j : fe_face_values.dof_indices())
{
const auto phi_j =
real_part.value(j, q_point) +
imag * imag_part.value(j, q_point);
fe_face_values[real_part].value(j, q_point) +
imag *
fe_face_values[imag_part].value(j, q_point);
const auto phi_j_T = tangential_part(phi_j, normal);

const auto temp =
Expand Down

0 comments on commit 59e7f8a

Please sign in to comment.