From e60a449c310ddbd7cd6707dda6c048c9f24cc0d8 Mon Sep 17 00:00:00 2001 From: Manaswinee Bezbaruah Date: Fri, 18 Mar 2022 11:34:39 -0500 Subject: [PATCH] add cartoon --- examples/step-81/doc/intro.dox | 24 ++- examples/step-81/doc/results.dox | 107 +++++++--- examples/step-81/step-81.cc | 329 +++++++++++++++++-------------- 3 files changed, 278 insertions(+), 182 deletions(-) diff --git a/examples/step-81/doc/intro.dox b/examples/step-81/doc/intro.dox index b9e5d19196..15bc1ba658 100644 --- a/examples/step-81/doc/intro.dox +++ b/examples/step-81/doc/intro.dox @@ -24,13 +24,13 @@ how to implement simple first-order absorbing boundary conditions and a more sophisticated "perfectly matched layer" for electromagnetic waves. -

Time-Harmonic Maxwell's Equations with interface conditions

+

Time-Harmonic Maxwell's Equations with interface conditions

We start the discussion with a short derivation of the governing equations and some pointers to literature. -

Derivation of time-harmonic Maxwell's equations

+

Derivation of time-harmonic Maxwell's equations

In two ($d=2$) or three ($d=3$) spatial dimensions, the time evolution of an electromagnetic @@ -105,7 +105,7 @@ write $\mathbf{E}(\mathbf{x})$, $\mathbf{H}(\mathbf{x})$, etc., when referring to the time-harmonic fields. -

Jump conditions on lower dimensional interfaces

+

Jump conditions on lower dimensional interfaces

Graphene is a two-dimensional carbon allotrope with a single atom layer that is arranged in a honeycomb lattice @cite Geim2004. Due to its @@ -133,7 +133,11 @@ field. This is best seen by visualizing Ampère's law, -@f[\text{TODO}@f] +@htmlonly +

+ Visualization of Ohm's law and Ampère's law leading to a jump condition over the interface +

+@endhtmlonly and then taking the limit of the upper and lower part of the line integral approaching the sheet. In contrast, the tangential part of the electric @@ -151,7 +155,7 @@ when approaching the interface from above or below the interface: $\mathbf{F}^\pm(\mathbf{x})=\lim_{\delta\to0,\delta>0}\mathbf{F}(\mathbf{x}\pm\delta\mathbf{\nu})$. -

Rescaling

+

Rescaling

We will be using a rescaled version of the Maxwell's equations described above. The rescaling has the following key differences:
@@ -210,7 +214,7 @@ i\varepsilon_r\hat{\mathbf{E}} + \nabla\times(\mu^{-1}\mathbf{H}) We will omit the hat in further discussion for ease of notation. -

Variational Statement

+

Variational Statement

Let $\Omega \subset \mathbb{R}^n$, $(n = 2,3)$ be a simply connected and bounded domain with Lipschitz-continuous and piecewise smooth boundary, $\partial\Omega$. @@ -310,7 +314,7 @@ A(\mathbf{E},\varphi) = F(\varphi) @f] -

Absorbing boundary conditions and perfectly matched layer

+

Absorbing boundary conditions and perfectly matched layer

Moreover, the above equations are supplemented by the Silver-Müller radiation condition, if the ambient (unbounded) medium is isotropic. This amounts to the @@ -326,7 +330,7 @@ x \not\in \Sigma In our case, we eliminate reflection from infinity by implementing a PML and avoid the explicit use of the last condition. -

Discretization Scheme

+

Discretization Scheme

The variational form is discretized on a non-uniform quadrilateral mesh with higher-order, curl-conforming Nédélec elements. This way the interface with a @@ -361,7 +365,7 @@ A_{ij} = F_i Using a skeleton similar to step-4, we have constructed a Maxwell class and we have used complex-valued FENedelec elements to solve our equations.
-

Perfectly Matched Layer

+

Perfectly Matched Layer

The SPP amplitude is negatively effected by the absorbing boundary condition and this causes the solution image to be distorted. In order to reduce the resonance and distortion in our solutions, we are implementing a Perfectly Matched Layer @@ -381,7 +385,7 @@ be a function of the radial distance $\rho$ from the origin $e_r$. The normal fi $\nu$ of $\Sigma$ is orthogonal to the radial direction $e_r$, which makes $\mathbf{J}_a \equiv 0$ and $\mathbf{M}_a \equiv 0$ within the PML.
- + Introduce a change of coordinates @f[ diff --git a/examples/step-81/doc/results.dox b/examples/step-81/doc/results.dox index dee3123006..ab14f3f31d 100644 --- a/examples/step-81/doc/results.dox +++ b/examples/step-81/doc/results.dox @@ -14,7 +14,7 @@ The following images are the outputs for the imaginary $E_x$ without the interfa @code # use absorbing boundary conditions? - set absorbing boundary condition boolean = false + set absorbing boundary condition = false # position of the dipole set dipole position = 0, 0 @@ -25,21 +25,24 @@ The following images are the outputs for the imaginary $E_x$ without the interfa # surface conductivity between material 1 and material 2 set sigma = 0, 0; 0, 0| 0, 0; 0, 0 @endcode - + Following are the output images: + + - -
+ Visualization of the solution of step-81 with no interface, Dirichlet boundary conditions and PML strength 0 +

Solution with no interface, Dirichlet boundary conditions and PML strength 0.

+ 		- Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 0 - - Visualization of the solution of step-81 with no interface, no absorbing boundary conditions and PML strength 0 + Visualization of the solution of step-81 with no interface, no absorbing boundary conditions and PML strength 0 +

Solution with no interface, absorbing boundary conditions and PML strength 0.

 - Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 4 + Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 4 +

Solution with no interface, absorbing boundary conditions and PML strength 4.
@@ -62,16 +65,19 @@ Once again, we will visualize the output with absorbing boundary conditions and - - + Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0 +

Solution with an interface, absorbing boundary conditions and PML strength 0.

+ + +
- Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 0 + Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0 +

Solution with an interface, Dirichlet boundary conditions and PML strength 0.

 - Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 0 - - Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 4 - + Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4 +

Solution with an interface, absorbing boundary conditions and PML strength 4.

+
@@ -79,26 +85,69 @@ Once again, we will visualize the output with absorbing boundary conditions and - - + +
- Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 0 + Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0 +

Solution with an interface, Dirichlet boundary conditions and PML strength 0.

 - Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 0 - - Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 4 + Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0 +

Solution with an interface, absorbing boundary conditions and PML strength 0.

+ 		+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4 +

Solution with an interface, absorbing boundary conditions and PML strength 4.
The SPP is confined near the interface that we created, however without absorbing boundary conditions, we don't observe a dissipation effect. On adding the absorbing boundary conditions, we observe distortion and resonance and we still don't notice any dissipation. As expected, the PML removes the distortion and resonance. The standing wave is also dissipating and getting absorbed within the PML, and as we increase the PML strength, the standing wave will dissipate more within the PML ring. +Here are some animations to demonstrate the effect of the PML + + + + + + + + +
+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0 +

Solution with an interface, Dirichlet boundary conditions and PML strength 0.

+ 		+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0 +

Solution with an interface, absorbing boundary conditions and PML strength 0.

+ 		+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4 +

Solution with an interface, absorbing boundary conditions and PML strength 4.

+
+ + + + + + + + + + +
+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0 +

Solution with an interface, Dirichlet boundary conditions and PML strength 0.

+ 		+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0 +

Solution with an interface, absorbing boundary conditions and PML strength 0.

+ 		+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4 +

Solution with an interface, absorbing boundary conditions and PML strength 4.

+
+

Notes

Real and Complex Matrices

-As is evident from the results, we are splitting our solution matrices into the real and the imaginary components. We started off using the $H^{curl}$ conforming Nédélec Elements, and we made two copies of the Finite Elements in order -to represent the real and the imaginary components of our input (FE_NedelecSZ was used instead of FE_Nedelec to avoid the sign conflicts issues present in traditional Nédélec elements). In the assembly, we create two vectors of dimension $dim$ that assist us in extracting the real and the imaginary components of our finite elements. +As is evident from the results, we are splitting our solution matrices into the real and the imaginary components. We started off using the $H^{curl}$ conforming Nédélec Elements, and we made two copies of the Finite Elements in order +to represent the real and the imaginary components of our input (FE_NedelecSZ was used instead of FE_Nedelec to avoid the sign conflicts issues present in traditional Nédélec elements). In the assembly, we create two vectors of dimension $dim$ that assist us in extracting the real and the imaginary components of our finite elements.

Rotations and Scaling

@@ -134,7 +183,7 @@ We will create a video demonstrating the wave in motion, which is essentially an template void Maxwell::output_results(unsigned int t) { - std::cout << "Running step:" << alpha << std::endl; + std::cout << "Running step:" << t << std::endl; DataOut<2> data_out; data_out.attach_dof_handler(dof_handler); Vector postprocessed; @@ -151,10 +200,10 @@ void Maxwell::output_results(unsigned int t) data_out.add_data_vector(postprocessed, {"E_x","E_y","null0","null1"}); data_out.build_patches(); const std::string filename = - "solution-" + Utilities::int_to_string(alpha) + ".vtk"; + "solution-" + Utilities::int_to_string(t) + ".vtk"; std::ofstream output(filename); data_out.write_vtk(output); - std::cout << "Done running step:" << alpha << std::endl; + std::cout << "Done running step:" << t << std::endl; @endcode @@ -168,8 +217,8 @@ for(int t = 0; t<=100; t++){ This would generate 100 solution .vtk files, which can be opened in a group on Paraview and then can be saved as an animation. We used FFMPEG to generate gifs.

Resulting videos

-Following are the resulting videos of our experiments. As we see, TODO +Following are the resulting videos of our experiments. As we see, TODO

Possibilities for Extension

-The current program doesn't allow for iterative solvers as the solutions will not converge with an iterative solver. One possible direction for future work is to implement an iterative solver and involve more preconditioners. An advantage of iterative solvers is the more efficient memory usage, and our current memory usage does not allow for a large number of DOFs. -Another possible direction would be to perform Local Mesh Refinement (instead of Global Mesh Refinement). This will also help us visualize more DOFs in a more memory and time efficient way. +The current program doesn't allow for iterative solvers as the solutions will not converge with an iterative solver. One possible direction for future work is to implement an iterative solver and involve more preconditioners. An advantage of iterative solvers is the more efficient memory usage, and our current memory usage does not allow for a large number of DOFs. +Another possible direction would be to perform Local Mesh Refinement (instead of Global Mesh Refinement). This will also help us visualize more DOFs in a more memory and time efficient way. diff --git a/examples/step-81/step-81.cc b/examples/step-81/step-81.cc index e318b4d5f1..0fbf4fd6d9 100644 --- a/examples/step-81/step-81.cc +++ b/examples/step-81/step-81.cc @@ -19,10 +19,10 @@ // @sect3{Include files} - // The set of include files is quite standard. The most notable incluse is - // the fe/fe_nedelec_sz.h file that allows us to use the FE_NedelecSZ elements. - // This is an implementation of the $H^{curl}$ conforming Nédélec Elements - // that resolves the sign conflict issues that arise from parametrization. +// The set of include files is quite standard. The most notable incluse is +// the fe/fe_nedelec_sz.h file that allows us to use the FE_NedelecSZ elements. +// This is an implementation of the $H^{curl}$ conforming Nédélec Elements +// that resolves the sign conflict issues that arise from parametrization. #include #include @@ -75,29 +75,31 @@ namespace Step81 using namespace dealii; using namespace std::complex_literals; - // @sect4{Parameters Class} - - // The Parameters class inherits ParameterAcceptor, and instantiates all the - // coefficients in our variational equations. - // These coefficients are passed through ParameterAcceptor and are editable - // through a .prm file - // More explanation on the use and inheritance from the ParameterAcceptor - // can be found in step-60. - - // epsilon is the Electric Permitivitty coefficient and it is a rank 2 tensor. Depending on the material, - // we assign the i^th diagonal element of the tensor to the material epsilon value - // (one of the private epsilon_1_ or epsilon_2_ variables). - // - // mu_inv is the inverese of the Magnetic Permiabillity coefficient and it is a complex number. - - // sigma is the Surface Conductivity coefficient between material left and material right - // and it is a rank 2 tensor. It is only changed if we are at the interface between two - // materials. If we are at an interface, we assign the i^th diagonal element of the - // tensor to the private sigma_ value. - - // J_a is the strength and orientation of the dipole. It is a rank 1 tensor that depends - // on the private dipole_position_, dipole_radius_, dipole_strength_, dipole_orientation_ - // variables. + // @sect4{Parameters Class} + + // The Parameters class inherits ParameterAcceptor, and instantiates all the + // coefficients in our variational equations. + // These coefficients are passed through ParameterAcceptor and are editable + // through a .prm file + // More explanation on the use and inheritance from the ParameterAcceptor + // can be found in step-60. + + // epsilon is the Electric Permitivitty coefficient and it is a rank 2 tensor. + // Depending on the material, we assign the i^th diagonal element of the + // tensor to the material epsilon value (one of the private epsilon_1_ or + // epsilon_2_ variables). + // + // mu_inv is the inverese of the Magnetic Permiabillity coefficient and it is + // a complex number. + + // sigma is the Surface Conductivity coefficient between material left and + // material right and it is a rank 2 tensor. It is only changed if we are at + // the interface between two materials. If we are at an interface, we assign + // the i^th diagonal element of the tensor to the private sigma_ value. + + // J_a is the strength and orientation of the dipole. It is a rank 1 tensor + // that depends on the private dipole_position_, dipole_radius_, + // dipole_strength_, dipole_orientation_ variables. template class Parameters : public ParameterAcceptor @@ -114,18 +116,16 @@ namespace Step81 using curl_type = Tensor<1, dim == 2 ? 1 : dim, rank0_type>; public: - rank2_type epsilon(const Point &x, - types::material_id material); + rank2_type epsilon(const Point &x, types::material_id material); - std::complex mu_inv(const Point &x, - types::material_id material); + std::complex mu_inv(const Point & x, + types::material_id material); - rank2_type sigma(const dealii::Point &x, - types::material_id left, - types::material_id right); + rank2_type sigma(const dealii::Point &x, + types::material_id left, + types::material_id right); - rank1_type J_a(const dealii::Point &point, - types::material_id id); + rank1_type J_a(const dealii::Point &point, types::material_id id); private: rank2_type epsilon_1; @@ -134,10 +134,10 @@ namespace Step81 std::complex mu_inv_2; rank2_type sigma_tensor; - double dipole_radius; - Point dipole_position; - Tensor<1, dim, double> dipole_orientation; - rank0_type dipole_strength; + double dipole_radius; + Point dipole_position; + Tensor<1, dim> dipole_orientation; + rank0_type dipole_strength; }; @@ -177,11 +177,9 @@ namespace Step81 add_parameter("dipole radius", dipole_radius, "radius of the dipole"); dipole_position = Point(0., 0.8); - add_parameter("dipole position", - dipole_position, - "position of the dipole"); + add_parameter("dipole position", dipole_position, "position of the dipole"); - dipole_orientation = Tensor<1, dim, double>{{0., 1.}}; + dipole_orientation = Tensor<1, dim>{{0., 1.}}; add_parameter("dipole orientation", dipole_orientation, "orientation of the dipole"); @@ -192,14 +190,14 @@ namespace Step81 template typename Parameters::rank2_type - Parameters::epsilon(const Point & /*x*/, + Parameters::epsilon(const Point & /*x*/, types::material_id material) { return (material == 1 ? epsilon_1 : epsilon_2); } template - std::complex Parameters::mu_inv(const Point & /*x*/, + std::complex Parameters::mu_inv(const Point & /*x*/, types::material_id material) { return (material == 1 ? mu_inv_1 : mu_inv_2); @@ -207,7 +205,7 @@ namespace Step81 template typename Parameters::rank2_type - Parameters::sigma(const dealii::Point & /*x*/, + Parameters::sigma(const dealii::Point & /*x*/, types::material_id left, types::material_id right) { @@ -216,7 +214,7 @@ namespace Step81 template typename Parameters::rank1_type - Parameters::J_a(const dealii::Point &point, + Parameters::J_a(const dealii::Point &point, types::material_id /*id*/) { rank1_type J_a; @@ -230,13 +228,23 @@ namespace Step81 return J_a; } - // @sect4{PerfectlyMatchedLayer Class} - // The PerfectlyMatchedLayer class inherits ParameterAcceptor, - // and it modifies our coefficients from Parameters. - // The radii and the strength of the PML is specified, and the - // coefficients will be modified using transformation - // matrices within the PML region. The radii and strength of - // the PML are editable through a .prm file + // @sect4{PerfectlyMatchedLayer Class} + // The PerfectlyMatchedLayer class inherits ParameterAcceptor, + // and it modifies our coefficients from Parameters. + // The radii and the strength of the PML is specified, and the + // coefficients will be modified using transformation + // matrices within the PML region. The radii and strength of + // the PML are editable through a .prm file + // The rotation function is the $T_{exer}$ mentioned in the + // perfectly matched layer section of the introduction. + // Moreover, the matrices A, B and C are defined as mentioned + // @f[ + // A = T_{e_xe_r}^{-1} + // \text{diag}\left(\frac{1}{\bar{d}^2},\frac{1}{d\bar{d}}\right)T_{e_xe_r},\qquad + // B = T_{e_xe_r}^{-1} \text{diag}\left(d,\bar{d}\right)T_{e_xe_r},\qquad + // C = T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}},\frac{1}{d}\right) + // T_{e_xe_r}.\qquad + // @f] template class PerfectlyMatchedLayer : public ParameterAcceptor @@ -256,20 +264,20 @@ namespace Step81 double outer_radius; double strength; - std::complex d_tensor(const Point point); + std::complex d(const Point point); - std::complex d_bar_tensor(const Point point); + std::complex d_bar(const Point point); - rank2_type T_exer(std::complex d_1, - std::complex d_2, - Point point); + rank2_type rotation(std::complex d_1, + std::complex d_2, + Point point); - rank2_type a_matrix(const Point point); + rank2_type a_matrix(const Point point); - rank2_type b_matrix(const Point point); + rank2_type b_matrix(const Point point); - rank2_type c_matrix(const Point point); + rank2_type c_matrix(const Point point); }; @@ -281,30 +289,30 @@ namespace Step81 add_parameter("inner radius", inner_radius, "inner radius of the PML shell"); - outer_radius = 15.; + outer_radius = 20.; add_parameter("outer radius", outer_radius, "outer radius of the PML shell"); strength = 8.; add_parameter("strength", strength, "strength of the PML"); - }; + } template typename std::complex - PerfectlyMatchedLayer::d_tensor(const Point point) + PerfectlyMatchedLayer::d(const Point point) { const auto radius = point.norm(); const double s = strength * ((radius - inner_radius) * (radius - inner_radius)) / ((outer_radius - inner_radius) * (outer_radius - inner_radius)); - return 1 + 1.0i * s; + return 1.0 + 1.0i * s; } template typename std::complex - PerfectlyMatchedLayer::d_bar_tensor(const Point point) + PerfectlyMatchedLayer::d_bar(const Point point) { const auto radius = point.norm(); const double s_bar = @@ -312,15 +320,15 @@ namespace Step81 ((radius - inner_radius) * (radius - inner_radius) * (radius - inner_radius)) / (radius * (outer_radius - inner_radius) * (outer_radius - inner_radius)); - return 1 + 1.0i * s_bar; + return 1.0 + 1.0i * s_bar; } template typename PerfectlyMatchedLayer::rank2_type - PerfectlyMatchedLayer::T_exer(std::complex d_1, - std::complex d_2, - Point point) + PerfectlyMatchedLayer::rotation(std::complex d_1, + std::complex d_2, + Point point) { rank2_type result; result[0][0] = point[0] * point[0] * d_1 + point[1] * point[1] * d_2; @@ -333,39 +341,43 @@ namespace Step81 template typename PerfectlyMatchedLayer::rank2_type - PerfectlyMatchedLayer::a_matrix(const Point point) + PerfectlyMatchedLayer::a_matrix(const Point point) { - const auto d = d_tensor(point); - const auto d_bar = d_bar_tensor(point); - return invert(T_exer(d * d, d * d_bar, point)) * - T_exer(d * d, d * d_bar, point); + const auto d = this->d(point); + const auto d_bar = this->d_bar(point); + return invert(rotation(d * d, d * d_bar, point)) * + rotation(d * d, d * d_bar, point); } template typename PerfectlyMatchedLayer::rank2_type - PerfectlyMatchedLayer::b_matrix(const Point point) + PerfectlyMatchedLayer::b_matrix(const Point point) { - const auto d = d_tensor(point); - const auto d_bar = d_bar_tensor(point); - return invert(T_exer(d, d_bar, point)) * T_exer(d, d_bar, point); + const auto d = this->d(point); + const auto d_bar = this->d_bar(point); + return invert(rotation(d, d_bar, point)) * rotation(d, d_bar, point); } template typename PerfectlyMatchedLayer::rank2_type - PerfectlyMatchedLayer::c_matrix(const Point point) + PerfectlyMatchedLayer::c_matrix(const Point point) { - const auto d = d_tensor(point); - const auto d_bar = d_bar_tensor(point); - return invert(T_exer(1. / d_bar, 1. / d, point)) * - T_exer(1. / d_bar, 1. / d, point); + const auto d = this->d(point); + const auto d_bar = this->d_bar(point); + return invert(rotation(1. / d_bar, 1. / d, point)) * + rotation(1. / d_bar, 1. / d, point); } // @sect4{Maxwell Class} // At this point we are ready to instantiate all the major functions of - // the finite element program and also a list of variables. + // the finite element program and also a list of variables. Most of these + // an exact copy of the functions in the tutorial programs. In addition, + // we instatiate the parameters and the perfectly matched layer. The + // default values of these parameters are set to show us a standing wave + // with absorbing boundary conditions and a PML. template class Maxwell : public ParameterAcceptor @@ -380,7 +392,7 @@ namespace Step81 unsigned int refinements; unsigned int fe_order; unsigned int quadrature_order; - bool absorbing_boundary; + bool absorbing_boundary; void parse_parameters_callback(); void make_grid(); @@ -435,11 +447,11 @@ namespace Step81 add_parameter("quadrature order", quadrature_order, "order of the quadrature"); - + absorbing_boundary = true; add_parameter("absorbing boundary condition", - absorbing_boundary, - "use absorbing boundary conditions?"); + absorbing_boundary, + "use absorbing boundary conditions?"); } @@ -459,24 +471,25 @@ namespace Step81 { GridGenerator::hyper_cube(triangulation, -scaling, scaling); triangulation.refine_global(refinements); - - if (!absorbing_boundary){ + + if (!absorbing_boundary) + { for (auto &face : triangulation.active_face_iterators()) - if (face->at_boundary()) - face->set_boundary_id(1); - }; - + if (face->at_boundary()) + face->set_boundary_id(1); + }; + for (auto &cell : triangulation.active_cell_iterators()) if (cell->center()[1] > 0.) cell->set_material_id(1); else cell->set_material_id(2); - + std::cout << "Number of active cells: " << triangulation.n_active_cells() << std::endl; } - + // Enumerate all the degrees of freedom and set up matrix and vector // objects to hold the system data. Enumerating is done by using // DoFHandler::distribute_dofs(). @@ -532,17 +545,20 @@ namespace Step81 } - // Assemble the stiffness matrix and the right-hand side: - //\f{align*}{ - // A_{ij} = \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 - // - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_i)_T) \cdot (\bar{\varphi}_j)_T\text{do}x - // - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_i)_T) \cdot (\nabla\times(\bar{\varphi}_j)_T)\text{d}x, - // \f} - // \f{align}{ - // F_i = 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 addition, we will be modifying the coefficients if the position of the cell is within the PML region. + // Assemble the stiffness matrix and the right-hand side: + //\f{align*}{ + // A_{ij} = \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 + // - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_i)_T) \cdot + // (\bar{\varphi}_j)_T\text{do}x + // - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_i)_T) \cdot + // (\nabla\times(\bar{\varphi}_j)_T)\text{d}x, \f} \f{align}{ + // F_i = 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 addition, we will be modifying the coefficients if the position of the + // cell is within the PML region. template void Maxwell::assemble_system() @@ -561,7 +577,7 @@ namespace Step81 update_quadrature_points | update_normal_vectors | update_JxW_values); - + const unsigned int dofs_per_cell = fe->dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); @@ -571,6 +587,22 @@ namespace Step81 Vector cell_rhs(dofs_per_cell); std::vector local_dof_indices(dofs_per_cell); + // This is assembling the interior of the domain on the left hand side. + // So we are assembling + // //\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 + // \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. + // \f} + // In doing so, we need test functions $\phi_i$ and $\phi_j$, and the curl + // of these test variables. We must be careful with the signs of the + // imaginary parts of these comples test variables. Moreover, we have a + // conditional that changes the parameters if the cell is in the PML region. for (const auto &cell : dof_handler.active_cell_iterators()) { fe_values.reinit(cell); @@ -587,9 +619,9 @@ namespace Step81 for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { - const Point &position = quadrature_points[q_point]; - const auto radius = position.norm(); - const auto inner_radius = perfectly_matched_layer.inner_radius; + const Point &position = quadrature_points[q_point]; + const auto radius = position.norm(); + const auto inner_radius = perfectly_matched_layer.inner_radius; auto mu_inv = parameters.mu_inv(position, id); auto epsilon = parameters.epsilon(position, id); @@ -599,7 +631,7 @@ namespace Step81 { auto A = perfectly_matched_layer.a_matrix(position); auto B = perfectly_matched_layer.b_matrix(position); - auto d = perfectly_matched_layer.d_tensor(position); + auto d = perfectly_matched_layer.d(position); mu_inv = mu_inv / d; epsilon = invert(A) * epsilon * invert(B); @@ -632,13 +664,23 @@ namespace Step81 } } + // Now we assemble the face and the boundary. The following loops will + // assemble + // //\f{align*}{ + // - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_i)_T) \cdot + // (\bar{\varphi}_j)_T\text{do}x \f} and \f{align}{ + // - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_i)_T) + // \cdot (\nabla\times(\bar{\varphi}_j)_T)\text{d}x, + // \f} + // respectively. The test variables and the PML are implemented + // similarly as the domain. for (const auto &face : cell->face_iterators()) { if (face->at_boundary()) { const auto id = face->boundary_id(); - if (id !=0) - { + if (id != 0) + { fe_face_values.reinit(cell, face); FEValuesViews::Vector real_part(fe_face_values, 0); FEValuesViews::Vector imag_part(fe_face_values, dim); @@ -646,10 +688,9 @@ namespace Step81 for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) { - const Point position = - quadrature_points[q_point]; - const auto radius = position.norm(); - const auto inner_radius = + const auto &position = quadrature_points[q_point]; + const auto radius = position.norm(); + const auto inner_radius = perfectly_matched_layer.inner_radius; auto mu_inv = parameters.mu_inv(position, id); @@ -659,18 +700,20 @@ namespace Step81 { auto A = perfectly_matched_layer.a_matrix(position); auto B = perfectly_matched_layer.b_matrix(position); - auto d = perfectly_matched_layer.d_tensor(position); + auto d = perfectly_matched_layer.d(position); mu_inv = mu_inv / d; epsilon = invert(A) * epsilon * invert(B); }; - const auto normal = fe_face_values.normal_vector(q_point); + const auto normal = + fe_face_values.normal_vector(q_point); for (unsigned int i = 0; i < dofs_per_cell; ++i) { - const auto phi_i = real_part.value(i, q_point) - - 1.0i * imag_part.value(i, q_point); + const auto phi_i = + real_part.value(i, q_point) - + 1.0i * imag_part.value(i, q_point); const auto phi_i_T = tangential_part(phi_i, normal); for (unsigned int j = 0; j < dofs_per_cell; ++j) @@ -686,13 +729,13 @@ namespace Step81 const auto sqrt_prod = prod; const auto temp = - -1.0i * - scalar_product((sqrt_prod * phi_j_T), phi_i_T); + -1.0i * scalar_product((sqrt_prod * phi_j_T), + phi_i_T); cell_matrix(i, j) += temp.real(); } /* j */ - } /* i */ - } /* q_point */ - } + } /* i */ + } /* q_point */ + } } else { @@ -712,10 +755,9 @@ namespace Step81 for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) { - const Point position = - quadrature_points[q_point]; - const auto radius = position.norm(); - const auto inner_radius = + const auto &position = quadrature_points[q_point]; + const auto radius = position.norm(); + const auto inner_radius = perfectly_matched_layer.inner_radius; auto sigma = parameters.sigma(position, id1, id2); @@ -757,7 +799,7 @@ namespace Step81 constraints.distribute_local_to_global( cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); } - }; + } // We use a direct solver from the SparseDirectUMFPACK to solve the system template @@ -768,17 +810,18 @@ namespace Step81 A_direct.vmult(solution, system_rhs); } -// The output is written into a vtk file with 4 components -template -void Maxwell::output_results() -{ + // The output is written into a vtk file with 4 components + template + void Maxwell::output_results() + { DataOut<2> data_out; data_out.attach_dof_handler(dof_handler); - data_out.add_data_vector(solution, {"real_Ex", "real_Ey", "imag_Ex", "imag_Ey"}); + data_out.add_data_vector(solution, + {"real_Ex", "real_Ey", "imag_Ex", "imag_Ey"}); data_out.build_patches(); std::ofstream output("solution.vtk"); data_out.write_vtk(output); -} + } template @@ -793,8 +836,8 @@ void Maxwell::output_results() } // namespace Step81 -// The following main function calls the class step-81(), initializes the ParameterAcceptor, -// and calls the run() function. +// The following main function calls the class step-81(), initializes the +// ParameterAcceptor, and calls the run() function. int main() { -- 2.39.5