From 604cf6141f8c94d2f49e3b4f1384931a62a815b2 Mon Sep 17 00:00:00 2001 From: kanschat Date: Mon, 15 Feb 2010 16:24:15 +0000 Subject: [PATCH] bundle cell and face data in loop git-svn-id: https://svn.dealii.org/trunk@20634 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/deal.II/mesh_worker_01.cc | 9 +- tests/deal.II/mesh_worker_02.cc | 279 +++++++++++++++++++++++ tests/deal.II/mesh_worker_02/cmp/generic | 29 +++ 3 files changed, 313 insertions(+), 4 deletions(-) create mode 100644 tests/deal.II/mesh_worker_02.cc create mode 100644 tests/deal.II/mesh_worker_02/cmp/generic diff --git a/tests/deal.II/mesh_worker_01.cc b/tests/deal.II/mesh_worker_01.cc index 6f83b2df05..37c57bb29f 100644 --- a/tests/deal.II/mesh_worker_01.cc +++ b/tests/deal.II/mesh_worker_01.cc @@ -157,13 +157,14 @@ test_simple(MGDoFHandler& mgdofs) { MappingQ1 mapping; - MeshWorker::IntegrationInfoBox info_box(dofs); + MeshWorker::IntegrationInfoBox info_box; info_box.initialize(integrator, fe, mapping); + MeshWorker::DoFInfo dof_info(dofs); - MeshWorker::integration_loop - ::CellInfo, typename Local::FaceInfo, dim> + MeshWorker::loop + , MeshWorker::IntegrationInfoBox > (dofs.begin_active(), dofs.end(), - info_box, + dof_info, info_box, std_cxx1x::bind (&Local::cell, local, _1, _2), std_cxx1x::bind (&Local::bdry, local, _1, _2), std_cxx1x::bind (&Local::face, local, _1, _2, _3, _4), diff --git a/tests/deal.II/mesh_worker_02.cc b/tests/deal.II/mesh_worker_02.cc new file mode 100644 index 0000000000..797ae07062 --- /dev/null +++ b/tests/deal.II/mesh_worker_02.cc @@ -0,0 +1,279 @@ +//---------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 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. +// +//---------------------------------------------------------------------- + +// Test consistency of assemblers MatrixSimple and MGMatrixSimple + +#include "../tests.h" +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +using namespace dealii; + + +// Use local matrices for Laplacian / interior penalty DG +template +class MatrixIntegrator : public Subscriptor +{ + public: + void cell(typename MeshWorker::IntegrationWorker::CellInfo& info) const; + void bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) const; + void face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, + typename MeshWorker::IntegrationWorker::FaceInfo& info2) const; +}; + + +template +void MatrixIntegrator::cell(typename MeshWorker::IntegrationWorker::CellInfo& info) const +{ + const FEValuesBase& fe = info.fe(); + FullMatrix& local_matrix = info.M1[0].matrix; + + for (unsigned int k=0; k +void MatrixIntegrator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) const +{ + const FEFaceValuesBase& fe = info.fe(); + FullMatrix& local_matrix = info.M1[0].matrix; + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty = 2. * deg * (deg+1) * info.face->measure() / info.cell->measure(); + + for (unsigned k=0;k +void MatrixIntegrator::face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, + typename MeshWorker::IntegrationWorker::FaceInfo& info2) const +{ + const FEFaceValuesBase& fe1 = info1.fe(); + const FEFaceValuesBase& fe2 = info2.fe(); + FullMatrix& matrix_v1u1 = info1.M1[0].matrix; + FullMatrix& matrix_v1u2 = info1.M2[0].matrix; + FullMatrix& matrix_v2u1 = info2.M2[0].matrix; + FullMatrix& matrix_v2u2 = info2.M1[0].matrix; + + const unsigned int deg = fe1.get_fe().tensor_degree(); + const double penalty1 = deg * (deg+1) * info1.face->measure() / info1.cell->measure(); + const double penalty2 = deg * (deg+1) * info2.face->measure() / info2.cell->measure(); + const double penalty = penalty1 + penalty2; + + for (unsigned k=0;k +void +assemble(const DoFHandler& dof_handler, SparseMatrix& matrix) +{ + const FiniteElement& fe = dof_handler.get_fe(); + MappingQ1 mapping; + + const MatrixIntegrator local; + MeshWorker::AssemblingIntegrator >, MatrixIntegrator > + integrator(local); + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; + integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points); + UpdateFlags update_flags = update_values | update_gradients; + integrator.add_update_flags(update_flags, true, true, true, true); + + integrator.initialize(matrix); + MeshWorker::IntegrationInfoBox info_box(dof_handler); + info_box.initialize(integrator, fe, mapping); + MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator); +} + + +template +void +assemble(const MGDoFHandler& dof_handler, + MGLevelObject > matrix, + MGLevelObject > dg_up, + MGLevelObject > dg_down) +{ + const FiniteElement& fe = dof_handler.get_fe(); + MappingQ1 mapping; + + const MatrixIntegrator local; + MeshWorker::AssemblingIntegrator >, MatrixIntegrator > + integrator(local); + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; + integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points); + UpdateFlags update_flags = update_values | update_gradients; + integrator.add_update_flags(update_flags, true, true, true, true); + + integrator.initialize(matrix); + integrator.initialize_fluxes(dg_up, dg_down); + MeshWorker::IntegrationInfoBox info_box(dof_handler); + info_box.initialize(integrator, fe, mapping); + MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator); +} + + +template +void +test_simple(MGDoFHandler& mgdofs) +{ + SparsityPattern pattern; + SparseMatrix matrix; + Vector v; + + const DoFHandler& dofs = mgdofs; + const FiniteElement& fe = dofs.get_fe(); + pattern.reinit (dofs.n_dofs(), dofs.n_dofs(), + (GeometryInfo::faces_per_cell + *GeometryInfo::max_children_per_face+1)*fe.dofs_per_cell); + DoFTools::make_flux_sparsity_pattern (dofs, pattern); + pattern.compress(); + matrix.reinit (pattern); + + assemble(dofs, matrix); + deallog << std::setprecision(3); + matrix.print(deallog); + + MGLevelObject mg_sparsity; + MGLevelObject mg_sparsity_dg_interface; + MGLevelObject > mg_matrix; + MGLevelObject > mg_matrix_dg_up; + MGLevelObject > mg_matrix_dg_down; + + const unsigned int n_levels = mgdofs.get_tria().n_levels(); + + mg_sparsity.resize(0, n_levels-1); + mg_sparsity_dg_interface.resize(0, n_levels-1); + mg_matrix.resize(0, n_levels-1); + mg_matrix_dg_up.resize(0, n_levels-1); + mg_matrix_dg_down.resize(0, n_levels-1); + + for (unsigned int level=mg_sparsity.get_minlevel(); + level<=mg_sparsity.get_maxlevel();++level) + { + CompressedSparsityPattern c_sparsity(mgdofs.n_dofs(level)); + CompressedSparsityPattern ci_sparsity; + if (level>0) + ci_sparsity.reinit(mgdofs.n_dofs(level-1), mgdofs.n_dofs(level)); + + MGTools::make_flux_sparsity_pattern(mgdofs, c_sparsity, level); + if (level>0) + MGTools::make_flux_sparsity_pattern_edge(mgdofs, ci_sparsity, level); + + mg_sparsity[level].copy_from(c_sparsity); + mg_matrix[level].reinit(mg_sparsity[level]); + if (level>0) + { + mg_sparsity_dg_interface[level].copy_from(ci_sparsity); + mg_matrix_dg_up[level].reinit(mg_sparsity_dg_interface[level]); + mg_matrix_dg_down[level].reinit(mg_sparsity_dg_interface[level]); + } + } + +} + + +template +void +test(const FiniteElement& fe) +{ + Triangulation tr; + GridGenerator::hyper_L(tr); + tr.begin()->set_refine_flag(); + tr.execute_coarsening_and_refinement(); +// tr.begin(2)->set_refine_flag(); +// tr.execute_coarsening_and_refinement(); +// tr.refine_global(1); + deallog << "Triangulation levels"; + for (unsigned int l=0;l::cell_iterator cell = tr.begin(); + cell != tr.end(); ++cell, ++cn) + cell->set_user_index(cn); + + MGDoFHandler dofs(tr); + dofs.distribute_dofs(fe); + deallog << "DoFHandler " << dofs.n_dofs() << " levels"; + for (unsigned int l=0;l > > fe2; + fe2.push_back(boost::shared_ptr >(new FE_DGP<2>(1))); +// fe2.push_back(boost::shared_ptr >(new FE_Q<2>(1))); + + for (unsigned int i=0;i