From 54211d02fbd678af0377b5ba5e7ed7a18a977e03 Mon Sep 17 00:00:00 2001 From: Joscha Gedicke Date: Wed, 9 Mar 2016 18:31:01 +0100 Subject: [PATCH] test assembler_simple_system_inhom_01 --- .../assembler_simple_system_inhom_01.cc | 691 ++++++++++++++++++ .../assembler_simple_system_inhom_01.output | 5 + 2 files changed, 696 insertions(+) create mode 100644 tests/integrators/assembler_simple_system_inhom_01.cc create mode 100644 tests/integrators/assembler_simple_system_inhom_01.output diff --git a/tests/integrators/assembler_simple_system_inhom_01.cc b/tests/integrators/assembler_simple_system_inhom_01.cc new file mode 100644 index 0000000000..60373cc371 --- /dev/null +++ b/tests/integrators/assembler_simple_system_inhom_01.cc @@ -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 +#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; + + +/* + * Right hand side + */ +template +class RightHandSide : public Function +{ +public: + RightHandSide () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + +template +double RightHandSide::value (const Point &p, + const unsigned int component) const +{ + Assert (component == 0, ExcNotImplemented()); + + double val = 1; // f = 1 + return val; +} + + +/* + * Boundary values + */ +template +class BoundaryValues : public Function +{ +public: + BoundaryValues () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + +template +double BoundaryValues::value (const Point &p, + const unsigned int component) const +{ + Assert (component == 0, ExcNotImplemented()); + + double val = 1; // u_D = 1 + return val; +} + + +/* + * Flux boundary values + */ +template +class FluxBoundaryValues : public Function +{ +public: + FluxBoundaryValues () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + +template +double FluxBoundaryValues::value (const Point &p, + const unsigned int component) const +{ + double val = 1; // g = 1 + return val; +} + + +/* + * Advection coefficient + */ +template +class Advection : public TensorFunction<1,dim> +{ +public: + Advection () : TensorFunction<1,dim> () {} + + virtual Tensor<1,dim> value (const Point &p) const; +}; + +template +Tensor<1,dim> +Advection::value (const Point &p) const +{ + Point value; // beta = (1,0)^t + value[0] = 1; + value[1] = 0; + + return value; +} + + +/* + * System integrator + */ +template +class SystemIntegrator : public MeshWorker::LocalIntegrator +{ +public: + void cell(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; +}; + +template +void SystemIntegrator::cell( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const +{ + // Matrix + const FEValuesBase &fe = info.fe_values(); + + FullMatrix &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 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 &b = dinfo.vector(0).block(0); + const RightHandSide right_hand_side; + + for (unsigned int k=0; k +void SystemIntegrator::boundary( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const +{ + const FEValuesBase &fe = info.fe_values(); + const unsigned int n_points = fe.n_quadrature_points; + const unsigned int n_dofs = fe.dofs_per_cell; + Vector &b = dinfo.vector(0).block(0); + const FluxBoundaryValues flux_bd; + + for (unsigned int k=0; k +void SystemIntegrator::face( + MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const +{ + FullMatrix &A11 = dinfo1.matrix(0,false).matrix; + FullMatrix &A12 = dinfo1.matrix(0,true).matrix; + FullMatrix &A21 = dinfo2.matrix(0,true).matrix; + FullMatrix &A22 = dinfo2.matrix(0,false).matrix; + const FEValuesBase &fe1 = info1.fe_values(0); + const FEValuesBase &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 = (Point)fe1.normal_vector(k); + for (unsigned int i=0; i +class MatrixIntegrator : public MeshWorker::LocalIntegrator +{ +public: + void cell(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; +}; + +template +void MatrixIntegrator::cell( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const +{ + const FEValuesBase &fe = info.fe_values(); + FullMatrix &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 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 +void MatrixIntegrator::boundary( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const +{ +} + + +template +void MatrixIntegrator::face( + MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const +{ + FullMatrix &A11 = dinfo1.matrix(0,false).matrix; + FullMatrix &A12 = dinfo1.matrix(0,true).matrix; + FullMatrix &A21 = dinfo2.matrix(0,true).matrix; + FullMatrix &A22 = dinfo2.matrix(0,false).matrix; + const FEValuesBase &fe1 = info1.fe_values(0); + const FEValuesBase &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 = (Point)fe1.normal_vector(k); + for (unsigned int i=0; i +class RHSIntegrator : public MeshWorker::LocalIntegrator +{ +public: + void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; +}; + + +template +void RHSIntegrator::cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const +{ + const FEValuesBase &fe = info.fe_values(); + const unsigned int n_points = fe.n_quadrature_points; + const unsigned int n_dofs = fe.dofs_per_cell; + Vector &b = dinfo.vector(0).block(0); + const RightHandSide right_hand_side; + + for (unsigned int k=0; k +void RHSIntegrator::boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const +{ + const FEValuesBase &fe = info.fe_values(); + const unsigned int n_points = fe.n_quadrature_points; + const unsigned int n_dofs = fe.dofs_per_cell; + Vector &b = dinfo.vector(0).block(0); + const FluxBoundaryValues flux_bd; + + for (unsigned int k=0; k +void RHSIntegrator::face(MeshWorker::DoFInfo &, + MeshWorker::DoFInfo &, + typename MeshWorker::IntegrationInfo &, + typename MeshWorker::IntegrationInfo &) const +{} + + +/* + * Main class + */ +template +class MeshWorkerConstraintMatrixTest +{ +public: + MeshWorkerConstraintMatrixTest(const FiniteElement &fe); + ~MeshWorkerConstraintMatrixTest(); + + void run(); + +private: + void setup_system(); + void createInhomConstraints(); + void assemble_system_MeshWorker(); + void assemble_MeshWorker(); + + Triangulation triangulation; + const MappingQ1 mapping; + DoFHandler dof_handler; + const FiniteElement &fe; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + SparseMatrix matrix; + SparseMatrix error_matrix; + + Vector system_rhs; + Vector rhs; + Vector error_rhs; + + ConstraintMatrix constraintsInhom; + ConstraintMatrix constraints; + ConstraintMatrix constraintsAll; +}; + + +template +MeshWorkerConstraintMatrixTest::MeshWorkerConstraintMatrixTest(const FiniteElement &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 +MeshWorkerConstraintMatrixTest::~MeshWorkerConstraintMatrixTest () +{ + dof_handler.clear (); +} + + +template +void MeshWorkerConstraintMatrixTest::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 +void MeshWorkerConstraintMatrixTest::assemble_system_MeshWorker() +{ + + MeshWorker::IntegrationInfoBox 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 dof_info(dof_handler); + MeshWorker::Assembler::SystemSimple, Vector > assembler; + + assembler.initialize(system_matrix, system_rhs); + assembler.initialize(constraintsAll); + + SystemIntegrator matrix_integrator; + MeshWorker::integration_loop ( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, matrix_integrator, assembler); + +} + +/* + * Assemble matrix and vector seperately + */ +template +void MeshWorkerConstraintMatrixTest::assemble_MeshWorker() +{ + + MeshWorker::IntegrationInfoBox 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 dof_info(dof_handler); + MeshWorker::Assembler::MatrixSimple > assemblerMatrix; + assemblerMatrix.initialize(constraints); + assemblerMatrix.initialize(matrix); + MeshWorker::Assembler::ResidualSimple > assemblerVector; + assemblerVector.initialize(constraints); + AnyData data; + data.add*>(&rhs, "RHS"); + assemblerVector.initialize(data); + + MatrixIntegrator matrix_integrator; + MeshWorker::integration_loop ( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, matrix_integrator, assemblerMatrix); + RHSIntegrator vector_integrator; + MeshWorker::integration_loop ( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, vector_integrator, assemblerVector); +} + + +template +void MeshWorkerConstraintMatrixTest::createInhomConstraints() +{ + this->constraintsInhom.clear(); + // boundary constraints + VectorTools::interpolate_boundary_values (this->dof_handler, + 0, + BoundaryValues(), + this->constraintsInhom); + // constraint of the origin + std::map > 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 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 +void MeshWorkerConstraintMatrixTest::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 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 index 0000000000..2d0ce2e623 --- /dev/null +++ b/tests/integrators/assembler_simple_system_inhom_01.output @@ -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 -- 2.39.5