]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fixed bug with 2D implementation. The tangents were being computed incorrectly on... 801/head
authorRoss Kynch <w0ss4g3@gmail.com>
Wed, 15 Apr 2015 11:39:13 +0000 (12:39 +0100)
committerRoss Kynch <w0ss4g3@gmail.com>
Thu, 23 Apr 2015 11:02:25 +0000 (12:02 +0100)
Have included a test to demonstrate this for a simple test polynomial.

doc/news/changes.h
include/deal.II/numerics/vector_tools.templates.h
tests/deal.II/nedelec_non_rect_2d.cc [new file with mode: 0644]
tests/deal.II/nedelec_non_rect_2d.output [new file with mode: 0644]

index 9d05c3c28195530af4585140754f44408266ef5f..8ac3f2417b4f390a6af768a27224563422a7b363 100644 (file)
@@ -468,6 +468,13 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> 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.
+  <br>
+  (Ross Kynch, 2015/04/23)
+  </li>
+  
   <li> Fixed: The TimerOutput class reported abnormally large cpu time when run
   with more than one process with MPI. Now this is fixed.
   <br>
index 919cd6db2ecfca5e85dfba4454a6e3b589f58486..2e51e18190a034d017731ad18e8aedf8cdbc7014 100644 (file)
@@ -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<dim>::faces_per_cell]
-            = { 1, 1, 0, 0 };
-
-          const double tol = 0.5 * cell->face (face)->diameter () / cell->get_fe ().degree;
-          const std::vector<Point<dim> > &
-          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<dim> shifted_reference_point_1
-                = reference_quadrature_points[q_point];
-              Point<dim> 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<dim> 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 (file)
index 0000000..a959747
--- /dev/null
@@ -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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+
+using namespace dealii;
+
+namespace polytest
+{
+  template<int dim>
+  class SimplePolynomial : public Function<dim>
+  {
+  public:
+    SimplePolynomial() : Function<dim>() {}
+    void vector_value_list (const std::vector<Point<dim> > &points,
+                            std::vector<Vector<double> >   &values) const;
+                          
+    void rhs_value_list (const std::vector<Point<dim> > &points,
+                          std::vector<Vector<double> >   &values) const;
+    
+  };
+  template<int dim>
+  void SimplePolynomial<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+                                                 std::vector<Vector<double> >   &values) const
+  {
+    Assert (dim == 2, ExcNotImplemented());
+    Assert (values.size() == points.size(),
+            ExcDimensionMismatch(values.size(), points.size()));
+    for (unsigned int i=0; i<points.size(); ++i)
+    {
+      const Point<dim> &p = points[i];
+      // non-zero curl-curl:
+      values[i][0] = 0.0;
+      values[i][1] = p[0]*p[0];
+
+    }
+    
+  }
+  
+  template<int dim>
+  void SimplePolynomial<dim>::rhs_value_list (const std::vector<Point<dim> > &points,
+                                                 std::vector<Vector<double> >   &values) const
+  {
+    Assert (dim == 2, ExcNotImplemented());
+    Assert (values.size() == points.size(),
+            ExcDimensionMismatch(values.size(), points.size()));
+    
+    for (unsigned int i=0; i<points.size(); ++i)
+    {
+      const Point<dim> &p = points[i];
+      // non-zero curl-curl:
+      values[i][0] = 0.0;
+      values[i][1] = -2.0 + p[0]*p[0];
+    }
+    
+  }
+  
+  template<int dim>
+  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<dim> tria;
+    DoFHandler<dim> dof_handler;
+    FE_Nedelec<dim> fe;
+    ConstraintMatrix constraints;
+    SparsityPattern sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    Vector<double> solution;
+    Vector<double> system_rhs;
+
+  };
+  
+  template<int dim>
+  polytest<dim>::polytest (unsigned int degree)
+  :
+  p_order(degree),
+  quad_order(2*degree + 3),
+  dof_handler(tria),
+  fe (degree)
+  {    
+  }
+  template <int dim>
+  polytest<dim>::~polytest ()
+  {
+    dof_handler.clear ();
+  }
+  
+  template<int dim>
+  void polytest<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs (fe);
+    solution.reinit (dof_handler.n_dofs());
+    system_rhs.reinit (dof_handler.n_dofs());
+    
+    
+    constraints.clear ();
+    SimplePolynomial<dim> 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<int dim>
+  void polytest<dim>::assemble_system()
+  {
+    const QGauss<dim> test_quad(quad_order);
+    FEValues<dim> fe_values_test (fe, test_quad,
+                                  update_values    |  update_gradients |
+                                  update_quadrature_points  |  update_JxW_values);
+                                  
+    const QGauss<dim> quadrature_formula(quad_order);
+    const unsigned int n_q_points = quadrature_formula.size();
+    
+    const QGauss<dim-1> face_quadrature_formula(quad_order);
+    const unsigned int n_face_q_points = face_quadrature_formula.size();
+    
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_quadrature_points  |  update_JxW_values);
+    
+    FEFaceValues<dim> 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<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+    Vector<double>       cell_rhs (dofs_per_cell);
+    
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+    
+    SimplePolynomial<dim>      right_hand_side;
+    std::vector<Vector<double> > rhs_value_list(n_q_points,
+                                                Vector<double>(fe.n_components()));
+
+    typename DoFHandler<dim>::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<n_q_points; ++q)
+      {
+        Tensor<1,dim> rhs_value;
+        for (unsigned int d=0; d<dim; ++d)
+        {
+          rhs_value[d] = rhs_value_list[q](d);
+        }
+        for (unsigned int j=0; j<dofs_per_cell; ++j)
+        {
+          for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+            cell_matrix(i,j)+= (
+                                fe_values[vec].curl(i,q)*fe_values[vec].curl(j,q)
+                                + fe_values[vec].value(i,q)*fe_values[vec].value(j,q)
+                               )*fe_values.JxW(q);
+          }
+          cell_rhs(j) += rhs_value*fe_values[vec].value(j,q)*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);
+    }
+  }  
+  
+  template<int dim>
+  void polytest<dim>::solve()
+  {
+    SparseDirectUMFPACK direct;
+    direct.initialize(system_matrix);
+    
+    direct.vmult (solution, system_rhs);
+    constraints.distribute (solution);    
+  }  
+  
+  template<int dim>
+  void polytest<dim>::output_error()
+  {
+    SimplePolynomial<dim> exact_solution;
+    Vector<double> diff_per_cell(tria.n_active_cells());
+    VectorTools::integrate_difference(dof_handler, solution, exact_solution,
+                                      diff_per_cell, QGauss<dim>(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<int dim>
+  void polytest<dim>::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<dim> 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 (file)
index 0000000..ff6c4e1
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.