From: Ross Kynch Date: Wed, 15 Apr 2015 11:39:13 +0000 (+0100) Subject: Fixed bug with 2D implementation. The tangents were being computed incorrectly on... X-Git-Tag: v8.3.0-rc1~217^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a2fd7d1f59cd51c391c1209e57930243c4624856;p=dealii.git Fixed bug with 2D implementation. The tangents were being computed incorrectly on geometries which dont align with the axis of the reference element. A way around this is to use the normal vectors instead, which come directly from the FEFaceValues object. Have included a test to demonstrate this for a simple test polynomial. --- diff --git a/doc/news/changes.h b/doc/news/changes.h index 9d05c3c281..8ac3f2417b 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -468,6 +468,13 @@ inconvenience this causes.

Specific improvements

    +
  1. Fixed: project_boundary_values_curl_conforming_l2() produced incorrect + results for non-uniform grids in 2D. Adjustment to the way 2D tangents to edges are + computed fixes this. +
    + (Ross Kynch, 2015/04/23) +
    2. +
  2. Fixed: The TimerOutput class reported abnormally large cpu time when run with more than one process with MPI. Now this is fixed.
    diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 919cd6db2e..2e51e18190 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -4161,44 +4161,12 @@ namespace VectorTools const FEValuesExtractors::Vector vec (first_vector_component); - // coordinate directions of the face. - const unsigned int - face_coordinate_direction[GeometryInfo::faces_per_cell] - = { 1, 1, 0, 0 }; - - const double tol = 0.5 * cell->face (face)->diameter () / cell->get_fe ().degree; - const std::vector > & - reference_quadrature_points = fe_face_values.get_quadrature_points (); - // Project the boundary function onto the shape functions for this edge // and set up a linear system of equations to get the values for the DoFs // associated with this edge. for (unsigned int q_point = 0; q_point < fe_face_values.n_quadrature_points; ++q_point) { - // Compute the tangent of the face - // at the quadrature point. - Point shifted_reference_point_1 - = reference_quadrature_points[q_point]; - Point shifted_reference_point_2 - = reference_quadrature_points[q_point]; - - shifted_reference_point_1 (face_coordinate_direction[face]) - += tol; - shifted_reference_point_2 (face_coordinate_direction[face]) - -= tol; - Tensor<1,dim> tangential - = (fe_face_values.get_mapping () - .transform_unit_to_real_cell (cell, - shifted_reference_point_1) - - - fe_face_values.get_mapping () - .transform_unit_to_real_cell (cell, - shifted_reference_point_2)) - / tol; - tangential - /= tangential.norm (); - // Compute the entires of the linear system // Note the system is symmetric so we could only compute the lower/upper triangle. // @@ -4207,22 +4175,37 @@ namespace VectorTools // // The RHS entries are: // \int_{edge} (tangential* boundary_value) * (tangential * edge_shape_function_i) dS. + // + // In 2D, tangential*vector is equivalent to cross_product(normal, vector), so we use this instead. + // This avoids possible issues with the computation of the tangent. + + // Store the normal at this quad point: + Point normal_at_q_point = fe_face_values.normal_vector(q_point); for (unsigned int j = 0; j < associated_edge_dofs; ++j) { const unsigned int j_face_idx = associated_edge_dof_to_face_dof[j]; + const unsigned int j_cell_idx = fe.face_to_cell_index (j_face_idx, face); + + Tensor<1,dim> phi_j = fe_face_values[vec].value (j_cell_idx, q_point); for (unsigned int i = 0; i < associated_edge_dofs; ++i) { const unsigned int i_face_idx = associated_edge_dof_to_face_dof[i]; + const unsigned int i_cell_idx = fe.face_to_cell_index (i_face_idx, face); + + Tensor<1,dim> phi_i = fe_face_values[vec].value (i_cell_idx, q_point); + + // Using n cross phi edge_matrix(i,j) += fe_face_values.JxW (q_point) - * (fe_face_values[vec].value (fe.face_to_cell_index (i_face_idx, face), q_point) * tangential) - * (fe_face_values[vec].value (fe.face_to_cell_index (j_face_idx, face), q_point) * tangential); + * ((phi_i[1]*normal_at_q_point(0) - phi_i[0]*normal_at_q_point(1)) + * (phi_j[1]*normal_at_q_point(0) - phi_j[0]*normal_at_q_point(1))); } + // Using n cross phi edge_rhs(j) += fe_face_values.JxW (q_point) - * (values[q_point] (first_vector_component) * tangential [0] - + values[q_point] (first_vector_component + 1) * tangential [1]) - * (fe_face_values[vec].value (fe.face_to_cell_index (j_face_idx, face), q_point) * tangential); + * ((values[q_point] (first_vector_component+1) * normal_at_q_point (0) + - values[q_point] (first_vector_component) * normal_at_q_point (1)) + * (phi_j[1]*normal_at_q_point(0) - phi_j[0]*normal_at_q_point(1))); } } diff --git a/tests/deal.II/nedelec_non_rect_2d.cc b/tests/deal.II/nedelec_non_rect_2d.cc new file mode 100644 index 0000000000..a95974753a --- /dev/null +++ b/tests/deal.II/nedelec_non_rect_2d.cc @@ -0,0 +1,319 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2015 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- +// +// By Ross Kynch +// +// Test to confirm that FE_Nedelec works on deformed elements in 2D when using +// Dirichlet boundary conditions (n x E = f, type). This is handled by the +// function project_boundary_values_curl_conforming_l2(). +// +// This test solves the real valued curl-curl equation in 2D: +// +// curl(curl(E)) + E = Js +// +// where the solution is: +// E = (0, x^2) +// +// so, Js = (0,-2) + E. +// +// i.e the solution should be "exact" at p=2. +// +// The domain is a distorted quad after 1 refinement. The distortion is +// performed by GridTools::distort_random. + +#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 + + +using namespace dealii; + +namespace polytest +{ + template + class SimplePolynomial : public Function + { + public: + SimplePolynomial() : Function() {} + void vector_value_list (const std::vector > &points, + std::vector > &values) const; + + void rhs_value_list (const std::vector > &points, + std::vector > &values) const; + + }; + template + void SimplePolynomial::vector_value_list (const std::vector > &points, + std::vector > &values) const + { + Assert (dim == 2, ExcNotImplemented()); + Assert (values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + for (unsigned int i=0; i &p = points[i]; + // non-zero curl-curl: + values[i][0] = 0.0; + values[i][1] = p[0]*p[0]; + + } + + } + + template + void SimplePolynomial::rhs_value_list (const std::vector > &points, + std::vector > &values) const + { + Assert (dim == 2, ExcNotImplemented()); + Assert (values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i=0; i &p = points[i]; + // non-zero curl-curl: + values[i][0] = 0.0; + values[i][1] = -2.0 + p[0]*p[0]; + } + + } + + template + class polytest + { + public: + polytest(unsigned int degree); + ~polytest (); + + void run(); + + private: + void setup_system(); + void assemble_system(); + void solve(); + + void output_error(); + + unsigned int p_order; + unsigned int quad_order; + + Triangulation tria; + DoFHandler dof_handler; + FE_Nedelec fe; + ConstraintMatrix constraints; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector solution; + Vector system_rhs; + + }; + + template + polytest::polytest (unsigned int degree) + : + p_order(degree), + quad_order(2*degree + 3), + dof_handler(tria), + fe (degree) + { + } + template + polytest::~polytest () + { + dof_handler.clear (); + } + + template + void polytest::setup_system() + { + dof_handler.distribute_dofs (fe); + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + + constraints.clear (); + SimplePolynomial boundary_function; + + VectorTools::project_boundary_values_curl_conforming_l2 (dof_handler, + 0, + boundary_function, + 0, + constraints); + constraints.close (); + DynamicSparsityPattern c_sparsity(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, + c_sparsity, + constraints,false); + + sparsity_pattern.copy_from(c_sparsity); + system_matrix.reinit (sparsity_pattern); + } + + template + void polytest::assemble_system() + { + const QGauss test_quad(quad_order); + FEValues fe_values_test (fe, test_quad, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const QGauss quadrature_formula(quad_order); + const unsigned int n_q_points = quadrature_formula.size(); + + const QGauss face_quadrature_formula(quad_order); + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + FEFaceValues fe_face_values(fe, face_quadrature_formula, + update_values | update_quadrature_points | + update_normal_vectors | update_JxW_values); + + const FEValuesExtractors::Vector vec(0); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + SimplePolynomial right_hand_side; + std::vector > rhs_value_list(n_q_points, + Vector(fe.n_components())); + + typename DoFHandler::active_cell_iterator cell, endc; + endc = dof_handler.end(); + cell = dof_handler.begin_active(); + for (; cell!=endc; ++cell) + { + fe_values_test.reinit(cell); + fe_values.reinit (cell); + cell_matrix=0; + cell_rhs=0; + + right_hand_side.rhs_value_list(fe_values.get_quadrature_points(), rhs_value_list); + for (unsigned int q=0; q rhs_value; + for (unsigned int d=0; dget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + } + } + + template + void polytest::solve() + { + SparseDirectUMFPACK direct; + direct.initialize(system_matrix); + + direct.vmult (solution, system_rhs); + constraints.distribute (solution); + } + + template + void polytest::output_error() + { + SimplePolynomial exact_solution; + Vector diff_per_cell(tria.n_active_cells()); + VectorTools::integrate_difference(dof_handler, solution, exact_solution, + diff_per_cell, QGauss(quad_order), VectorTools::L2_norm); + const double L2_error = diff_per_cell.l2_norm(); + + deallog << "p=" << p_order << " L2_error: " << L2_error << std::endl; + } + + template + void polytest::run() + { + GridGenerator::hyper_cube (tria, -1.2, 1.3); + // REFINE + tria.refine_global(1); + + // DISTORT ALL: + GridTools::distort_random (0.2, tria, false); + + setup_system(); + assemble_system(); + solve(); + output_error(); + } +} +int main () +{ + const unsigned int dim(2); + + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + + for (unsigned int p=0; p<3; ++p) + { + polytest::polytest poly(p); + poly.run(); + } +} diff --git a/tests/deal.II/nedelec_non_rect_2d.output b/tests/deal.II/nedelec_non_rect_2d.output new file mode 100644 index 0000000000..ff6c4e17f2 --- /dev/null +++ b/tests/deal.II/nedelec_non_rect_2d.output @@ -0,0 +1,4 @@ + +DEAL::p=0 L2_error: 0.512785 +DEAL::p=1 L2_error: 0.0649350 +DEAL::p=2 L2_error: 2.73094e-15