From: bangerth 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) +
  • Fixed: VectorTools::project_boundary_values_curl_conforming contained + a bug for some cases. This is now fixed. +
    + (Markus Bürg, 2014/03/10) +
  • +
  • Fixed: ParameterHandler will no longer output an error if the file to be read ends with "end" without a newline.
    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 +#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; +template +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 &B) const; + void setup_system (); + void assemble_system (); + void solve (); + void process_solution(const unsigned int cycle); + Triangulation triangulation; + DoFHandler dof_handler; + FE_Nedelec fe; + ConstraintMatrix constraints; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector solution; + Vector system_rhs; + + unsigned int p_order; + unsigned int quad_order; + + ConvergenceTable convergence_table; +}; +// EXACT SOLUTION CLASS +template +class ExactSolution : public Function +{ +public: + ExactSolution() : Function() {} + virtual double value (const Point &p, + const unsigned int component) const; + virtual void vector_value (const Point &p, + Vector &result) const; + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component) const; + virtual void vector_value_list (const std::vector > &points, + std::vector > &values) const; +private: + const double PI = dealii::numbers::PI; + const double bc_constant = 0.1; +}; +// RIGHT HAND SIDE CLASS +template +class RightHandSide : public Function +{ +public: + RightHandSide (); + virtual void vector_value (const Point &p, + Vector &values) const; + virtual void vector_value_list (const std::vector > &points, + std::vector > &value_list) const; +private: + const double PI = dealii::numbers::PI; + const double bc_constant = 0.1; +}; +// DEFINE EXACT SOLUTION MEMBERS +template +double ExactSolution::value(const Point &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 +void ExactSolution::vector_value(const Point &p, + Vector &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 +void ExactSolution::value_list (const std::vector > &points, + std::vector &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 &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 +void ExactSolution::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]; + 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 +RightHandSide::RightHandSide () : +Function (dim) +{} +template +inline +void RightHandSide::vector_value (const Point &p, + Vector &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 +void RightHandSide::vector_value_list (const std::vector > &points, + std::vector > &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::vector_value (points[p], value_list[p]); + } +} +// END RIGHT HAND SIDE MEMBERS + +template +MaxwellProblem::MaxwellProblem (const unsigned int order) +: +dof_handler (triangulation), +fe (order) +{ + p_order = order; + quad_order = p_order+2; +} +template +MaxwellProblem::~MaxwellProblem () +{ + dof_handler.clear (); +} + +template +double MaxwellProblem::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 +double MaxwellProblem::dotprod(const Tensor<1,dim> &A, const Vector &B) const +{ + double return_val = 0; + for(unsigned int k = 0; k < dim; k++) { + return_val += A[k]*B(k); + } + return return_val; +} + +template +void MaxwellProblem::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(), 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 +void MaxwellProblem::assemble_system () +{ + const QGauss quadrature_formula(quad_order); + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + FEValuesViews::Vector 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 cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + RightHandSide right_hand_side; + std::vector > rhs_values (n_q_points, + Vector(dim)); + Tensor<1,dim> value_i, value_j; + + typename DoFHandler::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_pointget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + } +} +template +void MaxwellProblem::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 +void MaxwellProblem::process_solution(const unsigned int cycle) +{ + const ExactSolution exact_solution; + Vector diff_per_cell(triangulation.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(); + + 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 +void MaxwellProblem::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