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