From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sat, 15 Mar 2014 02:13:08 +0000 (+0000)
Subject: Add a test by Ross Kynch for VectorTools::project_boundary_values_curl_conforming.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=66f898332f06a280008a123f0dfe2610599fcccc;p=dealii-svn.git

Add a test by Ross Kynch for VectorTools::project_boundary_values_curl_conforming.

git-svn-id: https://svn.dealii.org/trunk@32655 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 333194f0b3..258320b307 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -157,6 +157,12 @@ inconvenience this causes.
   (Wolfgang Bangerth, 2014/03/14)
   </li>
 
+  <li> Fixed: VectorTools::project_boundary_values_curl_conforming contained
+  a bug for some cases. This is now fixed.
+  <br>
+  (Markus B&uuml;rg, 2014/03/10)
+  </li>
+  
   <li> Fixed: ParameterHandler will no longer output an error if the file
   to be read ends with "end" without a newline.
   <br>
diff --git a/tests/deal.II/curl_curl_01.cc b/tests/deal.II/curl_curl_01.cc
new file mode 100644
index 0000000000..2718b9df20
--- /dev/null
+++ b/tests/deal.II/curl_curl_01.cc
@@ -0,0 +1,401 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2010 - 2013 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+/*
+ Test by Ross Kynch:
+  
+ Solves curl(curl(u)) + u = f in 2D
+ 
+ exact solution is given by:
+ u(0) = cos(pi*x)*sin(pi*y) + C
+ u(1) = -sin(pi*x)*cos(pi(y) + C
+ 
+ f(0) = (2*pi^2 + 1)*cos(pi*x)*sin(pi*x) + C
+ f(1) = -(2*pi^2 + 1)*sin(pi*x)*cos(pi(y) + C
+ 
+ where C is some constant. (To change C, edit bc_constant in ExactSolution and
+ RightHandSide classes.)
+
+ The test solves these equations using the Nedelec element of higher order and
+ computes the error.
+ 
+ */
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/convergence_table.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/compressed_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <fstream>
+#include <iostream>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/grid/grid_refinement.h>
+
+using namespace dealii;
+template <int dim>
+class MaxwellProblem
+{
+public:
+    MaxwellProblem (const unsigned int order);
+    ~MaxwellProblem ();
+    void run ();
+private:
+    double dotprod(const Tensor<1,dim> &A, const Tensor<1,dim> &B) const;
+    double dotprod(const Tensor<1,dim> &A, const Vector<double> &B) const;
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void process_solution(const unsigned int cycle);
+    Triangulation<dim>   triangulation;
+    DoFHandler<dim>      dof_handler;
+    FE_Nedelec<dim>            fe;
+    ConstraintMatrix     constraints;
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+    
+    unsigned int p_order;
+    unsigned int quad_order;
+    
+    ConvergenceTable	   convergence_table;
+};
+// EXACT SOLUTION CLASS
+template<int dim>
+class ExactSolution : public Function<dim>
+{
+public:
+    ExactSolution() : Function<dim>() {}
+    virtual double value (const Point<dim> &p,
+                          const unsigned int component) const;
+    virtual void vector_value (const Point<dim> &p,
+                               Vector<double> &result) const;
+    virtual void value_list (const std::vector<Point<dim> > &points,
+                             std::vector<double> &values,
+                             const unsigned int component) const;
+    virtual void vector_value_list (const std::vector<Point<dim> > &points,
+                                    std::vector<Vector<double> >   &values) const;
+private:
+    const double PI = dealii::numbers::PI;
+	const double bc_constant = 0.1;
+};
+// RIGHT HAND SIDE CLASS
+template <int dim>
+class RightHandSide :  public Function<dim>
+{
+public:
+    RightHandSide ();
+    virtual void vector_value (const Point<dim> &p,
+                               Vector<double>   &values) const;
+    virtual void vector_value_list (const std::vector<Point<dim> > &points,
+                                    std::vector<Vector<double> >   &value_list) const;
+private:
+    const double PI = dealii::numbers::PI;
+	const double bc_constant = 0.1;
+};
+// DEFINE EXACT SOLUTION MEMBERS
+template<int dim>
+double ExactSolution<dim>::value(const Point<dim> &p,
+                                 const unsigned int component) const
+{
+    Assert (dim >= 2, ExcNotImplemented());
+    AssertIndexRange(component, dim);
+    
+    double val = -1000;
+    switch(component) {
+        case 0:	val = cos(PI*p(0))*sin(PI*p(1)) + bc_constant;
+        case 1:	val = -sin(PI*p(0))*cos(PI*p(1)) + bc_constant;
+    }
+    return val;
+    
+}
+template<int dim>
+void ExactSolution<dim>::vector_value(const Point<dim> &p,
+                                      Vector<double> &result) const
+{
+    Assert(dim >= 2, ExcNotImplemented());
+    result(0) = cos(PI*p(0))*sin(PI*p(1)) + bc_constant;
+    result(1) = -sin(PI*p(0))*cos(PI*p(1)) + bc_constant;
+    
+}
+template <int dim>
+void ExactSolution<dim>::value_list (const std::vector<Point<dim> > &points,
+                                     std::vector<double> &values,
+                                     const unsigned int component) const
+{
+    Assert (values.size() == points.size(), ExcDimensionMismatch(values.size(), points.size()));
+    AssertIndexRange(component, dim);
+    for (unsigned int i=0; i<points.size(); ++i)
+    {
+        const Point<dim> &p = points[i];
+        switch(component)
+        {
+            case 0:
+                values[i] = cos(PI*p(0))*sin(PI*p(1)) + bc_constant;
+            case 1:
+                values[i] = -sin(PI*p(0))*cos(PI*p(1)) + bc_constant;
+        }
+    }
+}
+template <int dim>
+void ExactSolution<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];
+        values[i](0) = cos(PI*p(0))*sin(PI*p(1)) + bc_constant;
+        values[i](1) = -sin(PI*p(0))*cos(PI*p(1)) + bc_constant;
+    }
+}
+// END EXACT SOLUTION MEMBERS
+
+// DEFINE RIGHT HAND SIDE MEMBERS
+template <int dim>
+RightHandSide<dim>::RightHandSide () :
+Function<dim> (dim)
+{}
+template <int dim>
+inline
+void RightHandSide<dim>::vector_value (const Point<dim> &p,
+                                       Vector<double>   &values) const
+{
+    Assert (values.size() == dim, ExcDimensionMismatch (values.size(), dim));
+    Assert (dim >= 2, ExcNotImplemented());
+    
+    //2D solution
+    values(0) = (2*PI*PI + 1)*cos(PI*p(0))*sin(PI*p(1)) + bc_constant;
+    values(1) = -(2*PI*PI + 1)*sin(PI*p(0))*cos(PI*p(1)) + bc_constant;
+}
+template <int dim>
+void RightHandSide<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+                                            std::vector<Vector<double> >   &value_list) const
+{
+    Assert (value_list.size() == points.size(), ExcDimensionMismatch (value_list.size(), points.size()));
+    const unsigned int n_points = points.size();
+    for (unsigned int p=0; p<n_points; ++p)
+    {
+        RightHandSide<dim>::vector_value (points[p], value_list[p]);
+    }
+}
+// END RIGHT HAND SIDE MEMBERS
+
+template <int dim>
+MaxwellProblem<dim>::MaxwellProblem (const unsigned int order)
+:
+dof_handler (triangulation),
+fe (order)
+{
+    p_order = order;
+    quad_order = p_order+2;
+}
+template <int dim>
+MaxwellProblem<dim>::~MaxwellProblem ()
+{
+    dof_handler.clear ();
+}
+
+template<int dim>
+double MaxwellProblem<dim>::dotprod(const Tensor<1,dim> &A, const Tensor<1,dim> &B) const
+{
+    double return_val = 0;
+    for(unsigned int k = 0; k < dim; k++) {
+        return_val += A[k]*B[k];
+    }
+    return return_val;
+}
+
+template<int dim>
+double MaxwellProblem<dim>::dotprod(const Tensor<1,dim> &A, const Vector<double> &B) const
+{
+    double return_val = 0;
+    for(unsigned int k = 0; k < dim; k++) {
+        return_val += A[k]*B(k);
+    }
+    return return_val;
+}
+
+template <int dim>
+void MaxwellProblem<dim>::setup_system ()
+{
+    dof_handler.distribute_dofs (fe);
+    solution.reinit (dof_handler.n_dofs());
+    system_rhs.reinit (dof_handler.n_dofs());
+    constraints.clear ();
+    DoFTools::make_hanging_node_constraints (dof_handler,
+                                             constraints);
+    // FE_Nedelec boundary condition.
+    VectorTools::project_boundary_values_curl_conforming(dof_handler, 0, ExactSolution<dim>(), 0, constraints);
+    
+    constraints.close ();
+    CompressedSparsityPattern 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 MaxwellProblem<dim>::assemble_system ()
+{
+    const QGauss<dim>  quadrature_formula(quad_order);
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_quadrature_points  |  update_JxW_values);
+    FEValuesViews::Vector<dim> fe_views(fe_values, 0);
+    const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int   n_q_points    = quadrature_formula.size();
+    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);
+    
+    RightHandSide<dim>      right_hand_side;
+    std::vector<Vector<double> > rhs_values (n_q_points,
+                                             Vector<double>(dim));
+    Tensor<1,dim> value_i, value_j;
+    
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+    {
+        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_point=0; q_point<n_q_points; ++q_point)
+        {
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+            {
+                value_i[0] = fe_values.shape_value_component(i,q_point,0);
+                value_i[1] = fe_values.shape_value_component(i,q_point,1);
+                if (dim == 3) {
+                    value_i[2] = fe_values.shape_value_component(i,q_point,2);
+                }
+                for (unsigned int j=0; j<dofs_per_cell; ++j)
+                {
+                    value_j[0] = fe_values.shape_value_component(j,q_point,0);
+                    value_j[1] = fe_values.shape_value_component(j,q_point,1);
+                    if (dim == 3) {
+                        value_j[2] = fe_values.shape_value_component(j,q_point,2);
+                    }
+                    cell_matrix(i,j) += ( fe_views.curl(i,q_point)[0]*fe_views.curl(j,q_point)[0]
+                                         + dotprod(value_i,value_j) )*fe_values.JxW(q_point);
+                    
+                    
+                }
+                cell_rhs(i) += dotprod(value_i,rhs_values[q_point])*fe_values.JxW(q_point);
+            }
+        }
+        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 MaxwellProblem<dim>::solve ()
+{
+  /* CG:
+    SolverControl           solver_control (1000, 1e-8);
+    SolverCG<>              solver (solver_control);
+    PreconditionSSOR<> preconditioner;
+    preconditioner.initialize(system_matrix, 1.2);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+    constraints.distribute (solution);
+    */
+    
+    SparseDirectUMFPACK A_direct;
+    A_direct.initialize(system_matrix);
+  
+    A_direct.vmult (solution, system_rhs);
+    constraints.distribute (solution);
+    
+}
+template<int dim>
+void MaxwellProblem<dim>::process_solution(const unsigned int cycle)
+{
+    const ExactSolution<dim> exact_solution;
+    Vector<double> diff_per_cell(triangulation.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();
+    
+    convergence_table.add_value("cycle", cycle);
+    convergence_table.add_value("cells", triangulation.n_active_cells());
+    convergence_table.add_value("dofs", dof_handler.n_dofs());
+    convergence_table.add_value("L2 Error", L2_error);
+}
+
+template <int dim>
+void MaxwellProblem<dim>::run ()
+{
+    for (unsigned int cycle=0; cycle<3; ++cycle)
+    {
+        if (cycle == 0)
+        {
+            GridGenerator::hyper_cube (triangulation, -1, 1);
+            triangulation.refine_global (2);
+        }
+        else
+            triangulation.refine_global (1);
+        setup_system ();
+        assemble_system ();
+        solve ();
+        process_solution (cycle);
+    }
+    convergence_table.set_precision("L2 Error",8);
+    convergence_table.set_scientific("L2 Error",true);
+    convergence_table.write_text(deallog.get_file_stream());
+}
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  MaxwellProblem<2> maxwell(4);
+  maxwell.run ();
+}
diff --git a/tests/deal.II/curl_curl_01.output b/tests/deal.II/curl_curl_01.output
new file mode 100644
index 0000000000..1894cab461
--- /dev/null
+++ b/tests/deal.II/curl_curl_01.output
@@ -0,0 +1,5 @@
+
+cycle cells dofs     L2 Error    
+0     16    840   1.32580012e-04 
+1     64    3280  4.19612400e-06 
+2     256   12960 1.31545840e-07