From: kronbichler Date: Tue, 16 Apr 2013 15:56:15 +0000 (+0000) Subject: Add test for mass matrix evaluation in matrix-free for some elements. Avoid deprecate... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5137ceb97fe681be3815ec933fc0d638121c89b3;p=dealii-svn.git Add test for mass matrix evaluation in matrix-free for some elements. Avoid deprecated function. git-svn-id: https://svn.dealii.org/trunk@29305 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/matrix_free/compress_mapping.cc b/tests/matrix_free/compress_mapping.cc index 35a2c51191..d4e3c08091 100644 --- a/tests/matrix_free/compress_mapping.cc +++ b/tests/matrix_free/compress_mapping.cc @@ -121,62 +121,19 @@ void test_cube () -void create_parallelogram(Triangulation<2> &tria) -{ - const int dim = 2; - std::vector > points (4); - points[0] = Point (0, 0); - points[1] = Point (0, 1); - points[2] = Point (1 ,0.5); - points[3] = Point (1 ,1.5); - - std::vector > cells(1); - cells[0].vertices[0] = 0; - cells[0].vertices[1] = 2; - cells[0].vertices[2] = 1; - cells[0].vertices[3] = 3; - cells[0].material_id = 0; - - tria.create_triangulation (points, cells, SubCellData()); -} - - - -void create_parallelogram(Triangulation<3> &tria) -{ - const int dim = 3; - std::vector > points (8); - points[0] = Point (0,0,0); - points[1] = Point (0,1.,0.5); - points[2] = Point (0,0,1); - points[3] = Point (0,1.,1.5); - points[4] = Point (1.,0,1.); - points[5] = Point (1.,1.,1.5); - points[6] = Point (1.,0,2); - points[7] = Point (1.,1.,2.5); - - std::vector > cells(1); - cells[0].vertices[0] = 0; - cells[0].vertices[1] = 4; - cells[0].vertices[2] = 1; - cells[0].vertices[3] = 5; - cells[0].vertices[4] = 2; - cells[0].vertices[5] = 6; - cells[0].vertices[6] = 3; - cells[0].vertices[7] = 7; - cells[0].material_id = 0; - - tria.create_triangulation (points, cells, SubCellData()); -} - - - template void test_parallelogram () { deallog << "Parallelogram" << std::endl; Triangulation tria; - create_parallelogram(tria); + Point corners[dim]; + for (unsigned int d=0; d0) + corners[d][0] = 0.5+0.5*d; + } + GridGenerator::parallelepiped(tria, corners); tria.refine_global(5-dim); FE_Q fe (1); diff --git a/tests/matrix_free/estimate_condition_number_mass.cc b/tests/matrix_free/estimate_condition_number_mass.cc new file mode 100644 index 0000000000..0a78638d37 --- /dev/null +++ b/tests/matrix_free/estimate_condition_number_mass.cc @@ -0,0 +1,156 @@ +//------------------ estimate_condition_number_mass.cc --------------------- +// $Id$ +// Version: $Name$ +// +//------------------ estimate_condition_number_mass.cc --------------------- + + +// this function computes condition number estimates for the mass matrix at +// different polynomial degrees. The mesh uses a hypercube mesh with no +// hanging nodes and no other constraints + +#include "../tests.h" + +std::ofstream logfile("estimate_condition_number_mass/output"); + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +void +mass_operator (const MatrixFree &data, + Vector &dst, + const Vector &src, + const std::pair &cell_range) +{ + FEEvaluationGeneral fe_eval (data); + const unsigned int n_q_points = fe_eval.n_q_points; + + for(unsigned int cell=cell_range.first;cell +class MatrixFreeTest +{ + public: + typedef VectorizedArray vector_t; + + MatrixFreeTest(const MatrixFree &data_in): + data (data_in) + {}; + + void vmult (Vector &dst, + const Vector &src) const + { + dst = 0; + const std_cxx1x::function &, + Vector &, + const Vector &, + const std::pair&)> + wrap = mass_operator; + data.cell_loop (wrap, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void test (const FiniteElement &fe) +{ + typedef double number; + + Triangulation tria; + GridGenerator::hyper_cube (tria); + tria.refine_global(2); + + DoFHandler dof(tria); + dof.distribute_dofs(fe); + ConstraintMatrix constraints; + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + MatrixFree mf_data; + { + const QGauss<1> quad (fe_degree+1); + mf_data.reinit (dof, constraints, quad); + } + + MatrixFreeTest mf (mf_data); + Vector in (dof.n_dofs()), out (dof.n_dofs()); + + VectorTools::create_right_hand_side(dof, QGauss(fe_degree+1), + Functions::CosineFunction(), in); + + SolverControl control(10000, 1e-9*in.l2_norm()); + typename SolverCG<>::AdditionalData data; + data.compute_condition_number = true; + SolverCG<> solver(control,data); + solver.solve(mf, out, in, PreconditionIdentity()); + deallog << std::endl; +} + + +int main () +{ + deallog.attach(logfile); + deallog.depth_console(0); + + deallog << std::setprecision (2); + + { + deallog.threshold_double(1.e-9); + deallog.push("2d"); + test<2,1>(FE_Q<2>(1)); + test<2,1>(FE_DGQ<2>(1)); + test<2,2>(FE_Q<2>(2)); + test<2,2>(FE_DGQ<2>(2)); + test<2,4>(FE_Q<2>(4)); + test<2,4>(FE_Q<2>(QGaussLobatto<1>(5))); + test<2,4>(FE_DGQ<2>(4)); + test<2,4>(FE_Q_Hierarchical<2>(4)); + test<2,6>(FE_Q<2>(6)); + test<2,6>(FE_Q<2>(QGaussLobatto<1>(7))); + test<2,6>(FE_DGQ<2>(6)); + test<2,6>(FE_DGQArbitraryNodes<2>(QGaussLobatto<1>(7))); + test<2,10>(FE_Q<2>(10)); + test<2,10>(FE_Q<2>(QGaussLobatto<1>(11))); + test<2,16>(FE_Q<2>(QGaussLobatto<1>(17))); + deallog.pop(); + deallog.push("3d"); + test<3,1>(FE_Q<3>(1)); + test<3,2>(FE_Q<3>(2)); + test<3,5>(FE_Q<3>(5)); + test<3,5>(FE_Q<3>(QGaussLobatto<1>(6))); + deallog.pop(); + } +} + + diff --git a/tests/matrix_free/estimate_condition_number_mass/cmp/generic b/tests/matrix_free/estimate_condition_number_mass/cmp/generic new file mode 100644 index 0000000000..9a991fe035 --- /dev/null +++ b/tests/matrix_free/estimate_condition_number_mass/cmp/generic @@ -0,0 +1,90 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d:cg::Starting value 0.117 +DEAL:2d:cg::Convergence step 16 value 0 +DEAL:2d:cg::Condition number estimate: 15.4 +DEAL:2d:: +DEAL:2d::Testing FE_DGQ<2>(1) +DEAL:2d:cg::Starting value 0.0625 +DEAL:2d:cg::Convergence step 1 value 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(2) +DEAL:2d:cg::Starting value 0.0686 +DEAL:2d:cg::Convergence step 38 value 0 +DEAL:2d:cg::Condition number estimate: 29.1 +DEAL:2d:: +DEAL:2d::Testing FE_DGQ<2>(2) +DEAL:2d:cg::Starting value 0.0625 +DEAL:2d:cg::Convergence step 3 value 0 +DEAL:2d:cg::Condition number estimate: 6.95 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(4) +DEAL:2d:cg::Starting value 0.0367 +DEAL:2d:cg::Convergence step 92 value 0 +DEAL:2d:cg::Condition number estimate: 148. +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(QGaussLobatto(5)) +DEAL:2d:cg::Starting value 0.0355 +DEAL:2d:cg::Convergence step 62 value 0 +DEAL:2d:cg::Condition number estimate: 74.3 +DEAL:2d:: +DEAL:2d::Testing FE_DGQ<2>(4) +DEAL:2d:cg::Starting value 0.0353 +DEAL:2d:cg::Convergence step 7 value 0 +DEAL:2d:cg::Condition number estimate: 200. +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(6) +DEAL:2d:cg::Starting value 0.0310 +DEAL:2d:cg::Convergence step 209 value 0 +DEAL:2d:cg::Condition number estimate: 1.56e+03 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(QGaussLobatto(7)) +DEAL:2d:cg::Starting value 0.0241 +DEAL:2d:cg::Convergence step 86 value 0 +DEAL:2d:cg::Condition number estimate: 140. +DEAL:2d:: +DEAL:2d::Testing FE_DGQ<2>(6) +DEAL:2d:cg::Starting value 0.0305 +DEAL:2d:cg::Convergence step 15 value 0 +DEAL:2d:cg::Condition number estimate: 1.95e+03 +DEAL:2d:: +DEAL:2d::Testing FE_DGQArbitraryNodes<2>(QGaussLobatto(7)) +DEAL:2d:cg::Starting value 0.0240 +DEAL:2d:cg::Convergence step 12 value 0 +DEAL:2d:cg::Condition number estimate: 185. +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(10) +DEAL:2d:cg::Starting value 0.173 +DEAL:2d:cg::Convergence step 1305 value 0 +DEAL:2d:cg::Condition number estimate: 2.21e+06 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(QGaussLobatto(11)) +DEAL:2d:cg::Starting value 0.0148 +DEAL:2d:cg::Convergence step 132 value 0 +DEAL:2d:cg::Condition number estimate: 328. +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(QGaussLobatto(17)) +DEAL:2d:cg::Starting value 0.00937 +DEAL:2d:cg::Convergence step 200 value 0 +DEAL:2d:cg::Condition number estimate: 759. +DEAL:2d:: +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d:cg::Starting value 0.102 +DEAL:3d:cg::Convergence step 12 value 0 +DEAL:3d:cg::Condition number estimate: 57.8 +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(2) +DEAL:3d:cg::Starting value 0.0498 +DEAL:3d:cg::Convergence step 51 value 0 +DEAL:3d:cg::Condition number estimate: 161. +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(5) +DEAL:3d:cg::Starting value 0.0121 +DEAL:3d:cg::Convergence step 475 value 0 +DEAL:3d:cg::Condition number estimate: 8.96e+03 +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(QGaussLobatto(6)) +DEAL:3d:cg::Starting value 0.0137 +DEAL:3d:cg::Convergence step 167 value 0 +DEAL:3d:cg::Condition number estimate: 1.06e+03 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_05.cc b/tests/matrix_free/matrix_vector_05.cc index 838bf0790f..481c55899d 100644 --- a/tests/matrix_free/matrix_vector_05.cc +++ b/tests/matrix_free/matrix_vector_05.cc @@ -20,15 +20,18 @@ std::ofstream logfile("matrix_vector_05/output"); template void test () { - if (dim == 3) - return; Triangulation tria; - Tensor<2,dim> points; + Point points[dim]; points[0][0] = 0.25; points[0][1] = 0.123; points[1][0] = 0.09983712334; points[1][1] = 0.314159265358979; - GridGenerator::parallelogram (tria, points); + if (dim == 3) + { + //points[2][0] = 0.21; + //points[2][2] = 0.4123; + } + GridGenerator::parallelepiped (tria, points); typename Triangulation::active_cell_iterator cell = tria.begin_active (), endc = tria.end(); @@ -41,13 +44,12 @@ void test () if (cell->center().norm()<0.2) cell->set_refine_flag(); tria.execute_coarsening_and_refinement(); - if (dim < 3 || fe_degree < 2) + if (dim < 3) tria.refine_global(2); tria.begin(tria.n_levels()-1)->set_refine_flag(); tria.last()->set_refine_flag(); tria.execute_coarsening_and_refinement(); - cell = tria.begin_active (); - for (unsigned int i=0; i<5; ++i) + for (int i=0; i<7-2*fe_degree; ++i) { cell = tria.begin_active (); unsigned int counter = 0; @@ -56,6 +58,7 @@ void test () cell->set_refine_flag(); tria.execute_coarsening_and_refinement(); } + std::cout << dim << " " << tria.n_active_cells() << std::endl; FE_Q fe (fe_degree); DoFHandler dof (tria); diff --git a/tests/matrix_free/matrix_vector_05/cmp/generic b/tests/matrix_free/matrix_vector_05/cmp/generic index f2943c08e9..c8d1cd1e9f 100644 --- a/tests/matrix_free/matrix_vector_05/cmp/generic +++ b/tests/matrix_free/matrix_vector_05/cmp/generic @@ -5,3 +5,9 @@ DEAL:2d:: DEAL:2d::Testing FE_Q<2>(2) DEAL:2d::Norm of difference: 0 DEAL:2d:: +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(2) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_common.h b/tests/matrix_free/matrix_vector_common.h index f0518dc904..63f9199c3e 100644 --- a/tests/matrix_free/matrix_vector_common.h +++ b/tests/matrix_free/matrix_vector_common.h @@ -52,9 +52,6 @@ helmholtz_operator (const MatrixFree &data, for(unsigned int cell=cell_range.first;cell