From 6eed3e44ab81ad73d332bd9d60db8211671f8d6c Mon Sep 17 00:00:00 2001 From: kanschat Date: Tue, 9 Nov 2010 15:45:13 +0000 Subject: [PATCH] First test for integrators... failing due to program crashing in constructor of FE_Nedelec in 3D. Results being produced for 2D and 3D-P0 are good git-svn-id: https://svn.dealii.org/trunk@22645 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/integrators/Makefile | 10 ++ tests/integrators/cochain_01.cc | 273 ++++++++++++++++++++++++++++++++ tests/lib/test_grids.h | 47 +++--- 3 files changed, 305 insertions(+), 25 deletions(-) create mode 100644 tests/integrators/Makefile create mode 100644 tests/integrators/cochain_01.cc diff --git a/tests/integrators/Makefile b/tests/integrators/Makefile new file mode 100644 index 0000000000..6476e92010 --- /dev/null +++ b/tests/integrators/Makefile @@ -0,0 +1,10 @@ +############################################################ +# $Id$ +# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors +############################################################ + +include ../Makefile.paths +include $D/common/Make.global_options +include ../Makefile.rules +include Makefile.depend +include Makefile.tests diff --git a/tests/integrators/cochain_01.cc b/tests/integrators/cochain_01.cc new file mode 100644 index 0000000000..f410187d27 --- /dev/null +++ b/tests/integrators/cochain_01.cc @@ -0,0 +1,273 @@ +//-------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2005, 2006, 2010 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, whether differential operators produce a cochain complex on +// the standard Hilbert space sequence + +#include "../tests.h" +#include "../lib/test_grids.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 + +using namespace LocalIntegrators; + +const bool debugging = false; + +template +void cell_matrix( + MeshWorker::DoFInfo& dinfo, + typename MeshWorker::IntegrationInfo& info) +{ + unsigned int dm=0; // Matrix index + unsigned int de=0; // Element index + + L2::mass_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de)); + Laplace::cell_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de)); + Differential::gradient_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de), info.fe_values(de+1)); + + + ++de; + L2::mass_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de)); + Maxwell::curl_curl_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de)); + Maxwell::curl_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de), info.fe_values(de+1)); + + if (dim>2) + { + ++de; + L2::mass_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de)); + Elasticity::grad_div_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de)); + Differential::divergence_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de), info.fe_values(de+1)); + } + + ++de; + L2::mass_matrix(dinfo.matrix(dm++,false).matrix, info.fe_values(de)); +} + + + +template +void +test_cochain(const Triangulation& tr, const FiniteElement& fe) +{ + MappingQ1 mapping; + // Initialize DofHandler for a + // block system with local blocks + DoFHandler dof(tr); + dof.distribute_dofs(fe); + DoFRenumbering::component_wise(dof); + dof.initialize_local_block_info(); + + // Initialize hanging node constraints + ConstraintMatrix constraints; +// DoFTools::make_zero_boundary_constraints(dof, constraints); + DoFTools::make_hanging_node_constraints (dof, constraints); + constraints.close(); + // Setup sparsity pattern for the + // whole system + BlockSparsityPattern sparsity; + sparsity.reinit(dof.block_info().global().size(), + dof.block_info().global().size()); + BlockCompressedSparsityPattern c_sparsity(dof.block_info().global(), + dof.block_info().global()); + DoFTools::make_flux_sparsity_pattern(dof, c_sparsity, constraints, false); + sparsity.copy_from(c_sparsity); + + // Setup matrices + MatrixBlockVector > matrices; + + unsigned int d=0; + matrices.add(d,d,"mass-H1"); + matrices.add(d,d,"Laplacian"); + matrices.add(d+1,d,"grad"); + + ++d; + matrices.add(d,d,"mass-Hcurl"); + matrices.add(d,d,"curl-curl"); + matrices.add(d+1,d,"curl"); + + if (dim>2) + { + ++d; + matrices.add(d,d,"mass-Hdiv"); + matrices.add(d,d,"grad-div"); + matrices.add(d+1,d,"div"); + } + + ++d; + matrices.add(d,d,"mass-L2"); + + matrices.reinit(sparsity); + + if (debugging) + for (unsigned int i=0;i info_box; + UpdateFlags update_flags = update_values | update_gradients; + info_box.add_update_flags(update_flags, true, true, true, true); + info_box.initialize(fe, mapping, &dof.block_info()); + + MeshWorker::DoFInfo dof_info(dof.block_info()); + + MeshWorker::Assembler::MatrixLocalBlocksToGlobalBlocks > assembler; + assembler.initialize(&dof.block_info(), matrices); + assembler.initialize(constraints); + + MeshWorker::integration_loop( + dof.begin_active(), dof.end(), dof_info, info_box, + cell_matrix, 0, 0, assembler); + + for (unsigned int b=0;b source(dof.block_info().global()); + BlockVector result1(dof.block_info().global()); + BlockVector result2(dof.block_info().global()); + BlockVector aux(dof.block_info().global()); + for (unsigned int i=0;i > solver(control); + + for (unsigned int d=0;d0) + { + matrices.matrix(m+2).vmult(result2.block(d+1), aux.block(d)); + deallog << "d^2 " << result2.block(d+1).l2_norm() << std::endl; + matrices.matrix(m+1).vmult(result2.block(d), aux.block(d)); + deallog << "d*d^2 " << result2.block(d).l2_norm() << std::endl; + } + + matrices.matrix(m+2).vmult(result1.block(d+1), source.block(d)); + + solver.solve(matrices.matrix(m+3), aux.block(d+1), result1.block(d+1), PreconditionIdentity()); + + matrices.matrix(m+2).Tvmult(result1.block(d), aux.block(d+1)); + matrices.matrix(m+1).vmult(result2.block(d), source.block(d)); + + if (debugging) + deallog << "u " << source.block(d).l2_norm() << ' ' + << matrices.name(m+2) << ' ' << result1.block(d+1).l2_norm() + << ' ' << matrices.name(m+3) << ' ' << aux.block(d+1).l2_norm() + << ' ' << matrices.name(m+2) << "^T " << result1.block(d).l2_norm() + << ' ' << matrices.name(m+1) << ' ' << result2.block(d).l2_norm() << std::endl; + result2.block(d) -= result1.block(d); + deallog << "Difference d*d " << result2.block(d).l2_norm() << std::endl; + } +} + +void run (const Triangulation<2>& tr, unsigned int degree) +{ + std::ostringstream prefix; + prefix << "d2-p" << degree; + deallog.push(prefix.str()); + deallog << "Setup" << std::endl; + + FE_Q<2> h1(degree+1); + FE_Nedelec<2> hdiv(degree); + FE_DGQ<2> l2(degree); + + FESystem<2> fe(h1,1,hdiv,1,l2,1); + + test_cochain(tr, fe); + deallog.pop(); +} + + +void run (const Triangulation<3>& tr, unsigned int degree) +{ + std::ostringstream prefix; + prefix << "d3-p" << degree; + deallog.push(prefix.str()); + deallog << "Setup" << std::endl; + + FE_Q<3> h1(degree+1); + if (debugging) deallog << "H1" << std::endl; + FE_Nedelec<3> hcurl(degree); + if (debugging) deallog << "Hcurl" << std::endl; + FE_RaviartThomas<3> hdiv(degree); + if (debugging) deallog << "Hdiv" << std::endl; + FE_DGQ<3> l2(degree); + if (debugging) deallog << "L2" << std::endl; + + FESystem<3> fe(h1,1,hcurl,1,hdiv,1,l2,1); + + test_cochain(tr, fe); + deallog.pop(); +} + + +int main() +{ + const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output"); + std::ofstream logfile(logname.c_str()); + deallog.attach(logfile); +// deallog.depth_console(0); +// deallog.threshold_double(1.e-10); + + Triangulation<2> tr2; + TestGrids::hypercube(tr2, 1); + run(tr2, 0); + run(tr2, 1); + run(tr2, 2); + + Triangulation<3> tr3; + TestGrids::hypercube(tr3, 1); + run(tr3, 0); + run(tr3, 1); +} diff --git a/tests/lib/test_grids.h b/tests/lib/test_grids.h index dcd1564c3c..717435e10a 100644 --- a/tests/lib/test_grids.h +++ b/tests/lib/test_grids.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2006 by the deal.II authors +// Copyright (C) 2006, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -12,6 +12,7 @@ //---------------------------------------------------------------------- #include +#include /** * A set of test meshes for the deal.II test suite. @@ -34,7 +35,8 @@ * level * * - * @author Guido Kanschat, 2006 + * @author Guido Kanschat + * @date 2006, 2010 */ namespace TestGrids { @@ -55,7 +57,24 @@ namespace TestGrids template void hypercube(Triangulation& tr, unsigned int refinement = 0, - bool local = false); + bool local = false) + { + GridGenerator::hyper_cube(tr); + if (refinement && !local) + tr.refine_global(refinement); + if (refinement && local) + { + for (typename Triangulation::active_cell_iterator + cell = tr.begin_active(); cell != tr.end(); ++cell) + { + const Point& p = cell->center(); + bool negative = true; + for (unsigned int d=0;d= 0.)negative = false; + } + } + } + /** * Create a star-shaped mesh, * having more than the average @@ -91,26 +110,4 @@ namespace TestGrids */ template void laguna(Triangulation& tr); - - - template - void hypercube(Triangulation& tr, - unsigned int refinement, - bool local) - { - GridGenerator::hyper_cube(tr); - if (!local) - { - if (refinement > 0) - tr.refine_global(refinement); - } - else - { - for (unsigned int i=0;iset_refine_flag(); - tr.execute_coarsening_and_refinement(); - } - } - } } -- 2.39.5