From dc50b905cf9514a4f557f289ba8fc99edc2aede5 Mon Sep 17 00:00:00 2001 From: Daniel Garcia-Sanchez Date: Fri, 13 Sep 2019 10:59:50 +0200 Subject: [PATCH] Add maxwell_singularity tests --- tests/fe/fe_nedelec_singularity_01.cc | 798 ++++++++++++++++++ ...th_petsc_with_complex=true.mpirun=4.output | 103 +++ tests/fe/fe_nedelec_singularity_02.cc | 784 +++++++++++++++++ ...h_petsc_with_complex=false.mpirun=4.output | 103 +++ 4 files changed, 1788 insertions(+) create mode 100644 tests/fe/fe_nedelec_singularity_01.cc create mode 100644 tests/fe/fe_nedelec_singularity_01.with_petsc_with_mumps=true.with_petsc_with_complex=true.mpirun=4.output create mode 100644 tests/fe/fe_nedelec_singularity_02.cc create mode 100644 tests/fe/fe_nedelec_singularity_02.with_petsc_with_mumps=true.with_petsc_with_complex=false.mpirun=4.output diff --git a/tests/fe/fe_nedelec_singularity_01.cc b/tests/fe/fe_nedelec_singularity_01.cc new file mode 100644 index 0000000000..1541ab7d05 --- /dev/null +++ b/tests/fe/fe_nedelec_singularity_01.cc @@ -0,0 +1,798 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- +// +// By Daniel Garcia-Sanchez, CNRS +// +// Test a maxwell singularity in 3D. Maxwell singularities are common in sharp +// metallic edges such as the Fichera corner. Here we test the elements Nedelec +// and NedelecSZ using the L2 norm and the continuity of the solution. +// +// This test solves the complex valued curl-curl equation in 3D: +// +// curl((1/mu_r)curl(E)) +// -omega^2*epsilon_0*mu_0(epsilon_r-(i*sigma/(omega*epsilon_0)))E +// = RightHandSide +// +// The manufactured solution is: +// +// E_x = (((x^2 / sqrt(x^2 + y^2)) * (x^2 - (dimension_x / 2)^2) * +// (y^2 - (dimension_y / 2)^2) * (z^2 - (dimension_z / 2)^2)) / +// ((dimension_x / 2)^3 * (dimension_y / 2)^2 * (dimension_z / 2)^2)) +// E_y = ( ((y^2 / sqrt(x^2 + y^2)) * (x^2 - (dimension_x / 2)^2) * +// (y^2 - (dimension_y / 2)^2) * (z^2 - (dimension_z / 2)^2)) / +// ((dimension_x / 2)^2 * (dimension_y / 2)^3 * (dimension_z / 2)^2)) +// E_z = 10 * (x * (x^2 - (dimension_x / 2)^2) * (y^2 - (dimension_y / 2)^2) / +// ((dimension_x / 2)^2 * (dimension_y / 2)^2)) +// +// This manufactured solution has a singularity at x = y = 0 +// +// The right hand side can be calculated with a symbolic math package such as +// sympy. + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include "../tests.h" + + +using namespace dealii; + +namespace nedelec_singularity +{ + // For the sake of simplicity define the parameters as global variables. + static const double dimension_x = 0.04; + static const double dimension_y = 0.04; + static const double dimension_z = 0.04; + static const double epsilon_0 = 8.85418782e-12; + static const double mu_0 = 1.25663706e-06; + static const double epsilon_r = 1; + static const double mu_r = 1; + static const double sigma = 0.0001; + static const double omega = 6e9 * 2 * numbers::PI; + static unsigned int nb_probe_points = 100; + static unsigned int grid_level = 1; + static unsigned int coarse_mesh_divisions_z = 3; + + + + template + class ExactSolution : public Function + { + public: + ExactSolution(); + virtual double + value(const Point &p, const unsigned int component) const override; + }; + + + + template + ExactSolution::ExactSolution() + : Function(dim) + {} + + + + template + double + ExactSolution::value(const Point & p, + const unsigned int component) const + { + const double R_x = p[0]; + const double R_y = p[1]; + const double R_z = p[2]; + + switch (component) + { + case 0: + return 2 * std::pow(R_x, 2) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) / + (std::pow(dimension_x, 3) * std::pow(dimension_y, 2) * + std::pow(dimension_z, 2) * + std::sqrt(std::pow(R_x, 2) + std::pow(R_y, 2))); + break; + case 1: + return 2 * std::pow(R_x, 2) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) / + (std::pow(dimension_x, 3) * std::pow(dimension_y, 2) * + std::pow(dimension_z, 2) * + std::sqrt(std::pow(R_x, 2) + std::pow(R_y, 2))); + break; + case 2: + return 10 * R_x * (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) / + (std::pow(dimension_x, 2) * std::pow(dimension_y, 2)); + break; + default: + Assert(false, ExcNotImplemented()); + return 0; + } + } + + + + template + class RightHandSide : public Function> + { + public: + RightHandSide(); + virtual std::complex + value(const Point &p, const unsigned int component) const override; + }; + + + + template + RightHandSide::RightHandSide() + : Function>(dim) + {} + + + + template + std::complex + RightHandSide::value(const Point & p, + const unsigned int component) const + { + const double R_x = p[0]; + const double R_y = p[1]; + const double R_z = p[2]; + + const std::complex I(0, 1); + + switch (component) + { + case 0: + return 2. * R_x * + (-R_x * dimension_y * mu_0 * mu_r * omega * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) * + (epsilon_0 * epsilon_r * omega - I * sigma) - + 8 * R_x * dimension_y * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) + + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) * + (-3 * R_x * std::pow(R_y, 2) * dimension_y * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) - + 8 * R_x * dimension_y * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) + + R_x * dimension_y * (std::pow(R_x, 2) + std::pow(R_y, 2)) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (20 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) + + 3 * std::pow(R_y, 3) * dimension_x * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) + + 16 * R_y * dimension_x * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) * + (8 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) + + 2 * R_y * dimension_x * + (std::pow(R_x, 2) + std::pow(R_y, 2)) * + (std::pow(R_y, 2) * (-16 * std::pow(R_x, 2) + + 4 * std::pow(dimension_x, 2)) + + std::pow(R_y, 2) * (-16 * std::pow(R_y, 2) + + 4 * std::pow(dimension_y, 2)) - + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2))))) / + (std::pow(dimension_x, 3) * std::pow(dimension_y, 3) * + std::pow(dimension_z, 2) * mu_r * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 5.0 / 2.0)); + break; + case 1: + return 2. * R_y * + (-R_y * dimension_x * mu_0 * mu_r * omega * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) * + (epsilon_0 * epsilon_r * omega - I * sigma) - + 8 * R_y * dimension_x * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) + + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) * + (3 * std::pow(R_x, 3) * dimension_y * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) - + 3 * std::pow(R_x, 2) * R_y * dimension_x * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) + + 16 * R_x * dimension_y * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) * + (8 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) + + 2 * R_x * dimension_y * + (std::pow(R_x, 2) + std::pow(R_y, 2)) * + (std::pow(R_x, 2) * (-16 * std::pow(R_x, 2) + + 4 * std::pow(dimension_x, 2)) + + std::pow(R_x, 2) * (-16 * std::pow(R_y, 2) + + 4 * std::pow(dimension_y, 2)) - + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2))) - + 8 * R_y * dimension_x * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) + + R_y * dimension_x * (std::pow(R_x, 2) + std::pow(R_y, 2)) * + (20 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)))) / + (std::pow(dimension_x, 3) * std::pow(dimension_y, 3) * + std::pow(dimension_z, 2) * mu_r * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 5.0 / 2.0)); + break; + case 2: + return 2. * + (-5 * R_x * dimension_x * dimension_y * + std::pow(dimension_z, 2) * mu_0 * mu_r * omega * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (epsilon_0 * epsilon_r * omega - I * sigma) + + 8 * R_x * dimension_y * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (-std::pow(R_x, 2) * R_z * + std::sqrt(std::pow(R_x, 2) + std::pow(R_y, 2)) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) + + 2 * R_z * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), + 3.0 / 2.0) * + (8 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) - + 15 * dimension_x * std::pow(dimension_z, 2) * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2)) + + 8 * dimension_x * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (-5 * R_x * dimension_y * std::pow(dimension_z, 2) * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) - + std::pow(R_y, 3) * R_z * + std::sqrt(std::pow(R_x, 2) + std::pow(R_y, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) + + 2 * R_y * R_z * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), + 3.0 / 2.0) * + (8 * std::pow(R_y, 2) - std::pow(dimension_y, 2)))) / + (std::pow(dimension_x, 3) * std::pow(dimension_y, 3) * + std::pow(dimension_z, 2) * mu_r * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2)); + break; + default: + Assert(false, ExcNotImplemented()); + return 0; + } + } + + + + template + class NedelecSingularity + { + public: + NedelecSingularity(); + void + run(); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + output_results(); + + MPI_Comm mpi_communicator; + parallel::distributed::Triangulation triangulation; + const unsigned int fe_order; + const QGauss quadrature_formula; + FiniteElementT fe; + DoFHandler dof_handler; + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + AffineConstraints> constraints; + LinearAlgebraPETSc::MPI::SparseMatrix system_matrix; + LinearAlgebraPETSc::MPI::Vector locally_relevant_solution; + LinearAlgebraPETSc::MPI::Vector system_rhs; + }; + + + + template + NedelecSingularity::NedelecSingularity() + : mpi_communicator(MPI_COMM_WORLD) + , triangulation(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , fe_order(1) + , quadrature_formula(fe_order + 2) + , fe(fe_order) + , dof_handler(triangulation) + {} + + template + void + NedelecSingularity::setup_system() + { + dof_handler.distribute_dofs(fe); + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + locally_relevant_solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + const unsigned int first_vector_component = 0; + VectorTools::project_boundary_values_curl_conforming_l2( + dof_handler, + first_vector_component, + ZeroFunction>(dim), + 0, + constraints); + + constraints.close(); + + DynamicSparsityPattern dsp(locally_relevant_dofs); + + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern( + dsp, + dof_handler.compute_n_locally_owned_dofs_per_processor(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); + } + + + + template + void + NedelecSingularity::assemble_system() + { + FEValues fe_values(fe, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix> cell_matrix(dofs_per_cell, dofs_per_cell); + Vector> cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + const RightHandSide right_hand_side; + + std::vector>> rhs_values( + n_q_points, Vector>(dim)); + + const FEValuesExtractors::Vector electric_field(0); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit(cell); + + right_hand_side.vector_value_list(fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q = 0; q < n_q_points; ++q) + { + Tensor<1, dim, std::complex> rhs; + + for (unsigned int component = 0; component < dim; ++component) + { + // Convert vectors to tensors + rhs[component] = rhs_values[q][component]; + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const Tensor<1, dim> phi_i = + fe_values[electric_field].value(i, q); + const Tensor<1, dim> curl_phi_i = + fe_values[electric_field].curl(i, q); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + const Tensor<1, dim> phi_j = + fe_values[electric_field].value(j, q); + const Tensor<1, dim> curl_phi_j = + fe_values[electric_field].curl(j, q); + + std::complex matrix_sum = 0; + + matrix_sum += + std::pow(omega, 2) * + (-epsilon_0 * mu_0 * epsilon_r * phi_i * phi_j); + matrix_sum += omega * std::complex(0, 1) * + mu_0 * sigma * phi_i * phi_j; + matrix_sum += (1 / mu_r) * curl_phi_i * curl_phi_j; + + cell_matrix(i, j) += matrix_sum * fe_values.JxW(q); + } + + cell_rhs(i) += phi_i * rhs * fe_values.JxW(q); + } + } + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } + } + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + + + template + void + NedelecSingularity::solve() + { + LinearAlgebraPETSc::MPI::Vector completely_distributed_solution( + locally_owned_dofs, mpi_communicator); + + SolverControl solver_control; + PETScWrappers::SparseDirectMUMPS solver(solver_control, mpi_communicator); + solver.solve(system_matrix, completely_distributed_solution, system_rhs); + + constraints.distribute(completely_distributed_solution); + locally_relevant_solution = completely_distributed_solution; + } + + + + template + void + NedelecSingularity::output_results() + { + { + const ExactSolution exact_solution_function; + Vector difference_per_cell(triangulation.n_active_cells()); + + VectorTools::integrate_difference(dof_handler, + locally_relevant_solution, + exact_solution_function, + difference_per_cell, + QGauss(fe_order + 2), + VectorTools::L2_norm); + const double L2_error = + VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::L2_norm); + + deallog << " L2_error: " << L2_error << std::endl; + + // Check the continuity between between two adjacent elements. Nedelec + // enforces the continuity only on the tangencial component. Although, if + // the solution of the PDE is correct, the perpendicular component should + // also be continuous. An element boundary can be found at + // x = dimension_x/3. + const double delta = dimension_x / 1000.; + const Point point_a(dimension_x / 3. - delta, delta, delta); + const Point point_b(dimension_x / 3. + delta, delta, delta); + deallog << " Point_a = " << point_a << std::endl; + deallog << " Point_b = " << point_b << std::endl; + Vector> solution_a(dim); + Vector> solution_b(dim); + solution_a = 0; + solution_b = 0; + { + bool point_in_locally_owned_cell; + auto mapping = StaticMappingQ1::mapping; + // find the cell in which this point + // is, initialize a quadrature rule with + // it, and then a FEValues object + const std::pair::active_cell_iterator, + Point> + cell_point = GridTools::find_active_cell_around_point(mapping, + dof_handler, + point_a); + + point_in_locally_owned_cell = cell_point.first->is_locally_owned(); + if (point_in_locally_owned_cell) + { + VectorTools::point_value(dof_handler, + locally_relevant_solution, + point_a, + solution_a); + } + } + { + bool point_in_locally_owned_cell; + auto mapping = StaticMappingQ1::mapping; + // find the cell in which this point + // is, initialize a quadrature rule with + // it, and then a FEValues object + const std::pair::active_cell_iterator, + Point> + cell_point = GridTools::find_active_cell_around_point(mapping, + dof_handler, + point_b); + + point_in_locally_owned_cell = cell_point.first->is_locally_owned(); + if (point_in_locally_owned_cell) + { + VectorTools::point_value(dof_handler, + locally_relevant_solution, + point_b, + solution_b); + } + } + // Only one process has the solution_a or/and solution_b. This is a simple + // approach to send solution_a and solution_b to all the processes. + Utilities::MPI::sum(solution_a, mpi_communicator, solution_a); + Utilities::MPI::sum(solution_b, mpi_communicator, solution_b); + deallog << " Solution(point_a) : " << solution_a << std::endl; + deallog << " Solution(point_b) : " << solution_b << std::endl; + // Vector does not provide operator- + deallog << " Solution(point_b) - solution (point_a): " + << (solution_b -= solution_a) << std::endl; + } + + { + std::vector solution_names(1, "electric_field_x"); + if (dim >= 2) + { + solution_names.emplace_back("electric_field_y"); + } + if (dim == 3) + { + solution_names.emplace_back("electric_field_z"); + } + std::vector + interpretation(dim, DataComponentInterpretation::component_is_scalar); + + DataOut data_out; + data_out.add_data_vector(dof_handler, + locally_relevant_solution, + solution_names, + interpretation); + Vector subdomain(triangulation.n_active_cells()); + for (unsigned int i = 0; i < subdomain.size(); ++i) + subdomain(i) = triangulation.locally_owned_subdomain(); + data_out.add_data_vector(subdomain, "subdomain"); + + const RightHandSide rhs_function; + const ExactSolution exact_solution_function; + std::vector>> rhs( + dim, Vector>(triangulation.n_active_cells())); + std::vector>> exact_solution( + dim, Vector>(triangulation.n_active_cells())); + + // Loop over all the cells + for (const auto &cell : triangulation.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx) + { + rhs[dim_idx](cell->active_cell_index()) = + rhs_function.value(cell->center(), dim_idx); + exact_solution[dim_idx](cell->active_cell_index()) = + exact_solution_function.value(cell->center(), dim_idx); + } + } + // And on the cells that we are not interested in, set the respective + // value in the vector to a random value in order to make sure that if + // we were somehow wrong about our assumption that these elements + // would not appear in the output file, that we would find out by + // looking at the graphical output: + else + { + for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx) + { + rhs[dim_idx](cell->active_cell_index()) = -1e90; + } + } + } + + for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx) + { + data_out.add_data_vector(rhs[dim_idx], + "rhs_" + std::to_string(dim_idx)); + data_out.add_data_vector(exact_solution[dim_idx], + "exact_solution_" + std::to_string(dim_idx)); + } + + data_out.build_patches(2); + + unsigned int nb_number_positions; + if (std::is_same>::value) + { + data_out.write_vtu_in_parallel("result_nedelec.vtu", + mpi_communicator); + } + else if (std::is_same>::value) + { + data_out.write_vtu_in_parallel("result_nedelec_sz.vtu", + mpi_communicator); + } + else + { + Assert(false, ExcInternalError()); + } + } + } + + + + template + void + NedelecSingularity::run() + { + { + Point p0; + p0(0) = -dimension_x / 2; + p0(1) = -dimension_y / 2; + p0(2) = -dimension_z / 2; + Point p1; + p1(0) = dimension_x / 2; + p1(1) = dimension_y / 2; + p1(2) = dimension_z / 2; + double smallest_dimension = + std::min(dimension_z, std::min(dimension_x, dimension_y)); + std::vector divisions(dim); + divisions[0] = std::max(coarse_mesh_divisions_z, 1) * + int((p1(0) - p0(0)) / smallest_dimension); + divisions[1] = std::max(coarse_mesh_divisions_z, 1) * + int((p1(1) - p0(1)) / smallest_dimension); + divisions[2] = std::max(coarse_mesh_divisions_z, 1) * + int((p1(2) - p0(2)) / smallest_dimension); + GridGenerator::subdivided_hyper_rectangle(triangulation, + divisions, + p0, + p1); + } + + if (grid_level > 0) + { + triangulation.refine_global(grid_level); + } + + setup_system(); + deallog << " Number of active cells : " + << triangulation.n_active_cells() << std::endl; + deallog << " Number of degrees of freedom : " << dof_handler.n_dofs() + << std::endl; + + + assemble_system(); + solve(); + + output_results(); + } +} // namespace nedelec_singularity + +int +main(int argc, char *argv[]) +{ + try + { + const int dim = 3; + + dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, + argv, + 1); + + MPILogInitAll log; + + { + nedelec_singularity::NedelecSingularity> + nedelec_singularity_3d; + nedelec_singularity_3d.run(); + } + + { + nedelec_singularity::NedelecSingularity> + nedelec_singularity_3d; + nedelec_singularity_3d.run(); + } + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/fe/fe_nedelec_singularity_01.with_petsc_with_mumps=true.with_petsc_with_complex=true.mpirun=4.output b/tests/fe/fe_nedelec_singularity_01.with_petsc_with_mumps=true.with_petsc_with_complex=true.mpirun=4.output new file mode 100644 index 0000000000..5ec4400379 --- /dev/null +++ b/tests/fe/fe_nedelec_singularity_01.with_petsc_with_mumps=true.with_petsc_with_complex=true.mpirun=4.output @@ -0,0 +1,103 @@ + +DEAL:0:: Number of active cells : 209 +DEAL:0:: Number of degrees of freedom : 6084 +DEAL:0::Convergence step 1 value 0.00000 +DEAL:0:: L2_error: 0.00123995 +DEAL:0:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:0:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:0:: Solution(point_a) : (-0.449393,-8.77510e-05) (-0.0366468,-5.04085e-05) (0.0742523,-3.18423e-09) +DEAL:0:: +DEAL:0:: Solution(point_b) : (-0.475324,-9.14279e-05) (-0.0365900,-4.99325e-05) (0.0740284,-3.21944e-09) +DEAL:0:: +DEAL:0:: Solution(point_b) - solution (point_a): (-0.0259309,-3.67685e-06) (5.68532e-05,4.75991e-07) (-0.000223853,-3.52096e-11) +DEAL:0:: +DEAL:0:: Number of active cells : 209 +DEAL:0:: Number of degrees of freedom : 6084 +DEAL:0::Convergence step 1 value 0.00000 +DEAL:0:: L2_error: 0.00123995 +DEAL:0:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:0:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:0:: Solution(point_a) : (-0.449393,-8.77510e-05) (-0.0366468,-5.04085e-05) (0.0742523,-3.18440e-09) +DEAL:0:: +DEAL:0:: Solution(point_b) : (-0.475324,-9.14279e-05) (-0.0365900,-4.99325e-05) (0.0740284,-3.21962e-09) +DEAL:0:: +DEAL:0:: Solution(point_b) - solution (point_a): (-0.0259309,-3.67686e-06) (5.68532e-05,4.75991e-07) (-0.000223853,-3.52152e-11) +DEAL:0:: + +DEAL:1:: Number of active cells : 216 +DEAL:1:: Number of degrees of freedom : 6084 +DEAL:1::Convergence step 1 value 0.00000 +DEAL:1:: L2_error: 0.00123995 +DEAL:1:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:1:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:1:: Solution(point_a) : (-0.449393,-8.77510e-05) (-0.0366468,-5.04085e-05) (0.0742523,-3.18423e-09) +DEAL:1:: +DEAL:1:: Solution(point_b) : (-0.475324,-9.14279e-05) (-0.0365900,-4.99325e-05) (0.0740284,-3.21944e-09) +DEAL:1:: +DEAL:1:: Solution(point_b) - solution (point_a): (-0.0259309,-3.67685e-06) (5.68532e-05,4.75991e-07) (-0.000223853,-3.52096e-11) +DEAL:1:: +DEAL:1:: Number of active cells : 216 +DEAL:1:: Number of degrees of freedom : 6084 +DEAL:1::Convergence step 1 value 0.00000 +DEAL:1:: L2_error: 0.00123995 +DEAL:1:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:1:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:1:: Solution(point_a) : (-0.449393,-8.77510e-05) (-0.0366468,-5.04085e-05) (0.0742523,-3.18440e-09) +DEAL:1:: +DEAL:1:: Solution(point_b) : (-0.475324,-9.14279e-05) (-0.0365900,-4.99325e-05) (0.0740284,-3.21962e-09) +DEAL:1:: +DEAL:1:: Solution(point_b) - solution (point_a): (-0.0259309,-3.67686e-06) (5.68532e-05,4.75991e-07) (-0.000223853,-3.52152e-11) +DEAL:1:: + + +DEAL:2:: Number of active cells : 195 +DEAL:2:: Number of degrees of freedom : 6084 +DEAL:2::Convergence step 1 value 0.00000 +DEAL:2:: L2_error: 0.00123995 +DEAL:2:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:2:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:2:: Solution(point_a) : (-0.449393,-8.77510e-05) (-0.0366468,-5.04085e-05) (0.0742523,-3.18423e-09) +DEAL:2:: +DEAL:2:: Solution(point_b) : (-0.475324,-9.14279e-05) (-0.0365900,-4.99325e-05) (0.0740284,-3.21944e-09) +DEAL:2:: +DEAL:2:: Solution(point_b) - solution (point_a): (-0.0259309,-3.67685e-06) (5.68532e-05,4.75991e-07) (-0.000223853,-3.52096e-11) +DEAL:2:: +DEAL:2:: Number of active cells : 195 +DEAL:2:: Number of degrees of freedom : 6084 +DEAL:2::Convergence step 1 value 0.00000 +DEAL:2:: L2_error: 0.00123995 +DEAL:2:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:2:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:2:: Solution(point_a) : (-0.449393,-8.77510e-05) (-0.0366468,-5.04085e-05) (0.0742523,-3.18440e-09) +DEAL:2:: +DEAL:2:: Solution(point_b) : (-0.475324,-9.14279e-05) (-0.0365900,-4.99325e-05) (0.0740284,-3.21962e-09) +DEAL:2:: +DEAL:2:: Solution(point_b) - solution (point_a): (-0.0259309,-3.67686e-06) (5.68532e-05,4.75991e-07) (-0.000223853,-3.52152e-11) +DEAL:2:: + + +DEAL:3:: Number of active cells : 153 +DEAL:3:: Number of degrees of freedom : 6084 +DEAL:3::Convergence step 1 value 0.00000 +DEAL:3:: L2_error: 0.00123995 +DEAL:3:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:3:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:3:: Solution(point_a) : (-0.449393,-8.77510e-05) (-0.0366468,-5.04085e-05) (0.0742523,-3.18423e-09) +DEAL:3:: +DEAL:3:: Solution(point_b) : (-0.475324,-9.14279e-05) (-0.0365900,-4.99325e-05) (0.0740284,-3.21944e-09) +DEAL:3:: +DEAL:3:: Solution(point_b) - solution (point_a): (-0.0259309,-3.67685e-06) (5.68532e-05,4.75991e-07) (-0.000223853,-3.52096e-11) +DEAL:3:: +DEAL:3:: Number of active cells : 153 +DEAL:3:: Number of degrees of freedom : 6084 +DEAL:3::Convergence step 1 value 0.00000 +DEAL:3:: L2_error: 0.00123995 +DEAL:3:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:3:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:3:: Solution(point_a) : (-0.449393,-8.77510e-05) (-0.0366468,-5.04085e-05) (0.0742523,-3.18440e-09) +DEAL:3:: +DEAL:3:: Solution(point_b) : (-0.475324,-9.14279e-05) (-0.0365900,-4.99325e-05) (0.0740284,-3.21962e-09) +DEAL:3:: +DEAL:3:: Solution(point_b) - solution (point_a): (-0.0259309,-3.67686e-06) (5.68532e-05,4.75991e-07) (-0.000223853,-3.52152e-11) +DEAL:3:: + diff --git a/tests/fe/fe_nedelec_singularity_02.cc b/tests/fe/fe_nedelec_singularity_02.cc new file mode 100644 index 0000000000..d0144516a7 --- /dev/null +++ b/tests/fe/fe_nedelec_singularity_02.cc @@ -0,0 +1,784 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- +// +// By Daniel Garcia-Sanchez, CNRS +// +// Test a maxwell singularity in 3D. Maxwell singularities are common in sharp +// metallic edges such as the Fichera corner. Here we test the elements Nedelec +// and NedelecSZ using the L2 norm and the continuity of the solution. +// +// This test solves the real valued curl-curl equation in 3D: +// +// curl(curl(E))-omega^2*epsilon_0*mu_0*E = RightHandSide +// +// The manufactured solution is: +// +// E_x = (((x^2 / sqrt(x^2 + y^2)) * (x^2 - (dimension_x / 2)^2) * +// (y^2 - (dimension_y / 2)^2) * (z^2 - (dimension_z / 2)^2)) / +// ((dimension_x / 2)^3 * (dimension_y / 2)^2 * (dimension_z / 2)^2)) +// E_y = ( ((y^2 / sqrt(x^2 + y^2)) * (x^2 - (dimension_x / 2)^2) * +// (y^2 - (dimension_y / 2)^2) * (z^2 - (dimension_z / 2)^2)) / +// ((dimension_x / 2)^2 * (dimension_y / 2)^3 * (dimension_z / 2)^2)) +// E_z = 10 * (x * (x^2 - (dimension_x / 2)^2) * (y^2 - (dimension_y / 2)^2) / +// ((dimension_x / 2)^2 * (dimension_y / 2)^2)) +// +// This manufactured solution has a singularity at x = y = 0 +// +// The right hand side can be calculated with a symbolic math package such as +// sympy. + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include "../tests.h" + + +using namespace dealii; + +namespace nedelec_singularity +{ + // For the sake of simplicity define the parameters as global variables. + static const double dimension_x = 0.04; + static const double dimension_y = 0.04; + static const double dimension_z = 0.04; + static const double epsilon_0 = 8.85418782e-12; + static const double mu_0 = 1.25663706e-06; + static const double omega = 6e9 * 2 * numbers::PI; + static unsigned int nb_probe_points = 100; + static unsigned int grid_level = 1; + static unsigned int coarse_mesh_divisions_z = 3; + + + + template + class ExactSolution : public Function + { + public: + ExactSolution(); + virtual double + value(const Point &p, const unsigned int component) const override; + }; + + + + template + ExactSolution::ExactSolution() + : Function(dim) + {} + + + + template + double + ExactSolution::value(const Point & p, + const unsigned int component) const + { + const double R_x = p[0]; + const double R_y = p[1]; + const double R_z = p[2]; + + switch (component) + { + case 0: + return 2 * std::pow(R_x, 2) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) / + (std::pow(dimension_x, 3) * std::pow(dimension_y, 2) * + std::pow(dimension_z, 2) * + std::sqrt(std::pow(R_x, 2) + std::pow(R_y, 2))); + break; + case 1: + return 2 * std::pow(R_x, 2) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) / + (std::pow(dimension_x, 3) * std::pow(dimension_y, 2) * + std::pow(dimension_z, 2) * + std::sqrt(std::pow(R_x, 2) + std::pow(R_y, 2))); + break; + case 2: + return 10 * R_x * (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) / + (std::pow(dimension_x, 2) * std::pow(dimension_y, 2)); + break; + default: + Assert(false, ExcNotImplemented()); + return 0; + } + } + + + + template + class RightHandSide : public Function + { + public: + RightHandSide(); + virtual double + value(const Point &p, const unsigned int component) const override; + }; + + + + template + RightHandSide::RightHandSide() + : Function(dim) + {} + + + + template + double + RightHandSide::value(const Point & p, + const unsigned int component) const + { + const double R_x = p[0]; + const double R_y = p[1]; + const double R_z = p[2]; + + switch (component) + { + case 0: + return 2 * R_x * + (-3 * R_x * std::pow(R_y, 2) * dimension_y * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) + + R_x * dimension_y * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (-32 * std::pow(R_y, 2) - 32 * std::pow(R_z, 2) + + 8 * std::pow(dimension_y, 2) + + 8 * std::pow(dimension_z, 2) - + epsilon_0 * mu_0 * std::pow(omega, 2) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2))) + + R_x * dimension_y * (std::pow(R_x, 2) + std::pow(R_y, 2)) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (20 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) + + 3 * std::pow(R_y, 3) * dimension_x * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) + + 16 * R_y * dimension_x * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) * + (8 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) + + 2 * R_y * dimension_x * + (std::pow(R_x, 2) + std::pow(R_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) * + (std::pow(R_y, 2) * + (-16 * std::pow(R_x, 2) + 4 * std::pow(dimension_x, 2)) + + std::pow(R_y, 2) * + (-16 * std::pow(R_y, 2) + 4 * std::pow(dimension_y, 2)) - + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)))) / + (std::pow(dimension_x, 3) * std::pow(dimension_y, 3) * + std::pow(dimension_z, 2) * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 5.0 / 2.0)); + break; + case 1: + return 2 * R_y * + (3 * std::pow(R_x, 3) * dimension_y * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) - + 3 * std::pow(R_x, 2) * R_y * dimension_x * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) + + 16 * R_x * dimension_y * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) * + (8 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) + + 2 * R_x * dimension_y * + (std::pow(R_x, 2) + std::pow(R_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2)) * + (std::pow(R_x, 2) * + (-16 * std::pow(R_x, 2) + 4 * std::pow(dimension_x, 2)) + + std::pow(R_x, 2) * + (-16 * std::pow(R_y, 2) + 4 * std::pow(dimension_y, 2)) - + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2))) + + R_y * dimension_x * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 2) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (-32 * std::pow(R_x, 2) - 32 * std::pow(R_z, 2) + + 8 * std::pow(dimension_x, 2) + + 8 * std::pow(dimension_z, 2) - + epsilon_0 * mu_0 * std::pow(omega, 2) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2))) + + R_y * dimension_x * (std::pow(R_x, 2) + std::pow(R_y, 2)) * + (20 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) * + (4 * std::pow(R_z, 2) - std::pow(dimension_z, 2))) / + (std::pow(dimension_x, 3) * std::pow(dimension_y, 3) * + std::pow(dimension_z, 2) * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 5.0 / 2.0)); + break; + case 2: + return 2 * + (-8 * std::pow(R_x, 3) * R_z * dimension_y * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 5.0 / 2.0) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) + + 16 * R_x * R_z * dimension_y * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 7.0 / 2.0) * + (8 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) + + 5 * R_x * dimension_x * dimension_y * + std::pow(dimension_z, 2) * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 4) * + (-32 * std::pow(R_x, 2) - 96 * std::pow(R_y, 2) + + 8 * std::pow(dimension_x, 2) + + 24 * std::pow(dimension_y, 2) - + epsilon_0 * mu_0 * std::pow(omega, 2) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2))) - + 8 * std::pow(R_y, 3) * R_z * dimension_x * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 5.0 / 2.0) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (4 * std::pow(R_y, 2) - std::pow(dimension_y, 2)) + + 16 * R_y * R_z * dimension_x * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 7.0 / 2.0) * + (4 * std::pow(R_x, 2) - std::pow(dimension_x, 2)) * + (8 * std::pow(R_y, 2) - std::pow(dimension_y, 2))) / + (std::pow(dimension_x, 3) * std::pow(dimension_y, 3) * + std::pow(dimension_z, 2) * + std::pow(std::pow(R_x, 2) + std::pow(R_y, 2), 4)); + break; + default: + Assert(false, ExcNotImplemented()); + return 0; + } + } + + + + template + class NedelecSingularity + { + public: + NedelecSingularity(); + void + run(); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + output_results(); + + MPI_Comm mpi_communicator; + parallel::distributed::Triangulation triangulation; + const unsigned int fe_order; + const QGauss quadrature_formula; + FiniteElementT fe; + DoFHandler dof_handler; + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + AffineConstraints constraints; + LinearAlgebraPETSc::MPI::SparseMatrix system_matrix; + LinearAlgebraPETSc::MPI::Vector locally_relevant_solution; + LinearAlgebraPETSc::MPI::Vector system_rhs; + }; + + + + template + NedelecSingularity::NedelecSingularity() + : mpi_communicator(MPI_COMM_WORLD) + , triangulation(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , fe_order(1) + , quadrature_formula(fe_order + 2) + , fe(fe_order) + , dof_handler(triangulation) + {} + + template + void + NedelecSingularity::setup_system() + { + dof_handler.distribute_dofs(fe); + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + locally_relevant_solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + const unsigned int first_vector_component = 0; + VectorTools::project_boundary_values_curl_conforming_l2( + dof_handler, + first_vector_component, + ZeroFunction(dim), + 0, + constraints); + + constraints.close(); + + DynamicSparsityPattern dsp(locally_relevant_dofs); + + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern( + dsp, + dof_handler.compute_n_locally_owned_dofs_per_processor(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); + } + + + + template + void + NedelecSingularity::assemble_system() + { + FEValues fe_values(fe, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + const RightHandSide right_hand_side; + + std::vector> rhs_values(n_q_points, Vector(dim)); + + const FEValuesExtractors::Vector electric_field(0); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit(cell); + + right_hand_side.vector_value_list(fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q = 0; q < n_q_points; ++q) + { + Tensor<1, dim> rhs; + + for (unsigned int component = 0; component < dim; ++component) + { + // Convert vectors to tensors + rhs[component] = rhs_values[q][component]; + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const Tensor<1, dim> phi_i = + fe_values[electric_field].value(i, q); + const Tensor<1, dim> curl_phi_i = + fe_values[electric_field].curl(i, q); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + const Tensor<1, dim> phi_j = + fe_values[electric_field].value(j, q); + const Tensor<1, dim> curl_phi_j = + fe_values[electric_field].curl(j, q); + + double matrix_sum = 0; + + matrix_sum += std::pow(omega, 2) * + (-epsilon_0 * mu_0 * phi_i * phi_j); + matrix_sum += curl_phi_i * curl_phi_j; + + cell_matrix(i, j) += matrix_sum * fe_values.JxW(q); + } + + cell_rhs(i) += phi_i * rhs * fe_values.JxW(q); + } + } + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } + } + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + + + template + void + NedelecSingularity::solve() + { + LinearAlgebraPETSc::MPI::Vector completely_distributed_solution( + locally_owned_dofs, mpi_communicator); + + SolverControl solver_control; + PETScWrappers::SparseDirectMUMPS solver(solver_control, mpi_communicator); + solver.solve(system_matrix, completely_distributed_solution, system_rhs); + + constraints.distribute(completely_distributed_solution); + locally_relevant_solution = completely_distributed_solution; + } + + + + template + void + NedelecSingularity::output_results() + { + { + const ExactSolution exact_solution_function; + Vector difference_per_cell(triangulation.n_active_cells()); + + VectorTools::integrate_difference(dof_handler, + locally_relevant_solution, + exact_solution_function, + difference_per_cell, + QGauss(fe_order + 2), + VectorTools::L2_norm); + const double L2_error = + VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::L2_norm); + + deallog << " L2_error: " << L2_error << std::endl; + + // Check the continuity between between two adjacent elements. Nedelec + // enforces the continuity only on the tangencial component. Although, if + // the solution of the PDE is correct, the perpendicular component should + // also be continuous. An element boundary can be found at + // x = dimension_x/3. + const double delta = dimension_x / 1000.; + const Point point_a(dimension_x / 3. - delta, delta, delta); + const Point point_b(dimension_x / 3. + delta, delta, delta); + deallog << " Point_a = " << point_a << std::endl; + deallog << " Point_b = " << point_b << std::endl; + Vector solution_a(dim); + Vector solution_b(dim); + solution_a = 0; + solution_b = 0; + { + bool point_in_locally_owned_cell; + auto mapping = StaticMappingQ1::mapping; + // find the cell in which this point + // is, initialize a quadrature rule with + // it, and then a FEValues object + const std::pair::active_cell_iterator, + Point> + cell_point = GridTools::find_active_cell_around_point(mapping, + dof_handler, + point_a); + + point_in_locally_owned_cell = cell_point.first->is_locally_owned(); + if (point_in_locally_owned_cell) + { + VectorTools::point_value(dof_handler, + locally_relevant_solution, + point_a, + solution_a); + } + } + { + bool point_in_locally_owned_cell; + auto mapping = StaticMappingQ1::mapping; + // find the cell in which this point + // is, initialize a quadrature rule with + // it, and then a FEValues object + const std::pair::active_cell_iterator, + Point> + cell_point = GridTools::find_active_cell_around_point(mapping, + dof_handler, + point_b); + + point_in_locally_owned_cell = cell_point.first->is_locally_owned(); + if (point_in_locally_owned_cell) + { + VectorTools::point_value(dof_handler, + locally_relevant_solution, + point_b, + solution_b); + } + } + // Only one process has the solution_a or/and solution_b. This is a simple + // approach to send solution_a and solution_b to all the processes. + Utilities::MPI::sum(solution_a, mpi_communicator, solution_a); + Utilities::MPI::sum(solution_b, mpi_communicator, solution_b); + deallog << " Solution(point_a) : " << solution_a << std::endl; + deallog << " Solution(point_b) : " << solution_b << std::endl; + // Vector does not provide operator- + deallog << " Solution(point_b) - solution (point_a): " + << (solution_b -= solution_a) << std::endl; + } + + { + std::vector solution_names(1, "electric_field_x"); + if (dim >= 2) + { + solution_names.emplace_back("electric_field_y"); + } + if (dim == 3) + { + solution_names.emplace_back("electric_field_z"); + } + std::vector + interpretation(dim, DataComponentInterpretation::component_is_scalar); + + DataOut data_out; + data_out.add_data_vector(dof_handler, + locally_relevant_solution, + solution_names, + interpretation); + Vector subdomain(triangulation.n_active_cells()); + for (unsigned int i = 0; i < subdomain.size(); ++i) + subdomain(i) = triangulation.locally_owned_subdomain(); + data_out.add_data_vector(subdomain, "subdomain"); + + const RightHandSide rhs_function; + const ExactSolution exact_solution_function; + std::vector> rhs( + dim, Vector(triangulation.n_active_cells())); + std::vector> exact_solution( + dim, Vector(triangulation.n_active_cells())); + + // Loop over all the cells + for (const auto &cell : triangulation.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx) + { + rhs[dim_idx](cell->active_cell_index()) = + rhs_function.value(cell->center(), dim_idx); + exact_solution[dim_idx](cell->active_cell_index()) = + exact_solution_function.value(cell->center(), dim_idx); + } + } + // And on the cells that we are not interested in, set the respective + // value in the vector to a random value in order to make sure that if + // we were somehow wrong about our assumption that these elements + // would not appear in the output file, that we would find out by + // looking at the graphical output: + else + { + for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx) + { + rhs[dim_idx](cell->active_cell_index()) = -1e90; + } + } + } + + for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx) + { + data_out.add_data_vector(rhs[dim_idx], + "rhs_" + std::to_string(dim_idx)); + data_out.add_data_vector(exact_solution[dim_idx], + "exact_solution_" + std::to_string(dim_idx)); + } + + data_out.build_patches(2); + + unsigned int nb_number_positions; + if (std::is_same>::value) + { + data_out.write_vtu_in_parallel("result_nedelec.vtu", + mpi_communicator); + } + else if (std::is_same>::value) + { + data_out.write_vtu_in_parallel("result_nedelec_sz.vtu", + mpi_communicator); + } + else + { + Assert(false, ExcInternalError()); + } + } + } + + + + template + void + NedelecSingularity::run() + { + { + Point p0; + p0(0) = -dimension_x / 2; + p0(1) = -dimension_y / 2; + p0(2) = -dimension_z / 2; + Point p1; + p1(0) = dimension_x / 2; + p1(1) = dimension_y / 2; + p1(2) = dimension_z / 2; + double smallest_dimension = + std::min(dimension_z, std::min(dimension_x, dimension_y)); + std::vector divisions(dim); + divisions[0] = std::max(coarse_mesh_divisions_z, 1) * + int((p1(0) - p0(0)) / smallest_dimension); + divisions[1] = std::max(coarse_mesh_divisions_z, 1) * + int((p1(1) - p0(1)) / smallest_dimension); + divisions[2] = std::max(coarse_mesh_divisions_z, 1) * + int((p1(2) - p0(2)) / smallest_dimension); + GridGenerator::subdivided_hyper_rectangle(triangulation, + divisions, + p0, + p1); + } + + if (grid_level > 0) + { + triangulation.refine_global(grid_level); + } + + setup_system(); + deallog << " Number of active cells : " + << triangulation.n_active_cells() << std::endl; + deallog << " Number of degrees of freedom : " << dof_handler.n_dofs() + << std::endl; + + + assemble_system(); + solve(); + + output_results(); + } +} // namespace nedelec_singularity + +int +main(int argc, char *argv[]) +{ + try + { + const int dim = 3; + + dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, + argv, + 1); + + MPILogInitAll log; + + { + nedelec_singularity::NedelecSingularity> + nedelec_singularity_3d; + nedelec_singularity_3d.run(); + } + + { + nedelec_singularity::NedelecSingularity> + nedelec_singularity_3d; + nedelec_singularity_3d.run(); + } + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/fe/fe_nedelec_singularity_02.with_petsc_with_mumps=true.with_petsc_with_complex=false.mpirun=4.output b/tests/fe/fe_nedelec_singularity_02.with_petsc_with_mumps=true.with_petsc_with_complex=false.mpirun=4.output new file mode 100644 index 0000000000..9d4c834028 --- /dev/null +++ b/tests/fe/fe_nedelec_singularity_02.with_petsc_with_mumps=true.with_petsc_with_complex=false.mpirun=4.output @@ -0,0 +1,103 @@ + +DEAL:0:: Number of active cells : 209 +DEAL:0:: Number of degrees of freedom : 6084 +DEAL:0::Convergence step 1 value 0.00000 +DEAL:0:: L2_error: 0.00123995 +DEAL:0:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:0:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:0:: Solution(point_a) : -0.449393 -0.0366469 0.0742523 +DEAL:0:: +DEAL:0:: Solution(point_b) : -0.475324 -0.0365900 0.0740284 +DEAL:0:: +DEAL:0:: Solution(point_b) - solution (point_a): -0.0259309 5.68539e-05 -0.000223853 +DEAL:0:: +DEAL:0:: Number of active cells : 209 +DEAL:0:: Number of degrees of freedom : 6084 +DEAL:0::Convergence step 1 value 0.00000 +DEAL:0:: L2_error: 0.00123995 +DEAL:0:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:0:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:0:: Solution(point_a) : -0.449393 -0.0366469 0.0742523 +DEAL:0:: +DEAL:0:: Solution(point_b) : -0.475324 -0.0365900 0.0740284 +DEAL:0:: +DEAL:0:: Solution(point_b) - solution (point_a): -0.0259309 5.68539e-05 -0.000223853 +DEAL:0:: + +DEAL:1:: Number of active cells : 216 +DEAL:1:: Number of degrees of freedom : 6084 +DEAL:1::Convergence step 1 value 0.00000 +DEAL:1:: L2_error: 0.00123995 +DEAL:1:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:1:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:1:: Solution(point_a) : -0.449393 -0.0366469 0.0742523 +DEAL:1:: +DEAL:1:: Solution(point_b) : -0.475324 -0.0365900 0.0740284 +DEAL:1:: +DEAL:1:: Solution(point_b) - solution (point_a): -0.0259309 5.68539e-05 -0.000223853 +DEAL:1:: +DEAL:1:: Number of active cells : 216 +DEAL:1:: Number of degrees of freedom : 6084 +DEAL:1::Convergence step 1 value 0.00000 +DEAL:1:: L2_error: 0.00123995 +DEAL:1:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:1:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:1:: Solution(point_a) : -0.449393 -0.0366469 0.0742523 +DEAL:1:: +DEAL:1:: Solution(point_b) : -0.475324 -0.0365900 0.0740284 +DEAL:1:: +DEAL:1:: Solution(point_b) - solution (point_a): -0.0259309 5.68539e-05 -0.000223853 +DEAL:1:: + + +DEAL:2:: Number of active cells : 195 +DEAL:2:: Number of degrees of freedom : 6084 +DEAL:2::Convergence step 1 value 0.00000 +DEAL:2:: L2_error: 0.00123995 +DEAL:2:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:2:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:2:: Solution(point_a) : -0.449393 -0.0366469 0.0742523 +DEAL:2:: +DEAL:2:: Solution(point_b) : -0.475324 -0.0365900 0.0740284 +DEAL:2:: +DEAL:2:: Solution(point_b) - solution (point_a): -0.0259309 5.68539e-05 -0.000223853 +DEAL:2:: +DEAL:2:: Number of active cells : 195 +DEAL:2:: Number of degrees of freedom : 6084 +DEAL:2::Convergence step 1 value 0.00000 +DEAL:2:: L2_error: 0.00123995 +DEAL:2:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:2:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:2:: Solution(point_a) : -0.449393 -0.0366469 0.0742523 +DEAL:2:: +DEAL:2:: Solution(point_b) : -0.475324 -0.0365900 0.0740284 +DEAL:2:: +DEAL:2:: Solution(point_b) - solution (point_a): -0.0259309 5.68539e-05 -0.000223853 +DEAL:2:: + + +DEAL:3:: Number of active cells : 153 +DEAL:3:: Number of degrees of freedom : 6084 +DEAL:3::Convergence step 1 value 0.00000 +DEAL:3:: L2_error: 0.00123995 +DEAL:3:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:3:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:3:: Solution(point_a) : -0.449393 -0.0366469 0.0742523 +DEAL:3:: +DEAL:3:: Solution(point_b) : -0.475324 -0.0365900 0.0740284 +DEAL:3:: +DEAL:3:: Solution(point_b) - solution (point_a): -0.0259309 5.68539e-05 -0.000223853 +DEAL:3:: +DEAL:3:: Number of active cells : 153 +DEAL:3:: Number of degrees of freedom : 6084 +DEAL:3::Convergence step 1 value 0.00000 +DEAL:3:: L2_error: 0.00123995 +DEAL:3:: Point_a = 0.0132933 4.00000e-05 4.00000e-05 +DEAL:3:: Point_b = 0.0133733 4.00000e-05 4.00000e-05 +DEAL:3:: Solution(point_a) : -0.449393 -0.0366469 0.0742523 +DEAL:3:: +DEAL:3:: Solution(point_b) : -0.475324 -0.0365900 0.0740284 +DEAL:3:: +DEAL:3:: Solution(point_b) - solution (point_a): -0.0259309 5.68539e-05 -0.000223853 +DEAL:3:: + -- 2.39.5