From: Daniel Garcia-Sanchez Date: Thu, 23 Jul 2020 18:00:12 +0000 (+0200) Subject: Remove unused tests X-Git-Tag: v9.3.0-rc1~1236^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0da4494dc2a7134f71430693eaae0475dd26a5be;p=dealii.git Remove unused tests --- diff --git a/tests/fe/fe_nedelec_singularity_01.cc b/tests/fe/fe_nedelec_singularity_01.cc deleted file mode 100644 index 61bc4cf1af..0000000000 --- a/tests/fe/fe_nedelec_singularity_01.cc +++ /dev/null @@ -1,796 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2019 - 2020 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" - - - -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 std::complex - value(const Point &p, const unsigned int component) const override; - }; - - - - template - ExactSolution::ExactSolution() - : Function>(dim) - {} - - - - template - std::complex - 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, - Functions::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, - locally_owned_dofs, - 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 (const auto q : fe_values.quadrature_point_indices()) - { - 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 deleted file mode 100644 index 5ec4400379..0000000000 --- a/tests/fe/fe_nedelec_singularity_01.with_petsc_with_mumps=true.with_petsc_with_complex=true.mpirun=4.output +++ /dev/null @@ -1,103 +0,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.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 deleted file mode 100644 index 7540704502..0000000000 --- a/tests/fe/fe_nedelec_singularity_02.cc +++ /dev/null @@ -1,790 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2019 - 2020 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" - - - -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, - Functions::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, - locally_owned_dofs, - 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 (const auto q : fe_values.quadrature_point_indices()) - { - 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) : "; - for (const auto v : solution_a) - deallog << v << ' '; - deallog << std::endl << std::endl; - deallog << " Solution(point_b) : "; - for (const auto v : solution_b) - deallog << v << ' '; - deallog << std::endl << std::endl; - // Vector does not provide operator- - deallog << " Solution(point_b) - solution (point_a): "; - for (const auto v : (solution_b -= solution_a)) - deallog << v << ' '; - deallog << std::endl << 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 deleted file mode 100644 index 9d4c834028..0000000000 --- a/tests/fe/fe_nedelec_singularity_02.with_petsc_with_mumps=true.with_petsc_with_complex=false.mpirun=4.output +++ /dev/null @@ -1,103 +0,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: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:: -