From fa4cb5340815efde2a48884babb0a7698e1b9cb9 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Mon, 25 Aug 2014 08:48:03 -0400 Subject: [PATCH] indent with astyle --- examples/doxygen/theta_timestepping.cc | 12 +- examples/step-1/step-1.cc | 82 +- examples/step-13/step-13.cc | 82 +- examples/step-15/step-15.cc | 6 +- examples/step-16/step-16.cc | 16 +- examples/step-18/step-18.cc | 42 +- examples/step-22/step-22.cc | 22 +- examples/step-26/step-26.cc | 10 +- examples/step-3/step-3.cc | 70 +- examples/step-32/step-32.cc | 2 +- examples/step-33/step-33.cc | 24 +- examples/step-44/step-44.cc | 2 +- examples/step-51/step-51.cc | 36 +- examples/step-52/step-52.cc | 554 ++++----- examples/step-53/step-53.cc | 18 +- examples/step-6/step-6.cc | 8 +- examples/step-7/step-7.cc | 2 +- examples/step-9/step-9.cc | 10 +- include/deal.II/algorithms/any_data.h | 14 +- include/deal.II/algorithms/newton.templates.h | 18 +- include/deal.II/base/aligned_vector.h | 14 +- include/deal.II/base/data_out_base.h | 2 +- include/deal.II/base/event.h | 4 +- include/deal.II/base/function.h | 6 +- include/deal.II/base/function_lib.h | 188 +-- include/deal.II/base/function_parser.h | 466 ++++---- include/deal.II/base/iterator_range.h | 6 +- include/deal.II/base/named_data.h | 8 +- include/deal.II/base/parameter_handler.h | 4 +- include/deal.II/base/quadrature_lib.h | 4 +- include/deal.II/base/symmetric_tensor.h | 4 +- include/deal.II/base/table.h | 8 +- include/deal.II/base/tensor.h | 16 +- include/deal.II/base/tensor_base.h | 12 +- include/deal.II/base/time_stepping.h | 4 +- .../deal.II/base/time_stepping.templates.h | 18 +- include/deal.II/base/types.h | 6 +- include/deal.II/base/vectorization.h | 28 +- include/deal.II/base/work_stream.h | 390 +++---- include/deal.II/distributed/tria.h | 76 +- include/deal.II/dofs/dof_accessor.templates.h | 94 +- include/deal.II/dofs/dof_tools.h | 8 +- include/deal.II/fe/component_mask.h | 2 +- include/deal.II/fe/fe_dg_vector.templates.h | 2 +- include/deal.II/fe/fe_face.h | 4 +- include/deal.II/fe/fe_poly_face.templates.h | 70 +- include/deal.II/fe/fe_values.h | 2 +- include/deal.II/grid/grid_generator.h | 12 +- include/deal.II/grid/grid_out.h | 2 +- include/deal.II/grid/grid_tools.h | 36 +- include/deal.II/grid/manifold_lib.h | 42 +- include/deal.II/integrators/divergence.h | 24 +- include/deal.II/integrators/l2.h | 2 +- .../lac/chunk_sparse_matrix.templates.h | 2 +- .../deal.II/lac/constraint_matrix.templates.h | 2 +- include/deal.II/lac/lapack_full_matrix.h | 28 +- include/deal.II/lac/matrix_block.h | 12 +- include/deal.II/lac/petsc_full_matrix.h | 2 +- include/deal.II/lac/petsc_parallel_vector.h | 10 +- include/deal.II/lac/solver_control.h | 30 +- include/deal.II/lac/sparse_matrix.templates.h | 2 +- include/deal.II/lac/trilinos_precondition.h | 8 +- include/deal.II/lac/trilinos_sparse_matrix.h | 2 +- .../deal.II/lac/trilinos_sparsity_pattern.h | 2 +- include/deal.II/lac/trilinos_vector_base.h | 136 +-- .../deal.II/matrix_free/dof_info.templates.h | 2 +- include/deal.II/matrix_free/fe_evaluation.h | 148 +-- .../matrix_free/mapping_fe_evaluation.h | 24 +- include/deal.II/matrix_free/mapping_info.h | 2 +- include/deal.II/matrix_free/matrix_free.h | 24 +- include/deal.II/matrix_free/operators.h | 8 +- include/deal.II/matrix_free/shape_info.h | 16 +- .../matrix_free/shape_info.templates.h | 18 +- include/deal.II/meshworker/dof_info.h | 12 +- .../deal.II/meshworker/dof_info.templates.h | 2 +- include/deal.II/meshworker/functional.h | 26 +- include/deal.II/meshworker/integration_info.h | 24 +- include/deal.II/meshworker/local_integrator.h | 2 +- include/deal.II/meshworker/loop.h | 262 ++--- include/deal.II/meshworker/simple.h | 116 +- include/deal.II/meshworker/vector_selector.h | 10 +- .../meshworker/vector_selector.templates.h | 20 +- include/deal.II/multigrid/mg_transfer.h | 4 +- include/deal.II/numerics/data_out_faces.h | 2 +- include/deal.II/numerics/vector_tools.h | 14 +- .../deal.II/numerics/vector_tools.templates.h | 1016 ++++++++--------- source/algorithms/operator.cc | 6 +- source/base/data_out_base.cc | 4 +- source/base/function.cc | 76 +- source/base/function_lib.cc | 39 +- source/base/function_parser.cc | 136 +-- source/base/parameter_handler.cc | 184 +-- source/base/quadrature_lib.cc | 6 +- source/base/tensor_function.cc | 10 +- source/distributed/tria.cc | 210 ++-- source/dofs/dof_accessor_set.cc | 6 +- source/dofs/dof_handler_policy.cc | 6 +- source/dofs/dof_tools.cc | 42 +- source/dofs/dof_tools_sparsity.cc | 4 +- source/fe/fe.cc | 28 +- source/fe/fe_dgq.cc | 20 +- source/fe/fe_face.cc | 2 +- source/fe/fe_nedelec.cc | 44 +- source/fe/fe_poly_tensor.cc | 10 +- source/fe/fe_q_base.cc | 2 +- source/fe/fe_q_dg0.cc | 20 +- source/fe/fe_q_hierarchical.cc | 2 +- source/fe/fe_q_iso_q1.cc | 6 +- source/fe/fe_raviart_thomas.cc | 2 +- source/fe/fe_system.cc | 16 +- source/fe/fe_tools.cc | 102 +- source/fe/fe_values.cc | 6 +- source/grid/grid_generator.cc | 216 ++-- source/grid/grid_in.cc | 26 +- source/grid/grid_out.cc | 54 +- source/grid/grid_tools.cc | 14 +- source/grid/tria_accessor.cc | 14 +- source/hp/dof_handler.cc | 10 +- source/lac/lapack_full_matrix.cc | 14 +- source/lac/petsc_full_matrix.cc | 6 +- source/lac/petsc_parallel_vector.cc | 38 +- source/lac/trilinos_precondition.cc | 4 +- source/lac/trilinos_solver.cc | 8 +- source/lac/trilinos_sparse_matrix.cc | 2 +- source/lac/trilinos_sparsity_pattern.cc | 6 +- source/lac/trilinos_vector.cc | 2 +- source/lac/trilinos_vector_base.cc | 8 +- source/multigrid/mg_tools.cc | 12 +- source/numerics/data_out_dof_data.cc | 4 +- source/numerics/data_out_faces.cc | 2 +- source/numerics/data_out_rotation.cc | 4 +- source/numerics/derivative_approximation.cc | 944 +++++++-------- source/numerics/matrix_tools.cc | 176 +-- source/numerics/point_value_history.cc | 4 +- source/numerics/solution_transfer.cc | 10 +- 135 files changed, 3578 insertions(+), 3571 deletions(-) diff --git a/examples/doxygen/theta_timestepping.cc b/examples/doxygen/theta_timestepping.cc index caf390eba9..dc5b8fe624 100644 --- a/examples/doxygen/theta_timestepping.cc +++ b/examples/doxygen/theta_timestepping.cc @@ -47,8 +47,8 @@ class Implicit : public Operator > public: Implicit(const FullMatrix &matrix); void operator() (AnyData &out, const AnyData &in); - - private: + +private: SmartPointer, Implicit> matrix; FullMatrix m; }; @@ -93,7 +93,7 @@ Explicit::Explicit(const FullMatrix &M) void Explicit::operator() (AnyData &out, const AnyData &in) { - const double timestep = *in.read_ptr("Timestep"); + const double timestep = *in.read_ptr("Timestep"); if (this->notifications.test(Events::initial) || this->notifications.test(Events::new_timestep_size)) { m.equ(-timestep, *matrix); @@ -102,7 +102,7 @@ Explicit::operator() (AnyData &out, const AnyData &in) } this->notifications.clear(); m.vmult(*out.entry*>(0), - *in.read_ptr >("Previous iterate")); + *in.read_ptr >("Previous iterate")); } @@ -117,7 +117,7 @@ Implicit::Implicit(const FullMatrix &M) void Implicit::operator() (AnyData &out, const AnyData &in) { - const double timestep = *in.read_ptr("Timestep"); + const double timestep = *in.read_ptr("Timestep"); if (this->notifications.test(Events::initial) || this->notifications.test(Events::new_timestep_size)) { m.equ(timestep, *matrix); @@ -127,7 +127,7 @@ Implicit::operator() (AnyData &out, const AnyData &in) } this->notifications.clear(); m.vmult(*out.entry*>(0), - *in.read_ptr >("Previous time")); + *in.read_ptr >("Previous time")); } diff --git a/examples/step-1/step-1.cc b/examples/step-1/step-1.cc index 04e8fb1ca9..9abaf2acac 100644 --- a/examples/step-1/step-1.cc +++ b/examples/step-1/step-1.cc @@ -164,46 +164,46 @@ void second_grid () cell = triangulation.begin_active(), endc = triangulation.end(); for (; cell!=endc; ++cell) - { - // @note Writing a loop like this requires a lot of typing, but it - // is the only way of doing it in C++98 and C++03. However, if you - // have a C++11-compliant compiler, you can also use the C++11 - // range-based for loop style that requires significantly less - // typing. Take a look at @ref CPP11 "the deal.II C++11 page" to see - // how this works. - // - // Next, we want to loop over all vertices of the cells. Since we are - // in 2d, we know that each cell has exactly four vertices. However, - // instead of penning down a 4 in the loop bound, we make a first - // attempt at writing it in a dimension-independent way by which we - // find out about the number of vertices of a cell. Using the - // GeometryInfo class, we will later have an easier time getting the - // program to also run in 3d: we only have to change all occurrences - // of <2> to <3>, and do not - // have to audit our code for the hidden appearance of magic numbers - // like a 4 that needs to be replaced by an 8: - for (unsigned int v=0; - v < GeometryInfo<2>::vertices_per_cell; - ++v) - { - // If this cell is at the inner boundary, then at least one of its - // vertices must sit on the inner ring and therefore have a radial - // distance from the center of exactly 0.5, up to floating point - // accuracy. Compute this distance, and if we have found a vertex - // with this property flag this cell for later refinement. We can - // then also break the loop over all vertices and move on to the - // next cell. - const double distance_from_center - = center.distance (cell->vertex(v)); - - if (std::fabs(distance_from_center - inner_radius) < 1e-10) - { - cell->set_refine_flag (); - break; - } - } - } - + { + // @note Writing a loop like this requires a lot of typing, but it + // is the only way of doing it in C++98 and C++03. However, if you + // have a C++11-compliant compiler, you can also use the C++11 + // range-based for loop style that requires significantly less + // typing. Take a look at @ref CPP11 "the deal.II C++11 page" to see + // how this works. + // + // Next, we want to loop over all vertices of the cells. Since we are + // in 2d, we know that each cell has exactly four vertices. However, + // instead of penning down a 4 in the loop bound, we make a first + // attempt at writing it in a dimension-independent way by which we + // find out about the number of vertices of a cell. Using the + // GeometryInfo class, we will later have an easier time getting the + // program to also run in 3d: we only have to change all occurrences + // of <2> to <3>, and do not + // have to audit our code for the hidden appearance of magic numbers + // like a 4 that needs to be replaced by an 8: + for (unsigned int v=0; + v < GeometryInfo<2>::vertices_per_cell; + ++v) + { + // If this cell is at the inner boundary, then at least one of its + // vertices must sit on the inner ring and therefore have a radial + // distance from the center of exactly 0.5, up to floating point + // accuracy. Compute this distance, and if we have found a vertex + // with this property flag this cell for later refinement. We can + // then also break the loop over all vertices and move on to the + // next cell. + const double distance_from_center + = center.distance (cell->vertex(v)); + + if (std::fabs(distance_from_center - inner_radius) < 1e-10) + { + cell->set_refine_flag (); + break; + } + } + } + // Now that we have marked all the cells that we want refined, we let // the triangulation actually do this refinement. The function that does // so owes its long name to the fact that one can also mark cells for @@ -220,7 +220,7 @@ void second_grid () GridOut grid_out; grid_out.write_eps (triangulation, out); - std::cout << "Grid written to grid-2.eps" << std::endl; + std::cout << "Grid written to grid-2.eps" << std::endl; // At this point, all objects created in this function will be destroyed in // reverse order. Unfortunately, we defined the boundary object after the diff --git a/examples/step-13/step-13.cc b/examples/step-13/step-13.cc index 189f1f7077..b89d5fcc82 100644 --- a/examples/step-13/step-13.cc +++ b/examples/step-13/step-13.cc @@ -648,17 +648,17 @@ namespace Step13 // calls the virtual function assembling the right hand side. struct AssemblyScratchData { - AssemblyScratchData (const FiniteElement &fe, - const Quadrature &quadrature); - AssemblyScratchData (const AssemblyScratchData &scratch_data); + AssemblyScratchData (const FiniteElement &fe, + const Quadrature &quadrature); + AssemblyScratchData (const AssemblyScratchData &scratch_data); - FEValues fe_values; + FEValues fe_values; }; struct AssemblyCopyData { - FullMatrix cell_matrix; - std::vector local_dof_indices; + FullMatrix cell_matrix; + std::vector local_dof_indices; }; void @@ -666,8 +666,8 @@ namespace Step13 void local_assemble_matrix (const typename DoFHandler::active_cell_iterator &cell, - AssemblyScratchData &scratch_data, - AssemblyCopyData ©_data) const; + AssemblyScratchData &scratch_data, + AssemblyCopyData ©_data) const; void copy_local_to_global(const AssemblyCopyData ©_data, @@ -772,22 +772,22 @@ namespace Step13 Solver::assemble_linear_system (LinearSystem &linear_system) { Threads::Task<> rhs_task = Threads::new_task (&Solver::assemble_rhs, - *this, - linear_system.rhs); + *this, + linear_system.rhs); WorkStream::run(dof_handler.begin_active(), - dof_handler.end(), - std_cxx1x::bind(&Solver::local_assemble_matrix, - this, - std_cxx1x::_1, - std_cxx1x::_2, - std_cxx1x::_3), - std_cxx1x::bind(&Solver::copy_local_to_global, - this, - std_cxx1x::_1, - std_cxx1x::ref(linear_system)), - AssemblyScratchData(*fe, *quadrature), - AssemblyCopyData()); + dof_handler.end(), + std_cxx1x::bind(&Solver::local_assemble_matrix, + this, + std_cxx1x::_1, + std_cxx1x::_2, + std_cxx1x::_3), + std_cxx1x::bind(&Solver::copy_local_to_global, + this, + std_cxx1x::_1, + std_cxx1x::ref(linear_system)), + AssemblyScratchData(*fe, *quadrature), + AssemblyCopyData()); linear_system.hanging_node_constraints.condense (linear_system.matrix); // The syntax above using std_cxx1x::bind requires @@ -934,29 +934,29 @@ namespace Step13 template Solver::AssemblyScratchData:: AssemblyScratchData (const FiniteElement &fe, - const Quadrature &quadrature) - : - fe_values (fe, - quadrature, - update_gradients | update_JxW_values) + const Quadrature &quadrature) + : + fe_values (fe, + quadrature, + update_gradients | update_JxW_values) {} template Solver::AssemblyScratchData:: AssemblyScratchData (const AssemblyScratchData &scratch_data) - : - fe_values (scratch_data.fe_values.get_fe(), - scratch_data.fe_values.get_quadrature(), - update_gradients | update_JxW_values) + : + fe_values (scratch_data.fe_values.get_fe(), + scratch_data.fe_values.get_quadrature(), + update_gradients | update_JxW_values) {} template void Solver::local_assemble_matrix (const typename DoFHandler::active_cell_iterator &cell, - AssemblyScratchData &scratch_data, - AssemblyCopyData ©_data) const + AssemblyScratchData &scratch_data, + AssemblyCopyData ©_data) const { const unsigned int dofs_per_cell = fe->dofs_per_cell; const unsigned int n_q_points = quadrature->size(); @@ -971,8 +971,8 @@ namespace Step13 for (unsigned int i=0; iget_dof_indices (copy_data.local_dof_indices); } @@ -985,10 +985,10 @@ namespace Step13 LinearSystem &linear_system) const { for (unsigned int i=0; i &, - ConstraintMatrix &) + ConstraintMatrix &) = &DoFTools::make_hanging_node_constraints; // Start a side task then continue on the main thread Threads::Task<> side_task(std_cxx1x::bind(mhnc_p,std_cxx1x::cref(dof_handler), - std_cxx1x::ref(hanging_node_constraints))); + std_cxx1x::ref(hanging_node_constraints))); sparsity_pattern.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), diff --git a/examples/step-15/step-15.cc b/examples/step-15/step-15.cc index 99ff76f365..001cdb9d8f 100644 --- a/examples/step-15/step-15.cc +++ b/examples/step-15/step-15.cc @@ -550,9 +550,9 @@ namespace Step15 for (unsigned int i = 0; i < dofs_per_cell; ++i) cell_residual(i) -= (fe_values.shape_grad(i, q_point) - * coeff - * gradients[q_point] - * fe_values.JxW(q_point)); + * coeff + * gradients[q_point] + * fe_values.JxW(q_point)); } cell->get_dof_indices (local_dof_indices); diff --git a/examples/step-16/step-16.cc b/examples/step-16/step-16.cc index abe07378f7..58c246fb11 100644 --- a/examples/step-16/step-16.cc +++ b/examples/step-16/step-16.cc @@ -109,7 +109,7 @@ namespace Step16 { public: LaplaceIntegrator(); - virtual void cell(MeshWorker::DoFInfo &dinfo, MeshWorker::IntegrationInfo &info) const; + virtual void cell(MeshWorker::DoFInfo &dinfo, MeshWorker::IntegrationInfo &info) const; }; @@ -140,7 +140,7 @@ namespace Step16 // LocalIntegrators::Laplace::cell_matrix() or // LocalIntegrators::L2::L2(). Since we are assembling only a single // PDE, there is also only one of these objects with index zero. - + // In addition, we note that this integrator serves to compute the // matrices for the multilevel preconditioner as well as the matrix // and the right hand side for the global system. Since the @@ -269,12 +269,12 @@ namespace Step16 dof_handler.distribute_mg_dofs (fe); deallog << " Number of degrees of freedom: " - << dof_handler.n_dofs() - << " (by level: "; + << dof_handler.n_dofs() + << " (by level: "; for (unsigned int level=0; level filenames; + // Here we collect all filenames of the current timestep (same format as above) + std::vector filenames; for (unsigned int i=0; i > times_and_names; - times_and_names.push_back (std::pair (present_time, pvtu_master_filename)); - std::ofstream pvd_output ("solution.pvd"); - data_out.write_pvd_record (pvd_output, times_and_names); + // Finally, we write the paraview record, that references all .pvtu files and + // their respective time. Note that the variable times_and_names is declared + // static, so it will retain the entries from the pervious timesteps. + static std::vector > times_and_names; + times_and_names.push_back (std::pair (present_time, pvtu_master_filename)); + std::ofstream pvd_output ("solution.pvd"); + data_out.write_pvd_record (pvd_output, times_and_names); } - + } diff --git a/examples/step-22/step-22.cc b/examples/step-22/step-22.cc index d52c0bbde9..0088538e07 100644 --- a/examples/step-22/step-22.cc +++ b/examples/step-22/step-22.cc @@ -674,17 +674,17 @@ namespace Step22 } - // For the right-hand side we use the fact that the shape - // functions are only non-zero in one component (because our - // elements are primitive). Instead of multiplying the tensor - // representing the dim+1 values of shape function i with the - // whole right-hand side vector, we only look at the only - // non-zero component. The Function - // FiniteElement::system_to_component_index(i) will return - // which component this shape function lives in (0=x velocity, - // 1=y velocity, 2=pressure in 2d), which we use to pick out - // the correct component of the right-hand side vector to - // multiply with. + // For the right-hand side we use the fact that the shape + // functions are only non-zero in one component (because our + // elements are primitive). Instead of multiplying the tensor + // representing the dim+1 values of shape function i with the + // whole right-hand side vector, we only look at the only + // non-zero component. The Function + // FiniteElement::system_to_component_index(i) will return + // which component this shape function lives in (0=x velocity, + // 1=y velocity, 2=pressure in 2d), which we use to pick out + // the correct component of the right-hand side vector to + // multiply with. const unsigned int component_i = fe.system_to_component_index(i).first; diff --git a/examples/step-26/step-26.cc b/examples/step-26/step-26.cc index 72adfa1b85..5471d8e9de 100644 --- a/examples/step-26/step-26.cc +++ b/examples/step-26/step-26.cc @@ -134,7 +134,7 @@ namespace Step26 {} virtual double value (const Point &p, - const unsigned int component = 0) const; + const unsigned int component = 0) const; private: const double period; @@ -144,7 +144,7 @@ namespace Step26 template double RightHandSide::value (const Point &p, - const unsigned int component) const + const unsigned int component) const { Assert (component == 0, ExcInternalError()); Assert (dim == 2, ExcNotImplemented()); @@ -177,14 +177,14 @@ namespace Step26 { public: virtual double value (const Point &p, - const unsigned int component = 0) const; + const unsigned int component = 0) const; }; - + template double BoundaryValues::value (const Point &/*p*/, - const unsigned int component) const + const unsigned int component) const { Assert(component == 0, ExcInternalError()); return 0; diff --git a/examples/step-3/step-3.cc b/examples/step-3/step-3.cc index 0ca76a4654..49e7c75cbe 100644 --- a/examples/step-3/step-3.cc +++ b/examples/step-3/step-3.cc @@ -380,45 +380,45 @@ void Step3::assemble_system () // global right hand side to zero, before we fill them: cell_matrix = 0; cell_rhs = 0; - + // Now it is time to start integration over the cell, which we // do by looping over all quadrature points, which we will // number by q_index. for (unsigned int q_index=0; q_indexfe_values.shape_grad(i,q_index); this - // gradient is a 2-dimensional vector (in fact it is of type - // Tensor@<1,dim@>, with here dim=2) and the product of two - // such vectors is the scalar product, i.e. the product of - // the two shape_grad function calls is the dot - // product. This is in turn multiplied by the Jacobian - // determinant and the quadrature point weight (that one - // gets together by the call to FEValues::JxW() ). Finally, - // this is repeated for all shape functions $i$ and $j$: - for (unsigned int i=0; ife_values.shape_grad(i,q_index); this + // gradient is a 2-dimensional vector (in fact it is of type + // Tensor@<1,dim@>, with here dim=2) and the product of two + // such vectors is the scalar product, i.e. the product of + // the two shape_grad function calls is the dot + // product. This is in turn multiplied by the Jacobian + // determinant and the quadrature point weight (that one + // gets together by the call to FEValues::JxW() ). Finally, + // this is repeated for all shape functions $i$ and $j$: + for (unsigned int i=0; i (temperature[0], distributed_temperature_solution(i)); diff --git a/examples/step-33/step-33.cc b/examples/step-33/step-33.cc index 0f7d82b82d..c8b89e246e 100644 --- a/examples/step-33/step-33.cc +++ b/examples/step-33/step-33.cc @@ -2064,18 +2064,18 @@ namespace Step33 { switch (parameters.solver) { - // If the parameter file specified that a direct solver shall be used, - // then we'll get here. The process is straightforward, since deal.II - // provides a wrapper class to the Amesos direct solver within - // Trilinos. All we have to do is to create a solver control object - // (which is just a dummy object here, since we won't perform any - // iterations), and then create the direct solver object. When - // actually doing the solve, note that we don't pass a - // preconditioner. That wouldn't make much sense for a direct solver - // anyway. At the end we return the solver control statistics — - // which will tell that no iterations have been performed and that the - // final linear residual is zero, absent any better information that - // may be provided here: + // If the parameter file specified that a direct solver shall be used, + // then we'll get here. The process is straightforward, since deal.II + // provides a wrapper class to the Amesos direct solver within + // Trilinos. All we have to do is to create a solver control object + // (which is just a dummy object here, since we won't perform any + // iterations), and then create the direct solver object. When + // actually doing the solve, note that we don't pass a + // preconditioner. That wouldn't make much sense for a direct solver + // anyway. At the end we return the solver control statistics — + // which will tell that no iterations have been performed and that the + // final linear residual is zero, absent any better information that + // may be provided here: case Parameters::Solver::direct: { SolverControl solver_control (1,0); diff --git a/examples/step-44/step-44.cc b/examples/step-44/step-44.cc index 9ac78a6ba3..da8a710660 100644 --- a/examples/step-44/step-44.cc +++ b/examples/step-44/step-44.cc @@ -2728,7 +2728,7 @@ namespace Step44 // We've chosen by default a SSOR preconditioner as it appears to // provide the fastest solver convergence characteristics for this // problem on a single-thread machine. However, this might not be - // true for different problem sizes. + // true for different problem sizes. PreconditionSelector, Vector > preconditioner (parameters.preconditioner_type, parameters.preconditioner_relaxation); diff --git a/examples/step-51/step-51.cc b/examples/step-51/step-51.cc index 14be9e7875..da6e062820 100644 --- a/examples/step-51/step-51.cc +++ b/examples/step-51/step-51.cc @@ -106,27 +106,27 @@ namespace Step51 template <> const Point<1> SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] - = { Point<1>(-1.0 / 3.0), - Point<1>(0.0), - Point<1>(+1.0 / 3.0) - }; + = { Point<1>(-1.0 / 3.0), + Point<1>(0.0), + Point<1>(+1.0 / 3.0) + }; template <> const Point<2> SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] - = { Point<2>(-0.5, +0.5), - Point<2>(-0.5, -0.5), - Point<2>(+0.5, -0.5) - }; + = { Point<2>(-0.5, +0.5), + Point<2>(-0.5, -0.5), + Point<2>(+0.5, -0.5) + }; template <> const Point<3> SolutionBase<3>::source_centers[SolutionBase<3>::n_source_centers] - = { Point<3>(-0.5, +0.5, 0.25), - Point<3>(-0.6, -0.5, -0.125), - Point<3>(+0.5, -0.5, 0.5) - }; + = { Point<3>(-0.5, +0.5, 0.25), + Point<3>(-0.6, -0.5, -0.125), + Point<3>(+0.5, -0.5, 0.5) + }; template const double SolutionBase::width = 1./5.; @@ -388,7 +388,7 @@ namespace Step51 // Dirichlet boundary conditions, just as in a continuous Galerkin finite // element method. We can enforce the boundary conditions in an analogous // manner through the use of ConstrainMatrix constructs. In - // addition, hanging nodes are handled in the same way as for + // addition, hanging nodes are handled in the same way as for // continuous finite elements: For the face elements which // only define degrees of freedom on the face, this process sets the // solution on the refined to be the one from the coarse side. @@ -763,9 +763,9 @@ namespace Step51 for (unsigned int q=0; q convection - = scratch.convection_velocity.value(scratch.fe_values_local.quadrature_point(q)); + = scratch.convection_velocity.value(scratch.fe_values_local.quadrature_point(q)); const double JxW = scratch.fe_values_local.JxW(q); for (unsigned int k=0; k normal = scratch.fe_face_values.normal_vector(q); const Tensor<1,dim> convection - = scratch.convection_velocity.value(quadrature_point); + = scratch.convection_velocity.value(quadrature_point); // Here we compute the stabilization parameter discussed in the // introduction: since the diffusion is one and the diffusion @@ -903,7 +903,7 @@ namespace Step51 const unsigned int jj=scratch.fe_local_support_on_face[face][j]; scratch.ll_matrix(ii,jj) += tau_stab * scratch.u_phi[i] * scratch.u_phi[j] * JxW; } - + // When @p trace_reconstruct=true, we are solving for the local // solutions on an element by element basis. The local // right-hand-side is calculated by replacing the basis functions @p @@ -1209,7 +1209,7 @@ namespace Step51 component_interpretation (dim+1, DataComponentInterpretation::component_is_part_of_vector); component_interpretation[dim] - = DataComponentInterpretation::component_is_scalar; + = DataComponentInterpretation::component_is_scalar; data_out.add_data_vector (dof_handler_local, solution_local, names, component_interpretation); diff --git a/examples/step-52/step-52.cc b/examples/step-52/step-52.cc index 8011c4b5fc..749901a711 100644 --- a/examples/step-52/step-52.cc +++ b/examples/step-52/step-52.cc @@ -68,73 +68,73 @@ namespace Step52 // this class are not new and have been explained in previous tutorials. class Diffusion { - public: - Diffusion(); - - void run(); - - private: - void setup_system(); - - void assemble_system(); - - double get_source(double time,const Point<2> &point) const; - - // This function evaluates the diffusion equation $M^{-1}(f(t,y))$ at a given time and - // for a given y. - Vector evaluate_diffusion(const double time, const Vector &y) const; - - // Evaluate $\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1}$ or - // equivalently $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M$ at a given - // time, for a given $\tau$ and y. - Vector id_minus_tau_J_inverse(const double time, - const double tau, - const Vector &y); - - void output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const; - - // The next three functions are the driver for the explicit methods, the - // implicit methods, and the embedded explicit methods respectively. The - // driver function for embedded explicit methods returns the number of - // steps executed since this number is adapted. - void explicit_method(TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time); - - void implicit_method(TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time); - - unsigned int embedded_explicit_method(TimeStepping::runge_kutta_method method, - const unsigned int n_time_steps, - const double initial_time, - const double final_time); - - - unsigned int fe_degree; - - double diffusion_coefficient; - double absorption_xs; - - Triangulation<2> triangulation; - - FE_Q<2> fe; - - DoFHandler<2> dof_handler; - - ConstraintMatrix constraint_matrix; - - SparsityPattern sparsity_pattern; - - SparseMatrix system_matrix; - SparseMatrix mass_matrix; - SparseMatrix mass_minus_tau_Jacobian; - - SparseDirectUMFPACK inverse_mass_matrix; - - Vector solution; + public: + Diffusion(); + + void run(); + + private: + void setup_system(); + + void assemble_system(); + + double get_source(double time,const Point<2> &point) const; + + // This function evaluates the diffusion equation $M^{-1}(f(t,y))$ at a given time and + // for a given y. + Vector evaluate_diffusion(const double time, const Vector &y) const; + + // Evaluate $\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1}$ or + // equivalently $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M$ at a given + // time, for a given $\tau$ and y. + Vector id_minus_tau_J_inverse(const double time, + const double tau, + const Vector &y); + + void output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const; + + // The next three functions are the driver for the explicit methods, the + // implicit methods, and the embedded explicit methods respectively. The + // driver function for embedded explicit methods returns the number of + // steps executed since this number is adapted. + void explicit_method(TimeStepping::runge_kutta_method method, + const unsigned int n_time_steps, + const double initial_time, + const double final_time); + + void implicit_method(TimeStepping::runge_kutta_method method, + const unsigned int n_time_steps, + const double initial_time, + const double final_time); + + unsigned int embedded_explicit_method(TimeStepping::runge_kutta_method method, + const unsigned int n_time_steps, + const double initial_time, + const double final_time); + + + unsigned int fe_degree; + + double diffusion_coefficient; + double absorption_xs; + + Triangulation<2> triangulation; + + FE_Q<2> fe; + + DoFHandler<2> dof_handler; + + ConstraintMatrix constraint_matrix; + + SparsityPattern sparsity_pattern; + + SparseMatrix system_matrix; + SparseMatrix mass_matrix; + SparseMatrix mass_minus_tau_Jacobian; + + SparseDirectUMFPACK inverse_mass_matrix; + + Vector solution; }; @@ -142,11 +142,11 @@ namespace Step52 // We choose quadratic finite elements and we initialize the parameters. Diffusion::Diffusion() : - fe_degree(2), - diffusion_coefficient(1./30.), - absorption_xs(1.), - fe(fe_degree), - dof_handler(triangulation) + fe_degree(2), + diffusion_coefficient(1./30.), + absorption_xs(1.), + fe(fe_degree), + dof_handler(triangulation) {} @@ -174,8 +174,8 @@ namespace Step52 // @sect5{Diffusion::assemble_system} - // In this function, we compute - // $-\int D \nabla b_i \cdot \nabla b_j d\boldsymbol{r} - \int \Sigma_a b_i b_j d\boldsymbol{r}$ + // In this function, we compute + // $-\int D \nabla b_i \cdot \nabla b_j d\boldsymbol{r} - \int \Sigma_a b_i b_j d\boldsymbol{r}$ // and the mass matrix $\int b_i b_j d\boldsymbol{r}$. The mass matrix is then // inverted using a direct solver. void Diffusion::assemble_system() @@ -186,7 +186,7 @@ namespace Step52 const QGauss<2> quadrature_formula(fe_degree+1); FEValues<2> fe_values(fe, quadrature_formula, - update_values | update_gradients | update_JxW_values); + update_values | update_gradients | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; @@ -200,32 +200,32 @@ namespace Step52 typename DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); - - for (; cell!=endc; ++cell) - { - cell_matrix = 0.; - cell_mass_matrix = 0.; - - fe_values.reinit (cell); - - for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); - constraint_matrix.distribute_local_to_global(cell_matrix,local_dof_indices,system_matrix); - constraint_matrix.distribute_local_to_global(cell_mass_matrix,local_dof_indices,mass_matrix); - } + for (; cell!=endc; ++cell) + { + cell_matrix = 0.; + cell_mass_matrix = 0.; + + fe_values.reinit (cell); + + for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); + + constraint_matrix.distribute_local_to_global(cell_matrix,local_dof_indices,system_matrix); + constraint_matrix.distribute_local_to_global(cell_mass_matrix,local_dof_indices,mass_matrix); + } inverse_mass_matrix.initialize(mass_matrix); } @@ -246,8 +246,8 @@ namespace Step52 double source = 0.; source = intensity*(frequency*std::cos(frequency*time)*(b*x-x*x) + std::sin(frequency*time) * - (absorption_xs*(b*x-x*x)+2.*diffusion_coefficient)); - + (absorption_xs*(b*x-x*x)+2.*diffusion_coefficient)); + return source; } @@ -265,7 +265,7 @@ namespace Step52 const QGauss<2> quadrature_formula(fe_degree+1); FEValues<2> fe_values(fe, quadrature_formula, - update_values | update_quadrature_points | update_JxW_values); + update_values | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; @@ -278,26 +278,26 @@ namespace Step52 typename DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); - + for (; cell!=endc; ++cell) - { - cell_source = 0.; + { + cell_source = 0.; - fe_values.reinit (cell); + fe_values.reinit (cell); - for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); + cell->get_dof_indices(local_dof_indices); + + constraint_matrix.distribute_local_to_global(cell_source,local_dof_indices,tmp); + } - constraint_matrix.distribute_local_to_global(cell_source,local_dof_indices,tmp); - } - Vector value(dof_handler.n_dofs()); inverse_mass_matrix.vmult(value,tmp); @@ -312,10 +312,10 @@ namespace Step52 // is done in several steps: // - compute $M-\tau \frac{\partial f}{\partial y}$ // - inverse the matrix to get $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1}$ - // - compute $tmp=My$ - // - compute $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} tmp = \left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} My$. - Vector Diffusion::id_minus_tau_J_inverse(const double time, const double tau, - const Vector &y) + // - compute $tmp=My$ + // - compute $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} tmp = \left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} My$. + Vector Diffusion::id_minus_tau_J_inverse(const double time, const double tau, + const Vector &y) { Vector tmp(dof_handler.n_dofs()); Vector result(y); @@ -323,7 +323,7 @@ namespace Step52 mass_minus_tau_Jacobian.copy_from(mass_matrix); mass_minus_tau_Jacobian.add(-tau,system_matrix); - + inverse_mass_minus_tau_Jacobian.initialize(mass_minus_tau_Jacobian); mass_matrix.vmult(tmp,y); @@ -343,67 +343,67 @@ namespace Step52 std::string method_name; switch (method) - { + { case TimeStepping::FORWARD_EULER : - { - method_name = "forward_euler"; - break; - } + { + method_name = "forward_euler"; + break; + } case TimeStepping::RK_THIRD_ORDER : - { - method_name = "rk3"; - break; - } + { + method_name = "rk3"; + break; + } case TimeStepping::RK_CLASSIC_FOURTH_ORDER : - { - method_name = "rk4"; - break; - } + { + method_name = "rk4"; + break; + } case TimeStepping::BACKWARD_EULER : - { - method_name = "backward_euler"; - break; - } + { + method_name = "backward_euler"; + break; + } case TimeStepping::IMPLICIT_MIDPOINT : - { - method_name = "implicit_midpoint"; - break; - } + { + method_name = "implicit_midpoint"; + break; + } case TimeStepping::SDIRK_TWO_STAGES : - { - method_name = "sdirk"; - break; - } + { + method_name = "sdirk"; + break; + } case TimeStepping::HEUN_EULER : - { - method_name = "heun_euler"; - break; - } + { + method_name = "heun_euler"; + break; + } case TimeStepping::BOGACKI_SHAMPINE : - { - method_name = "bocacki_shampine"; - break; - } + { + method_name = "bocacki_shampine"; + break; + } case TimeStepping::DOPRI : - { - method_name = "dopri"; - break; - } + { + method_name = "dopri"; + break; + } case TimeStepping::FEHLBERG : - { - method_name = "fehlberg"; - break; - } + { + method_name = "fehlberg"; + break; + } case TimeStepping::CASH_KARP : - { - method_name = "cash_karp"; - break; - } + { + method_name = "cash_karp"; + break; + } default : - { - break; - } - } + { + break; + } + } DataOut<2> data_out; @@ -438,21 +438,21 @@ namespace Step52 TimeStepping::ExplicitRungeKutta > explicit_runge_kutta(method); output_results(0,method); for (unsigned int i=0; iDiffusion::implicit_method} // This function is equivalent to explicit_method but for implicit methods. - // When using implicit methods, we need to evaluate $M^{-1}(f(t,y))$ and + // When using implicit methods, we need to evaluate $M^{-1}(f(t,y))$ and // $\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1}$. void Diffusion::implicit_method(TimeStepping::runge_kutta_method method, const unsigned int n_time_steps, @@ -466,16 +466,16 @@ namespace Step52 TimeStepping::ImplicitRungeKutta > implicit_runge_kutta(method); output_results(0,method); for (unsigned int i=0; ifinal_time) - time_step = final_time-time; + { + // We choose the last time step such that the final time is exactly + // reached. + if (time+time_step>final_time) + time_step = final_time-time; - time = embedded_explicit_runge_kutta.evolve_one_time_step( - std_cxx1x::bind(&Diffusion::evaluate_diffusion,this,std_cxx1x::_1,std_cxx1x::_2), - time,time_step,solution); + time = embedded_explicit_runge_kutta.evolve_one_time_step( + std_cxx1x::bind(&Diffusion::evaluate_diffusion,this,std_cxx1x::_1,std_cxx1x::_2), + time,time_step,solution); - if ((n_steps+1)%10==0) - output_results(n_steps+1,method); + if ((n_steps+1)%10==0) + output_results(n_steps+1,method); - time_step = embedded_explicit_runge_kutta.get_status().delta_t_guess; - ++n_steps; - } + time_step = embedded_explicit_runge_kutta.get_status().delta_t_guess; + ++n_steps; + } return n_steps; } @@ -549,16 +549,16 @@ namespace Step52 typename Triangulation<2>::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end(); - + for (; cell!=endc; ++cell) for (unsigned int f=0; f::faces_per_cell; ++f) if (cell->face(f)->at_boundary()) - { - if ((cell->face(f)->center()[0]==0.) || (cell->face(f)->center()[0]==5.)) - cell->face(f)->set_boundary_indicator(1); - else - cell->face(f)->set_boundary_indicator(0); - } + { + if ((cell->face(f)->center()[0]==0.) || (cell->face(f)->center()[0]==5.)) + cell->face(f)->set_boundary_indicator(1); + else + cell->face(f)->set_boundary_indicator(0); + } setup_system(); @@ -572,90 +572,90 @@ namespace Step52 // Next, we solve the diffusion problem using different Runge-Kutta methods. std::cout << "Explicit methods:" << std::endl; explicit_method (TimeStepping::FORWARD_EULER, - n_time_steps, - initial_time, - final_time); + n_time_steps, + initial_time, + final_time); std::cout << "Forward Euler: error=" << solution.l2_norm() << std::endl; - + explicit_method (TimeStepping::RK_THIRD_ORDER, - n_time_steps, - initial_time, - final_time); + n_time_steps, + initial_time, + final_time); std::cout << "Third order Runge-Kutta: error=" << solution.l2_norm() << std::endl; explicit_method (TimeStepping::RK_CLASSIC_FOURTH_ORDER, - n_time_steps, - initial_time, - final_time); + n_time_steps, + initial_time, + final_time); std::cout << "Fourth order Runge-Kutta: error=" << solution.l2_norm() << std::endl; std::cout << std::endl; std::cout << "Implicit methods:" << std::endl; implicit_method (TimeStepping::BACKWARD_EULER, - n_time_steps, - initial_time, - final_time); + n_time_steps, + initial_time, + final_time); std::cout << "Backward Euler: error=" << solution.l2_norm() << std::endl; implicit_method (TimeStepping::IMPLICIT_MIDPOINT, - n_time_steps, - initial_time, - final_time); + n_time_steps, + initial_time, + final_time); std::cout << "Implicit Midpoint: error=" << solution.l2_norm() << std::endl; implicit_method (TimeStepping::CRANK_NICOLSON, - n_time_steps, - initial_time, - final_time); + n_time_steps, + initial_time, + final_time); std::cout << "Crank-Nicolson: error=" << solution.l2_norm() << std::endl; implicit_method (TimeStepping::SDIRK_TWO_STAGES, - n_time_steps, - initial_time, - final_time); + n_time_steps, + initial_time, + final_time); std::cout << "SDIRK: error=" << solution.l2_norm() << std::endl; std::cout << std::endl; - + std::cout << "Embedded explicit methods:" << std::endl; n_steps = embedded_explicit_method (TimeStepping::HEUN_EULER, - n_time_steps, - initial_time, - final_time); + n_time_steps, + initial_time, + final_time); std::cout << "Heun-Euler: error=" << solution.l2_norm() << std::endl; std::cout << " steps performed=" << n_steps << std::endl; - + n_steps = embedded_explicit_method (TimeStepping::BOGACKI_SHAMPINE, - n_time_steps, - initial_time, - final_time); + n_time_steps, + initial_time, + final_time); std::cout << "Bogacki-Shampine: error=" << solution.l2_norm() << std::endl; std::cout << " steps performed=" << n_steps << std::endl; n_steps = embedded_explicit_method (TimeStepping::DOPRI, - n_time_steps, - initial_time, - final_time); + n_time_steps, + initial_time, + final_time); std::cout << "Dopri: error=" << solution.l2_norm() << std::endl; std::cout << " steps performed=" << n_steps << std::endl; n_steps = embedded_explicit_method (TimeStepping::FEHLBERG, - n_time_steps, - initial_time, - final_time); + n_time_steps, + initial_time, + final_time); std::cout << "Fehlberg: error=" << solution.l2_norm() << std::endl; std::cout << " steps performed=" << n_steps << std::endl; - + n_steps = embedded_explicit_method (TimeStepping::CASH_KARP, - n_time_steps, - initial_time, - final_time); + n_time_steps, + initial_time, + final_time); std::cout << "Cash-Karp: error=" << solution.l2_norm() << std::endl; std::cout << " steps performed=" << n_steps << std::endl; } } - + // @sect3{The main() function} @@ -665,33 +665,33 @@ namespace Step52 int main () { try - { - Step52::Diffusion diffusion; - diffusion.run(); - } + { + Step52::Diffusion diffusion; + diffusion.run(); + } catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - }; + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; return 0; } diff --git a/examples/step-53/step-53.cc b/examples/step-53/step-53.cc index d30105ceae..e22d94c797 100644 --- a/examples/step-53/step-53.cc +++ b/examples/step-53/step-53.cc @@ -36,7 +36,7 @@ template class SphereGeometry : public Boundary { public: - SphereGeometry (const Point ¢er); + SphereGeometry (const Point ¢er); virtual Point get_new_point_on_line (const typename Triangulation::line_iterator &line) const; @@ -58,8 +58,8 @@ private: template SphereGeometry::SphereGeometry (const Point ¢er) -: -center (center) + : + center (center) {} @@ -100,7 +100,7 @@ SphereGeometry<3>::pull_back (const Point<3> &p) const const double r = relative_point.norm(); const double phi = std::atan2 (relative_point[1], relative_point[0]); const double theta = std::atan2 (relative_point[2], std::sqrt (relative_point[0]*relative_point[0] + - relative_point[1]*relative_point[1])); + relative_point[1]*relative_point[1])); std_cxx1x::array result; result[0] = r; @@ -238,7 +238,8 @@ SphereGeometry:: get_new_point_on_line (const typename Triangulation::line_iterator &line) const { std_cxx1x::array preimages[2] = { pull_back (line->vertex(0)), - pull_back (line->vertex(1)) }; + pull_back (line->vertex(1)) + }; return push_forward(average (preimages)); } @@ -252,7 +253,8 @@ get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) c std_cxx1x::array preimages[4] = { pull_back (quad->vertex(0)), pull_back (quad->vertex(1)), pull_back (quad->vertex(2)), - pull_back (quad->vertex(3)) }; + pull_back (quad->vertex(3)) + }; return push_forward(average(preimages)); } @@ -275,7 +277,7 @@ void make_grid () triangulation.refine_global(1); for (typename Triangulation::active_cell_iterator cell=triangulation.begin_active(); - cell!=triangulation.end(); ++cell) + cell!=triangulation.end(); ++cell) { if (cell->center().distance(center)< radius) cell->set_manifold_id(1); @@ -295,7 +297,7 @@ void make_grid () } - + // @sect3{The main function} diff --git a/examples/step-6/step-6.cc b/examples/step-6/step-6.cc index 07901f9bcb..b17903501e 100644 --- a/examples/step-6/step-6.cc +++ b/examples/step-6/step-6.cc @@ -444,10 +444,10 @@ void Step6::assemble_system () // @p cell_rhs into the global objects. cell->get_dof_indices (local_dof_indices); constraints.distribute_local_to_global (cell_matrix, - cell_rhs, - local_dof_indices, - system_matrix, - system_rhs); + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); } // Now we are done assembling the linear system. The constraint matrix took // care of applying the boundary conditions and also eliminated hanging node diff --git a/examples/step-7/step-7.cc b/examples/step-7/step-7.cc index c5c3ac2320..1139cae4a7 100644 --- a/examples/step-7/step-7.cc +++ b/examples/step-7/step-7.cc @@ -976,7 +976,7 @@ namespace Step7 template void HelmholtzProblem::run () { - const unsigned int n_cycles = (refinement_mode==global_refinement)?5:9; + const unsigned int n_cycles = (refinement_mode==global_refinement)?5:9; for (unsigned int cycle=0; cycle GradientEstimation::EstimateScratchData:: EstimateScratchData (const FiniteElement &fe, - const Vector &solution) + const Vector &solution) : fe_midpoint_value(fe, QMidpoint (), @@ -1033,8 +1033,8 @@ namespace Step9 SynchronousIterators begin_sync_it (IteratorTuple (dof_handler.begin_active(), error_per_cell.begin())), - end_sync_it (IteratorTuple (dof_handler.end(), - error_per_cell.end())); + end_sync_it (IteratorTuple (dof_handler.end(), + error_per_cell.end())); WorkStream::run (begin_sync_it, end_sync_it, @@ -1095,9 +1095,9 @@ namespace Step9 template void GradientEstimation::estimate_cell (const SynchronousIterators::active_cell_iterator, - Vector::iterator> > &cell, + Vector::iterator> > &cell, EstimateScratchData &scratch_data, - const EstimateCopyData &) + const EstimateCopyData &) { // We need space for the tensor Y, which is the sum of // outer products of the y-vectors. diff --git a/include/deal.II/algorithms/any_data.h b/include/deal.II/algorithms/any_data.h index 638b4ab577..065498dd66 100644 --- a/include/deal.II/algorithms/any_data.h +++ b/include/deal.II/algorithms/any_data.h @@ -176,8 +176,8 @@ public: /// List the contents to a stream template - void list (STREAM& os) const; - + void list (STREAM &os) const; + /// Conversion from old NamedData template AnyData(const NamedData &); @@ -472,14 +472,14 @@ AnyData::merge(const AnyData &other) template inline -void AnyData::list(STREAM& os) const +void AnyData::list(STREAM &os) const { - for (unsigned int i=0;i out; - VECTOR *p = &u; - out.add(p, "solution"); - p = Du; - out.add(p, "update"); - p = res; - out.add(p, "residual"); - *data_out << step; - *data_out << out; + NamedData out; + VECTOR *p = &u; + out.add(p, "solution"); + p = Du; + out.add(p, "update"); + p = res; + out.add(p, "residual"); + *data_out << step; + *data_out << out; } while (control.check(step++, resnorm) == SolverControl::iterate) diff --git a/include/deal.II/base/aligned_vector.h b/include/deal.II/base/aligned_vector.h index 0c7539838c..77ec5a95e1 100644 --- a/include/deal.II/base/aligned_vector.h +++ b/include/deal.II/base/aligned_vector.h @@ -327,7 +327,7 @@ namespace internal // never arrive here because they are non-trivial). if (std_cxx1x::is_trivial::value == true) - std::memcpy ((void*)(destination_+begin), source_+begin, + std::memcpy ((void *)(destination_+begin), source_+begin, (end-begin)*sizeof(T)); else if (copy_only_ == false) for (std::size_t i=begin; i::value == true && trivial_element) - std::memset ((void*)(destination_+begin), 0, (end-begin)*sizeof(T)); + std::memset ((void *)(destination_+begin), 0, (end-begin)*sizeof(T)); else // initialize memory and set for (std::size_t i=begin; i::AlignedVector (const AlignedVector &vec) template < class T > inline -AlignedVector& +AlignedVector & AlignedVector::operator = (const AlignedVector &vec) { resize(0); @@ -558,7 +558,7 @@ AlignedVector::reserve (const size_type size_alloc) if (_end_data != _data) { dealii::internal::AlignedVectorMove(new_data, new_data + old_size, - _data); + _data); #if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 _mm_free(new_data); #else @@ -807,7 +807,7 @@ typename AlignedVector::size_type AlignedVector::memory_consumption () const { size_type memory = sizeof(*this); - for (const T* t = _data ; t != _end_data; ++t) + for (const T *t = _data ; t != _end_data; ++t) memory += dealii::MemoryConsumption::memory_consumption(*t); memory += sizeof(T) * (_end_allocated-_end_data); return memory; @@ -829,7 +829,7 @@ bool operator == (const AlignedVector &lhs, if (lhs.size() != rhs.size()) return false; for (typename AlignedVector::const_iterator lit = lhs.begin(), - rit = rhs.begin(); lit != lhs.end(); ++lit, ++rit) + rit = rhs.begin(); lit != lhs.end(); ++lit, ++rit) if (*lit != *rit) return false; return true; diff --git a/include/deal.II/base/data_out_base.h b/include/deal.II/base/data_out_base.h index 4899c5f074..bd5721e7a8 100644 --- a/include/deal.II/base/data_out_base.h +++ b/include/deal.II/base/data_out_base.h @@ -1531,7 +1531,7 @@ namespace DataOutBase const EpsFlags &flags, std::ostream &out); - + /** * Write the given list of patches to the output stream in * GMV format. diff --git a/include/deal.II/base/event.h b/include/deal.II/base/event.h index 1cf7b27e73..b77c9450d7 100644 --- a/include/deal.II/base/event.h +++ b/include/deal.II/base/event.h @@ -151,12 +151,12 @@ namespace Algorithms * new. */ extern const Event initial; - + /** * The mesh has changed. */ extern const Event remesh; - + /** * The current derivative leads * to slow convergence of diff --git a/include/deal.II/base/function.h b/include/deal.II/base/function.h index ace10a1c22..e3c653fc4b 100644 --- a/include/deal.II/base/function.h +++ b/include/deal.II/base/function.h @@ -115,7 +115,7 @@ template class TensorFunction; */ template class Function : public FunctionTime, - public Subscriptor + public Subscriptor { public: /** @@ -231,7 +231,7 @@ public: * given point. */ virtual Tensor<1,dim, Number> gradient (const Point &p, - const unsigned int component = 0) const; + const unsigned int component = 0) const; /** * Return the gradient of all components of the function at the given @@ -353,7 +353,7 @@ public: std::vector > &values) const; virtual Tensor<1,dim, Number> gradient (const Point &p, - const unsigned int component = 0) const; + const unsigned int component = 0) const; virtual void vector_gradient (const Point &p, std::vector > &gradients) const; diff --git a/include/deal.II/base/function_lib.h b/include/deal.II/base/function_lib.h index 19f3824c64..e75bc47d72 100644 --- a/include/deal.II/base/function_lib.h +++ b/include/deal.II/base/function_lib.h @@ -1213,50 +1213,50 @@ namespace Functions template class InterpolatedTensorProductGridData : public Function { - public: - /** - * Constructor. - * @param coordinate_values An array of dim arrays. Each of the inner - * arrays contains the coordinate values $x_0,\ldotx, x_{K-1}$ and - * similarly for the other coordinate directions. These arrays - * need not have the same size. Obviously, we need dim such arrays - * for a dim-dimensional function object. The coordinate values - * within this array are assumed to be strictly ascending to allow - * for efficient lookup. - * @param data_values A dim-dimensional table of data at each of the - * mesh points defined by the coordinate arrays above. Note that the - * Table class has a number of conversion constructors that allow - * converting other data types into a table where you specify this - * argument. - */ - InterpolatedTensorProductGridData (const std_cxx1x::array,dim> &coordinate_values, - const Table &data_values); - - /** - * Compute the value of the function set by bilinear interpolation of the - * given data set. - * - * @param p The point at which the function is to be evaluated. - * @param component The vector component. Since this function is scalar, - * only zero is a valid argument here. - * @return The interpolated value at this point. If the point lies outside - * the set of coordinates, the function is extended by a constant. - */ - virtual - double - value (const Point &p, - const unsigned int component = 0) const; - - private: - /** - * The set of coordinate values in each of the coordinate directions. - */ - const std_cxx1x::array,dim> coordinate_values; - - /** - * The data that is to be interpolated. - */ - const Table data_values; + public: + /** + * Constructor. + * @param coordinate_values An array of dim arrays. Each of the inner + * arrays contains the coordinate values $x_0,\ldotx, x_{K-1}$ and + * similarly for the other coordinate directions. These arrays + * need not have the same size. Obviously, we need dim such arrays + * for a dim-dimensional function object. The coordinate values + * within this array are assumed to be strictly ascending to allow + * for efficient lookup. + * @param data_values A dim-dimensional table of data at each of the + * mesh points defined by the coordinate arrays above. Note that the + * Table class has a number of conversion constructors that allow + * converting other data types into a table where you specify this + * argument. + */ + InterpolatedTensorProductGridData (const std_cxx1x::array,dim> &coordinate_values, + const Table &data_values); + + /** + * Compute the value of the function set by bilinear interpolation of the + * given data set. + * + * @param p The point at which the function is to be evaluated. + * @param component The vector component. Since this function is scalar, + * only zero is a valid argument here. + * @return The interpolated value at this point. If the point lies outside + * the set of coordinates, the function is extended by a constant. + */ + virtual + double + value (const Point &p, + const unsigned int component = 0) const; + + private: + /** + * The set of coordinate values in each of the coordinate directions. + */ + const std_cxx1x::array,dim> coordinate_values; + + /** + * The data that is to be interpolated. + */ + const Table data_values; }; @@ -1296,56 +1296,56 @@ namespace Functions template class InterpolatedUniformGridData : public Function { - public: - /** - * Constructor - * @param interval_endpoints The left and right end points of the (uniformly - * subdivided) intervals in each of the coordinate directions. - * @param n_subdivisions The number of subintervals of the subintervals - * in each coordinate direction. A value of one for a coordinate - * means that the interval is considered as one subinterval consisting - * of the entire range. A value of two means that there are two subintervals - * each with one half of the range, etc. - * @param data_values A dim-dimensional table of data at each of the - * mesh points defined by the coordinate arrays above. Note that the - * Table class has a number of conversion constructors that allow - * converting other data types into a table where you specify this - * argument. - */ - InterpolatedUniformGridData (const std_cxx1x::array,dim> &interval_endpoints, - const std_cxx1x::array &n_subintervals, - const Table &data_values); - - /** - * Compute the value of the function set by bilinear interpolation of the - * given data set. - * - * @param p The point at which the function is to be evaluated. - * @param component The vector component. Since this function is scalar, - * only zero is a valid argument here. - * @return The interpolated value at this point. If the point lies outside - * the set of coordinates, the function is extended by a constant. - */ - virtual - double - value (const Point &p, - const unsigned int component = 0) const; - - private: - /** - * The set of interval endpoints in each of the coordinate directions. - */ - const std_cxx1x::array,dim> interval_endpoints; - - /** - * The number of subintervals in each of the coordinate directions. - */ - const std_cxx1x::array n_subintervals; - - /** - * The data that is to be interpolated. - */ - const Table data_values; + public: + /** + * Constructor + * @param interval_endpoints The left and right end points of the (uniformly + * subdivided) intervals in each of the coordinate directions. + * @param n_subdivisions The number of subintervals of the subintervals + * in each coordinate direction. A value of one for a coordinate + * means that the interval is considered as one subinterval consisting + * of the entire range. A value of two means that there are two subintervals + * each with one half of the range, etc. + * @param data_values A dim-dimensional table of data at each of the + * mesh points defined by the coordinate arrays above. Note that the + * Table class has a number of conversion constructors that allow + * converting other data types into a table where you specify this + * argument. + */ + InterpolatedUniformGridData (const std_cxx1x::array,dim> &interval_endpoints, + const std_cxx1x::array &n_subintervals, + const Table &data_values); + + /** + * Compute the value of the function set by bilinear interpolation of the + * given data set. + * + * @param p The point at which the function is to be evaluated. + * @param component The vector component. Since this function is scalar, + * only zero is a valid argument here. + * @return The interpolated value at this point. If the point lies outside + * the set of coordinates, the function is extended by a constant. + */ + virtual + double + value (const Point &p, + const unsigned int component = 0) const; + + private: + /** + * The set of interval endpoints in each of the coordinate directions. + */ + const std_cxx1x::array,dim> interval_endpoints; + + /** + * The number of subintervals in each of the coordinate directions. + */ + const std_cxx1x::array n_subintervals; + + /** + * The data that is to be interpolated. + */ + const Table data_values; }; } DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/base/function_parser.h b/include/deal.II/base/function_parser.h index 06c14811ee..bd0d90b0a1 100644 --- a/include/deal.II/base/function_parser.h +++ b/include/deal.II/base/function_parser.h @@ -180,232 +180,232 @@ template class Vector; template class FunctionParser : public Function { - public: - /** - * Constructor for Parsed functions. Its arguments are the same of the - * base class Function. The only difference is that this object needs to - * be initialized with initialize() method before you can use it. If an - * attempt to use this function is made before the initialize() method has - * been called, then an exception is thrown. - */ - FunctionParser (const unsigned int n_components = 1, - const double initial_time = 0.0); - - /** - * Destructor. Explicitly delete the FunctionParser objects (there is one - * for each component of the function). - */ - ~FunctionParser(); - - /** - * Type for the constant map. Used by the initialize() method. - */ - typedef std::map ConstMap; - - /** - * Iterator for the constants map. Used by the initialize() method. - */ - typedef ConstMap::iterator ConstMapIterator; - - /** - * Initialize the function. This methods accepts the following - * parameters: - * - * vars: a string with the variables that will be used by the - * expressions to be evaluated. Note that the variables can have any name - * (of course different from the function names defined above!), but the - * order IS important. The first variable will correspond to the first - * component of the point in which the function is evaluated, the second - * variable to the second component and so forth. If this function is also - * time dependent, then it is necessary to specify it by setting the - * time_dependent parameter to true. An exception is thrown if - * the number of variables specified here is different from dim (if this - * function is not time-dependent) or from dim+1 (if it is time- - * dependent). - * - * expressions: a list of strings containing the expressions that - * will be byte compiled by the internal parser (FunctionParser). Note - * that the size of this vector must match exactly the number of - * components of the FunctionParser, as declared in the constructor. If - * this is not the case, an exception is thrown. - * - * - * constants: a map of constants used to pass any necessary - * constant that we want to specify in our expressions (in the example - * above the number pi). An expression is valid if and only if it contains - * only defined variables and defined constants (other than the functions - * specified above). If a constant is given whose name is not valid (eg: - * constants["sin"] = 1.5;) an exception is thrown. - * - * time_dependent. If this is a time dependent function, then the - * last variable declared in vars is assumed to be the time - * variable, and this->get_time() is used to initialize it when evaluating - * the function. Naturally the number of variables parsed by the - * initialize() method in this case is dim+1. The value of this parameter - * defaults to false, i.e. do not consider time. - */ - void initialize (const std::string &vars, - const std::vector &expressions, - const ConstMap &constants, - const bool time_dependent = false); - - /** - * Same as above, but with an additional parameter: use_degrees. - * Parameter to decide if the trigonometric functions work in radians or - * degrees. The default for this parameter is false, i.e. use radians and - * not degrees. - * - * @note: this function is deprecated. Use the function without this - * argument instead (which has the default use_degrees=false). - */ - void initialize (const std::string &vars, - const std::vector &expressions, - const ConstMap &constants, - const bool time_dependent, - const bool use_degrees) DEAL_II_DEPRECATED; - - - /** - * Initialize the function. Same as above, but with an additional argument - * units - a map of units passed to FunctionParser via AddUnint. - * - * Can be used as "3cm". Have higher precedence in parsing, i.e. if cm=10 - * then 3/2cm is 3 /(2*10). - */ - void initialize (const std::string &vars, - const std::vector &expressions, - const ConstMap &constants, - const ConstMap &units, - const bool time_dependent = false, - const bool use_degrees = false) DEAL_II_DEPRECATED; - - /** - * Initialize the function. Same as above, but accepts a string rather - * than a vector of strings. If this is a vector valued function, its - * components are expected to be separated by a semicolon. An exception is - * thrown if this method is called and the number of components - * successfully parsed does not match the number of components of the base - * function. - */ - void initialize (const std::string &vars, - const std::string &expression, - const ConstMap &constants, - const bool time_dependent = false); - - /** - * Same as above, but with an additional parameter: use_degrees. - * Parameter to decide if the trigonometric functions work in radians or - * degrees. The default for this parameter is false, i.e. use radians and - * not degrees. - * - * @note: this function is deprecated. Use the function without this - * argument instead (which has the default use_degrees=false). - */ - void initialize (const std::string &vars, - const std::string &expression, - const ConstMap &constants, - const bool time_dependent, - const bool use_degrees) DEAL_II_DEPRECATED; - - /** - * Initialize the function. Same as above, but with units. - */ - void initialize (const std::string &vars, - const std::string &expression, - const ConstMap &constants, - const ConstMap &units, - const bool time_dependent = false, - const bool use_degrees = false) DEAL_II_DEPRECATED; - - /** - * A function that returns default names for variables, to be used in the - * first argument of the initialize() functions: it returns "x" in 1d, - * "x,y" in 2d, and "x,y,z" in 3d. - */ - static - std::string - default_variable_names (); - - /** - * Return the value of the function at the given point. Unless there is - * only one component (i.e. the function is scalar), you should state the - * component you want to have evaluated; it defaults to zero, i.e. the - * first component. - */ - virtual double value (const Point &p, - const unsigned int component = 0) const; - - /** - * Return all components of a vector-valued function at the given point @p - * p. - * - * values shall have the right size beforehand, i.e. - * #n_components. - */ - virtual void vector_value (const Point &p, - Vector &values) const; - - /** - * @addtogroup Exceptions - * @{ - */ - DeclException2 (ExcParseError, - int, char *, - << "Parsing Error at Column " << arg1 - << ". The parser said: " << arg2); - - DeclException2 (ExcInvalidExpressionSize, - int, int, - << "The number of components (" << arg1 - << ") is not equal to the number of expressions (" - << arg2 << ")."); - - //@} - - private: +public: + /** + * Constructor for Parsed functions. Its arguments are the same of the + * base class Function. The only difference is that this object needs to + * be initialized with initialize() method before you can use it. If an + * attempt to use this function is made before the initialize() method has + * been called, then an exception is thrown. + */ + FunctionParser (const unsigned int n_components = 1, + const double initial_time = 0.0); + + /** + * Destructor. Explicitly delete the FunctionParser objects (there is one + * for each component of the function). + */ + ~FunctionParser(); + + /** + * Type for the constant map. Used by the initialize() method. + */ + typedef std::map ConstMap; + + /** + * Iterator for the constants map. Used by the initialize() method. + */ + typedef ConstMap::iterator ConstMapIterator; + + /** + * Initialize the function. This methods accepts the following + * parameters: + * + * vars: a string with the variables that will be used by the + * expressions to be evaluated. Note that the variables can have any name + * (of course different from the function names defined above!), but the + * order IS important. The first variable will correspond to the first + * component of the point in which the function is evaluated, the second + * variable to the second component and so forth. If this function is also + * time dependent, then it is necessary to specify it by setting the + * time_dependent parameter to true. An exception is thrown if + * the number of variables specified here is different from dim (if this + * function is not time-dependent) or from dim+1 (if it is time- + * dependent). + * + * expressions: a list of strings containing the expressions that + * will be byte compiled by the internal parser (FunctionParser). Note + * that the size of this vector must match exactly the number of + * components of the FunctionParser, as declared in the constructor. If + * this is not the case, an exception is thrown. + * + * + * constants: a map of constants used to pass any necessary + * constant that we want to specify in our expressions (in the example + * above the number pi). An expression is valid if and only if it contains + * only defined variables and defined constants (other than the functions + * specified above). If a constant is given whose name is not valid (eg: + * constants["sin"] = 1.5;) an exception is thrown. + * + * time_dependent. If this is a time dependent function, then the + * last variable declared in vars is assumed to be the time + * variable, and this->get_time() is used to initialize it when evaluating + * the function. Naturally the number of variables parsed by the + * initialize() method in this case is dim+1. The value of this parameter + * defaults to false, i.e. do not consider time. + */ + void initialize (const std::string &vars, + const std::vector &expressions, + const ConstMap &constants, + const bool time_dependent = false); + + /** + * Same as above, but with an additional parameter: use_degrees. + * Parameter to decide if the trigonometric functions work in radians or + * degrees. The default for this parameter is false, i.e. use radians and + * not degrees. + * + * @note: this function is deprecated. Use the function without this + * argument instead (which has the default use_degrees=false). + */ + void initialize (const std::string &vars, + const std::vector &expressions, + const ConstMap &constants, + const bool time_dependent, + const bool use_degrees) DEAL_II_DEPRECATED; + + + /** + * Initialize the function. Same as above, but with an additional argument + * units - a map of units passed to FunctionParser via AddUnint. + * + * Can be used as "3cm". Have higher precedence in parsing, i.e. if cm=10 + * then 3/2cm is 3 /(2*10). + */ + void initialize (const std::string &vars, + const std::vector &expressions, + const ConstMap &constants, + const ConstMap &units, + const bool time_dependent = false, + const bool use_degrees = false) DEAL_II_DEPRECATED; + + /** + * Initialize the function. Same as above, but accepts a string rather + * than a vector of strings. If this is a vector valued function, its + * components are expected to be separated by a semicolon. An exception is + * thrown if this method is called and the number of components + * successfully parsed does not match the number of components of the base + * function. + */ + void initialize (const std::string &vars, + const std::string &expression, + const ConstMap &constants, + const bool time_dependent = false); + + /** + * Same as above, but with an additional parameter: use_degrees. + * Parameter to decide if the trigonometric functions work in radians or + * degrees. The default for this parameter is false, i.e. use radians and + * not degrees. + * + * @note: this function is deprecated. Use the function without this + * argument instead (which has the default use_degrees=false). + */ + void initialize (const std::string &vars, + const std::string &expression, + const ConstMap &constants, + const bool time_dependent, + const bool use_degrees) DEAL_II_DEPRECATED; + + /** + * Initialize the function. Same as above, but with units. + */ + void initialize (const std::string &vars, + const std::string &expression, + const ConstMap &constants, + const ConstMap &units, + const bool time_dependent = false, + const bool use_degrees = false) DEAL_II_DEPRECATED; + + /** + * A function that returns default names for variables, to be used in the + * first argument of the initialize() functions: it returns "x" in 1d, + * "x,y" in 2d, and "x,y,z" in 3d. + */ + static + std::string + default_variable_names (); + + /** + * Return the value of the function at the given point. Unless there is + * only one component (i.e. the function is scalar), you should state the + * component you want to have evaluated; it defaults to zero, i.e. the + * first component. + */ + virtual double value (const Point &p, + const unsigned int component = 0) const; + + /** + * Return all components of a vector-valued function at the given point @p + * p. + * + * values shall have the right size beforehand, i.e. + * #n_components. + */ + virtual void vector_value (const Point &p, + Vector &values) const; + + /** + * @addtogroup Exceptions + * @{ + */ + DeclException2 (ExcParseError, + int, char *, + << "Parsing Error at Column " << arg1 + << ". The parser said: " << arg2); + + DeclException2 (ExcInvalidExpressionSize, + int, int, + << "The number of components (" << arg1 + << ") is not equal to the number of expressions (" + << arg2 << ")."); + + //@} + +private: #ifdef DEAL_II_WITH_MUPARSER - /** - * place for the variables for each thread - */ - mutable Threads::ThreadLocalStorage > vars; - /** - * the muParser objects for each thread (and one for each component) - */ - mutable Threads::ThreadLocalStorage > fp; - - /** - * keep track of all the constants, required to initialize fp in each - * thread - */ - std::map constants; - /** - * variable names, required to initialize fp in each thread - */ - std::vector var_names; - /** - * the expressions, required to initialize fp in each thread - */ - std::vector expressions; - /** - * this function will initialize fp on the current thread - */ - void init_muparser() const; + /** + * place for the variables for each thread + */ + mutable Threads::ThreadLocalStorage > vars; + /** + * the muParser objects for each thread (and one for each component) + */ + mutable Threads::ThreadLocalStorage > fp; + + /** + * keep track of all the constants, required to initialize fp in each + * thread + */ + std::map constants; + /** + * variable names, required to initialize fp in each thread + */ + std::vector var_names; + /** + * the expressions, required to initialize fp in each thread + */ + std::vector expressions; + /** + * this function will initialize fp on the current thread + */ + void init_muparser() const; #endif - /** - * State of usability. This variable is checked every time the function is - * called for evaluation. It's set to true in the initialize() methods. - */ - bool initialized; - - /** - * Number of variables. If this is also a function of time, then the - * number of variables is dim+1, otherwise it is dim. In the case that - * this is a time dependent function, the time is supposed to be the last - * variable. If #n_vars is not identical to the number of the variables - * parsed by the initialize() method, then an exception is thrown. - */ - unsigned int n_vars; + /** + * State of usability. This variable is checked every time the function is + * called for evaluation. It's set to true in the initialize() methods. + */ + bool initialized; + + /** + * Number of variables. If this is also a function of time, then the + * number of variables is dim+1, otherwise it is dim. In the case that + * this is a time dependent function, the time is supposed to be the last + * variable. If #n_vars is not identical to the number of the variables + * parsed by the initialize() method, then an exception is thrown. + */ + unsigned int n_vars; }; @@ -414,16 +414,16 @@ std::string FunctionParser::default_variable_names () { switch (dim) - { - case 1: - return "x"; - case 2: - return "x,y"; - case 3: - return "x,y,z"; - default: - Assert (false, ExcNotImplemented()); - } + { + case 1: + return "x"; + case 2: + return "x,y"; + case 3: + return "x,y,z"; + default: + Assert (false, ExcNotImplemented()); + } return ""; } diff --git a/include/deal.II/base/iterator_range.h b/include/deal.II/base/iterator_range.h index c752cc2954..1533ea004a 100644 --- a/include/deal.II/base/iterator_range.h +++ b/include/deal.II/base/iterator_range.h @@ -124,7 +124,7 @@ public: * as discussed in the design sections of the IteratorRange class. */ class IteratorOverIterators : public std::iterator + typename Iterator::difference_type> { public: /** @@ -152,13 +152,13 @@ public: * Dereferencing operator. * @return The iterator within the collection currently pointed to. */ - const BaseIterator * operator-> () const; + const BaseIterator *operator-> () const; /** * Prefix increment operator. Move the current iterator to the next * element of the collection and return the new value. */ - IteratorOverIterators & operator ++ (); + IteratorOverIterators &operator ++ (); /** * Postfix increment operator. Move the current iterator to the next diff --git a/include/deal.II/base/named_data.h b/include/deal.II/base/named_data.h index 59e2a96b88..5357bd66d3 100644 --- a/include/deal.II/base/named_data.h +++ b/include/deal.II/base/named_data.h @@ -266,7 +266,7 @@ void NamedData::add(DATA &v, const std::string &n) { // see operator() below - + // Assert(!is_constant, ExcConstantObject()); names.push_back(n); data.push_back(v); @@ -279,7 +279,7 @@ void NamedData::add(const DATA &v, const std::string &n) { // see operator() below - + // Assert(!is_constant, ExcConstantObject()); DATA &aux = const_cast(v); data.push_back(aux); @@ -295,7 +295,7 @@ void NamedData::merge(NamedData &other) { // see operator() below - + // Assert(!is_constant, ExcConstantObject()); for (unsigned int i=0; i::merge(const NamedData &other) { // see operator() below - + // Assert(!is_constant, ExcConstantObject()); for (unsigned int i=0; iXML format, the output starts with one root element * ParameterHandler in order to get a valid XML document * and all subsections under it. - * + * * In Text format, the output contains the same information but * in a format so that the resulting file can be input into a latex * document such as a manual for the code for which this object handles @@ -1825,7 +1825,7 @@ public: * the higher subsection elements are printed. In XML format * this is required to get a valid XML document and output starts * with one root element ParameterHandler. - * + * * In most cases, you will not want to use this function directly, but * have it called recursively by the previous function. */ diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index 5bfe461aa0..0cdb84be60 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -451,10 +451,10 @@ class QSorted : public Quadrature public: /** The constructor takes an arbitrary quadrature formula. */ QSorted (const Quadrature); - + /** A rule to reorder pairs of points and weights.*/ bool operator()(const std::pair > &a, - const std::pair > &b); + const std::pair > &b); }; diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index 66814824fe..67710ad64e 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -1065,7 +1065,7 @@ namespace internal Accessor:: Accessor () : - tensor (*static_cast(0)), + tensor (*static_cast(0)), previous_indices () { Assert (false, ExcMessage ("You can't call the default constructor of this class.")); @@ -1106,7 +1106,7 @@ namespace internal Accessor:: Accessor () : - tensor (*static_cast(0)), + tensor (*static_cast(0)), previous_indices () { Assert (false, ExcMessage ("You can't call the default constructor of this class.")); diff --git a/include/deal.II/base/table.h b/include/deal.II/base/table.h index a861fa0089..ad244837a8 100644 --- a/include/deal.II/base/table.h +++ b/include/deal.II/base/table.h @@ -1852,7 +1852,7 @@ protected: * then called vector2d. */ typename AlignedVector::reference el (const unsigned int i, - const unsigned int j); + const unsigned int j); /** * Return the value of the @@ -1879,7 +1879,7 @@ protected: * then called vector2d. */ typename AlignedVector::const_reference el (const unsigned int i, - const unsigned int j) const; + const unsigned int j) const; }; @@ -2282,7 +2282,7 @@ namespace internal template void fill_Fortran_style (InputIterator, - TableBase &) + TableBase &) { Assert (false, ExcNotImplemented()); } @@ -2302,7 +2302,7 @@ TableBase::fill (InputIterator entries, if (C_style_indexing) for (typename AlignedVector::iterator p = values.begin(); - p != values.end(); ++p) + p != values.end(); ++p) *p = *entries++; else internal::Table::fill_Fortran_style (entries, *this); diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index c3c6dab702..05e71b4982 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -2102,10 +2102,10 @@ operator / (const Tensor &t, * @relates Tensor */ template - inline - Tensor > - operator * (const std::complex factor, - const Tensor &t) +inline +Tensor > +operator * (const std::complex factor, + const Tensor &t) { Tensor > tt; for (unsigned int d=0; d * @relates Tensor */ template - inline - Tensor > - operator * (const Tensor &t, - const std::complex factor) +inline +Tensor > +operator * (const Tensor &t, + const std::complex factor) { Tensor > tt; for (unsigned int d=0; d &t, template inline Tensor<1,dim,std::complex > - operator * (const Tensor<1,dim> &t, - const std::complex factor) -{ +operator * (const Tensor<1,dim> &t, + const std::complex factor) +{ Tensor<1,dim,std::complex > tt (false); for (unsigned int d=0; d > template inline Tensor<1,dim,std::complex > - operator * (const std::complex factor, - const Tensor<1,dim> &t) -{ +operator * (const std::complex factor, + const Tensor<1,dim> &t) +{ Tensor<1,dim,std::complex > tt (false); for (unsigned int d=0; d > &F, - std::vector > & J_inverse, + std::vector > &J_inverse, double t, double delta_t, VECTOR &y) = 0; @@ -147,7 +147,7 @@ namespace TimeStepping */ double evolve_one_time_step( std::vector > &F, - std::vector > & J_inverse, + std::vector > &J_inverse, double t, double delta_t, VECTOR &y); diff --git a/include/deal.II/base/time_stepping.templates.h b/include/deal.II/base/time_stepping.templates.h index b80037ec03..a63e818b7a 100644 --- a/include/deal.II/base/time_stepping.templates.h +++ b/include/deal.II/base/time_stepping.templates.h @@ -32,7 +32,7 @@ namespace TimeStepping template double RungeKutta::evolve_one_time_step( std::vector > &F, - std::vector > & J_inverse, + std::vector > &J_inverse, double t, double delta_t, VECTOR &y) @@ -161,7 +161,7 @@ namespace TimeStepping // Linear combinations of the stages. for (unsigned int i=0; in_stages; ++i) - y.sadd(1.,delta_t*this->b[i],f_stages[i]); + y.sadd(1.,delta_t *this->b[i],f_stages[i]); return (t+delta_t); } @@ -188,7 +188,7 @@ namespace TimeStepping { VECTOR Y(y); for (unsigned int j=0; ja[i][j],f_stages[j]); + Y.sadd(1.,delta_t *this->a[i][j],f_stages[j]); // Evaluate the function f at the point (t+c[i]*delta_t,Y). f_stages[i] = f(t+this->c[i]*delta_t,Y); } @@ -296,7 +296,7 @@ namespace TimeStepping { y = old_y; for (unsigned int i=0; in_stages; ++i) - y.sadd(1.,delta_t*this->b[i],f_stages[i]); + y.sadd(1.,delta_t *this->b[i],f_stages[i]); } return (t+delta_t); @@ -335,7 +335,7 @@ namespace TimeStepping { VECTOR old_y(z); for (unsigned int j=0; ja[i][j],f_stages[j]); + old_y.sadd(1.,delta_t *this->a[i][j],f_stages[j]); // Solve the nonlinear system using Newton's method const double new_t = t+this->c[i]*delta_t; @@ -705,8 +705,8 @@ namespace TimeStepping for (unsigned int i=0; in_stages; ++i) { - y.sadd(1.,delta_t*this->b1[i],f_stages[i]); - error.sadd(1.,delta_t*(b2[i]-b1[i]),f_stages[i]); + y.sadd(1.,delta_t *this->b1[i],f_stages[i]); + error.sadd(1.,delta_t *(b2[i]-b1[i]),f_stages[i]); } error_norm = error.l2_norm(); @@ -715,7 +715,7 @@ namespace TimeStepping { done = true; // Increase the guessed time step - double new_delta_t = delta_t*coarsen_param; + double new_delta_t = delta_t *coarsen_param; // Check that the guessed time step is smaller than the maximum time // step allowed. if (new_delta_t>max_delta_t) @@ -825,7 +825,7 @@ namespace TimeStepping { Y = y; for (unsigned int j = 0; j < i; ++j) - Y.sadd(1.0,delta_t*this->a[i][j],f_stages[j]); + Y.sadd(1.0,delta_t *this->a[i][j],f_stages[j]); f_stages[i] = f(t+this->c[i]*delta_t,Y); } } diff --git a/include/deal.II/base/types.h b/include/deal.II/base/types.h index fcc25fd692..f82d52e226 100644 --- a/include/deal.II/base/types.h +++ b/include/deal.II/base/types.h @@ -126,7 +126,7 @@ namespace types * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" */ typedef unsigned char boundary_id; - + /** * The type used to denote manifold indicators associated with every * object of the mesh. @@ -138,7 +138,7 @@ namespace types */ typedef unsigned int manifold_id; - + /** * @deprecated Old name for the typedef above. */ @@ -274,7 +274,7 @@ namespace numbers * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators" */ const types::manifold_id invalid_manifold_id = static_cast(-1); - + /** * A manifold_id we reserve for the default flat Cartesian manifold. * diff --git a/include/deal.II/base/vectorization.h b/include/deal.II/base/vectorization.h index ca843ed5dd..9ca30e57bc 100644 --- a/include/deal.II/base/vectorization.h +++ b/include/deal.II/base/vectorization.h @@ -177,7 +177,7 @@ public: * the given address. The memory need not be aligned by 64 bytes, as opposed * to casting a double address to VectorizedArray*. */ - void load (const double* ptr) + void load (const double *ptr) { data = _mm512_loadu_pd (ptr); } @@ -188,7 +188,7 @@ public: * 64 bytes, as opposed to casting a double address to * VectorizedArray*. */ - void store (double* ptr) const + void store (double *ptr) const { _mm512_storeu_pd (ptr, data); } @@ -369,7 +369,7 @@ public: * the given address. The memory need not be aligned by 64 bytes, as opposed * to casting a float address to VectorizedArray*. */ - void load (const float* ptr) + void load (const float *ptr) { data = _mm512_loadu_ps (ptr); } @@ -380,7 +380,7 @@ public: * 64 bytes, as opposed to casting a float address to * VectorizedArray*. */ - void store (float* ptr) const + void store (float *ptr) const { _mm512_storeu_ps (ptr, data); } @@ -573,7 +573,7 @@ public: * the given address. The memory need not be aligned by 32 bytes, as opposed * to casting a double address to VectorizedArray*. */ - void load (const double* ptr) + void load (const double *ptr) { data = _mm256_loadu_pd (ptr); } @@ -584,7 +584,7 @@ public: * 32 bytes, as opposed to casting a double address to * VectorizedArray*. */ - void store (double* ptr) const + void store (double *ptr) const { _mm256_storeu_pd (ptr, data); } @@ -772,7 +772,7 @@ public: * the given address. The memory need not be aligned by 32 bytes, as opposed * to casting a float address to VectorizedArray*. */ - void load (const float* ptr) + void load (const float *ptr) { data = _mm256_loadu_ps (ptr); } @@ -783,7 +783,7 @@ public: * 32 bytes, as opposed to casting a float address to * VectorizedArray*. */ - void store (float* ptr) const + void store (float *ptr) const { _mm256_storeu_ps (ptr, data); } @@ -979,7 +979,7 @@ public: * the given address. The memory need not be aligned by 16 bytes, as opposed * to casting a double address to VectorizedArray*. */ - void load (const double* ptr) + void load (const double *ptr) { data = _mm_loadu_pd (ptr); } @@ -990,7 +990,7 @@ public: * 16 bytes, as opposed to casting a double address to * VectorizedArray*. */ - void store (double* ptr) const + void store (double *ptr) const { _mm_storeu_pd (ptr, data); } @@ -1179,7 +1179,7 @@ public: * the given address. The memory need not be aligned by 16 bytes, as opposed * to casting a float address to VectorizedArray*. */ - void load (const float* ptr) + void load (const float *ptr) { data = _mm_loadu_ps (ptr); } @@ -1190,7 +1190,7 @@ public: * 16 bytes, as opposed to casting a float address to * VectorizedArray*. */ - void store (float* ptr) const + void store (float *ptr) const { _mm_storeu_ps (ptr, data); } @@ -1420,7 +1420,7 @@ public: * in the vectorized array, as opposed to casting a double address to * VectorizedArray*. */ - void load (const Number* ptr) + void load (const Number *ptr) { data = *ptr; } @@ -1431,7 +1431,7 @@ public: * the amount of bytes in the vectorized array, as opposed to casting a * double address to VectorizedArray*. */ - void store (Number* ptr) const + void store (Number *ptr) const { *ptr = data; } diff --git a/include/deal.II/base/work_stream.h b/include/deal.II/base/work_stream.h index fd30a26371..ccbd6b6565 100644 --- a/include/deal.II/base/work_stream.h +++ b/include/deal.II/base/work_stream.h @@ -349,15 +349,15 @@ namespace WorkStream // and have set the ring buffer to have exactly this size. so // if this function is called, we know that less than the // maximal number of items in currently in flight - // - // note that we need not lock access to this array since - // the current stage is run sequentially and we can therefore - // enter the following block only once at any given time. - // thus, there can be no race condition between checking that - // a flag is false and setting it to true. (there may be - // another thread where we release items and set 'false' - // flags to 'true', but that too does not produce any - // problems) + // + // note that we need not lock access to this array since + // the current stage is run sequentially and we can therefore + // enter the following block only once at any given time. + // thus, there can be no race condition between checking that + // a flag is false and setting it to true. (there may be + // another thread where we release items and set 'false' + // flags to 'true', but that too does not produce any + // problems) ItemType *current_item = 0; for (unsigned int i=0; i + typename ScratchData, + typename CopyData> class Worker { public: - /** - * Constructor. - */ - Worker (const std_cxx1x::function &worker, - const ScratchData &sample_scratch_data, - const CopyData &sample_copy_data) - : - worker (worker), - sample_scratch_data (sample_scratch_data), - sample_copy_data (sample_copy_data) - {} - - - /** - * The function that calls the worker function on a - * range of items denoted by the two arguments. - */ - void operator() (const tbb::blocked_range::const_iterator> &range) - { - // we need to find an unused scratch and corresponding copy - // data object in the list that - // corresponds to the current thread and then mark it as used. if - // we can't find one, create one - // - // as discussed in the discussion of the documentation of the - // IteratorRangeToItemStream::scratch_data variable, there is no - // need to synchronize access to this variable using a mutex - // as long as we have no yield-point in between. this means that - // we can't take an iterator into the list now and expect it to - // still be valid after calling the worker, but we at least do - // not have to lock the following section - ScratchData *scratch_data = 0; - CopyData *copy_data = 0; - { - typename ItemType::ScratchAndCopyDataList & - scratch_and_copy_data_list = data.get(); - - // see if there is an unused object. if so, grab it and mark - // it as used - for (typename ItemType::ScratchAndCopyDataList::iterator - p = scratch_and_copy_data_list.begin(); - p != scratch_and_copy_data_list.end(); ++p) - if (p->currently_in_use == false) - { - scratch_data = p->scratch_data.get(); - copy_data = p->copy_data.get(); - p->currently_in_use = true; - break; - } - - // if no element in the list was found, create one and mark it as used - if (scratch_data == 0) - { - Assert (copy_data==0, ExcInternalError()); - scratch_data = new ScratchData(sample_scratch_data); - copy_data = new CopyData(sample_copy_data); - - typename ItemType::ScratchAndCopyDataList::value_type - new_scratch_object (scratch_data, copy_data, true); - scratch_and_copy_data_list.push_back (new_scratch_object); - } - } - - // then call the worker and copier function on each - // element of the chunk we were given. since these - // functions are called on separate threads, nothing good - // can happen if they throw an exception and we are best - // off catching it and showing an error message - for (typename std::vector::const_iterator p=range.begin(); - p != range.end(); ++p) - { - try - { - worker (*p, - *scratch_data, - *copy_data); - } - catch (const std::exception &exc) - { - Threads::internal::handle_std_exception (exc); - } - catch (...) - { - Threads::internal::handle_unknown_exception (); - } - } - - // finally mark the scratch object as unused again. as above, there - // is no need to lock anything here since the object we work on - // is thread-local - { - typename ItemType::ScratchAndCopyDataList & - scratch_and_copy_data_list = data.get(); - - for (typename ItemType::ScratchAndCopyDataList::iterator p = - scratch_and_copy_data_list.begin(); p != scratch_and_copy_data_list.end(); - ++p) - if (p->scratch_data.get() == scratch_data) - { - Assert(p->currently_in_use == true, ExcInternalError()); - p->currently_in_use = false; - } - } - - } + /** + * Constructor. + */ + Worker (const std_cxx1x::function &worker, + const ScratchData &sample_scratch_data, + const CopyData &sample_copy_data) + : + worker (worker), + sample_scratch_data (sample_scratch_data), + sample_copy_data (sample_copy_data) + {} + + + /** + * The function that calls the worker function on a + * range of items denoted by the two arguments. + */ + void operator() (const tbb::blocked_range::const_iterator> &range) + { + // we need to find an unused scratch and corresponding copy + // data object in the list that + // corresponds to the current thread and then mark it as used. if + // we can't find one, create one + // + // as discussed in the discussion of the documentation of the + // IteratorRangeToItemStream::scratch_data variable, there is no + // need to synchronize access to this variable using a mutex + // as long as we have no yield-point in between. this means that + // we can't take an iterator into the list now and expect it to + // still be valid after calling the worker, but we at least do + // not have to lock the following section + ScratchData *scratch_data = 0; + CopyData *copy_data = 0; + { + typename ItemType::ScratchAndCopyDataList & + scratch_and_copy_data_list = data.get(); + + // see if there is an unused object. if so, grab it and mark + // it as used + for (typename ItemType::ScratchAndCopyDataList::iterator + p = scratch_and_copy_data_list.begin(); + p != scratch_and_copy_data_list.end(); ++p) + if (p->currently_in_use == false) + { + scratch_data = p->scratch_data.get(); + copy_data = p->copy_data.get(); + p->currently_in_use = true; + break; + } + + // if no element in the list was found, create one and mark it as used + if (scratch_data == 0) + { + Assert (copy_data==0, ExcInternalError()); + scratch_data = new ScratchData(sample_scratch_data); + copy_data = new CopyData(sample_copy_data); + + typename ItemType::ScratchAndCopyDataList::value_type + new_scratch_object (scratch_data, copy_data, true); + scratch_and_copy_data_list.push_back (new_scratch_object); + } + } + + // then call the worker and copier function on each + // element of the chunk we were given. since these + // functions are called on separate threads, nothing good + // can happen if they throw an exception and we are best + // off catching it and showing an error message + for (typename std::vector::const_iterator p=range.begin(); + p != range.end(); ++p) + { + try + { + worker (*p, + *scratch_data, + *copy_data); + } + catch (const std::exception &exc) + { + Threads::internal::handle_std_exception (exc); + } + catch (...) + { + Threads::internal::handle_unknown_exception (); + } + } + + // finally mark the scratch object as unused again. as above, there + // is no need to lock anything here since the object we work on + // is thread-local + { + typename ItemType::ScratchAndCopyDataList & + scratch_and_copy_data_list = data.get(); + + for (typename ItemType::ScratchAndCopyDataList::iterator p = + scratch_and_copy_data_list.begin(); p != scratch_and_copy_data_list.end(); + ++p) + if (p->scratch_data.get() == scratch_data) + { + Assert(p->currently_in_use == true, ExcInternalError()); + p->currently_in_use = false; + } + } + + } private: - typedef - typename Implementation3::IteratorRangeToItemStream::ItemType - ItemType; + typedef + typename Implementation3::IteratorRangeToItemStream::ItemType + ItemType; - typedef - typename ItemType::ScratchAndCopyDataList - ScratchAndCopyDataList; + typedef + typename ItemType::ScratchAndCopyDataList + ScratchAndCopyDataList; - Threads::ThreadLocalStorage data; + Threads::ThreadLocalStorage data; /** * Pointer to the function @@ -1312,12 +1312,12 @@ namespace WorkStream ScratchData &, CopyData &)> worker; - /** - * References to sample scratch and copy data for - * when we need them. - */ - const ScratchData &sample_scratch_data; - const CopyData &sample_copy_data; + /** + * References to sample scratch and copy data for + * when we need them. + */ + const ScratchData &sample_scratch_data; + const CopyData &sample_copy_data; }; } @@ -1470,11 +1470,11 @@ namespace WorkStream { // create the three stages of the pipeline internal::Implementation2::IteratorRangeToItemStream - iterator_range_to_item_stream (begin, end, - queue_length, - chunk_size, - sample_scratch_data, - sample_copy_data); + iterator_range_to_item_stream (begin, end, + queue_length, + chunk_size, + sample_scratch_data, + sample_copy_data); internal::Implementation2::Worker worker_filter (worker); internal::Implementation2::Copier copier_filter (copier); @@ -1575,65 +1575,65 @@ namespace WorkStream // loop over the various colors of what we're given for (unsigned int color=0; color 0) - { - if (static_cast& >(copier)) - { - // there is a copier function, so we have to go with - // the full three-stage design of the pipeline - internal::Implementation3::IteratorRangeToItemStream - iterator_range_to_item_stream (colored_iterators[color].begin(), - colored_iterators[color].end(), - queue_length, - chunk_size, - sample_scratch_data, - sample_copy_data); - - - internal::Implementation3::WorkerAndCopier - worker_and_copier_filter (worker, copier); - - // now create a pipeline from these stages - tbb::pipeline assembly_line; - assembly_line.add_filter (iterator_range_to_item_stream); - assembly_line.add_filter (worker_and_copier_filter); - - // and run it - assembly_line.run (queue_length); - - assembly_line.clear (); - } - else - { - // no copier function, we can implement things as a parallel for - Assert (static_cast& >(worker), - ExcMessage ("It makes no sense to call this function with " - "empty functions for both the worker and the " - "copier!")); - - typedef - internal::ParallelFor::Worker - ParallelForWorker; - - typedef - typename std::vector::const_iterator - RangeType; - - ParallelForWorker parallel_for_worker (worker, - sample_scratch_data, - sample_copy_data); - - tbb::parallel_for (tbb::blocked_range - (colored_iterators[color].begin(), - colored_iterators[color].end(), - /*grain_size=*/chunk_size), - std_cxx1x::bind (&ParallelForWorker::operator(), - std_cxx1x::ref(parallel_for_worker), - std_cxx1x::_1), - tbb::auto_partitioner()); - } - } + { + if (static_cast& >(copier)) + { + // there is a copier function, so we have to go with + // the full three-stage design of the pipeline + internal::Implementation3::IteratorRangeToItemStream + iterator_range_to_item_stream (colored_iterators[color].begin(), + colored_iterators[color].end(), + queue_length, + chunk_size, + sample_scratch_data, + sample_copy_data); + + + internal::Implementation3::WorkerAndCopier + worker_and_copier_filter (worker, copier); + + // now create a pipeline from these stages + tbb::pipeline assembly_line; + assembly_line.add_filter (iterator_range_to_item_stream); + assembly_line.add_filter (worker_and_copier_filter); + + // and run it + assembly_line.run (queue_length); + + assembly_line.clear (); + } + else + { + // no copier function, we can implement things as a parallel for + Assert (static_cast& >(worker), + ExcMessage ("It makes no sense to call this function with " + "empty functions for both the worker and the " + "copier!")); + + typedef + internal::ParallelFor::Worker + ParallelForWorker; + + typedef + typename std::vector::const_iterator + RangeType; + + ParallelForWorker parallel_for_worker (worker, + sample_scratch_data, + sample_copy_data); + + tbb::parallel_for (tbb::blocked_range + (colored_iterators[color].begin(), + colored_iterators[color].end(), + /*grain_size=*/chunk_size), + std_cxx1x::bind (&ParallelForWorker::operator(), + std_cxx1x::ref(parallel_for_worker), + std_cxx1x::_1), + tbb::auto_partitioner()); + } + } } #endif } diff --git a/include/deal.II/distributed/tria.h b/include/deal.II/distributed/tria.h index c67d68dcb6..740231ebe3 100644 --- a/include/deal.II/distributed/tria.h +++ b/include/deal.II/distributed/tria.h @@ -181,7 +181,7 @@ namespace parallel * * @note This class does not support anisotropic refinement, because * it relies on the p4est library that does not support this. Attempts - * to refine cells anisotropically will result in errors. + * to refine cells anisotropically will result in errors. * @note There is currently no support for distributing 1d triangulations. * * @@ -318,24 +318,24 @@ namespace parallel class Triangulation : public dealii::Triangulation { public: - /** - * A typedef that is used to to identify cell iterators. The - * concept of iterators is discussed at length in the - * @ref Iterators "iterators documentation module". - * - * The current typedef identifies cells in a triangulation. You - * can find the exact type it refers to in the base class's own - * typedef, but it should be TriaIterator >. The - * TriaIterator class works like a pointer that when you - * dereference it yields an object of type CellAccessor. - * CellAccessor is a class that identifies properties that - * are specific to cells in a triangulation, but it is derived - * (and consequently inherits) from TriaAccessor that describes - * what you can ask of more general objects (lines, faces, as - * well as cells) in a triangulation. - * - * @ingroup Iterators - */ + /** + * A typedef that is used to to identify cell iterators. The + * concept of iterators is discussed at length in the + * @ref Iterators "iterators documentation module". + * + * The current typedef identifies cells in a triangulation. You + * can find the exact type it refers to in the base class's own + * typedef, but it should be TriaIterator >. The + * TriaIterator class works like a pointer that when you + * dereference it yields an object of type CellAccessor. + * CellAccessor is a class that identifies properties that + * are specific to cells in a triangulation, but it is derived + * (and consequently inherits) from TriaAccessor that describes + * what you can ask of more general objects (lines, faces, as + * well as cells) in a triangulation. + * + * @ingroup Iterators + */ typedef typename dealii::Triangulation::cell_iterator cell_iterator; /** @@ -616,23 +616,23 @@ namespace parallel */ enum CellStatus { - /** - * The cell will not be refined or coarsened and might or might - * not move to a different processor. - */ + /** + * The cell will not be refined or coarsened and might or might + * not move to a different processor. + */ CELL_PERSIST, - /** - * The cell will be or was refined. - */ - CELL_REFINE, - /** - * The children of this cell will be or were coarsened into this cell. - */ - CELL_COARSEN, - /** - * Invalid status. Will not occur for the user. - */ - CELL_INVALID + /** + * The cell will be or was refined. + */ + CELL_REFINE, + /** + * The children of this cell will be or were coarsened into this cell. + */ + CELL_COARSEN, + /** + * Invalid status. Will not occur for the user. + */ + CELL_INVALID }; /** @@ -649,19 +649,19 @@ namespace parallel * Callers need to store the return value. It specifies an * offset of the position at which data can later be retrieved * during a call to notify_ready_to_unpack(). - * + * * The CellStatus argument in the callback function will tell you if the * given cell will be coarsened, refined, or will persist as is (this * can be different than the coarsen and refine flags set by you). If it * is - * + * * - CELL_PERIST: the cell won't be refined/coarsened, but might be * moved to a different processor * - CELL_REFINE: this cell will be refined into 4/8 cells, you can not * access the children (because they don't exist yet) * - CELL_COARSEN: the children of this cell will be coarsened into the * given cell (you can access the active children!) - * + * * When unpacking the data with notify_ready_to_unpack() you can access * the children of the cell if the status is CELL_REFINE but not for * CELL_COARSEN. As a consequence you need to handle coarsening while diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h index 3e928f8cfe..da7e7e25c6 100644 --- a/include/deal.II/dofs/dof_accessor.templates.h +++ b/include/deal.II/dofs/dof_accessor.templates.h @@ -81,11 +81,11 @@ DoFAccessor::DoFAccessor (const DoFAccessor inline types::global_dof_index DoFAccessor::dof_index (const unsigned int i, - const unsigned int fe_index) const + const unsigned int fe_index) const { // access the respective DoF return dealii::internal::DoFAccessor::Implementation::get_dof_index (*this->dof_handler, @@ -1467,7 +1467,7 @@ template inline types::global_dof_index DoFAccessor::mg_dof_index (const int level, - const unsigned int i) const + const unsigned int i) const { return this->dof_handler->template get_dof_index (level, this->present_index, 0, i); } @@ -1477,8 +1477,8 @@ template inline void DoFAccessor::set_dof_index (const unsigned int i, - const types::global_dof_index index, - const unsigned int fe_index) const + const types::global_dof_index index, + const unsigned int fe_index) const { // access the respective DoF dealii::internal::DoFAccessor::Implementation::set_dof_index (*this->dof_handler, @@ -1546,8 +1546,8 @@ template inline types::global_dof_index DoFAccessor::vertex_dof_index (const unsigned int vertex, - const unsigned int i, - const unsigned int fe_index) const + const unsigned int i, + const unsigned int fe_index) const { return dealii::internal::DoFAccessor::Implementation::get_vertex_dof_index @@ -1562,9 +1562,9 @@ template inline types::global_dof_index DoFAccessor::mg_vertex_dof_index (const int level, - const unsigned int vertex, - const unsigned int i, - const unsigned int fe_index) const + const unsigned int vertex, + const unsigned int i, + const unsigned int fe_index) const { Assert (this->dof_handler != 0, ExcInvalidObject ()); Assert (&this->dof_handler->get_fe () != 0, ExcInvalidObject ()); @@ -1578,9 +1578,9 @@ template inline void DoFAccessor::set_vertex_dof_index (const unsigned int vertex, - const unsigned int i, - const types::global_dof_index index, - const unsigned int fe_index) const + const unsigned int i, + const types::global_dof_index index, + const unsigned int fe_index) const { dealii::internal::DoFAccessor::Implementation::set_vertex_dof_index (*this->dof_handler, @@ -1595,10 +1595,10 @@ template inline void DoFAccessor::set_mg_vertex_dof_index (const int level, - const unsigned int vertex, - const unsigned int i, - const types::global_dof_index index, - const unsigned int fe_index) const + const unsigned int vertex, + const unsigned int i, + const types::global_dof_index index, + const unsigned int fe_index) const { Assert (this->dof_handler != 0, ExcInvalidObject ()); Assert (&this->dof_handler->get_fe () != 0, ExcInvalidObject ()); @@ -1612,8 +1612,8 @@ template inline void DoFAccessor::set_mg_dof_index (const int level, - const unsigned int i, - const types::global_dof_index index) const + const unsigned int i, + const types::global_dof_index index) const { this->dof_handler->template set_dof_index (level, this->present_index, 0, i, index); } @@ -1866,7 +1866,7 @@ template inline void DoFAccessor::get_dof_indices (std::vector &dof_indices, - const unsigned int fe_index) const + const unsigned int fe_index) const { Assert (this->dof_handler != 0, ExcNotInitialized()); Assert (&this->dof_handler->get_fe() != 0, ExcMessage ("DoFHandler not initialized")); @@ -1928,8 +1928,8 @@ DoFAccessor::get_dof_indices (std::vector inline void DoFAccessor::get_mg_dof_indices (const int level, - std::vector &dof_indices, - const unsigned int fe_index) const + std::vector &dof_indices, + const unsigned int fe_index) const { Assert (this->dof_handler != 0, ExcInvalidObject ()); Assert (&this->dof_handler->get_fe () != 0, ExcInvalidObject ()); @@ -1980,8 +1980,8 @@ void DoFAccessor::get_mg_dof_indices (const int template inline void DoFAccessor::set_mg_dof_indices (const int level, - const std::vector &dof_indices, - const unsigned int fe_index) + const std::vector &dof_indices, + const unsigned int fe_index) { Assert (this->dof_handler != 0, ExcInvalidObject ()); Assert (&this->dof_handler->get_fe () != 0, ExcInvalidObject ()); @@ -3309,7 +3309,7 @@ template inline void DoFCellAccessor::get_dof_values (const InputVector &values, - Vector &local_values) const + Vector &local_values) const { get_dof_values (values, local_values.begin(), local_values.end()); } @@ -3321,8 +3321,8 @@ template inline void DoFCellAccessor::get_dof_values (const InputVector &values, - ForwardIterator local_values_begin, - ForwardIterator local_values_end) const + ForwardIterator local_values_begin, + ForwardIterator local_values_end) const { Assert (this->is_artificial() == false, ExcMessage ("Can't ask for DoF indices on artificial cells.")); @@ -3351,9 +3351,9 @@ template inline void DoFCellAccessor::get_dof_values (const ConstraintMatrix &constraints, - const InputVector &values, - ForwardIterator local_values_begin, - ForwardIterator local_values_end) const + const InputVector &values, + ForwardIterator local_values_begin, + ForwardIterator local_values_end) const { Assert (this->is_artificial() == false, ExcMessage ("Can't ask for DoF indices on artificial cells.")); @@ -3382,7 +3382,7 @@ template inline void DoFCellAccessor::set_dof_values (const Vector &local_values, - OutputVector &values) const + OutputVector &values) const { Assert (this->is_artificial() == false, ExcMessage ("Can't ask for DoF indices on artificial cells.")); @@ -3412,12 +3412,12 @@ const FiniteElement & DoFCellAccessor::get_fe () const { Assert ((dynamic_cast*> - (this->dof_handler) != 0) - || - (this->has_children() == false), + (this->dof_handler) != 0) + || + (this->has_children() == false), ExcMessage ("In hp::DoFHandler objects, finite elements are only associated " - "with active cells. Consequently, you can not ask for the " - "active finite element on cells with children.")); + "with active cells. Consequently, you can not ask for the " + "active finite element on cells with children.")); return dealii::internal::DoFAccessor::get_fe (this->dof_handler->get_fe(), active_fe_index()); } @@ -3429,13 +3429,13 @@ unsigned int DoFCellAccessor::active_fe_index () const { Assert ((dynamic_cast*> - (this->dof_handler) != 0) - || - (this->has_children() == false), + (this->dof_handler) != 0) + || + (this->has_children() == false), ExcMessage ("You can not ask for the active_fe_index on a cell that has " "children because no degrees of freedom are assigned " "to this cell and, consequently, no finite element " - "is associated with it.")); + "is associated with it.")); return dealii::internal::DoFCellAccessor::Implementation::active_fe_index (*this); } @@ -3447,9 +3447,9 @@ void DoFCellAccessor::set_active_fe_index (const unsigned int i) { Assert ((dynamic_cast*> - (this->dof_handler) != 0) - || - (this->has_children() == false), + (this->dof_handler) != 0) + || + (this->has_children() == false), ExcMessage ("You can not set the active_fe_index on a cell that has " "children because no degrees of freedom will be assigned " "to this cell.")); diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h index 2ddf9d62b8..82c8ae6328 100644 --- a/include/deal.II/dofs/dof_tools.h +++ b/include/deal.II/dofs/dof_tools.h @@ -1184,7 +1184,7 @@ namespace DoFTools * @name Identifying subsets of degrees of freedom with particular properties * @{ */ - + /** * Extract the indices of the degrees of freedom belonging to * certain vector components of a vector-valued finite element. The @@ -1436,7 +1436,7 @@ namespace DoFTools extract_constant_modes (const DH &dof_handler, const ComponentMask &component_mask, std::vector > &constant_modes); - + /** * @} */ @@ -1887,7 +1887,7 @@ namespace DoFTools std::vector &dofs_per_component, std::vector target_component) DEAL_II_DEPRECATED; - + /** * For each active cell of a DoFHandler or hp::DoFHandler, extract * the active finite element index and fill the vector given as @@ -1992,7 +1992,7 @@ namespace DoFTools /** * @} */ - + /** * Create a mapping from degree of freedom indices to the index of * that degree of freedom on the boundary. After this operation, diff --git a/include/deal.II/fe/component_mask.h b/include/deal.II/fe/component_mask.h index 04099c9a14..9f4eaff2c3 100644 --- a/include/deal.II/fe/component_mask.h +++ b/include/deal.II/fe/component_mask.h @@ -102,7 +102,7 @@ public: * Set a particular entry in the mask to a value. */ void set (const unsigned int index, const bool value); - + /** * If this component mask has been initialized with a mask of * size greater than zero, then return the size of the mask diff --git a/include/deal.II/fe/fe_dg_vector.templates.h b/include/deal.II/fe/fe_dg_vector.templates.h index b404e45703..09bef38ae5 100644 --- a/include/deal.II/fe/fe_dg_vector.templates.h +++ b/include/deal.II/fe/fe_dg_vector.templates.h @@ -81,7 +81,7 @@ FE_DGVector::get_dpo_vector (const unsigned int deg) template bool FE_DGVector::has_support_on_face (const unsigned int, - const unsigned int) const + const unsigned int) const { return true; } diff --git a/include/deal.II/fe/fe_face.h b/include/deal.II/fe/fe_face.h index b94889963d..b3dc5c9915 100644 --- a/include/deal.II/fe/fe_face.h +++ b/include/deal.II/fe/fe_face.h @@ -236,12 +236,12 @@ protected: typename Mapping<1,spacedim>::InternalDataBase * get_face_data (const UpdateFlags, const Mapping<1,spacedim> &mapping, - const Quadrature<0>& quadrature) const ; + const Quadrature<0> &quadrature) const ; typename Mapping<1,spacedim>::InternalDataBase * get_subface_data (const UpdateFlags, const Mapping<1,spacedim> &mapping, - const Quadrature<0>& quadrature) const ; + const Quadrature<0> &quadrature) const ; virtual void fill_fe_values (const Mapping<1,spacedim> &mapping, diff --git a/include/deal.II/fe/fe_poly_face.templates.h b/include/deal.II/fe/fe_poly_face.templates.h index e84a3959e2..72488f6576 100644 --- a/include/deal.II/fe/fe_poly_face.templates.h +++ b/include/deal.II/fe/fe_poly_face.templates.h @@ -210,42 +210,42 @@ FE_PolyFace::fill_fe_face_values ( if (flags & update_values) for (unsigned int i=0; idofs_per_cell; ++k) + for (unsigned int k=0; kdofs_per_cell; ++k) data.shape_values(k,i) = 0.; - switch (dim) - { - case 3: - { - // Fill data for quad shape functions - if (this->dofs_per_quad !=0) - { - const unsigned int foffset = this->first_quad_index + this->dofs_per_quad * face; - for (unsigned int k=0; kdofs_per_quad; ++k) - data.shape_values(foffset+k,i) = fe_data.shape_values[k+this->first_face_quad_index][i]; - } - } - case 2: - { - // Fill data for line shape functions - if (this->dofs_per_line != 0) - { - const unsigned int foffset = this->first_line_index; - for (unsigned int line=0; line::lines_per_face; ++line) - { - for (unsigned int k=0; kdofs_per_line; ++k) - data.shape_values(foffset+GeometryInfo::face_to_cell_lines(face, line)*this->dofs_per_line+k,i) = fe_data.shape_values[k+(line*this->dofs_per_line)+this->first_face_line_index][i]; - } - } - } - case 1: - { - // Fill data for vertex shape functions - if (this->dofs_per_vertex != 0) - for (unsigned int lvertex=0; lvertex::vertices_per_face; ++lvertex) - data.shape_values(GeometryInfo::face_to_cell_vertices(face, lvertex),i) = fe_data.shape_values[lvertex][i]; - break; - } - } + switch (dim) + { + case 3: + { + // Fill data for quad shape functions + if (this->dofs_per_quad !=0) + { + const unsigned int foffset = this->first_quad_index + this->dofs_per_quad * face; + for (unsigned int k=0; kdofs_per_quad; ++k) + data.shape_values(foffset+k,i) = fe_data.shape_values[k+this->first_face_quad_index][i]; + } + } + case 2: + { + // Fill data for line shape functions + if (this->dofs_per_line != 0) + { + const unsigned int foffset = this->first_line_index; + for (unsigned int line=0; line::lines_per_face; ++line) + { + for (unsigned int k=0; kdofs_per_line; ++k) + data.shape_values(foffset+GeometryInfo::face_to_cell_lines(face, line)*this->dofs_per_line+k,i) = fe_data.shape_values[k+(line*this->dofs_per_line)+this->first_face_line_index][i]; + } + } + } + case 1: + { + // Fill data for vertex shape functions + if (this->dofs_per_vertex != 0) + for (unsigned int lvertex=0; lvertex::vertices_per_face; ++lvertex) + data.shape_values(GeometryInfo::face_to_cell_vertices(face, lvertex),i) = fe_data.shape_values[lvertex][i]; + break; + } + } } } diff --git a/include/deal.II/fe/fe_values.h b/include/deal.II/fe/fe_values.h index 0b1fc6d3f3..04f6f4c5c8 100644 --- a/include/deal.II/fe/fe_values.h +++ b/include/deal.II/fe/fe_values.h @@ -2207,7 +2207,7 @@ public: * @dealiiRequiresUpdateFlags{update_inverse_jacobians} */ const std::vector > &get_inverse_jacobians () const; - + /** * For a face, return the outward normal vector to the cell at the * ith quadrature point. diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index a89b4d88d5..a65d97cc86 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -82,7 +82,7 @@ namespace GridGenerator void hyper_cube (Triangulation &tria, const double left = 0., const double right= 1., - const bool colorize= false); + const bool colorize= false); /** * Same as hyper_cube(), but with the difference that not only one cell is @@ -797,7 +797,7 @@ namespace GridGenerator */ template void flatten_triangulation(const Triangulation &in_tria, - Triangulation &out_tria); + Triangulation &out_tria); /** * @} @@ -819,10 +819,10 @@ namespace GridGenerator template