DEAL::
DEAL::
-DEAL:: Number of active cells: 768
-DEAL:: Number of degrees of freedom: 3264
+DEAL:: Number of active cells: 192
+DEAL:: Number of degrees of freedom: 864
DEAL:: Number of constraints : 0
-DEAL::Reference matrix nonzeros: 49920, actually: 49766
-DEAL::Test matrix 1 nonzeros: 49920, actually: 42514
-DEAL::Matrix difference norm: 0
-DEAL::rhs difference norm: 0
-DEAL::Test matrix 2 nonzeros: 49920, actually: 42514
-DEAL::Matrix difference norm: 0
-DEAL::rhs difference norm: 0
-DEAL:cg::Starting value 28.
-DEAL:cg::Convergence step 48 value 3.5e-07
+DEAL:: Reference matrix nonzeros: 12672, actually: 8974
+DEAL:: Test matrix 1 nonzeros: 12672, actually: 8974
+DEAL:: Matrix difference norm: 0
+DEAL:: RHS difference norm: 0
+DEAL:: Test matrix 2 nonzeros: 12672, actually: 8974
+DEAL:: Matrix difference norm: 0
+DEAL:: RHS difference norm: 0
+DEAL:cg::Starting value 20.
+DEAL:cg::Convergence step 25 value 2.5e-07
DEAL::Distribute error: 0
DEAL::
DEAL::
-DEAL:: Number of active cells: 1173
-DEAL:: Number of degrees of freedom: 5732
-DEAL:: Number of constraints : 492
-DEAL::Reference matrix nonzeros: 98166, actually: 87212
-DEAL::Test matrix 1 nonzeros: 98166, actually: 78146
-DEAL::Matrix difference norm: 0
-DEAL::rhs difference norm: 0
-DEAL::Test matrix 2 nonzeros: 98166, actually: 78146
-DEAL::Matrix difference norm: 0
-DEAL::rhs difference norm: 0
-DEAL:cg::Starting value 31.
-DEAL:cg::Convergence step 84 value 4.4e-07
+DEAL:: Number of active cells: 309
+DEAL:: Number of degrees of freedom: 1575
+DEAL:: Number of constraints : 187
+DEAL:: Reference matrix nonzeros: 25901, actually: 17333
+DEAL:: Test matrix 1 nonzeros: 25901, actually: 17333
+DEAL:: Matrix difference norm: 0
+DEAL:: RHS difference norm: 0
+DEAL:: Test matrix 2 nonzeros: 25901, actually: 17333
+DEAL:: Matrix difference norm: 0
+DEAL:: RHS difference norm: 0
+DEAL:cg::Starting value 22.
+DEAL:cg::Convergence step 42 value 3.2e-07
DEAL::Distribute error: 0
DEAL::
DEAL::
-DEAL:: Number of active cells: 1644
-DEAL:: Number of degrees of freedom: 9764
-DEAL:: Number of constraints : 1434
-DEAL::Reference matrix nonzeros: 200326, actually: 161794
-DEAL::Test matrix 1 nonzeros: 200326, actually: 150182
-DEAL::Matrix difference norm: 0
-DEAL::rhs difference norm: 0
-DEAL::Test matrix 2 nonzeros: 200326, actually: 150182
-DEAL::Matrix difference norm: 0
-DEAL::rhs difference norm: 0
-DEAL:cg::Starting value 35.
-DEAL:cg::Convergence step 107 value 5.2e-07
+DEAL:: Number of active cells: 450
+DEAL:: Number of degrees of freedom: 2692
+DEAL:: Number of constraints : 475
+DEAL:: Reference matrix nonzeros: 51910, actually: 33908
+DEAL:: Test matrix 1 nonzeros: 51910, actually: 33908
+DEAL:: Matrix difference norm: 0
+DEAL:: RHS difference norm: 0
+DEAL:: Test matrix 2 nonzeros: 51910, actually: 33908
+DEAL:: Matrix difference norm: 0
+DEAL:: RHS difference norm: 0
+DEAL:cg::Starting value 25.
+DEAL:cg::Convergence step 61 value 3.8e-07
DEAL::Distribute error: 0
--- /dev/null
+//---------------- inhomogeneous_constraints_nonsymmetric.cc -------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2009 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------- inhomogeneous_constraints_nonsymmetric.cc -------------------
+
+
+// this function tests the correctness of the implementation of
+// inhomogeneous constraints on a nonsymmetric matrix that comes from a
+// discretization of the first derivative
+
+#include "../tests.h"
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <base/utilities.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_constraints.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/error_estimator.h>
+#include <lac/compressed_simple_sparsity_pattern.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+std::ofstream logfile("inhomogeneous_constraints_nonsymmetric/output");
+
+using namespace dealii;
+
+template <int dim>
+class AdvectionProblem
+{
+ public:
+ AdvectionProblem ();
+ ~AdvectionProblem ();
+
+ void run ();
+
+ private:
+ void setup_system ();
+ void test_equality ();
+ void assemble_reference ();
+ void assemble_test_1 ();
+ void assemble_test_2 ();
+
+ Triangulation<dim> triangulation;
+
+ DoFHandler<dim> dof_handler;
+ FE_Q<dim> fe;
+
+ ConstraintMatrix hanging_nodes_only;
+ ConstraintMatrix test_all_constraints;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> reference_matrix;
+ SparseMatrix<double> test_matrix;
+
+ Vector<double> reference_rhs;
+ Vector<double> test_rhs;
+};
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+ public:
+ RightHandSide () : Function<dim> () {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component) const;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ double product = 1;
+ for (unsigned int d=0; d<dim; ++d)
+ product *= (p[d]+1);
+ return product;
+}
+
+
+template <int dim>
+AdvectionProblem<dim>::AdvectionProblem ()
+ :
+ dof_handler (triangulation),
+ fe (2)
+{}
+
+
+template <int dim>
+AdvectionProblem<dim>::~AdvectionProblem ()
+{
+ dof_handler.clear ();
+}
+
+
+template <int dim>
+void AdvectionProblem<dim>::setup_system ()
+{
+ dof_handler.distribute_dofs (fe);
+
+ reference_rhs.reinit (dof_handler.n_dofs());
+ test_rhs.reinit (dof_handler.n_dofs());
+
+ hanging_nodes_only.clear ();
+ test_all_constraints.clear ();
+
+ // add boundary conditions as
+ // inhomogeneous constraints here. In
+ // contrast to step-27, we choose a
+ // constant function with value 1 here.
+ {
+ std::map<unsigned int,double> boundary_values;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ RightHandSide<dim>(),
+ boundary_values);
+ std::map<unsigned int,double>::const_iterator boundary_value =
+ boundary_values.begin();
+ for ( ; boundary_value !=boundary_values.end(); ++boundary_value)
+ {
+ test_all_constraints.add_line(boundary_value->first);
+ test_all_constraints.set_inhomogeneity (boundary_value->first,
+ boundary_value->second);
+ }
+ }
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ hanging_nodes_only);
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ test_all_constraints);
+ hanging_nodes_only.close ();
+ test_all_constraints.close ();
+
+ CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
+ dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, csp,
+ hanging_nodes_only, true);
+ sparsity_pattern.copy_from (csp);
+
+ reference_matrix.reinit (sparsity_pattern);
+ test_matrix.reinit (sparsity_pattern);
+}
+
+
+
+ // test whether we are equal with the
+ // standard matrix and right hand side
+template <int dim>
+void AdvectionProblem<dim>::test_equality ()
+{
+ // need to manually go through the
+ // matrix, since we can have different
+ // entries in constrained lines.
+ for (unsigned int i=0; i<reference_matrix.m(); ++i)
+ {
+ SparseMatrix<double>::const_iterator reference = reference_matrix.begin(i);
+ SparseMatrix<double>::iterator test = test_matrix.begin(i);
+ if (test_all_constraints.is_constrained(i) == false)
+ {
+ for ( ; test != test_matrix.end(i); ++test, ++reference)
+ test->value() -= reference->value();
+ }
+ else
+ for ( ; test != test_matrix.end(i); ++test)
+ test->value() = 0;
+ }
+
+ deallog << " Matrix difference norm: "
+ << test_matrix.frobenius_norm() << std::endl;
+ Assert (test_matrix.frobenius_norm() < 1e-13, ExcInternalError());
+
+ // same here -- Dirichlet lines will have
+ // nonzero rhs, whereas we will have zero
+ // rhs when using inhomogeneous
+ // constraints.
+ for (unsigned int i=0; i<reference_matrix.m(); ++i)
+ if (test_all_constraints.is_constrained(i) == false)
+ test_rhs(i) -= reference_rhs(i);
+ else
+ test_rhs(i) = 0;
+
+ deallog << " RHS difference norm: "
+ << test_rhs.l2_norm() << std::endl;
+
+ Assert (test_rhs.l2_norm() < 1e-14, ExcInternalError());
+}
+
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::assemble_reference ()
+{
+ reference_matrix = 0;
+ reference_rhs = 0;
+
+ QGauss<dim> quadrature_formula (3);
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ const RightHandSide<dim> rhs_function;
+ 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<unsigned int> local_dof_indices (dofs_per_cell);
+ std::vector<double> rhs_values (n_q_points);
+
+ 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);
+
+ rhs_function.value_list (fe_values.get_quadrature_points(),
+ rhs_values);
+
+ Tensor<1,dim> advection_direction;
+ advection_direction[0] = 1;
+ advection_direction[1] = 1;
+ advection_direction[dim-1] = -1;
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (fe_values.shape_value(i,q_point) *
+ advection_direction *
+ fe_values.shape_grad(j,q_point) *
+ fe_values.JxW(q_point));
+
+ cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+ rhs_values[q_point] *
+ fe_values.JxW(q_point));
+ }
+
+ local_dof_indices.resize (dofs_per_cell);
+ cell->get_dof_indices (local_dof_indices);
+
+ reference_matrix.add(local_dof_indices, cell_matrix);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ reference_rhs(local_dof_indices[i]) += cell_rhs(i);
+ }
+
+ hanging_nodes_only.condense (reference_matrix, reference_rhs);
+
+ // use some other rhs vector as dummy for
+ // application of Dirichlet conditions
+ std::map<unsigned int,double> boundary_values;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ RightHandSide<dim>(),
+ boundary_values);
+ MatrixTools::apply_boundary_values (boundary_values,
+ reference_matrix,
+ test_rhs,
+ reference_rhs);
+
+ deallog << " Reference matrix nonzeros: " << reference_matrix.n_nonzero_elements()
+ << ", actually: " << reference_matrix.n_actually_nonzero_elements ()
+ << std::endl;
+}
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::assemble_test_1 ()
+{
+ test_matrix = 0;
+ test_rhs = 0;
+
+
+ QGauss<dim> quadrature_formula (3);
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ const RightHandSide<dim> rhs_function;
+ 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<unsigned int> local_dof_indices (dofs_per_cell);
+ std::vector<double> rhs_values (n_q_points);
+
+ 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);
+
+ rhs_function.value_list (fe_values.get_quadrature_points(),
+ rhs_values);
+
+ Tensor<1,dim> advection_direction;
+ advection_direction[0] = 1;
+ advection_direction[1] = 1;
+ advection_direction[dim-1] = -1;
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (fe_values.shape_value(i,q_point) *
+ advection_direction *
+ fe_values.shape_grad(j,q_point) *
+ fe_values.JxW(q_point));
+
+ cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+ rhs_values[q_point] *
+ fe_values.JxW(q_point));
+ }
+
+ local_dof_indices.resize (dofs_per_cell);
+ cell->get_dof_indices (local_dof_indices);
+
+ test_matrix.add(local_dof_indices, cell_matrix);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ test_rhs(local_dof_indices[i]) += cell_rhs(i);
+
+ }
+
+ test_all_constraints.condense (test_matrix, test_rhs);
+ deallog << " Test matrix 1 nonzeros: " << test_matrix.n_nonzero_elements()
+ << ", actually: " << test_matrix.n_actually_nonzero_elements ()
+ << std::endl;
+
+ test_equality();
+}
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::assemble_test_2 ()
+{
+ test_matrix = 0;
+ test_rhs = 0;
+
+ QGauss<dim> quadrature_formula (3);
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ const RightHandSide<dim> rhs_function;
+ 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<unsigned int> local_dof_indices (dofs_per_cell);
+ std::vector<double> rhs_values (n_q_points);
+
+ 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);
+
+ rhs_function.value_list (fe_values.get_quadrature_points(),
+ rhs_values);
+
+ Tensor<1,dim> advection_direction;
+ advection_direction[0] = 1;
+ advection_direction[1] = 1;
+ advection_direction[dim-1] = -1;
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (fe_values.shape_value(i,q_point) *
+ advection_direction *
+ fe_values.shape_grad(j,q_point) *
+ fe_values.JxW(q_point));
+
+ cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+ rhs_values[q_point] *
+ fe_values.JxW(q_point));
+ }
+
+ local_dof_indices.resize (dofs_per_cell);
+ cell->get_dof_indices (local_dof_indices);
+
+ test_all_constraints.distribute_local_to_global (cell_matrix,
+ cell_rhs,
+ local_dof_indices,
+ test_matrix,
+ test_rhs);
+ }
+ deallog << " Test matrix 2 nonzeros: " << test_matrix.n_nonzero_elements()
+ << ", actually: " << test_matrix.n_actually_nonzero_elements ()
+ << std::endl;
+ test_equality();
+}
+
+
+template <int dim>
+void AdvectionProblem<dim>::run ()
+{
+ GridGenerator::hyper_cube (triangulation);
+ triangulation.refine_global (2);
+
+ // manually refine the first two cells to
+ // create some hanging nodes
+ {
+ typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active();
+ cell->set_refine_flag();
+ cell++;
+ cell->set_refine_flag();
+ }
+ triangulation.execute_coarsening_and_refinement();
+
+ setup_system ();
+
+ deallog << std::endl << std::endl
+ << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl
+ << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl
+ << " Number of constraints : "
+ << hanging_nodes_only.n_constraints()
+ << std::endl;
+
+ assemble_reference ();
+ assemble_test_1 ();
+ assemble_test_2 ();
+}
+
+
+
+int main ()
+{
+ deallog << std::setprecision (2);
+ logfile << std::setprecision (2);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-12);
+
+ {
+ AdvectionProblem<2> advection_problem;
+ advection_problem.run ();
+ }
+ {
+ AdvectionProblem<3> advection_problem;
+ advection_problem.run ();
+ }
+}