*
* This function tries to match all faces belonging to the first
* boundary with faces belonging to the second boundary with the help
- * of orthogonal_equality.
+ * of orthogonal_equality().
*
+ * The bitset that is returned together with the second face encodes the
+ * _relative_ orientation of the first face with respect to the second
+ * face, see the documentation of orthogonal_equality for further details.
+ *
* The @p offset is a vector tangential to the faces that is added to the
* location of vertices of the 'first' boundary when attempting to match
* them to the corresponding vertices of the 'second' boundary. This can
--- /dev/null
- std::ofstream logfile("work_stream_04/output");
+ // ---------------------------------------------------------------------
+ // $Id$
+ //
+ // Copyright (C) 2008 - 2013 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 WorkStream with empty functions
+
+ #include "../tests.h"
+ #include <iomanip>
+ #include <iomanip>
+ #include <fstream>
+ #include <cmath>
+
+ #include <deal.II/base/work_stream.h>
+
+
+ struct ScratchData
+ {};
+
+
+ void test ()
+ {
+ std::vector<unsigned int> v;
+ for (unsigned int i=0; i<20; ++i)
+ v.push_back (i);
+
+ WorkStream::run (v.begin(), v.end(),
+ std_cxx1x::function<void(const std::vector<unsigned int>::iterator,
+ ScratchData&,unsigned int&)>(),
+ std_cxx1x::function<void(const unsigned int&)>(),
+ ScratchData(),
+ 0U);
+ }
+
+
+
+
+ int main()
+ {
++ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test ();
+ }
--- /dev/null
- std::ofstream logfile("parameter_handler_2_read_from_string/output");
+ // ---------------------------------------------------------------------
+ // $Id$
+ //
+ // Copyright (C) 2003 - 2013 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.
+ //
+ // ---------------------------------------------------------------------
+
+
+
+ // like _02, but use the read_input_from_string function
+
+ #include "../tests.h"
+ #include <deal.II/base/logstream.h>
+ #include <deal.II/base/parameter_handler.h>
+ #include <fstream>
+
+
+ int main ()
+ {
- std::ifstream in("parameter_handler_2_read_from_string/prm");
++ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ ParameterHandler prm;
+ prm.enter_subsection ("Testing");
+ prm.declare_entry ("Function",
+ "a",
+ Patterns::List(Patterns::Selection("a|b|c|d|e|f|g|h")));
+ prm.leave_subsection ();
+
+ std::string input;
++ std::ifstream in(SOURCE_DIR "/parameter_handler_2_read_from_string/prm");
+ while (in)
+ {
+ std::string s;
+ std::getline (in, s);
+ input += s;
+ input += '\n';
+ }
+ prm.read_input_from_string(input.c_str());
+
+ std::string list;
+ prm.enter_subsection ("Testing");
+ list = prm.get ("Function");
+ prm.leave_subsection ();
+
+ deallog << list << std::endl;
+
+ return 0;
+ }
--- /dev/null
--- /dev/null
++
++DEAL::a, b, c
--- /dev/null
-std::ofstream logfile("step-10-high-order/output");
+ // ---------------------------------------------------------------------
+ // $Id$
+ //
+ // Copyright (C) 2005 - 2013 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.
+ //
+ // ---------------------------------------------------------------------
+
+
+
+ // step-10 with MappingQ(8) that shows pi with 14 digits already at the first
+ // iteration for the volume and MappingQ(20) for the boundary part that shows
+ // that things still work for very high order
+
+
+ #include "../tests.h"
+ #include <deal.II/base/logstream.h>
+ #include <fstream>
++std::ofstream logfile("output");
+
+
+ #include <deal.II/base/quadrature_lib.h>
+ #include <deal.II/base/convergence_table.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/grid/tria_boundary_lib.h>
+ #include <deal.II/grid/tria.h>
+ #include <deal.II/grid/grid_out.h>
+ #include <deal.II/hp/dof_handler.h>
+ #include <deal.II/dofs/dof_accessor.h>
+ #include <deal.II/fe/fe_q.h>
+ #include <deal.II/hp/fe_values.h>
+
+ #include <deal.II/fe/mapping_q.h>
+
+ #include <iomanip>
+ #include <fstream>
+ #include <cmath>
+
+ const long double pi = 3.141592653589793238462643;
+
+
+
+ template <int dim>
+ void compute_pi_by_area ()
+ {
+ deallog << "Computation of Pi by the area:" << std::endl
+ << "==============================" << std::endl;
+
+ const unsigned int degree = 10-dim;
+ deallog << "Degree = " << degree << std::endl;
+
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_ball (triangulation);
+
+ static const HyperBallBoundary<dim> boundary;
+ triangulation.set_boundary (0, boundary);
+
+ MappingQ<dim> mapping(degree);
+ const FE_Q<dim> dummy_fe (1);
+ const QGauss<dim> quadrature (degree+3);
+
+ DoFHandler<dim> dof_handler (triangulation);
+
+ FEValues<dim> x_fe_values (mapping, dummy_fe, quadrature,
+ update_JxW_values);
+
+ // in 3D, we obtain many digits only on a finer mesh
+ if (dim == 3)
+ triangulation.refine_global(1);
+ for (int refinement=0; refinement<4-dim;
+ ++refinement, triangulation.refine_global (1))
+ {
+ dof_handler.distribute_dofs (dummy_fe);
+
+ long double area = 0;
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ x_fe_values.reinit (cell);
+ const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+ for (unsigned int i=0; i<fe_values.n_quadrature_points; ++i)
+ area += fe_values.JxW (i);
+ };
+
+ // As a variation from step-10, no convergence table here because we
+ // should be in the regime of roundoff errors where comparing rates does
+ // not make sense as it is not necessarily stable between
+ // systems. Rather, check that we have at least 14 digits of pi correct
+ // on any level.
+ if (dim == 2)
+ {
+ deallog << "Evaluation of pi on " << triangulation.n_active_cells()
+ << " cells: " << area << std::endl;
+ // assert accuracy because numdiff might cut off digits from output
+ Assert(std::abs(area - pi) < 1e-14,
+ ExcMessage("Calculation not accurate"));
+ }
+ else
+ {
+ area *= 0.75;
+ deallog << "Evaluation of pi on in 3D " << triangulation.n_active_cells()
+ << " cells: " << area << std::endl;
+ //Assert(std::abs(area - pi) < 1e-12,
+ // ExcMessage("Calculation not accurate"));
+ }
+ };
+ deallog << std::endl;
+ }
+
+
+
+ template <int dim>
+ void compute_pi_by_perimeter ()
+ {
+ deallog << "Computation of Pi by the perimeter:" << std::endl
+ << "===================================" << std::endl;
+
+
+ const unsigned int degree = 20;
+ deallog << "Degree = " << degree << std::endl;
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_ball (triangulation);
+
+ static const HyperBallBoundary<dim> boundary;
+ triangulation.set_boundary (0, boundary);
+
+ const MappingQ<dim> mapping (degree);
+ const FE_Q<dim> fe (1);
+ const QGauss<dim-1> quadrature (degree);
+
+ DoFHandler<dim> dof_handler (triangulation);
+
+ FEFaceValues<dim> x_fe_face_values (mapping, fe, quadrature,
+ update_JxW_values);
+ for (unsigned int refinement=0; refinement<2;
+ ++refinement, triangulation.refine_global (1))
+ {
+ dof_handler.distribute_dofs (fe);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ long double perimeter = 0;
+ for (; cell!=endc; ++cell)
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ if (cell->face(face_no)->at_boundary())
+ {
+ x_fe_face_values.reinit (cell, face_no);
+ const FEFaceValues<dim> &fe_face_values
+ = x_fe_face_values.get_present_fe_values ();
+
+ for (unsigned int i=0; i<fe_face_values.n_quadrature_points; ++i)
+ perimeter += fe_face_values.JxW (i);
+ };
+ deallog << "Evaluation of pi on " << triangulation.n_active_cells()
+ << " cells: " << perimeter/2. << std::endl;
+ Assert(std::abs(perimeter/2. - pi) < 1e-14,
+ ExcMessage("Calculation not accurate"));
+ };
+
+ deallog << std::endl;
+ }
+
+
+ int main ()
+ {
+ deallog << std::setprecision(16);
+ logfile << std::setprecision(16);
+
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ compute_pi_by_area<2> ();
+ compute_pi_by_perimeter<2> ();
+
+ compute_pi_by_area<3> ();
+
+ return 0;
+ }
--- /dev/null
- std::ofstream logfile ("normal_flux_01/output");
+ // ---------------------------------------------------------------------
+ // $Id$
+ //
+ // Copyright (C) 2007 - 2013 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.
+ //
+ // ---------------------------------------------------------------------
+
+
+
+ // check the creation of normal flux boundary conditions for a finite
+ // element that consists of only a single set of vector components
+ // (i.e. it has dim components). Similar as the no-flux test in no_flux_01.cc
+
+ #include "../tests.h"
+ #include <deal.II/base/logstream.h>
+ #include <deal.II/base/function.h>
+ #include <deal.II/base/quadrature_lib.h>
+ #include <deal.II/lac/vector.h>
+ #include <deal.II/grid/grid_generator.h>
+ #include <deal.II/dofs/dof_handler.h>
+ #include <deal.II/lac/constraint_matrix.h>
+ #include <deal.II/fe/fe_q.h>
+ #include <deal.II/fe/fe_system.h>
+ #include <deal.II/fe/mapping_q1.h>
+ #include <deal.II/numerics/vector_tools.h>
+
+ #include <fstream>
+
+
+ template<int dim>
+ void test (const Triangulation<dim> &tr,
+ const FiniteElement<dim> &fe)
+ {
+ DoFHandler<dim> dof(tr);
+ dof.distribute_dofs(fe);
+
+ for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+ {
+ deallog << "FE=" << fe.get_name()
+ << ", case=" << i
+ << std::endl;
+
+ std::set<types::boundary_id> boundary_ids;
+ for (unsigned int j=0; j<=i; ++j)
+ boundary_ids.insert (j);
+
+ ConstraintMatrix cm;
+ VectorTools::compute_normal_flux_constraints (dof, 0, boundary_ids, cm);
+
+ cm.print (deallog.get_file_stream ());
+ }
+ }
+
+
+ template<int dim>
+ void test_hyper_cube()
+ {
+ Triangulation<dim> tr;
+ GridGenerator::hyper_cube(tr);
+
+ for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+ tr.begin_active()->face(i)->set_boundary_indicator (i);
+
+ tr.refine_global(2);
+
+ for (unsigned int degree=1; degree<4; ++degree)
+ {
+ FESystem<dim> fe (FE_Q<dim>(degree), dim);
+ test(tr, fe);
+ }
+ }
+
+
+ int main()
+ {
++ std::ofstream logfile ("output");
+ deallog << std::setprecision (2);
+ deallog << std::fixed;
+ deallog.attach(logfile);
+ deallog.depth_console (0);
+ deallog.threshold_double(1.e-12);
+
+ test_hyper_cube<2>();
+ test_hyper_cube<3>();
+ }
--- /dev/null
- std::ofstream logfile ("normal_flux_02/output");
+ // ---------------------------------------------------------------------
+ // $Id$
+ //
+ // Copyright (C) 2007 - 2013 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.
+ //
+ // ---------------------------------------------------------------------
+
+
+ // check the creation of normal flux boundary conditions for a finite
+ // element that consists of more than dim components and where
+ // therefore we have to pick the vector components from somewhere in
+ // the middle (opposite constraints as no_flux_02.cc does).
+
+
+ #include "../tests.h"
+ #include <deal.II/base/logstream.h>
+ #include <deal.II/base/function.h>
+ #include <deal.II/base/quadrature_lib.h>
+ #include <deal.II/lac/vector.h>
+ #include <deal.II/grid/grid_generator.h>
+ #include <deal.II/dofs/dof_handler.h>
+ #include <deal.II/lac/constraint_matrix.h>
+ #include <deal.II/fe/fe_q.h>
+ #include <deal.II/fe/fe_system.h>
+ #include <deal.II/fe/mapping_q1.h>
+ #include <deal.II/numerics/vector_tools.h>
+
+ #include <fstream>
+
+
+
+ template<int dim>
+ void test (const Triangulation<dim> &tr,
+ const FiniteElement<dim> &fe)
+ {
+ DoFHandler<dim> dof(tr);
+ dof.distribute_dofs(fe);
+
+ for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+ {
+ deallog << "FE=" << fe.get_name()
+ << ", case=" << i
+ << std::endl;
+
+ std::set<types::boundary_id> boundary_ids;
+ for (unsigned int j=0; j<=i; ++j)
+ boundary_ids.insert (j);
+
+ ConstraintMatrix cm;
+ VectorTools::compute_normal_flux_constraints (dof, 1, boundary_ids, cm);
+
+ cm.print (deallog.get_file_stream ());
+ }
+ }
+
+
+ template<int dim>
+ void test_hyper_cube()
+ {
+ Triangulation<dim> tr;
+ GridGenerator::hyper_cube(tr);
+
+ for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+ tr.begin_active()->face(i)->set_boundary_indicator (i);
+
+ tr.refine_global(2);
+
+ for (unsigned int degree=1; degree<4; ++degree)
+ {
+ FESystem<dim> fe (FE_Q<dim>(degree+1), 1,
+ FE_Q<dim>(degree), dim,
+ FE_Q<dim>(degree+1), 1);
+ test(tr, fe);
+ }
+ }
+
+
+ int main()
+ {
++ std::ofstream logfile ("output");
+ deallog << std::setprecision (2);
+ deallog << std::fixed;
+ deallog.attach(logfile);
+ deallog.depth_console (0);
+ deallog.threshold_double(1.e-12);
+
+ test_hyper_cube<2>();
+ test_hyper_cube<3>();
+ }
--- /dev/null
- std::ofstream logfile ("normal_flux_03/output");
+ // ---------------------------------------------------------------------
+ // $Id$
+ //
+ // Copyright (C) 2007 - 2013 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.
+ //
+ // ---------------------------------------------------------------------
+
+
+ // like normal_flux_01 but check on a hyper_sphere geometry
+
+ #include "../tests.h"
+ #include <deal.II/base/logstream.h>
+ #include <deal.II/base/function.h>
+ #include <deal.II/base/quadrature_lib.h>
+ #include <deal.II/lac/vector.h>
+ #include <deal.II/grid/grid_generator.h>
+ #include <deal.II/grid/tria_boundary_lib.h>
+ #include <deal.II/dofs/dof_handler.h>
+ #include <deal.II/lac/constraint_matrix.h>
+ #include <deal.II/fe/fe_q.h>
+ #include <deal.II/fe/fe_system.h>
+ #include <deal.II/fe/mapping_q1.h>
+ #include <deal.II/numerics/vector_tools.h>
+
+ #include <fstream>
+
+
+ template<int dim>
+ void test (const Triangulation<dim> &tr,
+ const FiniteElement<dim> &fe)
+ {
+ DoFHandler<dim> dof(tr);
+ dof.distribute_dofs(fe);
+
+ deallog << "FE=" << fe.get_name()
+ << std::endl;
+
+ std::set<types::boundary_id> boundary_ids;
+ boundary_ids.insert (0);
+
+ ConstraintMatrix cm;
+ VectorTools::compute_normal_flux_constraints (dof, 0, boundary_ids, cm);
+
+ cm.print (deallog.get_file_stream ());
+ }
+
+
+
+ template<int dim>
+ void test_hyper_sphere()
+ {
+ Triangulation<dim> tr;
+ GridGenerator::hyper_ball(tr);
+
+ static const HyperBallBoundary<dim> boundary;
+ tr.set_boundary (0, boundary);
+
+ tr.refine_global(2);
+
+ for (unsigned int degree=1; degree<4; ++degree)
+ {
+ FESystem<dim> fe (FE_Q<dim>(degree), dim);
+ test(tr, fe);
+ }
+ }
+
+
+ int main()
+ {
++ std::ofstream logfile ("output");
+ deallog << std::setprecision (2);
+
+ deallog.attach(logfile);
+ deallog.depth_console (0);
+ deallog.threshold_double(1.e-12);
+
+ test_hyper_sphere<2>();
+ test_hyper_sphere<3>();
+ }
--- /dev/null
- std::ofstream logfile ("normal_flux_hp_01/output");
+ // ---------------------------------------------------------------------
+ // $Id$
+ //
+ // Copyright (C) 2013 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.
+ //
+ // ---------------------------------------------------------------------
+
+
+ #include "../tests.h"
+ #include <deal.II/base/logstream.h>
+ #include <deal.II/base/function.h>
+ #include <deal.II/base/quadrature_lib.h>
+ #include <deal.II/lac/vector.h>
+ #include <deal.II/grid/grid_generator.h>
+ #include <deal.II/hp/dof_handler.h>
+ #include <deal.II/lac/constraint_matrix.h>
+ #include <deal.II/fe/fe_q.h>
+ #include <deal.II/fe/fe_system.h>
+ #include <deal.II/fe/mapping_q1.h>
+ #include <deal.II/hp/fe_collection.h>
+ #include <deal.II/numerics/vector_tools.h>
+
+ #include <fstream>
+
+
+ template<int dim>
+ void test (const Triangulation<dim> &tr,
+ const hp::FECollection<dim> &fe)
+ {
+ hp::DoFHandler<dim> dof(tr);
+ for (typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof.begin_active(); cell != dof.end(); ++cell)
+ cell->set_active_fe_index (cell->index() % 2);
+ dof.distribute_dofs(fe);
+
+ for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+ {
+ deallog << "FE=" << fe[0].get_name()
+ << ", case=" << i
+ << std::endl;
+
+ std::set<types::boundary_id> boundary_ids;
+ for (unsigned int j=0; j<=i; ++j)
+ boundary_ids.insert (j);
+
+ ConstraintMatrix cm;
+ VectorTools::compute_normal_flux_constraints (dof, 0, boundary_ids, cm);
+
+ cm.print (deallog.get_file_stream ());
+ }
+ }
+
+
+ template<int dim>
+ void test_hyper_cube()
+ {
+ Triangulation<dim> tr;
+ GridGenerator::hyper_cube(tr);
+
+ for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+ tr.begin_active()->face(i)->set_boundary_indicator (i);
+
+ tr.refine_global(2);
+
+ for (unsigned int degree=1; degree<2; ++degree)
+ {
+ hp::FECollection<dim> fe;
+ for (unsigned int deg=degree; deg<degree+2; ++deg)
+ fe.push_back (FESystem<dim> (FE_Q<dim>(deg), dim));
+ test(tr, fe);
+ }
+ }
+
+
+ int main()
+ {
++ std::ofstream logfile ("output");
+ deallog << std::setprecision (2);
+ deallog << std::fixed;
+ deallog.attach(logfile);
+ deallog.depth_console (0);
+ deallog.threshold_double(1.e-12);
+
+ test_hyper_cube<2>();
+ test_hyper_cube<3>();
+ }
--- /dev/null
- initlog(__FILE__);
+ // ---------------------------------------------------------------------
+ // $Id$
+ //
+ // Copyright (C) 1998 - 2013 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 the cell matrices generated in FETools and the local renumbering vector.
+
+ #include "../tests.h"
+ #include <iostream>
+ #include <iomanip>
+ #include <fstream>
+
+ #include <deal.II/base/logstream.h>
+
+ #include <deal.II/fe/fe_q.h>
+ #include <deal.II/fe/fe_dgq.h>
+ #include <deal.II/fe/fe_dgp.h>
+ #include <deal.II/fe/fe_raviart_thomas.h>
+ #include <deal.II/fe/fe_nedelec.h>
+ #include <deal.II/fe/fe_system.h>
+
+ #include <deal.II/fe/fe_tools.h>
+
+ void logvec (const std::vector<types::global_dof_index>& v, const std::vector<types::global_dof_index>& w)
+ {
+ deallog << '[';
+ for (unsigned int i=0; i<w.size(); ++i)
+ deallog << ' ' << w[i];
+ deallog << " ]" << std::endl << '[';;
+ for (unsigned int i=0; i<v.size(); ++i)
+ deallog << ' ' << v[i];
+ deallog << " ]" << std::endl;
+ }
+
+
+ template<int dim>
+ void test_renumbering(const FiniteElement<dim> &fe)
+ {
+ std::vector<types::global_dof_index> v(fe.dofs_per_cell);
+ std::vector<types::global_dof_index> w(fe.n_blocks());
+
+ deallog << fe.get_name() << std::endl;
+
+ FETools::compute_block_renumbering(fe, v, w, false);
+ logvec(v,w);
+ FETools::compute_block_renumbering(fe, v, w, true);
+ logvec(v,w);
+ }
+
+
+ template<int dim>
+ void test_renumbering()
+ {
+ deallog.push("Renumber");
+
+ FE_Q<dim> q1(1);
+ FE_Q<dim> q3(3);
+ FE_DGQ<dim> dg1(1);
+ FE_RaviartThomas<dim> rt1(1);
+ FE_Nedelec<dim> n1(1);
+ test_renumbering(q1);
+ test_renumbering(q3);
+ test_renumbering(dg1);
+ test_renumbering(rt1);
+ test_renumbering(n1);
+
+ FESystem<dim> q1_3(q1,3);
+ test_renumbering(q1_3);
+ FESystem<dim> q3_3(q3,3);
+ test_renumbering(q3_3);
+ FESystem<dim> q3_2_q1_3(q3, 2, q1, 3);
+ test_renumbering(q3_2_q1_3);
+ test_renumbering(FESystem<dim>(rt1,1,dg1,1));
+ test_renumbering(FESystem<dim>(rt1,1,q1,1));
+ test_renumbering(FESystem<dim>(rt1,2,q1,2));
+ test_renumbering(FESystem<dim>(n1,1,q1,1));
+
+ deallog.pop();
+ }
+
+
+ int main()
+ {
++ initlog();
+ // test_renumbering<1>();
+ test_renumbering<2>();
+ test_renumbering<3>();
+ }
--- /dev/null
- const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+ // ---------------------------------------------------------------------
+ // $Id$
+ //
+ // Copyright (C) 2012 - 2013 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.
+ //
+ // ---------------------------------------------------------------------
+
+
+ /**
+ * @file Test whether
+ */
+
+ #include "../tests.h"
+ #include <deal.II/base/function.h>
+ #include <deal.II/base/function_lib.h>
+ #include <deal.II/grid/tria.h>
+ #include <deal.II/grid/grid_generator.h>
+
+ #include <deal.II/dofs/dof_handler.h>
+
+ #include <deal.II/fe/fe_dgp.h>
+ #include <deal.II/fe/fe_q.h>
+ #include <deal.II/fe/fe_raviart_thomas.h>
+ #include <deal.II/fe/fe_system.h>
+
+ #include <deal.II/integrators/laplace.h>
+
+ #include <deal.II/meshworker/dof_info.h>
+ #include <deal.II/meshworker/integration_info.h>
+ #include <deal.II/meshworker/local_integrator.h>
+ #include <deal.II/meshworker/simple.h>
+ #include <deal.II/meshworker/loop.h>
+
+ #include <deal.II/multigrid/sparse_matrix_collection.h>
+
+ using namespace dealii;
+ using namespace LocalIntegrators;
+
+
+ template <int dim>
+ class LaplaceMatrix : public MeshWorker::LocalIntegrator<dim>
+ {
+ public:
+ LaplaceMatrix();
+ virtual void cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const;
+ virtual void boundary(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const;
+ virtual void face(MeshWorker::DoFInfo<dim>& dinfo1, MeshWorker::DoFInfo<dim>& dinfo2,
+ MeshWorker::IntegrationInfo<dim>& info1, MeshWorker::IntegrationInfo<dim>& info2) const;
+ };
+
+
+ template <int dim>
+ LaplaceMatrix<dim>::LaplaceMatrix()
+ {}
+
+
+ template <int dim>
+ void LaplaceMatrix<dim>::cell(MeshWorker::DoFInfo<dim>& dinfo, MeshWorker::IntegrationInfo<dim>& info) const
+ {
+ Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values());
+ }
+
+
+ template <int dim>
+ void LaplaceMatrix<dim>::boundary(MeshWorker::DoFInfo<dim>& dinfo,
+ typename MeshWorker::IntegrationInfo<dim>& info) const
+ {
+ const unsigned int deg = info.fe_values().get_fe().tensor_degree();
+ Laplace::nitsche_matrix(dinfo.matrix(0,false).matrix, info.fe_values(),
+ Laplace::compute_penalty(dinfo, dinfo, deg, deg));
+ }
+
+
+ template <int dim>
+ void LaplaceMatrix<dim>::face(
+ MeshWorker::DoFInfo<dim>& dinfo1, MeshWorker::DoFInfo<dim>& dinfo2,
+ MeshWorker::IntegrationInfo<dim>& info1, MeshWorker::IntegrationInfo<dim>& info2) const
+ {
+ if (info1.fe_values().get_fe().conforms(FiniteElementData<dim>::H1))
+ return;
+
+ const unsigned int deg = info1.fe_values().get_fe().tensor_degree();
+
+ if (info1.fe_values().get_fe().conforms(FiniteElementData<dim>::Hdiv) &&
+ !info1.fe_values().get_fe().conforms(FiniteElementData<dim>::Hcurl))
+ Laplace::ip_tangential_matrix(dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix,
+ dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix,
+ info1.fe_values(), info2.fe_values(),
+ Laplace::compute_penalty(dinfo1, dinfo2, deg, deg));
+ else
+ Laplace::ip_matrix(dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix,
+ dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix,
+ info1.fe_values(), info2.fe_values(),
+ Laplace::compute_penalty(dinfo1, dinfo2, deg, deg));
+ }
+
+
+ template <int dim>
+ void assemble_mg_matrix(DoFHandler<dim>& dof_handler,
+ MeshWorker::LocalIntegrator<dim>& matrix_integrator, mg::SparseMatrixCollection<double>& mg)
+ {
+ MGConstrainedDoFs mg_constraints;
+ mg_constraints.clear();
+ mg_constraints.initialize(dof_handler);
+
+ mg.set_zero();
+
+ MappingQ1<dim> mapping;
+
+ MeshWorker::IntegrationInfoBox<dim> info_box;
+ UpdateFlags update_flags = update_values | update_gradients | update_hessians;
+ info_box.add_update_flags_all(update_flags);
+ info_box.initialize(dof_handler.get_fe(), mapping);
+
+ MeshWorker::DoFInfo<dim> dof_info(dof_handler);
+
+ MeshWorker::Assembler::MGMatrixSimple<SparseMatrix<double> > assembler;
+ assembler.initialize(mg_constraints);
+ assembler.initialize(mg.matrix);
+ assembler.initialize_interfaces(mg.matrix_in, mg.matrix_out);
+ assembler.initialize_fluxes(mg.matrix_up, mg.matrix_down);
+
+ MeshWorker::integration_loop<dim, dim> (
+ dof_handler.begin_mg(), dof_handler.end_mg(),
+ dof_info, info_box, matrix_integrator, assembler);
+
+ const unsigned int nlevels = dof_handler.get_tria().n_levels();
+ for (unsigned int level=0;level<nlevels;++level)
+ {
+ for(unsigned int i=0; i<dof_handler.n_dofs(level); ++i)
+ if(mg.matrix[level].diag_element(i)==0)
+ mg.matrix[level].set(i,i,1.);
+ }
+ }
+
+
+ template <int dim>
+ void test(FiniteElement<dim> &fe)
+ {
+ deallog << fe.get_name() << std::endl;
+
+ Triangulation<dim> tr;
+ std::vector<unsigned int> repititions (dim, 1);
+ repititions[0] = 2;
+ const Point<dim> p1 = (dim == 1 ? Point<dim> (-1.) : (dim == 2 ? Point<dim>(-1.,-1.) : Point<dim> (-1.,-1.,-1.)));
+ const Point<dim> p2 = (dim == 1 ? Point<dim> (1.) : (dim == 2 ? Point<dim>(1.,1.) : Point<dim> (1.,1.,1.)));
+ GridGenerator::subdivided_hyper_rectangle(tr, repititions, p1, p2);
+ tr.begin_active()->set_refine_flag();
+ tr.execute_coarsening_and_refinement();
+
+ LaplaceMatrix<dim> matrix_integrator;
+
+ MGDoFHandler<dim> dof(tr);
+ dof.distribute_dofs(fe);
+ mg::SparseMatrixCollection<double> mg;
+ mg.resize(0,tr.n_levels()-1);
+ mg.reinit(dof);
+ assemble_mg_matrix(dof, matrix_integrator, mg);
+
+ for (unsigned int level=0;level<tr.n_levels();++level)
+ {
+ const unsigned int prec = 1;
+ const unsigned int wd = 5;
+
+ deallog << "Level " << level << std::endl << "mg" << std::endl;
+ mg.matrix[level].print_formatted(deallog.get_file_stream(), prec, false, wd, "0.");
+ if (level>0)
+ {
+ deallog << "in" << std::endl;
+ mg.matrix_in[level].print_formatted(deallog.get_file_stream(), prec, false, wd, "0.");
+ deallog << "out" << std::endl;
+ mg.matrix_out[level].print_formatted(deallog.get_file_stream(), prec, false, wd, "0.");
+ deallog << "up" << std::endl;
+ mg.matrix_up[level].print_formatted(deallog.get_file_stream(), prec, false, wd, "0.");
+ deallog << "down" << std::endl;
+ mg.matrix_down[level].print_formatted(deallog.get_file_stream(), prec, false, wd, "0.");
+ }
+ }
+ }
+
+
+ int main()
+ {
++ const std::string logname("output");
+ std::ofstream logfile(logname.c_str());
+ deallog.attach(logfile);
+ deallog.depth_console (0);
+
+ FE_DGP<2> p0(0);
+ FE_DGP<2> p1(1);
+ FE_RaviartThomas<2> rt0(0);
+ FE_RaviartThomas<2> rt1(1);
+ FE_Q<2> q2(2);
+
+ FESystem<2> sys1(p0, 2, p1, 1);
+ FESystem<2> sys2(p0, 2, rt0, 1);
+ FESystem<2> sys3(rt0, 1, p0, 2);
+ FESystem<2> sys4(p1, 2, q2, 2);
+ FESystem<2> sys5(q2, 2, p1, 2);
+
+ test(rt0);
+ test(rt1);
+ // test(sys2);
+ // test(sys3);
+ // test(sys4);
+ // test(sys5);
+ }
--- /dev/null
- initlog(__FILE__);
+ // ---------------------------------------------------------------------
+ // $Id$
+ //
+ // Copyright (C) 2000 - 2013 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.
+ //
+ // ---------------------------------------------------------------------
+
+ #include "../tests.h"
+ #include <deal.II/base/logstream.h>
+ #include <deal.II/multigrid/sparse_matrix_collection.h>
+
+ #include <fstream>
+ #include <iomanip>
+ #include <iomanip>
+ #include <algorithm>
+
+ using namespace std;
+
+
+
+ int main()
+ {
++ initlog();
+
+ mg::SparseMatrixCollection<float> smc;
+ smc.resize(0, 5);
+
+ deallog << "matrix " << smc.matrix.min_level() << '-' << smc.matrix.max_level() << std::endl;
+ deallog << "matrix_in " << smc.matrix_in.min_level() << '-' << smc.matrix_in.max_level() << std::endl;
+ deallog << "matrix_out " << smc.matrix_out.min_level() << '-' << smc.matrix_out.max_level() << std::endl;
+ deallog << "matrix_up " << smc.matrix_up.min_level() << '-' << smc.matrix_up.max_level() << std::endl;
+ deallog << "matrix_down " << smc.matrix_down.min_level() << '-' << smc.matrix_down.max_level() << std::endl;
+ }