Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-81.h
Go to the documentation of this file.
1
737,
738 *   types::material_id material)
739 *   {
740 *   return (material == 1 ? epsilon_1 : epsilon_2);
741 *   }
742 *  
743 *   template <int dim>
744 *   std::complex<double> Parameters<dim>::mu_inv(const Point<dim> & /*x*/,
745 *   types::material_id material)
746 *   {
747 *   return (material == 1 ? mu_inv_1 : mu_inv_2);
748 *   }
749 *  
750 *   template <int dim>
751 *   typename Parameters<dim>::rank2_type
752 *   Parameters<dim>::sigma(const Point<dim> & /*x*/,
753 *   types::material_id left,
754 *   types::material_id right)
755 *   {
756 *   return (left == right ? rank2_type() : sigma_tensor);
757 *   }
758 *  
759 *   template <int dim>
760 *   typename Parameters<dim>::rank1_type
761 *   Parameters<dim>::J_a(const Point<dim> &point, types::material_id /*id*/)
762 *   {
763 *   rank1_type J_a;
764 *   const auto distance = (dipole_position - point).norm() / dipole_radius;
765 *   if (distance > 1.)
766 *   return J_a;
767 *   double scale = std::cos(distance * M_PI / 2.) *
768 *   std::cos(distance * M_PI / 2.) / (M_PI / 2. - 2. / M_PI) /
769 *   dipole_radius / dipole_radius;
770 *   J_a = dipole_strength * dipole_orientation * scale;
771 *   return J_a;
772 *   }
773 *  
774 * @endcode
775 *
776 *
777 * <a name="step_81-PerfectlyMatchedLayerClass"></a>
778 * <h4>PerfectlyMatchedLayer Class</h4>
779 * The PerfectlyMatchedLayer class inherits ParameterAcceptor as well. It
780 * implements the transformation matrices used to modify the permittivity
781 * and permeability tensors supplied from the Parameters class. The
782 * actual transformation of the material tensors will be done in the
783 * assembly loop. The radii and the strength of the PML is specified, and
784 * the coefficients will be modified using transformation matrices within
785 * the PML region. The radii and strength of the PML are editable through
786 * a .prm file. The rotation function @f$T_{exer}@f$ is the same as
787 * introduced in the perfectly matched layer section of the introduction.
788 * Similarly, the matrices A, B and C are defined as follows
789 * @f[
790 * A = T_{e_xe_r}^{-1}
791 * \text{diag}\left(\frac{1}{\bar{d}^2},\frac{1}{d\bar{d}}\right)T_{e_xe_r},\qquad
792 * B = T_{e_xe_r}^{-1} \text{diag}\left(d,\bar{d}\right)T_{e_xe_r},\qquad
793 * C = T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}},\frac{1}{d}\right)
794 * T_{e_xe_r}.\qquad
795 * @f]
796 *
797
798 *
799 *
800 * @code
801 *   template <int dim>
802 *   class PerfectlyMatchedLayer : public ParameterAcceptor
803 *   {
804 *   public:
805 *   static_assert(dim == 2,
806 *   "The perfectly matched layer is only implemented in 2d.");
807 *  
808 *   Parameters<dim> parameters;
809 *  
810 *   using rank1_type = Tensor<1, dim, std::complex<double>>;
811 *  
812 *   using rank2_type = Tensor<2, dim, std::complex<double>>;
813 *  
814 *   PerfectlyMatchedLayer();
815 *  
816 *   std::complex<double> d(const Point<dim> point);
817 *  
818 *   std::complex<double> d_bar(const Point<dim> point);
819 *  
820 *  
821 *   rank2_type rotation(std::complex<double> d_1,
822 *   std::complex<double> d_2,
823 *   Point<dim> point);
824 *  
825 *   rank2_type a_matrix(const Point<dim> point);
826 *  
827 *   rank2_type b_matrix(const Point<dim> point);
828 *  
829 *   rank2_type c_matrix(const Point<dim> point);
830 *  
831 *   private:
832 *   double inner_radius;
833 *   double outer_radius;
834 *   double strength;
835 *   };
836 *  
837 *  
838 *   template <int dim>
839 *   PerfectlyMatchedLayer<dim>::PerfectlyMatchedLayer()
840 *   : ParameterAcceptor("PerfectlyMatchedLayer")
841 *   {
842 *   inner_radius = 12.;
843 *   add_parameter("inner radius",
844 *   inner_radius,
845 *   "inner radius of the PML shell");
846 *   outer_radius = 20.;
847 *   add_parameter("outer radius",
848 *   outer_radius,
849 *   "outer radius of the PML shell");
850 *   strength = 8.;
851 *   add_parameter("strength", strength, "strength of the PML");
852 *   }
853 *  
854 *  
855 *   template <int dim>
856 *   typename std::complex<double>
857 *   PerfectlyMatchedLayer<dim>::d(const Point<dim> point)
858 *   {
859 *   const auto radius = point.norm();
860 *   if (radius > inner_radius)
861 *   {
862 *   const double s =
863 *   strength * ((radius - inner_radius) * (radius - inner_radius)) /
864 *   ((outer_radius - inner_radius) * (outer_radius - inner_radius));
865 *   return {1.0, s};
866 *   }
867 *   else
868 *   {
869 *   return 1.0;
870 *   }
871 *   }
872 *  
873 *  
874 *   template <int dim>
875 *   typename std::complex<double>
876 *   PerfectlyMatchedLayer<dim>::d_bar(const Point<dim> point)
877 *   {
878 *   const auto radius = point.norm();
879 *   if (radius > inner_radius)
880 *   {
881 *   const double s_bar =
882 *   strength / 3. *
883 *   ((radius - inner_radius) * (radius - inner_radius) *
884 *   (radius - inner_radius)) /
885 *   (radius * (outer_radius - inner_radius) *
886 *   (outer_radius - inner_radius));
887 *   return {1.0, s_bar};
888 *   }
889 *   else
890 *   {
891 *   return 1.0;
892 *   }
893 *   }
894 *  
895 *  
896 *   template <int dim>
897 *   typename PerfectlyMatchedLayer<dim>::rank2_type
898 *   PerfectlyMatchedLayer<dim>::rotation(std::complex<double> d_1,
899 *   std::complex<double> d_2,
900 *   Point<dim> point)
901 *   {
902 *   rank2_type result;
903 *   result[0][0] = point[0] * point[0] * d_1 + point[1] * point[1] * d_2;
904 *   result[0][1] = point[0] * point[1] * (d_1 - d_2);
905 *   result[1][0] = point[0] * point[1] * (d_1 - d_2);
906 *   result[1][1] = point[1] * point[1] * d_1 + point[0] * point[0] * d_2;
907 *   return result;
908 *   }
909 *  
910 *  
911 *   template <int dim>
912 *   typename PerfectlyMatchedLayer<dim>::rank2_type
913 *   PerfectlyMatchedLayer<dim>::a_matrix(const Point<dim> point)
914 *   {
915 *   const auto d = this->d(point);
916 *   const auto d_bar = this->d_bar(point);
917 *   return invert(rotation(d * d, d * d_bar, point)) *
918 *   rotation(d * d, d * d_bar, point);
919 *   }
920 *  
921 *  
922 *   template <int dim>
923 *   typename PerfectlyMatchedLayer<dim>::rank2_type
924 *   PerfectlyMatchedLayer<dim>::b_matrix(const Point<dim> point)
925 *   {
926 *   const auto d = this->d(point);
927 *   const auto d_bar = this->d_bar(point);
928 *   return invert(rotation(d, d_bar, point)) * rotation(d, d_bar, point);
929 *   }
930 *  
931 *  
932 *   template <int dim>
933 *   typename PerfectlyMatchedLayer<dim>::rank2_type
934 *   PerfectlyMatchedLayer<dim>::c_matrix(const Point<dim> point)
935 *   {
936 *   const auto d = this->d(point);
937 *   const auto d_bar = this->d_bar(point);
938 *   return invert(rotation(1. / d_bar, 1. / d, point)) *
939 *   rotation(1. / d_bar, 1. / d, point);
940 *   }
941 *  
942 *  
943 * @endcode
944 *
945 *
946 * <a name="step_81-MaxwellClass"></a>
947 * <h4>Maxwell Class</h4>
948 * At this point we are ready to declare all the major building blocks of
949 * the finite element program which consists of the usual setup and
950 * assembly routines. Most of the structure has already been introduced
951 * in previous tutorial programs. The Maxwell class also holds private
952 * instances of the Parameters and PerfectlyMatchedLayers classes
953 * introduced above. The default values of these parameters are set to
954 * show us a standing wave with absorbing boundary conditions and a PML.
955 *
956
957 *
958 *
959 * @code
960 *   template <int dim>
961 *   class Maxwell : public ParameterAcceptor
962 *   {
963 *   public:
964 *   Maxwell();
965 *   void run();
966 *  
967 *   private:
968 *   /* run time parameters */
969 *   double scaling;
970 *   unsigned int refinements;
971 *   unsigned int fe_order;
972 *   unsigned int quadrature_order;
973 *   bool absorbing_boundary;
974 *  
975 *   void parse_parameters_callback();
976 *   void make_grid();
977 *   void setup_system();
978 *   void assemble_system();
979 *   void solve();
980 *   void output_results();
981 *  
982 *   Parameters<dim> parameters;
983 *   PerfectlyMatchedLayer<dim> perfectly_matched_layer;
984 *  
986 *   DoFHandler<dim> dof_handler;
987 *  
988 *   std::unique_ptr<FiniteElement<dim>> fe;
989 *  
990 *   AffineConstraints<double> constraints;
991 *   SparsityPattern sparsity_pattern;
992 *   SparseMatrix<double> system_matrix;
993 *   Vector<double> solution;
994 *   Vector<double> system_rhs;
995 *   };
996 *  
997 * @endcode
998 *
999 *
1000 * <a name="step_81-ClassTemplateDefinitionsandImplementation"></a>
1001 * <h3>Class Template Definitions and Implementation</h3>
1002 *
1003
1004 *
1005 *
1006 * <a name="step_81-TheConstructor"></a>
1007 * <h4>The Constructor</h4>
1008 * The Constructor simply consists of default initialization a number of
1009 * discretization parameters (such as the domain size, mesh refinement,
1010 * and the order of finite elements and quadrature) and declaring a
1011 * corresponding entry via ParameterAcceptor::add_parameter(). All of
1012 * these can be modified by editing the .prm file. Absorbing boundary
1013 * conditions can be controlled with the absorbing_boundary boolean. If
1014 * absorbing boundary conditions are disabled we simply enforce
1015 * homogeneous Dirichlet conditions on the tangential component of the
1016 * electric field. In the context of time-harmonic Maxwell's equations
1017 * these are also known as perfectly conducting boundary conditions.
1018 *
1019
1020 *
1021 *
1022 * @code
1023 *   template <int dim>
1024 *   Maxwell<dim>::Maxwell()
1025 *   : ParameterAcceptor("Maxwell")
1026 *   , dof_handler(triangulation)
1027 *   {
1028 *   ParameterAcceptor::parse_parameters_call_back.connect(
1029 *   [&]() { parse_parameters_callback(); });
1030 *  
1031 *   scaling = 20;
1032 *   add_parameter("scaling", scaling, "scale of the hypercube geometry");
1033 *  
1034 *   refinements = 8;
1035 *   add_parameter("refinements",
1036 *   refinements,
1037 *   "number of refinements of the geometry");
1038 *  
1039 *   fe_order = 0;
1040 *   add_parameter("fe order", fe_order, "order of the finite element space");
1041 *  
1042 *   quadrature_order = 1;
1043 *   add_parameter("quadrature order",
1044 *   quadrature_order,
1045 *   "order of the quadrature");
1046 *  
1047 *   absorbing_boundary = true;
1048 *   add_parameter("absorbing boundary condition",
1049 *   absorbing_boundary,
1050 *   "use absorbing boundary conditions?");
1051 *   }
1052 *  
1053 *  
1054 *   template <int dim>
1055 *   void Maxwell<dim>::parse_parameters_callback()
1056 *   {
1057 *   fe = std::make_unique<FESystem<dim>>(FE_NedelecSZ<dim>(fe_order), 2);
1058 *   }
1059 *  
1060 * @endcode
1061 *
1062 * The Maxwell::make_grid() routine creates the mesh for the
1063 * computational domain which in our case is a scaled square domain.
1064 * Additionally, a material interface is introduced by setting the
1065 * material id of the upper half (@f$y>0@f$) to 1 and of the lower half
1066 * (@f$y<0@f$) of the computational domain to 2.
1067 * We are using a block decomposition into real and imaginary matrices
1068 * for the solution matrices. More details on this are available
1069 * under the Results section.
1070 *
1071
1072 *
1073 *
1074 * @code
1075 *   template <int dim>
1076 *   void Maxwell<dim>::make_grid()
1077 *   {
1078 *   GridGenerator::hyper_cube(triangulation, -scaling, scaling);
1079 *   triangulation.refine_global(refinements);
1080 *  
1081 *   if (!absorbing_boundary)
1082 *   {
1083 *   for (auto &face : triangulation.active_face_iterators())
1084 *   if (face->at_boundary())
1085 *   face->set_boundary_id(1);
1086 *   };
1087 *  
1088 *   for (auto &cell : triangulation.active_cell_iterators())
1089 *   if (cell->center()[1] > 0.)
1090 *   cell->set_material_id(1);
1091 *   else
1092 *   cell->set_material_id(2);
1093 *  
1094 *  
1095 *   std::cout << "Number of active cells: " << triangulation.n_active_cells()
1096 *   << std::endl;
1097 *   }
1098 *  
1099 * @endcode
1100 *
1101 * The Maxwell::setup_system() routine follows the usual routine of
1102 * enumerating all the degrees of freedom and setting up the matrix and
1103 * vector objects to hold the system data. Enumerating is done by using
1104 * DoFHandler::distribute_dofs().
1105 *
1106
1107 *
1108 *
1109 * @code
1110 *   template <int dim>
1111 *   void Maxwell<dim>::setup_system()
1112 *   {
1113 *   dof_handler.distribute_dofs(*fe);
1114 *   std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
1115 *   << std::endl;
1116 *  
1117 *   solution.reinit(dof_handler.n_dofs());
1118 *   system_rhs.reinit(dof_handler.n_dofs());
1119 *  
1120 *   constraints.clear();
1121 *  
1122 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1123 *  
1124 *   VectorTools::project_boundary_values_curl_conforming_l2(
1125 *   dof_handler,
1126 *   0, /* real part */
1127 *   Functions::ZeroFunction<dim>(2 * dim),
1128 *   0, /* boundary id */
1129 *   constraints);
1130 *   VectorTools::project_boundary_values_curl_conforming_l2(
1131 *   dof_handler,
1132 *   dim, /* imaginary part */
1133 *   Functions::ZeroFunction<dim>(2 * dim),
1134 *   0, /* boundary id */
1135 *   constraints);
1136 *  
1137 *   constraints.close();
1138 *  
1139 *   DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
1140 *   DoFTools::make_sparsity_pattern(dof_handler,
1141 *   dsp,
1142 *   constraints,
1143 *   /* keep_constrained_dofs = */ true);
1144 *   sparsity_pattern.copy_from(dsp);
1145 *   system_matrix.reinit(sparsity_pattern);
1146 *   }
1147 *  
1148 * @endcode
1149 *
1150 * This is a helper function that takes the tangential component of a tensor.
1151 *
1152 * @code
1153 *   template <int dim>
1154 *   DEAL_II_ALWAYS_INLINE inline Tensor<1, dim, std::complex<double>>
1155 *   tangential_part(const Tensor<1, dim, std::complex<double>> &tensor,
1156 *   const Tensor<1, dim> &normal)
1157 *   {
1158 *   auto result = tensor;
1159 *   result[0] = normal[1] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
1160 *   result[1] = -normal[0] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
1161 *   return result;
1162 *   }
1163 *  
1164 *  
1165 * @endcode
1166 *
1167 * Assemble the stiffness matrix and the right-hand side:
1168 * \f{align*}{
1169 * A_{ij} = \int_\Omega (\mu_r^{-1}\nabla \times \varphi_j) \cdot
1170 * (\nabla\times\bar{\varphi}_i)\text{d}x
1171 * - \int_\Omega \varepsilon_r\varphi_j \cdot \bar{\varphi}_i\text{d}x
1172 * - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_j)_T) \cdot
1173 * (\bar{\varphi}_i)_T\text{do}x
1174 * - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_j)_T) \cdot
1175 * (\nabla\times(\bar{\varphi}_i)_T)\text{d}x, \f} \f{align}{
1176 * F_i = i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x - \int_\Omega
1177 * \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i}) \text{d}x.
1178 * \f}
1179 * In addition, we will be modifying the coefficients if the position of the
1180 * cell is within the PML region.
1181 *
1182
1183 *
1184 *
1185 * @code
1186 *   template <int dim>
1187 *   void Maxwell<dim>::assemble_system()
1188 *   {
1189 *   const QGauss<dim> quadrature_formula(quadrature_order);
1190 *   const QGauss<dim - 1> face_quadrature_formula(quadrature_order);
1191 *  
1192 *   FEValues<dim, dim> fe_values(*fe,
1193 *   quadrature_formula,
1194 *   update_values | update_gradients |
1195 *   update_quadrature_points |
1196 *   update_JxW_values);
1197 *   FEFaceValues<dim, dim> fe_face_values(*fe,
1198 *   face_quadrature_formula,
1199 *   update_values | update_gradients |
1200 *   update_quadrature_points |
1201 *   update_normal_vectors |
1202 *   update_JxW_values);
1203 *  
1204 *   const unsigned int dofs_per_cell = fe->dofs_per_cell;
1205 *  
1206 *   const unsigned int n_q_points = quadrature_formula.size();
1207 *   const unsigned int n_face_q_points = face_quadrature_formula.size();
1208 *  
1209 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1210 *   Vector<double> cell_rhs(dofs_per_cell);
1211 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1212 *  
1213 * @endcode
1214 *
1215 * Next, let us assemble on the interior of the domain on the left hand
1216 * side. So we are computing
1217 * \f{align*}{
1218 * \int_\Omega (\mu_r^{-1}\nabla \times \varphi_i) \cdot
1219 * (\nabla\times\bar{\varphi}_j)\text{d}x
1220 * -
1221 * \int_\Omega \varepsilon_r\varphi_i \cdot \bar{\varphi}_j\text{d}x
1222 * \f}
1223 * and
1224 * \f{align}{
1225 * i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x
1226 * - \int_\Omega \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i})
1227 * \text{d}x.
1228 * \f}
1229 * In doing so, we need test functions @f$\varphi_i@f$ and @f$\varphi_j@f$, and the
1230 * curl of these test variables. We must be careful with the signs of the
1231 * imaginary parts of these complex test variables. Moreover, we have a
1232 * conditional that changes the parameters if the cell is in the PML region.
1233 *
1234 * @code
1235 *   const FEValuesExtractors::Vector real_part(0);
1236 *   const FEValuesExtractors::Vector imag_part(dim);
1237 *   for (const auto &cell : dof_handler.active_cell_iterators())
1238 *   {
1239 *   fe_values.reinit(cell);
1240 *  
1241 *   cell_matrix = 0.;
1242 *   cell_rhs = 0.;
1243 *  
1244 *   cell->get_dof_indices(local_dof_indices);
1245 *   const auto id = cell->material_id();
1246 *  
1247 *   const auto &quadrature_points = fe_values.get_quadrature_points();
1248 *  
1249 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1250 *   {
1251 *   const Point<dim> &position = quadrature_points[q_point];
1252 *  
1253 *   auto mu_inv = parameters.mu_inv(position, id);
1254 *   auto epsilon = parameters.epsilon(position, id);
1255 *   const auto J_a = parameters.J_a(position, id);
1256 *  
1257 *   const auto A = perfectly_matched_layer.a_matrix(position);
1258 *   const auto B = perfectly_matched_layer.b_matrix(position);
1259 *   const auto d = perfectly_matched_layer.d(position);
1260 *  
1261 *   mu_inv = mu_inv / d;
1262 *   epsilon = invert(A) * epsilon * invert(B);
1263 *  
1264 *   for (const auto i : fe_values.dof_indices())
1265 *   {
1266 *   constexpr std::complex<double> imag{0., 1.};
1267 *  
1268 *   const auto phi_i =
1269 *   fe_values[real_part].value(i, q_point) -
1270 *   imag * fe_values[imag_part].value(i, q_point);
1271 *   const auto curl_phi_i =
1272 *   fe_values[real_part].curl(i, q_point) -
1273 *   imag * fe_values[imag_part].curl(i, q_point);
1274 *  
1275 *   const auto rhs_value =
1276 *   (imag * scalar_product(J_a, phi_i)) * fe_values.JxW(q_point);
1277 *   cell_rhs(i) += rhs_value.real();
1278 *  
1279 *   for (const auto j : fe_values.dof_indices())
1280 *   {
1281 *   const auto phi_j =
1282 *   fe_values[real_part].value(j, q_point) +
1283 *   imag * fe_values[imag_part].value(j, q_point);
1284 *   const auto curl_phi_j =
1285 *   fe_values[real_part].curl(j, q_point) +
1286 *   imag * fe_values[imag_part].curl(j, q_point);
1287 *  
1288 *   const auto temp =
1289 *   (scalar_product(mu_inv * curl_phi_j, curl_phi_i) -
1290 *   scalar_product(epsilon * phi_j, phi_i)) *
1291 *   fe_values.JxW(q_point);
1292 *   cell_matrix(i, j) += temp.real();
1293 *   }
1294 *   }
1295 *   }
1296 *  
1297 * @endcode
1298 *
1299 * Now we assemble the face and the boundary. The following loops will
1300 * assemble
1301 * \f{align*}{
1302 * - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_i)_T) \cdot
1303 * (\bar{\varphi}_j)_T\text{do}x \f} and \f{align}{
1304 * - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_i)_T)
1305 * \cdot (\nabla\times(\bar{\varphi}_j)_T)\text{d}x,
1306 * \f}
1307 * respectively. The test variables and the PML are implemented
1308 * similarly as the domain.
1309 *
1310 * @code
1311 *   const FEValuesExtractors::Vector real_part(0);
1312 *   const FEValuesExtractors::Vector imag_part(dim);
1313 *   for (const auto &face : cell->face_iterators())
1314 *   {
1315 *   if (face->at_boundary())
1316 *   {
1317 *   const auto id = face->boundary_id();
1318 *   if (id != 0)
1319 *   {
1320 *   fe_face_values.reinit(cell, face);
1321 *  
1322 *   for (unsigned int q_point = 0; q_point < n_face_q_points;
1323 *   ++q_point)
1324 *   {
1325 *   const auto &position = quadrature_points[q_point];
1326 *  
1327 *   auto mu_inv = parameters.mu_inv(position, id);
1328 *   auto epsilon = parameters.epsilon(position, id);
1329 *  
1330 *   const auto A =
1331 *   perfectly_matched_layer.a_matrix(position);
1332 *   const auto B =
1333 *   perfectly_matched_layer.b_matrix(position);
1334 *   const auto d = perfectly_matched_layer.d(position);
1335 *  
1336 *   mu_inv = mu_inv / d;
1337 *   epsilon = invert(A) * epsilon * invert(B);
1338 *  
1339 *   const auto normal =
1340 *   fe_face_values.normal_vector(q_point);
1341 *  
1342 *   for (const auto i : fe_face_values.dof_indices())
1343 *   {
1344 *   constexpr std::complex<double> imag{0., 1.};
1345 *  
1346 *   const auto phi_i =
1347 *   fe_face_values[real_part].value(i, q_point) -
1348 *   imag *
1349 *   fe_face_values[imag_part].value(i, q_point);
1350 *   const auto phi_i_T = tangential_part(phi_i, normal);
1351 *  
1352 *   for (const auto j : fe_face_values.dof_indices())
1353 *   {
1354 *   const auto phi_j =
1355 *   fe_face_values[real_part].value(j, q_point) +
1356 *   imag *
1357 *   fe_face_values[imag_part].value(j, q_point);
1358 *   const auto phi_j_T =
1359 *   tangential_part(phi_j, normal) *
1360 *   fe_face_values.JxW(q_point);
1361 *  
1362 *   const auto prod = mu_inv * epsilon;
1363 *   const auto sqrt_prod = prod;
1364 *  
1365 *   const auto temp =
1366 *   -imag * scalar_product((sqrt_prod * phi_j_T),
1367 *   phi_i_T);
1368 *   cell_matrix(i, j) += temp.real();
1369 *   } /* j */
1370 *   } /* i */
1371 *   } /* q_point */
1372 *   }
1373 *   }
1374 *   else
1375 *   {
1376 * @endcode
1377 *
1378 * We are on an interior face:
1379 *
1380 * @code
1381 *   const auto face_index = cell->face_iterator_to_index(face);
1382 *  
1383 *   const auto id1 = cell->material_id();
1384 *   const auto id2 = cell->neighbor(face_index)->material_id();
1385 *  
1386 *   if (id1 == id2)
1387 *   continue; /* skip this face */
1388 *  
1389 *   fe_face_values.reinit(cell, face);
1390 *  
1391 *   for (unsigned int q_point = 0; q_point < n_face_q_points;
1392 *   ++q_point)
1393 *   {
1394 *   const auto &position = quadrature_points[q_point];
1395 *  
1396 *   auto sigma = parameters.sigma(position, id1, id2);
1397 *  
1398 *   const auto B = perfectly_matched_layer.b_matrix(position);
1399 *   const auto C = perfectly_matched_layer.c_matrix(position);
1400 *   sigma = invert(C) * sigma * invert(B);
1401 *  
1402 *   const auto normal = fe_face_values.normal_vector(q_point);
1403 *  
1404 *   for (const auto i : fe_face_values.dof_indices())
1405 *   {
1406 *   constexpr std::complex<double> imag{0., 1.};
1407 *  
1408 *   const auto phi_i =
1409 *   fe_face_values[real_part].value(i, q_point) -
1410 *   imag * fe_face_values[imag_part].value(i, q_point);
1411 *   const auto phi_i_T = tangential_part(phi_i, normal);
1412 *  
1413 *   for (const auto j : fe_face_values.dof_indices())
1414 *   {
1415 *   const auto phi_j =
1416 *   fe_face_values[real_part].value(j, q_point) +
1417 *   imag *
1418 *   fe_face_values[imag_part].value(j, q_point);
1419 *   const auto phi_j_T = tangential_part(phi_j, normal);
1420 *  
1421 *   const auto temp =
1422 *   -imag *
1423 *   scalar_product((sigma * phi_j_T), phi_i_T) *
1424 *   fe_face_values.JxW(q_point);
1425 *   cell_matrix(i, j) += temp.real();
1426 *   } /* j */
1427 *   } /* i */
1428 *   } /* q_point */
1429 *   }
1430 *   }
1431 *  
1432 *   constraints.distribute_local_to_global(
1433 *   cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1434 *   }
1435 *   }
1436 *  
1437 * @endcode
1438 *
1439 * We use a direct solver from the SparseDirectUMFPACK to solve the system
1440 *
1441 * @code
1442 *   template <int dim>
1443 *   void Maxwell<dim>::solve()
1444 *   {
1445 *   SparseDirectUMFPACK A_direct;
1446 *   A_direct.initialize(system_matrix);
1447 *   A_direct.vmult(solution, system_rhs);
1448 *   }
1449 *  
1450 * @endcode
1451 *
1452 * The output is written into a vtk file with 4 components
1453 *
1454 * @code
1455 *   template <int dim>
1456 *   void Maxwell<dim>::output_results()
1457 *   {
1458 *   DataOut<2> data_out;
1459 *   data_out.attach_dof_handler(dof_handler);
1460 *   data_out.add_data_vector(solution,
1461 *   {"real_Ex", "real_Ey", "imag_Ex", "imag_Ey"});
1462 *   data_out.build_patches();
1463 *   std::ofstream output("solution.vtk");
1464 *   data_out.write_vtk(output);
1465 *   }
1466 *  
1467 *  
1468 *   template <int dim>
1469 *   void Maxwell<dim>::run()
1470 *   {
1471 *   make_grid();
1472 *   setup_system();
1473 *   assemble_system();
1474 *   solve();
1475 *   output_results();
1476 *   }
1477 *  
1478 *   } // namespace Step81
1479 *  
1480 * @endcode
1481 *
1482 * The following main function calls the class @ref step_81 "step-81"(), initializes the
1483 * ParameterAcceptor, and calls the run() function.
1484 *
1485
1486 *
1487 *
1488 * @code
1489 *   int main()
1490 *   {
1491 *   try
1492 *   {
1493 *   using namespace dealii;
1494 *  
1495 *   Step81::Maxwell<2> maxwell_2d;
1496 *   ParameterAcceptor::initialize("parameters.prm");
1497 *   maxwell_2d.run();
1498 *   }
1499 *   catch (std::exception &exc)
1500 *   {
1501 *   std::cerr << std::endl
1502 *   << std::endl
1503 *   << "----------------------------------------------------"
1504 *   << std::endl;
1505 *   std::cerr << "Exception on processing: " << std::endl
1506 *   << exc.what() << std::endl
1507 *   << "Aborting!" << std::endl
1508 *   << "----------------------------------------------------"
1509 *   << std::endl;
1510 *   return 1;
1511 *   }
1512 *   catch (...)
1513 *   {
1514 *   std::cerr << std::endl
1515 *   << std::endl
1516 *   << "----------------------------------------------------"
1517 *   << std::endl;
1518 *   std::cerr << "Unknown exception!" << std::endl
1519 *   << "Aborting!" << std::endl
1520 *   << "----------------------------------------------------"
1521 *   << std::endl;
1522 *   return 1;
1523 *   }
1524 *   return 0;
1525 *   }
1526 * @endcode
1527<a name="step_81-Results"></a><h1>Results</h1>
1528
1529
1530The solution is written to a .vtk file with four components. These are the
1531real and imaginary parts of the @f$E_x@f$ and @f$E_y@f$ solution waves. With the
1532current setup, the output should read
1533
1534@code
1535Number of active cells: 4096
1536Number of degrees of freedom: 16640
1537Program ended with exit code: 0
1538@endcode
1539
1540<a name="step_81-AbsorbingboundaryconditionsandthePML"></a><h3> Absorbing boundary conditions and the PML </h3>
1541
1542
1543The following images are the outputs for the imaginary @f$E_x@f$ without the
1544interface and with the dipole centered at @f$(0,0)@f$. In order to remove the
1545interface, the surface conductivity is set to 0. First, we turn off the
1546absorbing boundary conditions and the PML. Second, we want to see the
1547effect of the PML when absorbing boundary conditions apply. So we set
1548absorbing boundary conditions to true and leave the PML strength to 0.
1549Lastly, we increase the strength of the PML to 4. Change the following in
1550the .prm file :
1551
1552@code
1553# use absorbing boundary conditions?
1554 set absorbing boundary condition = false
1555
1556# position of the dipole
1557 set dipole position = 0, 0
1558
1559# strength of the PML
1560 set strength = 0
1561
1562# surface conductivity between material 1 and material 2
1563 set sigma = 0, 0; 0, 0| 0, 0; 0, 0
1564@endcode
1565
1566Following are the output images:
1567
1568<table width="80%" align="center">
1569 <tr>
1570 <td align="center">
1571 <img src="https://www.dealii.org/images/steps/developer/step-81-nointerface_noabs_PML0.png" alt="Visualization of the solution of step-81 with no interface, Dirichlet boundary conditions and PML strength 0" height="210"/>
1572 <p> Solution with no interface, Dirichlet boundary conditions and PML strength 0.</p>
1573 </td>
1574 <td></td>
1575 <td align="center">
1576 <img src="https://www.dealii.org/images/steps/developer/step-81-nointerface_abs_PML0.png" alt="Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 0" height="210">
1577 <p> Solution with no interface, absorbing boundary conditions and PML strength 0.</p>
1578 </td>
1579 <td></td>
1580 <td align="center">
1581 <img src="https://www.dealii.org/images/steps/developer/step-81-nointerface_abs_PML4.png" alt="Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 4" height="210">
1582 <p> Solution with no interface, absorbing boundary conditions and PML strength 4.</p>
1583 </td>
1584 </tr>
1585</table>
1586
1587We observe that with absorbing boundary conditions and in absence of the
1588PML, there is a lot of distortion and resonance (the real parts will not be
1589generated without a PML). This is, as we stipulated, due to reflection from
1590infinity. As we see, a much more coherent image is generated with an
1591appropriate PML.
1592
1593<a name="step_81-SurfacePlasmonPolariton"></a><h3> Surface Plasmon Polariton </h3>
1594
1595Now, let's generate a standing wave by adding an interface at the center.
1596In order to observe this effect, we offset the center of the dipole to @f$(0,
15970.8)@f$ and set the surface conductivity back to @f$(0.001, 0.2)@f$:
1598
1599@code
1600# position of the dipole
1601 set dipole position = 0, 0.8
1602
1603# surface conductivity between material 1 and material 2
1604 set sigma = 0.001, 0.2; 0, 0| 0, 0; 0.001, 0.2
1605@endcode
1606
1607Once again, we will visualize the output with absorbing boundary conditions
1608and PML strength 0 and with absorbing boundary conditions and PML strength
16094. The following tables are the imaginary part of @f$E_x@f$ and the real part
1610of @f$E_x@f$.
1611
1612<table width="80%" align="center">
1613 <tr>
1614 <td align="center">
1615 <img src="https://www.dealii.org/images/steps/developer/step-81-imagEx_noabs_PML0.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
1616 <p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
1617 </td>
1618 <td></td>
1619 <td align="center">
1620 <img src="https://www.dealii.org/images/steps/developer/step-81-imagEx_abs_PML0.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
1621 <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
1622 </td>
1623 <td></td>
1624 <td align="center">
1625 <img src="https://www.dealii.org/images/steps/developer/step-81-imagEx_abs_PML4.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4" height="210">
1626 <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
1627 </td>
1628 </tr>
1629</table>
1630
1631
1632<table width="80%" align="center">
1633 <tr>
1634 <td align="center">
1635 <img src="https://www.dealii.org/images/steps/developer/step-81-realEx_noabs_PML0.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
1636 <p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
1637 </td>
1638 <td></td>
1639 <td align="center">
1640 <img src="https://www.dealii.org/images/steps/developer/step-81-realEx_abs_PML0.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
1641 <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
1642 </td>
1643 <td></td>
1644 <td align="center">
1645 <img src="https://www.dealii.org/images/steps/developer/step-81-realEx_abs_PML4.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4" height="210">
1646 <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
1647 </td>
1648 </tr>
1649</table>
1650
1651The SPP is confined near the interface that we created, however without
1652absorbing boundary conditions, we don't observe a dissipation effect. On
1653adding the absorbing boundary conditions, we observe distortion and
1654resonance and we still don't notice any dissipation. As expected, the PML
1655removes the distortion and resonance. The standing wave is also dissipating
1656and getting absorbed within the PML, and as we increase the PML strength,
1657the standing wave will dissipate more within the PML ring.
1658
1659Here are some animations to demonstrate the effect of the PML
1660<table width="80%" align="center">
1661 <tr>
1662 <td align="center">
1663 <img src="https://www.dealii.org/images/steps/developer/step-81-dirichlet_Ex.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
1664 <p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
1665 </td>
1666 <td></td>
1667 <td align="center">
1668 <img src="https://www.dealii.org/images/steps/developer/step-81-absorbing_Ex.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
1669 <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
1670 </td>
1671 <td></td>
1672 <td align="center">
1673 <img src="https://www.dealii.org/images/steps/developer/step-81-perfectly_matched_layer_Ex.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4" height="210">
1674 <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
1675 </td>
1676 </tr>
1677</table>
1678
1679
1680<table width="80%" align="center">
1681 <tr>
1682 <td align="center">
1683 <img src="https://www.dealii.org/images/steps/developer/step-81-dirichlet_Ey.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
1684 <p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
1685 </td>
1686 <td></td>
1687 <td align="center">
1688 <img src="https://www.dealii.org/images/steps/developer/step-81-absorbing_Ey.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
1689 <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
1690 </td>
1691 <td></td>
1692 <td align="center">
1693 <img src="https://www.dealii.org/images/steps/developer/step-81-perfectly_matched_layer_Ey.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4" height="210">
1694 <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
1695 </td>
1696 </tr>
1697</table>
1698
1699<a name="step_81-Notes"></a><h3> Notes </h3>
1700
1701
1702<a name="step_81-RealandComplexMatrices"></a><h4> Real and Complex Matrices </h4>
1703
1704As is evident from the results, we are splitting our solution matrices into
1705the real and the imaginary components. We started off using the @f$H^{curl}@f$
1706conforming Nédélec Elements, and we made two copies of the Finite Elements
1707in order to represent the real and the imaginary components of our input
1708(FE_NedelecSZ was used instead of FE_Nedelec to avoid the sign conflicts
1709issues present in traditional Nédélec elements). In the assembly, we create
1710two vectors of dimension @f$dim@f$ that assist us in extracting the real and
1711the imaginary components of our finite elements.
1712
1713
1714<a name="step_81-RotationsandScaling"></a><h4> Rotations and Scaling </h4>
1715
1716As we see in our assembly, our finite element is rotated and scaled as
1717follows:
1718
1719@code
1720const auto phi_i = real_part.value(i, q_point) - 1.0i * imag_part.value(i, q_point);
1721@endcode
1722
1723This @f$\phi_i@f$ variable doesn't need to be scaled in this way, we may choose
1724any arbitrary scaling constants @f$a@f$ and @f$b@f$. If we choose this scaling, the
1725@f$\phi_j@f$ must also be modified with the same scaling, as follows:
1726
1727@code
1728const auto phi_i = a*real_part.value(i, q_point) -
1729 bi * imag_part.value(i, q_point);
1730
1731const auto phi_j = a*real_part.value(i, q_point) +
1732 bi * imag_part.value(i, q_point);
1733@endcode
1734
1735Moreover, the cell_rhs need not be the real part of the rhs_value. Say if
1736we modify to take the imaginary part of the computed rhs_value, we must
1737also modify the cell_matrix accordingly to take the imaginary part of temp.
1738However, making these changes to both sides of the equation will not affect
1739our solution, and we will still be able to generate the surface plasmon
1740polariton.
1741
1742@code
1743cell_rhs(i) += rhs_value.imag();
1744
1745cell_matrix(i) += temp.imag();
1746@endcode
1747
1748<a name="step_81-Postprocessing"></a><h4> Postprocessing </h4>
1749
1750We will create a video demonstrating the wave in motion, which is
1751essentially an implementation of @f$e^{-i\omega t}(Re(E) + i*Im(E))@f$ as we
1752increment time. This is done by slightly changing the output function to
1753generate a series of .vtk files, which will represent out solution wave as
1754we increment time. Introduce an input variable @f$t@f$ in the output_results()
1755class as output_results(unsigned int t). Then change the class itself to
1756the following:
1757
1758@code
1759template <int dim>
1760void Maxwell<dim>::output_results(unsigned int t)
1761{
1762 std::cout << "Running step:" << t << std::endl;
1763 DataOut<2> data_out;
1764 data_out.attach_dof_handler(dof_handler);
1765 Vector<double> postprocessed;
1766 postprocessed.reinit(solution);
1767 for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
1768 {
1769 if (i % 4 == 0)
1770 {
1771 postprocessed[i] = std::cos(2 * M_PI * 0.04 * t) * solution[i] -
1772 std::sin(2 * M_PI * 0.04 * t) * solution[i + 1];
1773 }
1774 else if (i % 4 == 2)
1775 {
1776 postprocessed[i] = std::cos(2 * M_PI * 0.04 * t) * solution[i] -
1777 std::sin(2 * M_PI * 0.04 * t) * solution[i + 1];
1778 }
1779 }
1780 data_out.add_data_vector(postprocessed, {"E_x", "E_y", "null0", "null1"});
1781 data_out.build_patches();
1782 const std::string filename =
1783 "solution-" + Utilities::int_to_string(t) + ".vtk";
1784 std::ofstream output(filename);
1785 data_out.write_vtk(output);
1786 std::cout << "Done running step:" << t << std::endl;
1787}
1788@endcode
1789
1790Finally, in the run() function, replace output_results() with
1791@code
1792for (int t = 0; t <= 100; t++)
1793 {
1794 output_results(t);
1795 }
1796@endcode
1797
1798This would generate 100 solution .vtk files, which can be opened in a group
1799on Paraview and then can be saved as an animation. We used FFMPEG to
1800generate gifs.
1801
1802<a name="step_81-PossibilitiesforExtension"></a><h3> Possibilities for Extension </h3>
1803
1804
1805The example step could be extended in a number of different directions.
1806<ul>
1807 <li>
1808 The current program uses a direct solver to solve the linear system.
1809 This is efficient for two spatial dimensions where scattering problems
1810 up to a few millions degrees of freedom can be solved. In 3D, however,
1811 the increased stencil size of the Nedelec element pose a severe
1812 limiting factor on the problem size that can be computed. As an
1813 alternative, the idea to use iterative solvers can be entertained.
1814 This, however requires specialized preconditioners. For example, just
1815 using an iterative Krylov space solver (such as SolverGMRES) on above
1816 problem will requires many thousands of iterations to converge.
1817 Unfortunately, time-harmonic Maxwell's equations lack the usual notion
1818 of local smoothing properties, which renders the usual suspects, such
1819 as a geometric multigrid (see the Multigrid class), largely useless. A
1820 possible extension would be to implement an additive Schwarz preconditioner
1821 (based on domain decomposition, see for example
1822 @cite Gopalakrishnan2003), or a sweeping preconditioner (see for
1823 example @cite Ying2012).
1824 </li>
1825 <li>
1826 Another possible extension of the current program is to introduce local
1827 mesh refinement (either based on a residual estimator, or based on the
1828 dual weighted residual method, see @ref step_14 "step-14"). This is in particular of
1829 interest to counter the increased computational cost caused by the
1830 scale separation between the SPP and the dipole.
1831 </li>
1832</ul>
1833 *
1834 *
1835<a name="step_81-PlainProg"></a>
1836<h1> The plain program</h1>
1837@include "step-81.cc"
1838*/
void add_parameter(const std::string &entry, ParameterType &parameter, const std::string &documentation="", ParameterHandler &prm_=prm, const Patterns::PatternBase &pattern= *Patterns::Tools::Convert< ParameterType >::to_pattern())
Definition point.h:111
numbers::NumberTraits< Number >::real_type norm() const
Point< 3 > center
__global__ void set(Number *val, const Number s, const size_type N)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:564
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)