]> https://gitweb.dealii.org/ - dealii.git/commitdiff
test assembler_simple_system_inhom_01
authorJoscha Gedicke <joscha.gedicke@iwr.uni-heidelberg.de>
Wed, 9 Mar 2016 17:31:01 +0000 (18:31 +0100)
committerJoscha Gedicke <joscha.gedicke@iwr.uni-heidelberg.de>
Wed, 9 Mar 2016 19:13:54 +0000 (20:13 +0100)
tests/integrators/assembler_simple_system_inhom_01.cc [new file with mode: 0644]
tests/integrators/assembler_simple_system_inhom_01.output [new file with mode: 0644]

diff --git a/tests/integrators/assembler_simple_system_inhom_01.cc b/tests/integrators/assembler_simple_system_inhom_01.cc
new file mode 100644 (file)
index 0000000..60373cc
--- /dev/null
@@ -0,0 +1,691 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2012 - 2016 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 SystemSimple with mixed homogeneous and inhomogeneous constraints.
+ * Compare SystemSimple with matrix and rhs computed from
+ * MatrixSimple, ResidualSimple with homogeneous constraints
+ * and global constraints.condense(matrix,rhs) operation to
+ * incorporate the inhomogeneous constraints.
+ * Tests are run for
+ * -> FE_Q<2>(1), FE_DGQ<2>(1)
+ * -> mesh: adaptive with hanging nodes
+ * -> homogeneous hanging nodes constraints
+ * -> inhomogeneous boundary constraints plus constraint of all degrees
+ *    at the origin
+ * -> matrix is the (nonsymmetric) advection matrix plus the jump of the
+ *    normal flux times the average of the values along inner edges.
+ * -> vector is standard rhs plus neumann boundary integral.
+ *
+ * Output: norm of the difference of the matrices and vectors
+ */
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.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/compressed_sparsity_pattern.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/dofs/dof_accessor.h>
+
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/meshworker/dof_info.h>
+#include <deal.II/meshworker/integration_info.h>
+#include <deal.II/meshworker/simple.h>
+#include <deal.II/meshworker/loop.h>
+
+#include <iostream>
+#include <fstream>
+
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <deal.II/base/tensor_function.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/numerics/data_out.h>
+
+
+using namespace dealii;
+
+
+/*
+ * Right hand side
+ */
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide () : Function<dim>() {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const;
+};
+
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim> &p,
+                                  const unsigned int component) const
+{
+  Assert (component == 0, ExcNotImplemented());
+
+  double val = 1; // f = 1
+  return val;
+}
+
+
+/*
+ * Boundary values
+ */
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues () : Function<dim>() {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const;
+};
+
+template <int dim>
+double BoundaryValues<dim>::value (const Point<dim> &p,
+                                   const unsigned int component) const
+{
+  Assert (component == 0, ExcNotImplemented());
+
+  double val = 1; // u_D = 1
+  return val;
+}
+
+
+/*
+ * Flux boundary values
+ */
+template <int dim>
+class FluxBoundaryValues : public Function<dim>
+{
+public:
+  FluxBoundaryValues () : Function<dim>() {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const;
+};
+
+template <int dim>
+double FluxBoundaryValues<dim>::value (const Point<dim> &p,
+                                       const unsigned int component) const
+{
+  double val = 1; // g = 1
+  return val;
+}
+
+
+/*
+ * Advection coefficient
+ */
+template <int dim>
+class Advection : public TensorFunction<1,dim>
+{
+public:
+  Advection () : TensorFunction<1,dim> () {}
+
+  virtual Tensor<1,dim> value (const Point<dim> &p) const;
+};
+
+template <int dim>
+Tensor<1,dim>
+Advection<dim>::value (const Point<dim> &p) const
+{
+  Point<dim> value; // beta = (1,0)^t
+  value[0] = 1;
+  value[1] = 0;
+
+  return value;
+}
+
+
+/*
+ * System integrator
+ */
+template <int dim>
+class SystemIntegrator : public MeshWorker::LocalIntegrator<dim>
+{
+public:
+  void cell(MeshWorker::DoFInfo<dim> &dinfo,
+            typename MeshWorker::IntegrationInfo<dim> &info) const;
+  void boundary(MeshWorker::DoFInfo<dim> &dinfo,
+                typename MeshWorker::IntegrationInfo<dim> &info) const;
+  void face(MeshWorker::DoFInfo<dim> &dinfo1,
+            MeshWorker::DoFInfo<dim> &dinfo2,
+            typename MeshWorker::IntegrationInfo<dim> &info1,
+            typename MeshWorker::IntegrationInfo<dim> &info2) const;
+};
+
+template <int dim>
+void SystemIntegrator<dim>::cell(
+  MeshWorker::DoFInfo<dim> &dinfo,
+  typename MeshWorker::IntegrationInfo<dim> &info) const
+{
+  // Matrix
+  const FEValuesBase<dim> &fe = info.fe_values();
+
+  FullMatrix<double> &local_matrix = dinfo.matrix(0).matrix;
+  const unsigned int n_points = fe.n_quadrature_points;
+  const unsigned int n_dofs = fe.dofs_per_cell;
+  const Advection<dim> advection;
+
+  for (unsigned int k = 0; k < n_points; ++k)
+    {
+      const Tensor<1,dim> advection_directions = advection.value(fe.quadrature_point(k));
+      for (unsigned int i = 0; i < n_dofs; ++i)
+        for (unsigned int j = 0; j < n_dofs; ++j)
+          local_matrix(i, j) += ( fe.shape_grad(i,k)* fe.shape_grad(j,k)
+                                  + fe.shape_value(i,k)*(advection_directions*fe.shape_grad(j,k))
+                                ) * fe.JxW(k);
+    }
+
+  // RHS
+  Vector<double> &b = dinfo.vector(0).block(0);
+  const RightHandSide<dim>  right_hand_side;
+
+  for (unsigned int k=0; k<n_points; ++k)
+    {
+      const double fval = right_hand_side.value(fe.quadrature_point(k));
+      for (unsigned int i=0; i<n_dofs; ++i)
+        b(i) += fe.JxW(k)*(fe.shape_value(i,k) * fval);
+    }
+}
+
+
+template <int dim>
+void SystemIntegrator<dim>::boundary(
+  MeshWorker::DoFInfo<dim> &dinfo,
+  typename MeshWorker::IntegrationInfo<dim> &info) const
+{
+  const FEValuesBase<dim> &fe = info.fe_values();
+  const unsigned int n_points = fe.n_quadrature_points;
+  const unsigned int n_dofs = fe.dofs_per_cell;
+  Vector<double> &b = dinfo.vector(0).block(0);
+  const FluxBoundaryValues<dim>  flux_bd;
+
+  for (unsigned int k=0; k<n_points; ++k)
+    {
+      const double flux_bd_value = flux_bd.value(fe.quadrature_point(k));
+      for (unsigned int i=0; i<n_dofs; ++i)
+        b(i) += fe.JxW(k)*(fe.shape_value(i,k) * flux_bd_value);
+    }
+}
+
+
+template <int dim>
+void SystemIntegrator<dim>::face(
+  MeshWorker::DoFInfo<dim> &dinfo1,
+  MeshWorker::DoFInfo<dim> &dinfo2,
+  typename MeshWorker::IntegrationInfo<dim> &info1,
+  typename MeshWorker::IntegrationInfo<dim> &info2) const
+{
+  FullMatrix<double> &A11 = dinfo1.matrix(0,false).matrix;
+  FullMatrix<double> &A12 = dinfo1.matrix(0,true).matrix;
+  FullMatrix<double> &A21 = dinfo2.matrix(0,true).matrix;
+  FullMatrix<double> &A22 = dinfo2.matrix(0,false).matrix;
+  const FEValuesBase<dim> &fe1 = info1.fe_values(0);
+  const FEValuesBase<dim> &fe2 = info2.fe_values(0);
+  const unsigned int n_points = fe1.n_quadrature_points;
+  const unsigned int n_dofs = fe1.dofs_per_cell;
+  double h_e;
+  if (dim==3) h_e = std::sqrt(dinfo1.face->measure());
+  else h_e = dinfo1.face->measure();
+
+  for (unsigned int k=0; k<n_points; ++k)
+    {
+      const double dx = fe1.JxW(k);
+      const Point<dim> &n = (Point<dim>)fe1.normal_vector(k);
+      for (unsigned int i=0; i<n_dofs; ++i)
+        {
+          // average
+          const double value1i = 0.5*fe1.shape_value(i,k);
+          const double value2i = 0.5*fe2.shape_value(i,k);
+
+          for (unsigned int j=0; j<n_dofs; ++j)
+            {
+              // normal jump
+              const double grad1j = -n * fe1.shape_grad(j,k);
+              const double grad2j = n * fe2.shape_grad(j,k);
+
+              A11(i,j) += fe1.JxW(k)*( value1i*grad1j );
+              A12(i,j) += fe1.JxW(k)*( value1i*grad2j );
+              A21(i,j) += fe1.JxW(k)*( value2i*grad1j );
+              A22(i,j) += fe1.JxW(k)*( value2i*grad2j );
+            }
+        }
+    }
+}
+
+
+/*
+ * Matrix integrator
+ */
+template <int dim>
+class MatrixIntegrator : public MeshWorker::LocalIntegrator<dim>
+{
+public:
+  void cell(MeshWorker::DoFInfo<dim> &dinfo,
+            typename MeshWorker::IntegrationInfo<dim> &info) const;
+  void boundary(MeshWorker::DoFInfo<dim> &dinfo,
+                typename MeshWorker::IntegrationInfo<dim> &info) const;
+  void face(MeshWorker::DoFInfo<dim> &dinfo1,
+            MeshWorker::DoFInfo<dim> &dinfo2,
+            typename MeshWorker::IntegrationInfo<dim> &info1,
+            typename MeshWorker::IntegrationInfo<dim> &info2) const;
+};
+
+template <int dim>
+void MatrixIntegrator<dim>::cell(
+  MeshWorker::DoFInfo<dim> &dinfo,
+  typename MeshWorker::IntegrationInfo<dim> &info) const
+{
+  const FEValuesBase<dim> &fe = info.fe_values();
+  FullMatrix<double> &local_matrix = dinfo.matrix(0).matrix;
+
+  const unsigned int n_points = fe.n_quadrature_points;
+  const unsigned int n_dofs = fe.dofs_per_cell;
+  const Advection<dim> advection;
+
+  for (unsigned int k = 0; k < n_points; ++k)
+    {
+      const Tensor<1,dim> advection_directions = advection.value(fe.quadrature_point(k));
+      for (unsigned int i = 0; i < n_dofs; ++i)
+        for (unsigned int j = 0; j < n_dofs; ++j)
+          local_matrix(i, j) += ( fe.shape_grad(i,k)* fe.shape_grad(j,k)
+                                  + fe.shape_value(i,k)*(advection_directions*fe.shape_grad(j,k))
+                                ) * fe.JxW(k);
+    }
+}
+
+
+template <int dim>
+void MatrixIntegrator<dim>::boundary(
+  MeshWorker::DoFInfo<dim> &dinfo,
+  typename MeshWorker::IntegrationInfo<dim> &info) const
+{
+}
+
+
+template <int dim>
+void MatrixIntegrator<dim>::face(
+  MeshWorker::DoFInfo<dim> &dinfo1,
+  MeshWorker::DoFInfo<dim> &dinfo2,
+  typename MeshWorker::IntegrationInfo<dim> &info1,
+  typename MeshWorker::IntegrationInfo<dim> &info2) const
+{
+  FullMatrix<double> &A11 = dinfo1.matrix(0,false).matrix;
+  FullMatrix<double> &A12 = dinfo1.matrix(0,true).matrix;
+  FullMatrix<double> &A21 = dinfo2.matrix(0,true).matrix;
+  FullMatrix<double> &A22 = dinfo2.matrix(0,false).matrix;
+  const FEValuesBase<dim> &fe1 = info1.fe_values(0);
+  const FEValuesBase<dim> &fe2 = info2.fe_values(0);
+  const unsigned int n_points = fe1.n_quadrature_points;
+  const unsigned int n_dofs = fe1.dofs_per_cell;
+  double h_e;
+  if (dim==3) h_e = std::sqrt(dinfo1.face->measure());
+  else h_e = dinfo1.face->measure();
+
+  for (unsigned int k=0; k<n_points; ++k)
+    {
+      const double dx = fe1.JxW(k);
+      const Point<dim> &n = (Point<dim>)fe1.normal_vector(k);
+      for (unsigned int i=0; i<n_dofs; ++i)
+        {
+          // average
+          const double value1i = 0.5*fe1.shape_value(i,k);
+          const double value2i = 0.5*fe2.shape_value(i,k);
+
+          for (unsigned int j=0; j<n_dofs; ++j)
+            {
+              // normal jump
+              const double grad1j = -n * fe1.shape_grad(j,k);
+              const double grad2j = n * fe2.shape_grad(j,k);
+
+              A11(i,j) += fe1.JxW(k)*( value1i*grad1j );
+              A12(i,j) += fe1.JxW(k)*( value1i*grad2j );
+              A21(i,j) += fe1.JxW(k)*( value2i*grad1j );
+              A22(i,j) += fe1.JxW(k)*( value2i*grad2j );
+            }
+        }
+    }
+}
+
+
+/*
+ * Vector integrator
+ */
+template <int dim>
+class RHSIntegrator : public MeshWorker::LocalIntegrator<dim>
+{
+public:
+  void cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
+  void boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const;
+  void face(MeshWorker::DoFInfo<dim> &dinfo1,
+            MeshWorker::DoFInfo<dim> &dinfo2,
+            typename MeshWorker::IntegrationInfo<dim> &info1,
+            typename MeshWorker::IntegrationInfo<dim> &info2) const;
+};
+
+
+template <int dim>
+void RHSIntegrator<dim>::cell(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const
+{
+  const FEValuesBase<dim> &fe = info.fe_values();
+  const unsigned int n_points = fe.n_quadrature_points;
+  const unsigned int n_dofs = fe.dofs_per_cell;
+  Vector<double> &b = dinfo.vector(0).block(0);
+  const RightHandSide<dim>  right_hand_side;
+
+  for (unsigned int k=0; k<n_points; ++k)
+    {
+      const double fval = right_hand_side.value(fe.quadrature_point(k));
+      for (unsigned int i=0; i<n_dofs; ++i)
+        b(i) += fe.JxW(k)*(fe.shape_value(i,k) * fval);
+    }
+}
+
+
+template <int dim>
+void RHSIntegrator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo, typename MeshWorker::IntegrationInfo<dim> &info) const
+{
+  const FEValuesBase<dim> &fe = info.fe_values();
+  const unsigned int n_points = fe.n_quadrature_points;
+  const unsigned int n_dofs = fe.dofs_per_cell;
+  Vector<double> &b = dinfo.vector(0).block(0);
+  const FluxBoundaryValues<dim>  flux_bd;
+
+  for (unsigned int k=0; k<n_points; ++k)
+    {
+      const double flux_bd_value = flux_bd.value(fe.quadrature_point(k));
+      for (unsigned int i=0; i<n_dofs; ++i)
+        b(i) += fe.JxW(k)*(fe.shape_value(i,k) * flux_bd_value);
+    }
+}
+
+
+template <int dim>
+void RHSIntegrator<dim>::face(MeshWorker::DoFInfo<dim> &,
+                              MeshWorker::DoFInfo<dim> &,
+                              typename MeshWorker::IntegrationInfo<dim> &,
+                              typename MeshWorker::IntegrationInfo<dim> &) const
+{}
+
+
+/*
+ * Main class
+ */
+template<int dim>
+class MeshWorkerConstraintMatrixTest
+{
+public:
+  MeshWorkerConstraintMatrixTest(const FiniteElement<dim> &fe);
+  ~MeshWorkerConstraintMatrixTest();
+
+  void run();
+
+private:
+  void setup_system();
+  void createInhomConstraints();
+  void assemble_system_MeshWorker();
+  void assemble_MeshWorker();
+
+  Triangulation<dim> triangulation;
+  const MappingQ1<dim> mapping;
+  DoFHandler<dim> dof_handler;
+  const FiniteElement<dim> &fe;
+
+  SparsityPattern sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+  SparseMatrix<double> matrix;
+  SparseMatrix<double> error_matrix;
+
+  Vector<double> system_rhs;
+  Vector<double> rhs;
+  Vector<double> error_rhs;
+
+  ConstraintMatrix constraintsInhom;
+  ConstraintMatrix constraints;
+  ConstraintMatrix constraintsAll;
+};
+
+
+template<int dim>
+MeshWorkerConstraintMatrixTest<dim>::MeshWorkerConstraintMatrixTest(const FiniteElement<dim> &fe) :
+  mapping(),
+  dof_handler(triangulation),
+  fe(fe)
+{
+  GridGenerator::hyper_cube(this->triangulation, -1, 1);
+
+  //refine with hanging node
+  this->triangulation.refine_global (1);
+  this->triangulation.begin_active()->set_refine_flag ();
+  this->triangulation.execute_coarsening_and_refinement ();
+  this->triangulation.refine_global (1);
+}
+
+
+template <int dim>
+MeshWorkerConstraintMatrixTest<dim>::~MeshWorkerConstraintMatrixTest ()
+{
+  dof_handler.clear ();
+}
+
+
+template<int dim>
+void MeshWorkerConstraintMatrixTest<dim>::setup_system()
+{
+  dof_handler.distribute_dofs(fe);
+
+  system_rhs.reinit(dof_handler.n_dofs());
+  rhs.reinit(dof_handler.n_dofs());
+  error_rhs.reinit(dof_handler.n_dofs());
+
+  constraints.clear();
+  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+  constraints.close();
+
+  DynamicSparsityPattern c_sparsity(dof_handler.n_dofs());
+
+  std::string str("FE_DGQ");
+  std::string fe_name = fe.get_name();
+  std::size_t found = fe_name.find(str);
+  if ( found!=std::string::npos)
+    DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity, constraints,
+                                         false);
+  else
+    DoFTools::make_sparsity_pattern(dof_handler, c_sparsity, constraints,
+                                    false);
+
+  sparsity_pattern.copy_from(c_sparsity);
+
+  system_matrix.reinit(sparsity_pattern);
+  matrix.reinit(sparsity_pattern);
+  error_matrix.reinit(sparsity_pattern);
+}
+
+/*
+ * Assemble with SystemSimple
+ */
+template<int dim>
+void MeshWorkerConstraintMatrixTest<dim>::assemble_system_MeshWorker()
+{
+
+  MeshWorker::IntegrationInfoBox<dim> info_box;
+
+  const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
+  info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points,
+                                       n_gauss_points);
+
+  UpdateFlags update_flags = update_quadrature_points |
+                             update_values |
+                             update_gradients;
+
+  info_box.add_update_flags_all(update_flags);
+  info_box.initialize(fe, mapping);
+
+  MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+  MeshWorker::Assembler::SystemSimple<SparseMatrix<double>, Vector<double> > assembler;
+
+  assembler.initialize(system_matrix, system_rhs);
+  assembler.initialize(constraintsAll);
+
+  SystemIntegrator<dim> matrix_integrator;
+  MeshWorker::integration_loop<dim, dim> (
+    dof_handler.begin_active(), dof_handler.end(),
+    dof_info, info_box, matrix_integrator, assembler);
+
+}
+
+/*
+ * Assemble matrix and vector seperately
+ */
+template<int dim>
+void MeshWorkerConstraintMatrixTest<dim>::assemble_MeshWorker()
+{
+
+  MeshWorker::IntegrationInfoBox<dim> info_box;
+
+  const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
+  info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points,
+                                       n_gauss_points);
+
+  UpdateFlags update_flags = update_quadrature_points |
+                             update_values |
+                             update_gradients;
+
+  info_box.add_update_flags_all(update_flags);
+  info_box.initialize(fe, mapping);
+
+  MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+  MeshWorker::Assembler::MatrixSimple<SparseMatrix<double> > assemblerMatrix;
+  assemblerMatrix.initialize(constraints);
+  assemblerMatrix.initialize(matrix);
+  MeshWorker::Assembler::ResidualSimple<Vector<double> > assemblerVector;
+  assemblerVector.initialize(constraints);
+  AnyData data;
+  data.add<Vector<double>*>(&rhs, "RHS");
+  assemblerVector.initialize(data);
+
+  MatrixIntegrator<dim> matrix_integrator;
+  MeshWorker::integration_loop<dim, dim> (
+    dof_handler.begin_active(), dof_handler.end(),
+    dof_info, info_box, matrix_integrator, assemblerMatrix);
+  RHSIntegrator<dim> vector_integrator;
+  MeshWorker::integration_loop<dim, dim> (
+    dof_handler.begin_active(), dof_handler.end(),
+    dof_info, info_box, vector_integrator, assemblerVector);
+}
+
+
+template <int dim>
+void MeshWorkerConstraintMatrixTest<dim>::createInhomConstraints()
+{
+  this->constraintsInhom.clear();
+  // boundary constraints
+  VectorTools::interpolate_boundary_values (this->dof_handler,
+                                            0,
+                                            BoundaryValues<dim>(),
+                                            this->constraintsInhom);
+  // constraint of the origin
+  std::map<types::global_dof_index, Point<dim> > support_points;
+  DoFTools::map_dofs_to_support_points(this->mapping,this->dof_handler,support_points);
+  for (unsigned int dof_index=0; dof_index < this->dof_handler.n_dofs(); ++dof_index)
+    {
+      if (!this->constraints.is_constrained(dof_index) && !this->constraintsInhom.is_constrained(dof_index))
+        {
+          const Point<dim> p = support_points[dof_index];
+
+          if  (std::sqrt(p.square()) < 1e-6)
+            {
+              this->constraintsInhom.add_line (dof_index);
+              this->constraintsInhom.set_inhomogeneity (dof_index, 2);
+            }
+        }
+    }
+  this->constraintsInhom.close();
+}
+
+
+template<int dim>
+void MeshWorkerConstraintMatrixTest<dim>::run()
+{
+  setup_system();
+  createInhomConstraints();
+  constraintsAll.clear();
+  constraintsAll.merge(constraints);
+  constraintsAll.merge(constraintsInhom);
+  constraintsAll.close();
+
+  assemble_system_MeshWorker();
+
+  assemble_MeshWorker();
+  constraintsInhom.condense(matrix,rhs);
+
+  // make diagonal entries and constraints component equal
+  constraintsAll.distribute(system_rhs);
+  constraintsAll.distribute(rhs);
+  for (unsigned int i=0; i<error_matrix.m(); ++i)
+    if (constraintsAll.is_constrained(i))
+      {
+        system_matrix.diag_element(i) = 1;
+        matrix.diag_element(i) = 1;
+      }
+
+  // evaluate difference
+  error_matrix.copy_from(system_matrix);
+  error_matrix.add(-1.0,matrix);
+  error_rhs = system_rhs;
+  error_rhs -= rhs;
+  deallog << "difference rhs "<< error_rhs.l2_norm() << std::endl;
+  deallog << "difference matrix "<< error_matrix.frobenius_norm() << std::endl;
+
+}
+
+
+int main()
+{
+  initlog();
+  deallog.threshold_double(1.e-10);
+
+  FE_Q<2> fe(1);
+  deallog.push(fe.get_name());
+  MeshWorkerConstraintMatrixTest<2> test(fe);
+  test.run();
+  deallog.pop ();
+
+  FE_DGQ<2> fe2(1);
+  deallog.push(fe2.get_name());
+  MeshWorkerConstraintMatrixTest<2> test2(fe2);
+  test2.run();
+  deallog.pop ();
+
+}
+
diff --git a/tests/integrators/assembler_simple_system_inhom_01.output b/tests/integrators/assembler_simple_system_inhom_01.output
new file mode 100644 (file)
index 0000000..2d0ce2e
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL:FE_Q<2>(1)::difference rhs 0
+DEAL:FE_Q<2>(1)::difference matrix 0
+DEAL:FE_DGQ<2>(1)::difference rhs 0
+DEAL:FE_DGQ<2>(1)::difference matrix 0

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.