From: bangerth Date: Thu, 23 Aug 2012 03:58:24 +0000 (+0000) Subject: Merge up from the upstream branch. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=94b73a2b7d0c7cce71244a46d18c06675a01c00d;p=dealii-svn.git Merge up from the upstream branch. git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_3@26096 0785d39b-7218-0410-832d-ea1e28bc413d --- 94b73a2b7d0c7cce71244a46d18c06675a01c00d diff --cc deal.II/contrib/boost-1.49.0/include/boost/fusion/tuple/detail/preprocessed/tuple.hpp index 8f9b94a5cd,8f9b94a5cd..e1cd851940 --- a/deal.II/contrib/boost-1.49.0/include/boost/fusion/tuple/detail/preprocessed/tuple.hpp +++ b/deal.II/contrib/boost-1.49.0/include/boost/fusion/tuple/detail/preprocessed/tuple.hpp @@@ -18,4 -18,4 +18,4 @@@ #include #else #error "FUSION_MAX_VECTOR_SIZE out of bounds for preprocessed headers" --#endif ++#endif diff --cc deal.II/examples/step-42/step-42.cc index 8fbb364edc,3a9a898ef8..54bdc04736 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@@ -433,1053 -476,921 +476,920 @@@ namespace Step4 } template - void Obstacle::vector_value (const Point &p, - Vector &values) const + void PlasticityContactProblem::make_grid () { - for (unsigned int c=0; cn_components; ++c) - values(c) = Obstacle::value (p, c); + std::vector repet(3); + repet[0] = 1;//20; + repet[1] = 1; + repet[2] = 1; + + Point p1 (0,0,0); + Point p2 (1.0, 1.0, 1.0); + GridGenerator::subdivided_hyper_rectangle (triangulation, repet, p1, p2); + + Triangulation<3>::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + + /* boundary_indicators: + _______ + / 9 /| + /______ / | + 8| | 8| + | 8 | / + |_______|/ + 6 + */ + + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + { + if (cell->face (face)->center ()[2] == p2(2)) + cell->face (face)->set_boundary_indicator (9); + if (cell->face (face)->center ()[0] == p1(0) || + cell->face (face)->center ()[0] == p2(0) || + cell->face (face)->center ()[1] == p1(1) || + cell->face (face)->center ()[1] == p2(1)) + cell->face (face)->set_boundary_indicator (8); + if (cell->face (face)->center ()[2] == p1(2)) + cell->face (face)->set_boundary_indicator (6); + } + + triangulation.refine_global (n_refinements_global); + + // Lokale Verfeinerung des Gitters + for (int step=0; step::faces_per_cell; ++face) + { + // if (cell->face (face)->at_boundary() + // && cell->face (face)->boundary_indicator () == 9) + // { + // cell->set_refine_flag (); + // break; + // } + // else if (cell->level () == n_refinements + n_refinements_local - 1) + // { + // cell->set_refine_flag (); + // break; + // } + + // if (cell->face (face)->at_boundary() + // && cell->face (face)->boundary_indicator () == 9) + // { + // if (cell->face (face)->vertex (0)(0) <= 0.7 && + // cell->face (face)->vertex (1)(0) >= 0.3 && + // cell->face (face)->vertex (0)(1) <= 0.875 && + // cell->face (face)->vertex (2)(1) >= 0.125) + // { + // cell->set_refine_flag (); + // break; + // } + // } + + if (step == 0 && + cell->center ()(2) < n_refinements_local*9.0/64.0) + { + cell->set_refine_flag (); + break; + } + }; + triangulation.execute_coarsening_and_refinement (); + }; } - } + template + void PlasticityContactProblem::setup_system () + { + // setup dofs + { + dof_handler.distribute_dofs (fe); - // @sect3{Implementation of the Step4 class} - - // Next for the implementation of the class - // template that makes use of the functions - // above. As before, we will write everything - - template - Step4::Step4 (int _n_refinements_global, int _n_refinements_local) - : - n_refinements_global (_n_refinements_global), - n_refinements_local (_n_refinements_local), - mpi_communicator (MPI_COMM_WORLD), - triangulation (mpi_communicator), - fe (FE_Q(1), dim), - dof_handler (triangulation), - pcout (std::cout, - (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), - sigma_0 (400), - gamma (1.e-2), - e_modul (2.0e5), - nu (0.3) - { - // double _E, double _nu, double _sigma_0, double _gamma - plast_lin_hard.reset (new ConstitutiveLaw (e_modul, nu, sigma_0, gamma, mpi_communicator, pcout)); - } + locally_owned_dofs = dof_handler.locally_owned_dofs (); + locally_relevant_dofs.clear(); + DoFTools::extract_locally_relevant_dofs (dof_handler, + locally_relevant_dofs); + } - template - void Step4::make_grid () - { - std::vector repet(3); - repet[0] = 1;//20; - repet[1] = 1; - repet[2] = 1; - - Point p1 (0,0,0); - Point p2 (1.0, 1.0, 1.0); - GridGenerator::subdivided_hyper_rectangle (triangulation, repet, p1, p2); - - Triangulation<3>::active_cell_iterator - cell = triangulation.begin_active(), - endc = triangulation.end(); - - /* boundary_indicators: - _______ - / 9 /| - /______ / | - 8| | 8| - | 8 | / - |_______|/ - 6 - */ - - for (; cell!=endc; ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) - { - if (cell->face (face)->center ()[2] == p2(2)) - cell->face (face)->set_boundary_indicator (9); - if (cell->face (face)->center ()[0] == p1(0) || - cell->face (face)->center ()[0] == p2(0) || - cell->face (face)->center ()[1] == p1(1) || - cell->face (face)->center ()[1] == p2(1)) - cell->face (face)->set_boundary_indicator (8); - if (cell->face (face)->center ()[2] == p1(2)) - cell->face (face)->set_boundary_indicator (6); - } + // setup hanging nodes and dirichlet constraints + { + constraints_hanging_nodes.reinit (locally_relevant_dofs); + DoFTools::make_hanging_node_constraints (dof_handler, + constraints_hanging_nodes); + constraints_hanging_nodes.close (); + + pcout << "Number of active cells: " + << triangulation.n_active_cells() + << std::endl + << "Total number of cells: " + << triangulation.n_cells() + << std::endl + << "Number of degrees of freedom: " + << dof_handler.n_dofs () + << std::endl; + + dirichlet_constraints (); + } - triangulation.refine_global (n_refinements_global); + // Initialzation for matrices and vectors + { + solution.reinit (locally_relevant_dofs, mpi_communicator); + system_rhs_newton.reinit (locally_owned_dofs, mpi_communicator); + old_solution.reinit (system_rhs_newton); + resid_vector.reinit (system_rhs_newton); + diag_mass_matrix_vector.reinit (system_rhs_newton); + active_set.set_size (locally_relevant_dofs.size ()); + } - // Lokale Verfeinerung des Gitters - for (int step=0; step::faces_per_cell; ++face) - { - // if (cell->face (face)->at_boundary() - // && cell->face (face)->boundary_indicator () == 9) - // { - // cell->set_refine_flag (); - // break; - // } - // else if (cell->level () == n_refinements + n_refinements_local - 1) - // { - // cell->set_refine_flag (); - // break; - // } - - // if (cell->face (face)->at_boundary() - // && cell->face (face)->boundary_indicator () == 9) - // { - // if (cell->face (face)->vertex (0)(0) <= 0.7 && - // cell->face (face)->vertex (1)(0) >= 0.3 && - // cell->face (face)->vertex (0)(1) <= 0.875 && - // cell->face (face)->vertex (2)(1) >= 0.125) - // { - // cell->set_refine_flag (); - // break; - // } - // } - - if (step == 0 && - cell->center ()(2) < n_refinements_local*9.0/64.0) - { - cell->set_refine_flag (); - break; - } - }; - triangulation.execute_coarsening_and_refinement (); - }; - } + TrilinosWrappers::SparsityPattern sp (locally_owned_dofs, + mpi_communicator); - template - void Step4::setup_system () - { - // setup dofs - { - dof_handler.distribute_dofs (fe); + DoFTools::make_sparsity_pattern (dof_handler, sp, constraints_dirichlet_hanging_nodes, false, + Utilities::MPI::this_mpi_process(mpi_communicator)); + + sp.compress(); + + system_matrix_newton.reinit (sp); - locally_owned_dofs = dof_handler.locally_owned_dofs (); - locally_relevant_dofs.clear(); - DoFTools::extract_locally_relevant_dofs (dof_handler, - locally_relevant_dofs); + TrilinosWrappers::SparseMatrix mass_matrix; + mass_matrix.reinit (sp); + assemble_mass_matrix_diagonal (mass_matrix); + const unsigned int + start = (system_rhs_newton.local_range().first), + end = (system_rhs_newton.local_range().second); + for (unsigned int j=start; j + void PlasticityContactProblem::assemble_nl_system (TrilinosWrappers::MPI::Vector &u) { - // constraints_hanging_nodes.clear (); - constraints_hanging_nodes.reinit (locally_relevant_dofs); - DoFTools::make_hanging_node_constraints (dof_handler, - constraints_hanging_nodes); - constraints_hanging_nodes.close (); - - pcout << "Number of active cells: " - << triangulation.n_active_cells() - << std::endl - << "Total number of cells: " - << triangulation.n_cells() - << std::endl - << "Number of degrees of freedom: " - << dof_handler.n_dofs () - << std::endl; - - dirichlet_constraints (); + QGauss quadrature_formula(2); + QGauss face_quadrature_formula(2); + + FEValues fe_values (fe, quadrature_formula, + UpdateFlags(update_values | + update_gradients | + update_q_points | + update_JxW_values)); + + FEFaceValues fe_values_face (fe, face_quadrature_formula, + update_values | update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size (); + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + const EquationData::RightHandSide right_hand_side; + std::vector > right_hand_side_values (n_q_points, + Vector(dim)); + std::vector > right_hand_side_values_face (n_face_q_points, + Vector(dim)); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + - std::vector local_dof_indices (dofs_per_cell); - ++ std::vector local_dof_indices (dofs_per_cell); + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + const FEValuesExtractors::Vector displacement (0); + + TrilinosWrappers::MPI::Vector test_rhs(solution); + const double kappa = 1.0; + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + fe_values.reinit (cell); + cell_matrix = 0; + cell_rhs = 0; + + right_hand_side.vector_value_list (fe_values.get_quadrature_points(), + right_hand_side_values); + + std::vector > strain_tensor (n_q_points); + fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor); + + for (unsigned int q_point=0; q_point stress_strain_tensor_linearized; + SymmetricTensor<4,dim> stress_strain_tensor; + SymmetricTensor<2,dim> stress_tensor; + + plast_lin_hard->linearized_plast_linear_hardening (stress_strain_tensor_linearized, + stress_strain_tensor, + strain_tensor[q_point]); + + // if (q_point == 0) + // std::cout<< stress_strain_tensor_linearized <get_strain(fe_values, i, q_point); + + for (unsigned int j=0; jget_strain(fe_values, j, q_point) * + fe_values.JxW (q_point)); + } + + // the linearized part a(v^i;v^i,v) of the rhs + cell_rhs(i) += (stress_tensor * + strain_tensor[q_point] * + fe_values.JxW (q_point)); + + // the residual part a(v^i;v) of the rhs + cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor * + plast_lin_hard->get_strain(fe_values, i, q_point) * + fe_values.JxW (q_point)); + + // the residual part F(v) of the rhs + Tensor<1,dim> rhs_values; + rhs_values = 0; + cell_rhs(i) += (fe_values[displacement].value (i, q_point) * + rhs_values * + fe_values.JxW (q_point)); + } + } + + for (unsigned int face=0; face::faces_per_cell; ++face) + { + if (cell->face (face)->at_boundary() + && cell->face (face)->boundary_indicator () == 9) + { + fe_values_face.reinit (cell, face); + + right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(), + right_hand_side_values_face); + + for (unsigned int q_point=0; q_point rhs_values; + rhs_values = 0; + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, cell_rhs, + local_dof_indices, + system_matrix_newton, system_rhs_newton, true); + }; + + system_matrix_newton.compress (); + system_rhs_newton.compress (Add); } - // Initialzation for matrices and vectors + template + void PlasticityContactProblem::residual_nl_system (TrilinosWrappers::MPI::Vector &u) { - solution.reinit (locally_relevant_dofs, mpi_communicator); - system_rhs_newton.reinit (locally_owned_dofs, mpi_communicator); - old_solution.reinit (system_rhs_newton); - resid_vector.reinit (system_rhs_newton); - diag_mass_matrix_vector.reinit (system_rhs_newton); - active_set.set_size (locally_relevant_dofs.size ()); + QGauss quadrature_formula(2); + QGauss face_quadrature_formula(2); + + FEValues fe_values (fe, quadrature_formula, + UpdateFlags(update_values | + update_gradients | + update_q_points | + update_JxW_values)); + + FEFaceValues fe_values_face (fe, face_quadrature_formula, + update_values | update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size (); + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + const EquationData::RightHandSide right_hand_side; + std::vector > right_hand_side_values (n_q_points, + Vector(dim)); + std::vector > right_hand_side_values_face (n_face_q_points, + Vector(dim)); + + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const FEValuesExtractors::Vector displacement (0); + + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + unsigned int elast_points = 0; + unsigned int plast_points = 0; + double yield = 0; + unsigned int cell_number = 0; + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + fe_values.reinit (cell); + cell_rhs = 0; + + right_hand_side.vector_value_list (fe_values.get_quadrature_points(), + right_hand_side_values); + + std::vector > strain_tensor (n_q_points); + fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor); + + for (unsigned int q_point=0; q_point stress_strain_tensor; + SymmetricTensor<2,dim> stress_tensor; + + plast_lin_hard->plast_linear_hardening (stress_strain_tensor, strain_tensor[q_point], + elast_points, plast_points, yield); + + for (unsigned int i=0; iget_strain(fe_values, i, q_point) * + fe_values.JxW (q_point)); + + Tensor<1,dim> rhs_values; + rhs_values = 0; + cell_rhs(i) += ((fe_values[displacement].value (i, q_point) * + rhs_values) * + fe_values.JxW (q_point)); + }; + }; + + for (unsigned int face=0; face::faces_per_cell; ++face) + { + if (cell->face (face)->at_boundary() + && cell->face (face)->boundary_indicator () == 9) + { + fe_values_face.reinit (cell, face); + + right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(), + right_hand_side_values_face); + + for (unsigned int q_point=0; q_point rhs_values; + rhs_values = 0; + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_rhs, + local_dof_indices, + system_rhs_newton); + + cell_number += 1; + }; + + system_rhs_newton.compress (); + + unsigned int sum_elast_points = Utilities::MPI::sum(elast_points, mpi_communicator); + unsigned int sum_plast_points = Utilities::MPI::sum(plast_points, mpi_communicator); + pcout<< "Elast-Points = " << sum_elast_points < + void PlasticityContactProblem::dirichlet_constraints () + { + /* boundary_indicators: + _______ + / 9 /| + /______ / | + 8| | 8| + | 8 | / + |_______|/ + 6 + */ + + constraints_dirichlet_hanging_nodes.reinit (locally_relevant_dofs); + constraints_dirichlet_hanging_nodes.merge (constraints_hanging_nodes); + + std::vector component_mask (dim, true); + component_mask[0] = true; + component_mask[1] = true; + component_mask[2] = true; + VectorTools::interpolate_boundary_values (dof_handler, + 6, + EquationData::BoundaryValues(), + constraints_dirichlet_hanging_nodes, + component_mask); + + component_mask[0] = true; + component_mask[1] = true; + component_mask[2] = false; + VectorTools::interpolate_boundary_values (dof_handler, + 8, + EquationData::BoundaryValues(), + constraints_dirichlet_hanging_nodes, + component_mask); + constraints_dirichlet_hanging_nodes.close (); + } - const EquationData::RightHandSide right_hand_side; - std::vector > right_hand_side_values (n_q_points, - Vector(dim)); - std::vector > right_hand_side_values_face (n_face_q_points, - Vector(dim)); + template + void PlasticityContactProblem::solve () + { + pcout << "Solving ..." << std::endl; + Timer t; - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); + TrilinosWrappers::MPI::Vector distributed_solution (system_rhs_newton); + distributed_solution = solution; - std::vector local_dof_indices (dofs_per_cell); + constraints_hanging_nodes.set_zero (distributed_solution); - typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + // Solving iterative - const FEValuesExtractors::Vector displacement (0); + MPI_Barrier (mpi_communicator); + t.restart(); - TrilinosWrappers::MPI::Vector test_rhs(solution); - const double kappa = 1.0; - for (; cell!=endc; ++cell) - if (cell->is_locally_owned()) - { - fe_values.reinit (cell); - cell_matrix = 0; - cell_rhs = 0; - - right_hand_side.vector_value_list (fe_values.get_quadrature_points(), - right_hand_side_values); - - std::vector > strain_tensor (n_q_points); - fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor); - - for (unsigned int q_point=0; q_point stress_strain_tensor_linearized; - SymmetricTensor<4,dim> stress_strain_tensor; - SymmetricTensor<2,dim> stress_tensor; - - plast_lin_hard->linearized_plast_linear_hardening (stress_strain_tensor_linearized, - stress_strain_tensor, - strain_tensor[q_point]); - - // if (q_point == 0) - // std::cout<< stress_strain_tensor_linearized <get_strain(fe_values, i, q_point); - - for (unsigned int j=0; jget_strain(fe_values, j, q_point) * - fe_values.JxW (q_point)); - } - - // the linearized part a(v^i;v^i,v) of the rhs - cell_rhs(i) += (stress_tensor * - strain_tensor[q_point] * - fe_values.JxW (q_point)); - - // the residual part a(v^i;v) of the rhs - cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor * - plast_lin_hard->get_strain(fe_values, i, q_point) * - fe_values.JxW (q_point)); - - // the residual part F(v) of the rhs - Tensor<1,dim> rhs_values; - rhs_values = 0; - cell_rhs(i) += (fe_values[displacement].value (i, q_point) * - rhs_values * - fe_values.JxW (q_point)); - } - } - - for (unsigned int face=0; face::faces_per_cell; ++face) - { - if (cell->face (face)->at_boundary() - && cell->face (face)->boundary_indicator () == 9) - { - fe_values_face.reinit (cell, face); - - right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(), - right_hand_side_values_face); - - for (unsigned int q_point=0; q_point rhs_values; - rhs_values = 0; - for (unsigned int i=0; iget_dof_indices (local_dof_indices); - constraints.distribute_local_to_global (cell_matrix, cell_rhs, - local_dof_indices, - system_matrix_newton, system_rhs_newton, true); - }; + preconditioner_u.initialize (system_matrix_newton, additional_data); - system_matrix_newton.compress (); - system_rhs_newton.compress (Add); - } + MPI_Barrier (mpi_communicator); + t.stop(); + if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + run_time[6] += t.wall_time(); - template - void Step4::residual_nl_system (TrilinosWrappers::MPI::Vector &u, - Vector &sigma_eff_vector) - { - QGauss quadrature_formula(2); - QGauss face_quadrature_formula(2); - - FEValues fe_values (fe, quadrature_formula, - UpdateFlags(update_values | - update_gradients | - update_q_points | - update_JxW_values)); - - FEFaceValues fe_values_face (fe, face_quadrature_formula, - update_values | update_quadrature_points | - update_JxW_values); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size (); - const unsigned int n_face_q_points = face_quadrature_formula.size(); - - const EquationData::RightHandSide right_hand_side; - std::vector > right_hand_side_values (n_q_points, - Vector(dim)); - std::vector > right_hand_side_values_face (n_face_q_points, - Vector(dim)); - - Vector cell_rhs (dofs_per_cell); - Vector cell_sigma_eff (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - const FEValuesExtractors::Vector displacement (0); - - typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - unsigned int elast_points = 0; - unsigned int plast_points = 0; - double sigma_eff = 0; - double yield = 0; - unsigned int cell_number = 0; - for (; cell!=endc; ++cell) - if (cell->is_locally_owned()) - { - fe_values.reinit (cell); - cell_rhs = 0; - - right_hand_side.vector_value_list (fe_values.get_quadrature_points(), - right_hand_side_values); - - std::vector > strain_tensor (n_q_points); - fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor); - - for (unsigned int q_point=0; q_point stress_strain_tensor; - SymmetricTensor<2,dim> stress_tensor; - - plast_lin_hard->plast_linear_hardening (stress_strain_tensor, strain_tensor[q_point], - elast_points, plast_points, sigma_eff, yield); - - // sigma_eff_vector (cell_number) += sigma_eff; - sigma_eff_vector (cell_number) += yield; - - /* if (q_point == 0) - std::cout<< stress_strain_tensor <get_strain(fe_values, i, q_point) * - fe_values.JxW (q_point)); - - Tensor<1,dim> rhs_values; - rhs_values = 0; - cell_rhs(i) += ((fe_values[displacement].value (i, q_point) * - rhs_values) * - fe_values.JxW (q_point)); - }; - }; - - for (unsigned int face=0; face::faces_per_cell; ++face) - { - if (cell->face (face)->at_boundary() - && cell->face (face)->boundary_indicator () == 9) - { - fe_values_face.reinit (cell, face); - - right_hand_side.vector_value_list (fe_values_face.get_quadrature_points(), - right_hand_side_values_face); - - for (unsigned int q_point=0; q_point rhs_values; - rhs_values = 0; - for (unsigned int i=0; iget_dof_indices (local_dof_indices); - constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_rhs, - local_dof_indices, - system_rhs_newton); - - sigma_eff_vector(cell_number) /= n_q_points; - cell_number += 1; - }; + MPI_Barrier (mpi_communicator); + t.restart(); - system_rhs_newton.compress (); + // ReductionControl reduction_control (10000, 1e-15, 1e-4); + // SolverCG + // solver (reduction_control, mpi_communicator); + // solver.solve (system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u); - unsigned int sum_elast_points = Utilities::MPI::sum(elast_points, mpi_communicator); - unsigned int sum_plast_points = Utilities::MPI::sum(plast_points, mpi_communicator); - pcout<< "Elast-Points = " << sum_elast_points < - void Step4::projection_active_set () - { - clock_t start_proj, end_proj; - - const EquationData::Obstacle obstacle; - std::vector vertex_touched (dof_handler.n_dofs (), false); - - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - TrilinosWrappers::MPI::Vector distributed_solution (system_rhs_newton); - distributed_solution = solution; - TrilinosWrappers::MPI::Vector lambda (solution); - lambda = resid_vector; - TrilinosWrappers::MPI::Vector diag_mass_matrix_vector_relevant (solution); - diag_mass_matrix_vector_relevant = diag_mass_matrix_vector; - - constraints.reinit(locally_relevant_dofs); - active_set.clear (); - IndexSet active_set_locally_owned; - active_set_locally_owned.set_size (locally_owned_dofs.size ()); - const double c = 100.0*e_modul; - - for (; cell!=endc; ++cell) - if (cell->is_locally_owned()) - for (unsigned int face=0; face::faces_per_cell; ++face) - if (cell->face (face)->at_boundary() - && cell->face (face)->boundary_indicator () == 9) - for (unsigned int v=0; v::vertices_per_cell; ++v) - { - unsigned int index_z = cell->face (face)->vertex_dof_index (v,2); - - if (vertex_touched[cell->face (face)->vertex_index(v)] == false) - vertex_touched[cell->face (face)->vertex_index(v)] = true; - else - continue; - - // the local row where - Point point (cell->face (face)->vertex (v)[0],/* + solution (index_x),*/ - cell->face (face)->vertex (v)[1], - cell->face (face)->vertex (v)[2]); - - double obstacle_value = obstacle.value (point, 2); - double solution_index_z = solution (index_z); - double gap = obstacle_value - point (2); - - if (lambda (index_z) + - c*diag_mass_matrix_vector_relevant (index_z)*(solution_index_z - gap) > 0) - { - constraints.add_line (index_z); - constraints.set_inhomogeneity (index_z, gap); - - distributed_solution (index_z) = gap; - - if (locally_relevant_dofs.is_element (index_z)) - active_set.add_index (index_z); - - if (locally_owned_dofs.is_element (index_z)) - active_set_locally_owned.add_index (index_z); - - // std::cout<< index_z << ", " - // << "Error: " << lambda (index_z) + - // diag_mass_matrix_vector_relevant (index_z)*c*(solution_index_z - gap) - // << ", " << lambda (index_z) - // << ", " << diag_mass_matrix_vector_relevant (index_z) - // << ", " << obstacle_value - // << ", " << solution_index_z - // < - void Step4::dirichlet_constraints () - { - /* boundary_indicators: - _______ - / 9 /| - /______ / | - 8| | 8| - | 8 | / - |_______|/ - 6 - */ - - constraints_dirichlet_hanging_nodes.reinit (locally_relevant_dofs); - constraints_dirichlet_hanging_nodes.merge (constraints_hanging_nodes); - - std::vector component_mask (dim, true); - component_mask[0] = true; - component_mask[1] = true; - component_mask[2] = true; - VectorTools::interpolate_boundary_values (dof_handler, - 6, - EquationData::BoundaryValues(), - constraints_dirichlet_hanging_nodes, - component_mask); - - component_mask[0] = true; - component_mask[1] = true; - component_mask[2] = false; - VectorTools::interpolate_boundary_values (dof_handler, - 8, - EquationData::BoundaryValues(), - constraints_dirichlet_hanging_nodes, - component_mask); - constraints_dirichlet_hanging_nodes.close (); - } + number_iterations += solver_control.last_step(); - template - void Step4::solve () - { - pcout << "Solving ..." << std::endl; - Timer t; - - TrilinosWrappers::MPI::Vector distributed_solution (system_rhs_newton); - distributed_solution = solution; - - constraints_hanging_nodes.set_zero (distributed_solution); - - // Solving iterative - - MPI_Barrier (mpi_communicator); - t.restart(); - - preconditioner_u.initialize (system_matrix_newton, additional_data); - - MPI_Barrier (mpi_communicator); - t.stop(); - if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) - run_time[6] += t.wall_time(); - - MPI_Barrier (mpi_communicator); - t.restart(); - - // ReductionControl reduction_control (10000, 1e-15, 1e-4); - // SolverCG - // solver (reduction_control, mpi_communicator); - // solver.solve (system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u); - - PrimitiveVectorMemory mem; - TrilinosWrappers::MPI::Vector tmp (system_rhs_newton); - const double solver_tolerance = 1e-4 * - system_matrix_newton.residual (tmp, distributed_solution, system_rhs_newton); - SolverControl solver_control (system_matrix_newton.m(), solver_tolerance); - SolverFGMRES - solver(solver_control, mem, - SolverFGMRES:: - AdditionalData(30, true)); - solver.solve(system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u); - - pcout << "Initial error: " << solver_control.initial_value() < - void Step4::solve_newton () - { - double resid=0; - double resid_old=100000; - TrilinosWrappers::MPI::Vector res (system_rhs_newton); - TrilinosWrappers::MPI::Vector tmp_vector (system_rhs_newton); - Timer t; - - std::vector > constant_modes; - std::vector components (dim,true); - components[dim] = false; - DoFTools::extract_constant_modes (dof_handler, components, - constant_modes); - - additional_data.elliptic = true; - additional_data.n_cycles = 1; - additional_data.w_cycle = false; - additional_data.output_details = false; - additional_data.smoother_sweeps = 2; - additional_data.aggregation_threshold = 1e-2; - - IndexSet active_set_old (active_set); - Vector sigma_eff_vector; - sigma_eff_vector.reinit (triangulation.n_active_cells()); - unsigned int j = 0; - unsigned int number_assemble_system = 0; - for (; j<=100;j++) - { - pcout<< " " <(i)); - old_solution = tmp_vector; - old_solution.sadd(1-a,a, distributed_solution); - - MPI_Barrier (mpi_communicator); - t.restart(); - system_rhs_newton = 0; - sigma_eff_vector = 0; - solution = old_solution; - residual_nl_system (solution, sigma_eff_vector); - res = system_rhs_newton; - - const unsigned int - start_res = (res.local_range().first), - end_res = (res.local_range().second); - for (unsigned int n=start_res; n(i)); + old_solution = tmp_vector; + old_solution.sadd(1-a,a, distributed_solution); + + MPI_Barrier (mpi_communicator); + t.restart(); + system_rhs_newton = 0; + solution = old_solution; + residual_nl_system (solution); + res = system_rhs_newton; + + const unsigned int + start_res = (res.local_range().first), + end_res = (res.local_range().second); + for (unsigned int n=start_res; n::type_dof_data, - data_component_interpretation); - data_out.add_data_vector (lambda, std::vector(dim, "Residual"), - DataOut::type_dof_data, - data_component_interpretation); - data_out.add_data_vector (active_set, std::vector(dim, "ActiveSet"), - DataOut::type_dof_data, - data_component_interpretation); + DataOut data_out; - Vector subdomain (triangulation.n_active_cells()); - for (unsigned int i=0; i + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_out.add_data_vector (solution, std::vector(dim, "Displacement"), + DataOut::type_dof_data, + data_component_interpretation); + data_out.add_data_vector (lambda, std::vector(dim, "Residual"), + DataOut::type_dof_data, + data_component_interpretation); + data_out.add_data_vector (active_set, std::vector(dim, "ActiveSet"), + DataOut::type_dof_data, + data_component_interpretation); - const std::string filename = (title + "-" + - Utilities::int_to_string - (triangulation.locally_owned_subdomain(), 4)); + Vector subdomain (triangulation.n_active_cells()); + for (unsigned int i=0; i filenames; - for (unsigned int i=0; - i - void Step4::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const - { - pcout<< "Moving mesh." < vertex_touched (triangulation.n_vertices(), - false); - - for (typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active (); - cell != dof_handler.end(); ++cell) - if (cell->is_locally_owned()) - for (unsigned int v=0; v::vertices_per_cell; ++v) - { - if (vertex_touched[cell->vertex_index(v)] == false) - { - vertex_touched[cell->vertex_index(v)] = true; - - Point vertex_displacement; - for (unsigned int d=0; dvertex_dof_index(v,d)) != 0) - vertex_displacement[d] - = _complete_displacement(cell->vertex_dof_index(v,d)); - } - - cell->vertex(v) += vertex_displacement; - } - } - } + if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + { + std::vector filenames; + for (unsigned int i=0; + i - void Step4::output_results (TrilinosWrappers::MPI::Vector vector, - const std::string& title) const - { - DataOut data_out; + TrilinosWrappers::MPI::Vector tmp (solution); + tmp *= -1; + move_mesh (tmp); + } - data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (vector, "vector_to_plot"); + template + void PlasticityContactProblem::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const + { + pcout<< "Moving mesh." < vertex_touched (triangulation.n_vertices(), + false); + + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (); + cell != dof_handler.end(); ++cell) + if (cell->is_locally_owned()) + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + if (vertex_touched[cell->vertex_index(v)] == false) + { + vertex_touched[cell->vertex_index(v)] = true; + + Point vertex_displacement; + for (unsigned int d=0; dvertex_dof_index(v,d)) != 0) + vertex_displacement[d] + = _complete_displacement(cell->vertex_dof_index(v,d)); + } + + cell->vertex(v) += vertex_displacement; + } + } + } - data_out.build_patches (); + template + void PlasticityContactProblem::output_results (TrilinosWrappers::MPI::Vector vector, + const std::string& title) const + { + DataOut data_out; - std::ofstream output_vtk (dim == 2 ? - (title + ".vtk").c_str () : - (title + ".vtk").c_str ()); - data_out.write_vtk (output_vtk); - } + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (vector, "vector_to_plot"); - template - void Step4::output_results (Vector vector, const std::string& title) const - { - DataOut data_out; + data_out.build_patches (); - data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (vector, "vector_to_plot"); + std::ofstream output_vtk (dim == 2 ? + (title + ".vtk").c_str () : + (title + ".vtk").c_str ()); + data_out.write_vtk (output_vtk); + } - data_out.build_patches (); + template + void PlasticityContactProblem::run () + { + pcout << "Solving problem in " << dim << " space dimensions." << std::endl; - + - std::ofstream output_vtk (dim == 2 ? - (title + ".vtk").c_str () : - (title + ".vtk").c_str ()); - data_out.write_vtk (output_vtk); - } + run_time.resize (8); - + - template - void Step4::run () - { - pcout << "Solving problem in " << dim << " space dimensions." << std::endl; + clock_t start, end; - + - run_time.resize (8); + start = clock(); + make_grid(); + // mesh_surface (); - clock_t start, end; + setup_system (); - + - start = clock(); - make_grid(); - // mesh_surface (); - - setup_system (); - - end = clock(); - run_time[0] = (double)(end-start)/CLOCKS_PER_SEC; + end = clock(); + run_time[0] = (double)(end-start)/CLOCKS_PER_SEC; - solve_newton (); + solve_newton (); + } } - // @sect3{The main function} - // And this is the main function. It also - // looks mostly like in step-3, but if you - // look at the code below, note how we first - // create a variable of type - // Step4@<2@> (forcing - // the compiler to compile the class template - // with dim replaced by - // 2) and run a 2d simulation, - // and then we do the whole thing over in 3d. - // - // In practice, this is probably not what you - // would do very frequently (you probably - // either want to solve a 2d problem, or one - // in 3d, but not both at the same - // time). However, it demonstrates the - // mechanism by which we can simply change - // which dimension we want in a single place, - // and thereby force the compiler to - // recompile the dimension independent class - // templates for the dimension we - // request. The emphasis here lies on the - // fact that we only need to change a single - // place. This makes it rather trivial to - // debug the program in 2d where computations - // are fast, and then switch a single place - // to a 3 to run the much more computing - // intensive program in 3d for `real' - // computations. - // - // Each of the two blocks is enclosed in - // braces to make sure that the - // laplace_problem_2d variable - // goes out of scope (and releases the memory - // it holds) before we move on to allocate - // memory for the 3d case. Without the - // additional braces, the - // laplace_problem_2d variable - // would only be destroyed at the end of the - // function, i.e. after running the 3d - // problem, and would needlessly hog memory - // while the 3d run could actually use it. - // - // Finally, the first line of the function is - // used to suppress some output. Remember - // that in the previous example, we had the - // output from the linear solvers about the - // starting residual and the number of the - // iteration where convergence was - // detected. This can be suppressed through - // the deallog.depth_console(0) - // call. - // - // The rationale here is the following: the - // deallog (i.e. deal-log, not de-allog) - // variable represents a stream to which some - // parts of the library write output. It - // redirects this output to the console and - // if required to a file. The output is - // nested in a way so that each function can - // use a prefix string (separated by colons) - // for each line of output; if it calls - // another function, that may also use its - // prefix which is then printed after the one - // of the calling function. Since output from - // functions which are nested deep below is - // usually not as important as top-level - // output, you can give the deallog variable - // a maximal depth of nested output for - // output to console and file. The depth zero - // which we gave here means that no output is - // written. By changing it you can get more - // information about the innards of the - // library. int main (int argc, char *argv[]) { + using namespace dealii; + using namespace Step42; + deallog.depth_console (0); clock_t start, end; diff --cc deal.II/include/deal.II/base/types.h index 80997dd2c1,be5ecda538..9d6bdabbbb --- a/deal.II/include/deal.II/base/types.h +++ b/deal.II/include/deal.II/base/types.h @@@ -63,60 -69,20 +69,60 @@@ namespace type * the @ref distributed module for * more information. */ - const unsigned int artificial_subdomain_id = static_cast(-2); + const unsigned int artificial_subdomain_id = static_cast(-2); - /** - * The type used to denote global dof - * indices. - */ + +#ifdef DEAL_II_USE_LARGE_INDEX_TYPE + /** + * The type used for global indices of + * degrees of freedom. While in sequential + * computations the 4 billion indices of + * 32-bit unsigned integers is plenty, + * parallel computations using the + * parallel::distributed::Triangulation + * class can overflow this number and we + * need a bigger index space. + * + * The data type always indicates an + * unsigned integer type. + */ + typedef unsigned long long int global_dof_index; + + /** + * An identifier that denotes the MPI type + * associated with types::global_dof_index. + */ +# define DEAL_II_DOF_INDEX_MPI_TYPE MPI_UNSIGNED_LONG_LONG +#else + /** + * The type used for global indices of + * degrees of freedom. While in sequential + * computations the 4 billion indices of + * 32-bit unsigned integers is plenty, + * parallel computations using the + * parallel::distributed::Triangulation + * class can overflow this number and we + * need a bigger index space. + * + * The data type always indicates an + * unsigned integer type. + */ typedef unsigned int global_dof_index; - /** - * An invalid value for indices of degrees - * of freedom. - */ + /** + * An identifier that denotes the MPI type + * associated with types::global_dof_index. + */ +# define DEAL_II_DOF_INDEX_MPI_TYPE MPI_UNSIGNED +#endif + + /** + * An invalid value for indices of degrees + * of freedom. + */ const global_dof_index invalid_dof_index = static_cast(-1); + /** * The type used to denote boundary indicators associated with every * piece of the boundary and, in the case of meshes that describe @@@ -140,10 -118,8 +158,10 @@@ * Invalid material_id which we need in several places as a default value. * We assume that all material_ids lie in the range [0, invalid_material_id). */ - const material_id_t invalid_material_id = static_cast(-1); + const material_id invalid_material_id = static_cast(-1); + + } diff --cc deal.II/include/deal.II/dofs/dof_handler.h index 6ae454508e,b236cdea75..d352b8d2f5 --- a/deal.II/include/deal.II/dofs/dof_handler.h +++ b/deal.II/include/deal.II/dofs/dof_handler.h @@@ -318,9 -318,7 +318,9 @@@ class DoFHandler : public Subscripto * object to the finite element. */ virtual void distribute_dofs (const FiniteElement &fe, - const unsigned int offset = 0); + const types::global_dof_index offset = 0); - ++ + virtual void distribute_mg_dofs (const FiniteElement& fe, const types::global_dof_index offset = 0); /** * After distribute_dofs() with @@@ -1042,9 -1040,7 +1042,9 @@@ * DoFs which are constrained by * hanging nodes, see @ref constraints. */ - unsigned int n_dofs () const; + types::global_dof_index n_dofs () const; - ++ + types::global_dof_index n_dofs (const unsigned int level) const; /** * Return the number of degrees of freedom @@@ -1076,8 -1072,8 +1076,8 @@@ * indicators under * consideration. */ - unsigned int + types::global_dof_index - n_boundary_dofs (const std::set &boundary_indicators) const; + n_boundary_dofs (const std::set &boundary_indicators) const; /** * Access to an object informing @@@ -1378,39 -1374,11 +1378,39 @@@ */ DoFHandler & operator = (const DoFHandler &); + class MGVertexDoFs { + private: + unsigned int coarsest_level; + unsigned int finest_level; + types::global_dof_index* indices; + types::global_dof_index* indices_offset; - ++ + public: + DeclException0 (ExcNoMemory); + MGVertexDoFs (); + ~MGVertexDoFs (); + unsigned int get_coarsest_level () const; + unsigned int get_finest_level () const; + types::global_dof_index get_index (const unsigned int level, const unsigned int dof_number) const; + void init (const unsigned int coarsest_level, const unsigned int finest_level, const unsigned int dofs_per_vertex); + void set_index (const unsigned int level, const unsigned int dof_number, const types::global_dof_index index); + }; - ++ + void clear_mg_space (); + /** * Free all used memory. */ void clear_space (); - + + void reserve_space (); - ++ + template + types::global_dof_index get_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index) const; + + template + void set_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index) const; + /** * Space to store the DoF numbers * for the different @@@ -1419,9 -1387,7 +1419,9 @@@ * Triangulation objects. */ std::vector*> levels; - + + std::vector*> mg_levels; + /** * Space to store DoF numbers of * faces. They are not stored in @@@ -1431,13 -1397,7 +1431,13 @@@ * array. */ dealii::internal::DoFHandler::DoFFaces *faces; - + + dealii::internal::DoFHandler::DoFFaces* mg_faces; - ++ + std::vector mg_vertex_dofs; - ++ + std::vector mg_used_dofs; + /** * Make accessor objects friends. */ @@@ -1457,9 -1417,9 +1457,9 @@@ #ifndef DOXYGEN -template <> unsigned int DoFHandler<1>::n_boundary_dofs () const; -template <> unsigned int DoFHandler<1>::n_boundary_dofs (const FunctionMap &) const; -template <> unsigned int DoFHandler<1>::n_boundary_dofs (const std::set &) const; +template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs () const; +template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs (const FunctionMap &) const; - template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set &) const; ++template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set &) const; /* ----------------------- Inline functions ---------------------------------- */ diff --cc deal.II/include/deal.II/dofs/dof_tools.h index fd7eda09b9,22a9d90963..b6f7b2a531 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@@ -493,10 -493,10 +493,10 @@@ namespace DoFTool template void make_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity_pattern, - const ConstraintMatrix &constraints = ConstraintMatrix(), - const bool keep_constrained_dofs = true, - const types::subdomain_id subdomain_id = types::invalid_subdomain_id); + SparsityPattern &sparsity_pattern, + const ConstraintMatrix &constraints = ConstraintMatrix(), + const bool keep_constrained_dofs = true, - const types::subdomain_id_t subdomain_id = types::invalid_subdomain_id); ++ const types::subdomain_id subdomain_id = types::invalid_subdomain_id); /** * Locate non-zero entries for @@@ -654,11 -654,11 +654,11 @@@ template void make_sparsity_pattern (const DH &dof, - const Table<2, Coupling> &coupling, - SparsityPattern &sparsity_pattern, - const ConstraintMatrix &constraints = ConstraintMatrix(), - const bool keep_constrained_dofs = true, - const types::subdomain_id subdomain_id = types::invalid_subdomain_id); + const Table<2, Coupling> &coupling, + SparsityPattern &sparsity_pattern, + const ConstraintMatrix &constraints = ConstraintMatrix(), + const bool keep_constrained_dofs = true, - const types::subdomain_id_t subdomain_id = types::invalid_subdomain_id); ++ const types::subdomain_id subdomain_id = types::invalid_subdomain_id); /** * @deprecated This is the old @@@ -812,10 -812,10 +812,10 @@@ template void make_flux_sparsity_pattern (const DH &dof_handler, - SparsityPattern &sparsity_pattern, - const ConstraintMatrix &constraints, - const bool keep_constrained_dofs = true, - const types::subdomain_id subdomain_id = numbers::invalid_unsigned_int); + SparsityPattern &sparsity_pattern, + const ConstraintMatrix &constraints, + const bool keep_constrained_dofs = true, - const types::subdomain_id_t subdomain_id = numbers::invalid_unsigned_int); ++ const types::subdomain_id subdomain_id = numbers::invalid_unsigned_int); /** * This function does the same as @@@ -1335,44 -1335,8 +1335,44 @@@ extract_boundary_dofs (const DH &dof_handler, const std::vector &component_mask, std::vector &selected_dofs, - const std::set &boundary_indicators = std::set()); + const std::set &boundary_indicators = std::set()); + /** + * This function does the same as the previous one but it + * returns its result as an IndexSet rather than a std::vector@. + * Thus, it can also be called for DoFHandler objects that are + * defined on parallel::distributed::Triangulation objects. + * + * @note If the DoFHandler object is indeed defined on a + * parallel::distributed::Triangulation, then the @p selected_dofs + * index set will contain only those degrees of freedom on the + * boundary that belong to the locally relevant set (see + * @ref GlossLocallyRelevantDof "locally relevant DoFs"). + * + * @param dof_handler The object that describes which degrees of freedom + * live on which cell + * @param component_mask A mask denoting the vector components of the + * finite element that should be considered (see also + * @ref GlossComponentMask). + * @param selected_dofs The IndexSet object that is returned and that + * will contain the indices of degrees of freedom that are + * located on the boundary (and correspond to the selected + * vector components and boundary indicators, depending on + * the values of the @p component_mask and @p boundary_indicators + * arguments). + * @param boundary_indicators If empty, this function extracts the + * indices of the degrees of freedom for all parts of the boundary. + * If it is a non-empty list, then the function only considers + * boundary faces with the boundary indicators listed in this + * argument. + */ + template + void + extract_boundary_dofs (const DH &dof_handler, + const std::vector &component_mask, + IndexSet &selected_dofs, + const std::set &boundary_indicators = std::set()); + /** * This function does the same as the previous one but it * returns its result as an IndexSet rather than a std::vector@. @@@ -1442,9 -1406,9 +1442,9 @@@ template void extract_dofs_with_support_on_boundary (const DH &dof_handler, - const std::vector &component_mask, - std::vector &selected_dofs, - const std::set &boundary_indicators = std::set()); + const std::vector &component_mask, + std::vector &selected_dofs, - const std::set &boundary_indicators = std::set()); ++ const std::set &boundary_indicators = std::set()); /** * @name Hanging Nodes @@@ -1494,8 -1458,8 +1494,8 @@@ template void extract_subdomain_dofs (const DH &dof_handler, - const types::subdomain_id_t subdomain_id, - const types::subdomain_id subdomain_id, - std::vector &selected_dofs); ++ const types::subdomain_id subdomain_id, + std::vector &selected_dofs); /** @@@ -1625,7 -1589,7 +1625,7 @@@ template void get_subdomain_association (const DH &dof_handler, - std::vector &subdomain); - std::vector &subdomain); ++ std::vector &subdomain); /** * Count how many degrees of freedom are @@@ -1669,7 -1633,7 +1669,7 @@@ template unsigned int count_dofs_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id_t subdomain); - const types::subdomain_id subdomain); ++ const types::subdomain_id subdomain); /** * Count how many degrees of freedom are @@@ -1706,8 -1670,8 +1706,8 @@@ template void count_dofs_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id_t subdomain, - const types::subdomain_id subdomain, - std::vector &n_dofs_on_subdomain); ++ const types::subdomain_id subdomain, + std::vector &n_dofs_on_subdomain); /** * Return a set of indices that denotes @@@ -1748,7 -1712,8 +1748,7 @@@ template IndexSet dof_indices_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id_t subdomain); - const types::subdomain_id subdomain); - ++ const types::subdomain_id subdomain); // @} /** * @name Dof indices for patches @@@ -2454,8 -2419,8 +2454,8 @@@ template void map_dof_to_boundary_indices (const DH &dof_handler, - const std::set &boundary_indicators, - const std::set &boundary_indicators, - std::vector &mapping); ++ const std::set &boundary_indicators, + std::vector &mapping); /** * Return a list of support diff --cc deal.II/include/deal.II/hp/dof_handler.h index 3876c81c8d,c15ffc2a42..81e9ee9864 --- a/deal.II/include/deal.II/hp/dof_handler.h +++ b/deal.II/include/deal.II/hp/dof_handler.h @@@ -863,8 -862,8 +863,8 @@@ namespace h * indicators under * consideration. */ - unsigned int + types::global_dof_index - n_boundary_dofs (const std::set &boundary_indicators) const; + n_boundary_dofs (const std::set &boundary_indicators) const; /** * Return the number of @@@ -1098,26 -1097,12 +1098,26 @@@ * copy such an object. */ DoFHandler & operator = (const DoFHandler &); - + + class MGVertexDoFs { + public: + MGVertexDoFs (); + ~MGVertexDoFs (); + types::global_dof_index get_index (const unsigned int level, const unsigned int dof_number) const; + void set_index (const unsigned int level, const unsigned int dof_number, const types::global_dof_index index); + }; + /** * Free all used memory. */ void clear_space (); - + + template + types::global_dof_index get_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index) const; - ++ + template + void set_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index) const; + /** * Create default tables for * the active_fe_indices in @@@ -1291,9 -1277,7 +1291,9 @@@ * the present class. */ std::vector vertex_dofs_offsets; - + + std::vector mg_vertex_dofs; + /** * Array to store the * information, if a cell on @@@ -1401,18 -1385,7 +1401,37 @@@ return *tria; } + template + inline - types::global_dof_index DoFHandler::MGVertexDoFs::get_index (const unsigned int level, const unsigned int dof_number) const { ++ DoFHandler::MGVertexDoFs::MGVertexDoFs() ++ { ++ Assert (false, ExcNotImplemented ()); ++ } ++ ++ template ++ inline ++ DoFHandler::MGVertexDoFs::~MGVertexDoFs() ++ { ++ Assert (false, ExcNotImplemented ()); ++ } ++ ++ template ++ inline ++ types::global_dof_index DoFHandler::MGVertexDoFs::get_index (const unsigned int, ++ const unsigned int) const ++ { + Assert (false, ExcNotImplemented ()); + return invalid_dof_index; + } - + + template + inline - void DoFHandler::MGVertexDoFs::set_index (const unsigned int level, const unsigned int dof_number, const types::global_dof_index index) { ++ void DoFHandler::MGVertexDoFs::set_index (const unsigned int, ++ const unsigned int, ++ types::global_dof_index) ++ { + Assert (false, ExcNotImplemented ()); + } #endif diff --cc deal.II/include/deal.II/multigrid/mg_transfer.templates.h index bcd44fdf6e,94b82212dc..7e51e549cd --- a/deal.II/include/deal.II/multigrid/mg_transfer.templates.h +++ b/deal.II/include/deal.II/multigrid/mg_transfer.templates.h @@@ -356,16 -215,11 +356,11 @@@ MGTransferPrebuilt::memory_cons { std::size_t result = sizeof(*this); result += sizeof(unsigned int) * sizes.size(); - #ifdef DEAL_PREFER_MATRIX_EZ - std::vector > >::const_iterator m; - const std::vector > >::const_iterator end = prolongation_matrices.end(); - for (m = prolongation_matrices.begin(); m != end ; ++m) - result += *m->memory_consumption(); - #else + for (unsigned int i=0;imemory_consumption() - + prolongation_sparsities[i]->memory_consumption(); + + prolongation_sparsities[i]->memory_consumption(); - #endif + return result; } diff --cc deal.II/include/deal.II/numerics/vector_tools.h index 0b962ce803,f6a1387560..0c85f89126 --- a/deal.II/include/deal.II/numerics/vector_tools.h +++ b/deal.II/include/deal.II/numerics/vector_tools.h @@@ -767,9 -767,9 +767,9 @@@ namespace VectorTool void interpolate_boundary_values (const Mapping &mapping, const DH &dof, - const types::boundary_id_t boundary_component, + const types::boundary_id boundary_component, const Function &boundary_function, - std::map &boundary_values, + std::map &boundary_values, const std::vector &component_mask = std::vector()); /** @@@ -781,9 -781,9 +781,9 @@@ template void interpolate_boundary_values (const DH &dof, - const types::boundary_id_t boundary_component, + const types::boundary_id boundary_component, const Function &boundary_function, - std::map &boundary_values, + std::map &boundary_values, const std::vector &component_mask = std::vector()); diff --cc deal.II/include/deal.II/numerics/vector_tools.templates.h index 09c182635e,2df4c30f05..61a33e0d60 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@@ -1941,9 -1942,9 +1942,9 @@@ namespace VectorTool void interpolate_boundary_values (const Mapping &mapping, const DH &dof, - const types::boundary_id_t boundary_component, + const types::boundary_id boundary_component, const Function &boundary_function, - std::map &boundary_values, + std::map &boundary_values, const std::vector &component_mask) { typename FunctionMap::type function_map; @@@ -1957,9 -1958,9 +1958,9 @@@ template void interpolate_boundary_values (const DH &dof, - const types::boundary_id_t boundary_component, + const types::boundary_id boundary_component, const Function &boundary_function, - std::map &boundary_values, + std::map &boundary_values, const std::vector &component_mask) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); @@@ -2140,8 -2141,8 +2141,8 @@@ else AssertDimension (dof.get_fe().n_components(), component_mapping.size()); - std::vector dof_to_boundary_mapping; + std::vector dof_to_boundary_mapping; - std::set selected_boundary_components; + std::set selected_boundary_components; for (typename FunctionMap::type::const_iterator i=boundary_functions.begin(); i!=boundary_functions.end(); ++i) selected_boundary_components.insert (i->first); diff --cc deal.II/source/dofs/dof_handler.cc index dc84a2fbd2,9e466cd2d0..4be24de818 --- a/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/source/dofs/dof_handler.cc @@@ -257,427 -257,99 +257,427 @@@ namespace interna * was selected when calling * @p distribute_dofs the last time. */ - template - static - void reserve_space (DoFHandler<1,spacedim> &dof_handler) - { - dof_handler.vertex_dofs - .resize(dof_handler.tria->n_vertices() * - dof_handler.selected_fe->dofs_per_vertex, - DoFHandler<1,spacedim>::invalid_dof_index); - - for (unsigned int i=0; in_levels(); ++i) - { - dof_handler.levels - .push_back (new internal::DoFHandler::DoFLevel<1>); - - dof_handler.levels.back()->lines.dofs - .resize (dof_handler.tria->n_raw_lines(i) * - dof_handler.selected_fe->dofs_per_line, - DoFHandler<1,spacedim>::invalid_dof_index); - - dof_handler.levels.back()->cell_dof_indices_cache - .resize (dof_handler.tria->n_raw_lines(i) * - dof_handler.selected_fe->dofs_per_cell, - DoFHandler<1,spacedim>::invalid_dof_index); - } - } - - - template - static - void reserve_space (DoFHandler<2,spacedim> &dof_handler) - { - dof_handler.vertex_dofs - .resize(dof_handler.tria->n_vertices() * - dof_handler.selected_fe->dofs_per_vertex, - DoFHandler<2,spacedim>::invalid_dof_index); - - for (unsigned int i=0; in_levels(); ++i) - { - dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<2>); - - dof_handler.levels.back()->quads.dofs - .resize (dof_handler.tria->n_raw_quads(i) * - dof_handler.selected_fe->dofs_per_quad, - DoFHandler<2,spacedim>::invalid_dof_index); - - dof_handler.levels.back()->cell_dof_indices_cache - .resize (dof_handler.tria->n_raw_quads(i) * - dof_handler.selected_fe->dofs_per_cell, - DoFHandler<2,spacedim>::invalid_dof_index); - } - - dof_handler.faces = new internal::DoFHandler::DoFFaces<2>; - dof_handler.faces->lines.dofs - .resize (dof_handler.tria->n_raw_lines() * - dof_handler.selected_fe->dofs_per_line, - DoFHandler<2,spacedim>::invalid_dof_index); - } - - - template - static - void reserve_space (DoFHandler<3,spacedim> &dof_handler) - { - dof_handler.vertex_dofs - .resize(dof_handler.tria->n_vertices() * - dof_handler.selected_fe->dofs_per_vertex, - DoFHandler<3,spacedim>::invalid_dof_index); - - for (unsigned int i=0; in_levels(); ++i) - { - dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<3>); - - dof_handler.levels.back()->hexes.dofs - .resize (dof_handler.tria->n_raw_hexs(i) * - dof_handler.selected_fe->dofs_per_hex, - DoFHandler<3,spacedim>::invalid_dof_index); - - dof_handler.levels.back()->cell_dof_indices_cache - .resize (dof_handler.tria->n_raw_hexs(i) * - dof_handler.selected_fe->dofs_per_cell, - DoFHandler<3,spacedim>::invalid_dof_index); - } - dof_handler.faces = new internal::DoFHandler::DoFFaces<3>; - - dof_handler.faces->lines.dofs - .resize (dof_handler.tria->n_raw_lines() * - dof_handler.selected_fe->dofs_per_line, - DoFHandler<3,spacedim>::invalid_dof_index); - dof_handler.faces->quads.dofs - .resize (dof_handler.tria->n_raw_quads() * - dof_handler.selected_fe->dofs_per_quad, - DoFHandler<3,spacedim>::invalid_dof_index); - } + template + static + void reserve_space (DoFHandler<1,spacedim> &dof_handler) + { + dof_handler.vertex_dofs + .resize(dof_handler.tria->n_vertices() * + dof_handler.selected_fe->dofs_per_vertex, + DoFHandler<1,spacedim>::invalid_dof_index); + + for (unsigned int i=0; in_levels(); ++i) + { + dof_handler.levels + .push_back (new internal::DoFHandler::DoFLevel<1>); + + dof_handler.levels.back()->lines.dofs + .resize (dof_handler.tria->n_raw_lines(i) * + dof_handler.selected_fe->dofs_per_line, + DoFHandler<1,spacedim>::invalid_dof_index); + + dof_handler.levels.back()->cell_dof_indices_cache + .resize (dof_handler.tria->n_raw_lines(i) * + dof_handler.selected_fe->dofs_per_cell, + DoFHandler<1,spacedim>::invalid_dof_index); + } + } + + + template + static + void reserve_space (DoFHandler<2,spacedim> &dof_handler) + { + dof_handler.vertex_dofs + .resize(dof_handler.tria->n_vertices() * + dof_handler.selected_fe->dofs_per_vertex, + DoFHandler<2,spacedim>::invalid_dof_index); + + for (unsigned int i=0; in_levels(); ++i) + { + dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<2>); + + dof_handler.levels.back()->quads.dofs + .resize (dof_handler.tria->n_raw_quads(i) * + dof_handler.selected_fe->dofs_per_quad, + DoFHandler<2,spacedim>::invalid_dof_index); + + dof_handler.levels.back()->cell_dof_indices_cache + .resize (dof_handler.tria->n_raw_quads(i) * + dof_handler.selected_fe->dofs_per_cell, + DoFHandler<2,spacedim>::invalid_dof_index); + } + + dof_handler.faces = new internal::DoFHandler::DoFFaces<2>; + dof_handler.faces->lines.dofs + .resize (dof_handler.tria->n_raw_lines() * + dof_handler.selected_fe->dofs_per_line, + DoFHandler<2,spacedim>::invalid_dof_index); + } + + + template + static + void reserve_space (DoFHandler<3,spacedim> &dof_handler) + { + dof_handler.vertex_dofs + .resize(dof_handler.tria->n_vertices() * + dof_handler.selected_fe->dofs_per_vertex, + DoFHandler<3,spacedim>::invalid_dof_index); + + for (unsigned int i=0; in_levels(); ++i) + { + dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<3>); + + dof_handler.levels.back()->hexes.dofs + .resize (dof_handler.tria->n_raw_hexs(i) * + dof_handler.selected_fe->dofs_per_hex, + DoFHandler<3,spacedim>::invalid_dof_index); + + dof_handler.levels.back()->cell_dof_indices_cache + .resize (dof_handler.tria->n_raw_hexs(i) * + dof_handler.selected_fe->dofs_per_cell, + DoFHandler<3,spacedim>::invalid_dof_index); + } + dof_handler.faces = new internal::DoFHandler::DoFFaces<3>; + + dof_handler.faces->lines.dofs + .resize (dof_handler.tria->n_raw_lines() * + dof_handler.selected_fe->dofs_per_line, + DoFHandler<3,spacedim>::invalid_dof_index); + dof_handler.faces->quads.dofs + .resize (dof_handler.tria->n_raw_quads() * + dof_handler.selected_fe->dofs_per_quad, + DoFHandler<3,spacedim>::invalid_dof_index); + } - ++ + template + static + void reserve_space_mg (DoFHandler<1, spacedim>& dof_handler) { + Assert (dof_handler.get_tria ().n_levels () > 0, ExcMessage ("Invalid triangulation")); + dof_handler.clear_mg_space (); - ++ + const dealii::Triangulation<1, spacedim>& tria = dof_handler.get_tria (); + const unsigned int& dofs_per_line = dof_handler.get_fe ().dofs_per_line; + const unsigned int& n_levels = tria.n_levels (); - ++ + for (unsigned int i = 0; i < n_levels; ++i) { + dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<1>); + dof_handler.mg_levels.back ()->lines.dofs = std::vector (tria.n_raw_lines (i) * dofs_per_line, DoFHandler<1>::invalid_dof_index); + } - ++ + const unsigned int& n_vertices = tria.n_vertices (); - ++ + dof_handler.mg_vertex_dofs.resize (n_vertices); - ++ + std::vector max_level (n_vertices, 0); + std::vector min_level (n_vertices, n_levels); - ++ + for (typename dealii::Triangulation<1, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell) { + const unsigned int level = cell->level (); - ++ + for (unsigned int vertex = 0; vertex < GeometryInfo<1>::vertices_per_cell; ++vertex) { + const unsigned int vertex_index = cell->vertex_index (vertex); - ++ + if (min_level[vertex_index] > level) + min_level[vertex_index] = level; - ++ + if (max_level[vertex_index] < level) + max_level[vertex_index] = level; + } + } - ++ + for (unsigned int vertex = 0; vertex < n_vertices; ++vertex) + if (tria.vertex_used (vertex)) { + Assert (min_level[vertex] < n_levels, ExcInternalError ()); + Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ()); + dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], dof_handler.get_fe ().dofs_per_vertex); + } - ++ + else { + Assert (min_level[vertex] == n_levels, ExcInternalError ()); + Assert (max_level[vertex] == 0, ExcInternalError ()); + dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0); + } + } - ++ + template + static + void reserve_space_mg (DoFHandler<2, spacedim>& dof_handler) { + Assert (dof_handler.get_tria ().n_levels () > 0, ExcMessage ("Invalid triangulation")); + dof_handler.clear_mg_space (); - ++ + const dealii::FiniteElement<2, spacedim>& fe = dof_handler.get_fe (); + const dealii::Triangulation<2, spacedim>& tria = dof_handler.get_tria (); + const unsigned int& n_levels = tria.n_levels (); - ++ + for (unsigned int i = 0; i < n_levels; ++i) { + dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<2>); + dof_handler.mg_levels.back ()->quads.dofs = std::vector (tria.n_raw_quads (i) * fe.dofs_per_quad, DoFHandler<2>::invalid_dof_index); + } - ++ + dof_handler.mg_faces = new internal::DoFHandler::DoFFaces<2>; + dof_handler.mg_faces->lines.dofs = std::vector (tria.n_raw_lines () * fe.dofs_per_line, DoFHandler<2>::invalid_dof_index); - ++ + const unsigned int& n_vertices = tria.n_vertices (); - ++ + dof_handler.mg_vertex_dofs.resize (n_vertices); - ++ + std::vector max_level (n_vertices, 0); + std::vector min_level (n_vertices, n_levels); - ++ + for (typename dealii::Triangulation<2, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell) { + const unsigned int level = cell->level (); - ++ + for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_cell; ++vertex) { + const unsigned int vertex_index = cell->vertex_index (vertex); - ++ + if (min_level[vertex_index] > level) + min_level[vertex_index] = level; - ++ + if (max_level[vertex_index] < level) + max_level[vertex_index] = level; + } + } - ++ + for (unsigned int vertex = 0; vertex < n_vertices; ++vertex) + if (tria.vertex_used (vertex)) { + Assert (min_level[vertex] < n_levels, ExcInternalError ()); + Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ()); + dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], fe.dofs_per_vertex); + } - ++ + else { + Assert (min_level[vertex] == n_levels, ExcInternalError ()); + Assert (max_level[vertex] == 0, ExcInternalError ()); + dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0); + } + } - ++ + template + static + void reserve_space_mg (DoFHandler<3, spacedim>& dof_handler) { + Assert (dof_handler.get_tria ().n_levels () > 0, ExcMessage ("Invalid triangulation")); + dof_handler.clear_mg_space (); - ++ + const dealii::FiniteElement<3, spacedim>& fe = dof_handler.get_fe (); + const dealii::Triangulation<3, spacedim>& tria = dof_handler.get_tria (); + const unsigned int& n_levels = tria.n_levels (); - ++ + for (unsigned int i = 0; i < n_levels; ++i) { + dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<3>); + dof_handler.mg_levels.back ()->hexes.dofs = std::vector (tria.n_raw_hexs (i) * fe.dofs_per_hex, DoFHandler<3>::invalid_dof_index); + } - ++ + dof_handler.mg_faces = new internal::DoFHandler::DoFFaces<3>; + dof_handler.mg_faces->lines.dofs = std::vector (tria.n_raw_lines () * fe.dofs_per_line, DoFHandler<3>::invalid_dof_index); + dof_handler.mg_faces->quads.dofs = std::vector (tria.n_raw_quads () * fe.dofs_per_quad, DoFHandler<3>::invalid_dof_index); - ++ + const unsigned int& n_vertices = tria.n_vertices (); - ++ + dof_handler.mg_vertex_dofs.resize (n_vertices); - ++ + std::vector max_level (n_vertices, 0); + std::vector min_level (n_vertices, n_levels); - ++ + for (typename dealii::Triangulation<3, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell) { + const unsigned int level = cell->level (); - ++ + for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_cell; ++vertex) { + const unsigned int vertex_index = cell->vertex_index (vertex); - ++ + if (min_level[vertex_index] > level) + min_level[vertex_index] = level; - ++ + if (max_level[vertex_index] < level) + max_level[vertex_index] = level; + } + } - ++ + for (unsigned int vertex = 0; vertex < n_vertices; ++vertex) + if (tria.vertex_used (vertex)) { + Assert (min_level[vertex] < n_levels, ExcInternalError ()); + Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ()); + dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], fe.dofs_per_vertex); + } - ++ + else { + Assert (min_level[vertex] == n_levels, ExcInternalError ()); + Assert (max_level[vertex] == 0, ExcInternalError ()); + dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0); + } + } - ++ + template + static + types::global_dof_index distribute_dofs_on_cell (typename DoFHandler<1, spacedim>::cell_iterator& cell, types::global_dof_index next_free_dof) { + const FiniteElement<1, spacedim>& fe = cell->get_fe (); - ++ + if (fe.dofs_per_vertex > 0) + for (unsigned int vertex = 0; vertex < GeometryInfo<1>::vertices_per_cell; ++vertex) { + typename DoFHandler<1, spacedim>::cell_iterator neighbor = cell->neighbor (vertex); - ++ + if (neighbor.state () == IteratorState::valid) + if (neighbor->user_flag_set () && (neighbor->level () == cell->level ())) { + if (vertex == 0) + for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof) + cell->set_mg_vertex_dof_index (cell->level (), 0, dof, neighbor->mg_vertex_dof_index (cell->level (), 1, dof)); - ++ + else + for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof) + cell->set_mg_vertex_dof_index (cell->level (), 1, dof, neighbor->mg_vertex_dof_index (cell->level (), 0, dof)); - ++ + continue; + } - ++ + for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof) + cell->set_mg_vertex_dof_index (cell->level (), vertex, dof, next_free_dof++); + } - ++ + if (fe.dofs_per_line > 0) + for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof) + cell->set_mg_dof_index (cell->level (), dof, next_free_dof++); - ++ + cell->set_user_flag (); + return next_free_dof; + } - ++ + template + static + types::global_dof_index distribute_dofs_on_cell (typename DoFHandler<2, spacedim>::cell_iterator& cell, types::global_dof_index next_free_dof) { + const FiniteElement<2, spacedim>& fe = cell->get_fe (); - ++ + if (fe.dofs_per_vertex > 0) + for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_cell; ++vertex) + if (cell->mg_vertex_dof_index (cell->level (), vertex, 0) == DoFHandler<2>::invalid_dof_index) + for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof) + cell->set_mg_vertex_dof_index (cell->level (), vertex, dof, next_free_dof++); - ++ + if (fe.dofs_per_line > 0) + for (unsigned int face = 0; face < GeometryInfo<2>::faces_per_cell; ++face) { + typename DoFHandler<2, spacedim>::line_iterator line = cell->line (face); - ++ + if (line->mg_dof_index (cell->level (), 0) == DoFHandler<2>::invalid_dof_index) + for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof) + line->set_mg_dof_index (cell->level (), dof, next_free_dof++); + } - ++ + if (fe.dofs_per_quad > 0) + for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof) + cell->set_mg_dof_index (cell->level (), dof, next_free_dof++); - ++ + cell->set_user_flag (); + return next_free_dof; + } - ++ + template + static + types::global_dof_index distribute_dofs_on_cell (typename DoFHandler<3, spacedim>::cell_iterator& cell, types::global_dof_index next_free_dof) { + const FiniteElement<3, spacedim>& fe = cell->get_fe (); - ++ + if (fe.dofs_per_vertex > 0) + for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_cell; ++vertex) + if (cell->mg_vertex_dof_index (cell->level (), vertex, 0) == DoFHandler<3>::invalid_dof_index) + for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof) + cell->set_mg_vertex_dof_index (cell->level (), vertex, dof, next_free_dof++); - ++ + if (fe.dofs_per_line > 0) + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell; ++line) { + typename DoFHandler<3, spacedim>::line_iterator line_it = cell->line (line); - ++ + if (line_it->mg_dof_index (cell->level (), 0) == DoFHandler<3>::invalid_dof_index) + for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof) + line_it->set_mg_dof_index (cell->level (), dof, next_free_dof++); + } - ++ + if (fe.dofs_per_quad > 0) + for (unsigned int face = 0; face < GeometryInfo<3>::quads_per_cell; ++face) { + typename DoFHandler<3, spacedim>::quad_iterator quad = cell->quad (face); - ++ + if (quad->mg_dof_index (cell->level (), 0) == DoFHandler<3>::invalid_dof_index) + for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof) + quad->set_mg_dof_index (cell->level (), dof, next_free_dof++); + } - ++ + if (fe.dofs_per_hex > 0) + for (unsigned int dof = 0; dof < fe.dofs_per_hex; ++dof) + cell->set_mg_dof_index (cell->level (), dof, next_free_dof++); - ++ + cell->set_user_flag (); + return next_free_dof; + } - ++ + template + static + types::global_dof_index get_dof_index (const DoFHandler<1, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<1>& mg_level, internal::DoFHandler::DoFFaces<1>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>) { + return mg_level.lines.get_dof_index (dof_handler, obj_index, fe_index, local_index); + } - ++ + template + static + types::global_dof_index get_dof_index (const DoFHandler<2, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<2>&, internal::DoFHandler::DoFFaces<2>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>) { + return mg_faces.lines.get_dof_index (dof_handler, obj_index, fe_index, local_index); + } - ++ + template + static + types::global_dof_index get_dof_index (const DoFHandler<2, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<2>& mg_level, internal::DoFHandler::DoFFaces<2>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<2>) { + return mg_level.quads.get_dof_index (dof_handler, obj_index, fe_index, local_index); + } - ++ + template + static + types::global_dof_index get_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>&, internal::DoFHandler::DoFFaces<3>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>) { + return mg_faces.lines.get_dof_index (dof_handler, obj_index, fe_index, local_index); + } - ++ + template + static + types::global_dof_index get_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>&, internal::DoFHandler::DoFFaces<3>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<2>) { + return mg_faces.quads.get_dof_index (dof_handler, obj_index, fe_index, local_index); + } - ++ + template + static + types::global_dof_index get_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>& mg_level, internal::DoFHandler::DoFFaces<3>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<3>) { + return mg_level.hexes.get_dof_index (dof_handler, obj_index, fe_index, local_index); + } - ++ + template + static + void set_dof_index (const DoFHandler<1, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<1>& mg_level, internal::DoFHandler::DoFFaces<1>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>) { + mg_level.lines.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index); + } - ++ + template + static + void set_dof_index (const DoFHandler<2, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<2>&, internal::DoFHandler::DoFFaces<2>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>) { + mg_faces.lines.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index); + } - ++ + template + static + void set_dof_index (const DoFHandler<2, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<2>& mg_level, internal::DoFHandler::DoFFaces<2>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<2>) { + mg_level.quads.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index); + } - ++ + template + static + void set_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>&, internal::DoFHandler::DoFFaces<3>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>) { + mg_faces.lines.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index); + } - ++ + template + static + void set_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>&, internal::DoFHandler::DoFFaces<3>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<2>) { + mg_faces.quads.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index); + } - ++ + template + static + void set_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>& mg_level, internal::DoFHandler::DoFFaces<3>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<3>) { + mg_level.hexes.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index); + } }; } } @@@ -1819,15 -1489,15 +1819,15 @@@ types::global_dof_index DoFHandler<1>:: template <> - types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set &boundary_indicators) const -unsigned int DoFHandler<1>::n_boundary_dofs (const std::set &boundary_indicators) const ++types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set &boundary_indicators) const { // check that only boundary // indicators 0 and 1 are allowed // in 1d - for (std::set::const_iterator i=boundary_indicators.begin(); + for (std::set::const_iterator i=boundary_indicators.begin(); i!=boundary_indicators.end(); ++i) Assert ((*i == 0) || (*i == 1), - ExcInvalidBoundaryIndicator()); + ExcInvalidBoundaryIndicator()); return boundary_indicators.size()*get_fe().dofs_per_vertex; } @@@ -1858,15 -1528,15 +1858,15 @@@ types::global_dof_index DoFHandler<1,2> template <> - types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set &boundary_indicators) const -unsigned int DoFHandler<1,2>::n_boundary_dofs (const std::set &boundary_indicators) const ++types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set &boundary_indicators) const { // check that only boundary // indicators 0 and 1 are allowed // in 1d - for (std::set::const_iterator i=boundary_indicators.begin(); + for (std::set::const_iterator i=boundary_indicators.begin(); i!=boundary_indicators.end(); ++i) Assert ((*i == 0) || (*i == 1), - ExcInvalidBoundaryIndicator()); + ExcInvalidBoundaryIndicator()); return boundary_indicators.size()*get_fe().dofs_per_vertex; } @@@ -1932,11 -1602,11 +1932,11 @@@ DoFHandler::n_boundary_do template -unsigned int +types::global_dof_index - DoFHandler::n_boundary_dofs (const std::set &boundary_indicators) const + DoFHandler::n_boundary_dofs (const std::set &boundary_indicators) const { Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(), - ExcInvalidBoundaryIndicator()); + ExcInvalidBoundaryIndicator()); std::set boundary_dofs; @@@ -1964,25 -1634,16 +1964,25 @@@ std::size_ DoFHandler::memory_consumption () const { std::size_t mem = (MemoryConsumption::memory_consumption (tria) + - MemoryConsumption::memory_consumption (selected_fe) + - MemoryConsumption::memory_consumption (block_info_object) + - MemoryConsumption::memory_consumption (levels) + - MemoryConsumption::memory_consumption (*faces) + - MemoryConsumption::memory_consumption (faces) + - sizeof (number_cache) + - MemoryConsumption::memory_consumption (vertex_dofs)); + MemoryConsumption::memory_consumption (selected_fe) + + MemoryConsumption::memory_consumption (block_info_object) + + MemoryConsumption::memory_consumption (levels) + + MemoryConsumption::memory_consumption (*faces) + + MemoryConsumption::memory_consumption (faces) + + sizeof (number_cache) + + MemoryConsumption::memory_consumption (vertex_dofs)); for (unsigned int i=0; imemory_consumption (); - ++ + mem += MemoryConsumption::memory_consumption (*mg_faces); - ++ + for (unsigned int i = 0; i < mg_vertex_dofs.size (); ++i) + mem += sizeof (MGVertexDoFs) + (1 + mg_vertex_dofs[i].get_finest_level () - mg_vertex_dofs[i].get_coarsest_level ()) * sizeof (types::global_dof_index); - ++ + mem += MemoryConsumption::memory_consumption (mg_used_dofs); return mem; } @@@ -2027,53 -1688,6 +2027,53 @@@ distribute_dofs (const FiniteElement +void DoFHandler::distribute_mg_dofs (const FiniteElement& fe, const types::global_dof_index offset) { + Assert ((dynamic_cast*> (&*tria) == 0), ExcMessage ("Invalid triangulation")); + distribute_dofs (fe); + reserve_space (); - ++ + const unsigned int& n_levels = (*tria).n_levels (); - ++ + mg_used_dofs.resize (n_levels, 0); - ++ + std::vector user_flags; - ++ + (*tria).save_user_flags (user_flags); + const_cast&>(*tria).clear_user_flags (); - ++ + for (unsigned int level = 0; level < n_levels; ++level) { + types::global_dof_index next_free_dof = offset; - ++ + for (cell_iterator cell = begin (level); cell != end (level); ++cell) + next_free_dof = internal::DoFHandler::Implementation::distribute_dofs_on_cell (cell, next_free_dof); - ++ + mg_used_dofs[level] = next_free_dof; + } - ++ + const_cast&>(*tria).load_user_flags (user_flags); + block_info_object.initialize (*this, true, true); +} + +template +void DoFHandler::reserve_space () { + internal::DoFHandler::Implementation::reserve_space_mg (*this); +} + +template +void DoFHandler::clear_mg_space () { + for (unsigned int i = 0; i < mg_levels.size (); ++i) + delete mg_levels[i]; - ++ + mg_levels.clear (); + delete mg_faces; + mg_faces = NULL; - ++ + std::vector tmp; - ++ + std::swap (mg_vertex_dofs, tmp); +} + template void DoFHandler::initialize_local_block_info () @@@ -2196,72 -1809,6 +2196,72 @@@ void DoFHandler::clear_sp std::swap (vertex_dofs, tmp); } +template +template +types::global_dof_index DoFHandler::get_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index) const { + return internal::DoFHandler::Implementation::get_dof_index (*this, *this->mg_levels[obj_level], *this->mg_faces, obj_index, fe_index, local_index, internal::int2type ()); +} + +template +template +void DoFHandler::set_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index) const { + internal::DoFHandler::Implementation::set_dof_index (*this, *this->mg_levels[obj_level], *this->mg_faces, obj_index, fe_index, local_index, global_index, internal::int2type ()); +} + +template +DoFHandler::MGVertexDoFs::MGVertexDoFs (): coarsest_level (numbers::invalid_unsigned_int), finest_level (0), indices (0), indices_offset (0) { +} + +template +DoFHandler::MGVertexDoFs::~MGVertexDoFs () { + delete[] indices; + delete[] indices_offset; +} + +template +void DoFHandler::MGVertexDoFs::init (const unsigned int cl, const unsigned int fl, const unsigned int dofs_per_vertex) { + if (indices != 0) { + delete[] indices; + indices = 0; + } - ++ + if (indices_offset != 0) { + delete[] indices_offset; + indices_offset = 0; + } - ++ + coarsest_level = cl; + finest_level = fl; - ++ + if (cl > fl) + return; - ++ + const unsigned int n_levels = finest_level - coarsest_level + 1; + const unsigned int n_indices = n_levels * dofs_per_vertex; - ++ + indices = new unsigned int[n_indices]; + Assert (indices != 0, ExcNoMemory ()); - ++ + for (unsigned int i = 0; i < n_indices; ++i) + indices[i] = DoFHandler::invalid_dof_index; - ++ + indices_offset = new unsigned int[n_levels]; + Assert (indices != 0, ExcNoMemory ()); - ++ + for (unsigned int i = 0; i < n_levels; ++i) + indices_offset[i] = i * dofs_per_vertex; +} + +template +unsigned int DoFHandler::MGVertexDoFs::get_coarsest_level () const { + return coarsest_level; +} + +template +unsigned int DoFHandler::MGVertexDoFs::get_finest_level () const { + return finest_level; +} + /*-------------- Explicit Instantiations -------------------------------*/ #include "dof_handler.inst" diff --cc deal.II/source/dofs/dof_handler_policy.cc index 4e12ba94c2,5c8185b46f..ea53851a2d --- a/deal.II/source/dofs/dof_handler_policy.cc +++ b/deal.II/source/dofs/dof_handler_policy.cc @@@ -223,9 -224,9 +223,9 @@@ namespace interna */ template static - unsigned int - distribute_dofs (const unsigned int offset, + types::global_dof_index + distribute_dofs (const types::global_dof_index offset, - const types::subdomain_id_t subdomain_id, + const types::subdomain_id subdomain_id, DoFHandler &dof_handler) { const dealii::Triangulation & tria @@@ -911,9 -912,9 +911,9 @@@ // mark all own cells, that miss some // dof_data and collect the neighbors // that are going to send stuff to us - std::set senders; + std::set senders; { - std::vector local_dof_indices; + std::vector local_dof_indices; typename DoFHandler::active_cell_iterator cell, endc = dof_handler.end(); diff --cc deal.II/source/dofs/dof_renumbering.cc index fad2b36a3a,99bec3b823..f74aba9bd0 --- a/deal.II/source/dofs/dof_renumbering.cc +++ b/deal.II/source/dofs/dof_renumbering.cc @@@ -1038,9 -1037,24 +1038,24 @@@ namespace DoFRenumberin + template + void + sort_selected_dofs_back (DH& dof_handler, + const std::vector& selected_dofs, + const unsigned int level) + { + std::vector renumbering(dof_handler.n_dofs(level), + DH::invalid_dof_index); + compute_sort_selected_dofs_back(renumbering, dof_handler, selected_dofs, level); + + dof_handler.renumber_dofs(level, renumbering); + } + + + template void - compute_sort_selected_dofs_back (std::vector& new_indices, + compute_sort_selected_dofs_back (std::vector& new_indices, const DH& dof_handler, const std::vector& selected_dofs) { diff --cc deal.II/source/dofs/dof_tools.cc index 85c08f4f0e,e23dd691a8..d42c9c7f96 --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@@ -57,10 -57,10 +57,10 @@@ namespace DoFTool template void make_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity, - const ConstraintMatrix &constraints, - const bool keep_constrained_dofs, - const types::subdomain_id subdomain_id) + SparsityPattern &sparsity, + const ConstraintMatrix &constraints, + const bool keep_constrained_dofs, - const types::subdomain_id_t subdomain_id) ++ const types::subdomain_id subdomain_id) { const unsigned int n_dofs = dof.n_dofs(); @@@ -119,11 -119,11 +119,11 @@@ template void make_sparsity_pattern (const DH &dof, - const Table<2,Coupling> &couplings, - SparsityPattern &sparsity, - const ConstraintMatrix &constraints, - const bool keep_constrained_dofs, - const types::subdomain_id subdomain_id) + const Table<2,Coupling> &couplings, + SparsityPattern &sparsity, + const ConstraintMatrix &constraints, + const bool keep_constrained_dofs, - const types::subdomain_id_t subdomain_id) ++ const types::subdomain_id subdomain_id) { const unsigned int n_dofs = dof.n_dofs(); @@@ -518,10 -518,10 +518,10 @@@ template void make_flux_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity, - const ConstraintMatrix &constraints, - const bool keep_constrained_dofs, + SparsityPattern &sparsity, + const ConstraintMatrix &constraints, + const bool keep_constrained_dofs, - const types::subdomain_id_t subdomain_id) + const types::subdomain_id subdomain_id) { const unsigned int n_dofs = dof.n_dofs(); @@@ -4071,9 -4071,9 +4071,9 @@@ template void extract_dofs_with_support_on_boundary (const DH &dof_handler, - const std::vector &component_select, - std::vector &selected_dofs, - const std::set &boundary_indicators) + const std::vector &component_select, + std::vector &selected_dofs, - const std::set &boundary_indicators) ++ const std::set &boundary_indicators) { AssertDimension (component_select.size(), n_components(dof_handler)); Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(), @@@ -4306,11 -4306,11 +4306,11 @@@ template void extract_subdomain_dofs (const DH &dof_handler, - const types::subdomain_id_t subdomain_id, - const types::subdomain_id subdomain_id, - std::vector &selected_dofs) ++ const types::subdomain_id subdomain_id, + std::vector &selected_dofs) { Assert(selected_dofs.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); + ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); // preset all values by false std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); @@@ -4491,7 -4491,7 +4491,7 @@@ template void get_subdomain_association (const DH &dof_handler, - std::vector &subdomain_association) - std::vector &subdomain_association) ++ std::vector &subdomain_association) { // if the Triangulation is distributed, the // only thing we can usefully ask is for @@@ -4529,14 -4529,14 +4529,14 @@@ endc = dof_handler.end(); for (; cell!=endc; ++cell) { - Assert (cell->is_artificial() == false, - ExcMessage ("You can't call this function for meshes that " - "have artificial cells.")); + Assert (cell->is_artificial() == false, + ExcMessage ("You can't call this function for meshes that " + "have artificial cells.")); - const types::subdomain_id_t subdomain_id = cell->subdomain_id(); - const types::subdomain_id subdomain_id = cell->subdomain_id(); - const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; - local_dof_indices.resize (dofs_per_cell); - cell->get_dof_indices (local_dof_indices); ++ const types::subdomain_id subdomain_id = cell->subdomain_id(); + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + local_dof_indices.resize (dofs_per_cell); + cell->get_dof_indices (local_dof_indices); // set subdomain ids. if dofs // already have their values @@@ -4572,9 -4572,9 +4572,9 @@@ template unsigned int count_dofs_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id_t subdomain) - const types::subdomain_id subdomain) ++ const types::subdomain_id subdomain) { - std::vector subdomain_association (dof_handler.n_dofs()); + std::vector subdomain_association (dof_handler.n_dofs()); get_subdomain_association (dof_handler, subdomain_association); return std::count (subdomain_association.begin(), @@@ -4587,7 -4587,7 +4587,7 @@@ template IndexSet dof_indices_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id_t subdomain) - const types::subdomain_id subdomain) ++ const types::subdomain_id subdomain) { // If we have a distributed::Triangulation only @@@ -4648,12 -4648,12 +4648,12 @@@ template void count_dofs_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id_t subdomain, - const types::subdomain_id subdomain, - std::vector &n_dofs_on_subdomain) ++ const types::subdomain_id subdomain, + std::vector &n_dofs_on_subdomain) { Assert (n_dofs_on_subdomain.size() == dof_handler.get_fe().n_components(), - ExcDimensionMismatch (n_dofs_on_subdomain.size(), - dof_handler.get_fe().n_components())); + ExcDimensionMismatch (n_dofs_on_subdomain.size(), + dof_handler.get_fe().n_components())); std::fill (n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0); // in debug mode, make sure that there are @@@ -4875,15 -4875,15 +4875,15 @@@ #if defined(DEAL_II_USE_P4EST) && defined(DEAL_II_COMPILER_SUPPORTS_MPI) const unsigned int dim = DH::dimension; const unsigned int spacedim = DH::space_dimension; - + if (const parallel::distributed::Triangulation * tria - = (dynamic_cast*> - (&dof_handler.get_tria()))) + = (dynamic_cast*> + (&dof_handler.get_tria()))) { - std::vector local_dof_count = dofs_per_component; + std::vector local_dof_count = dofs_per_component; - MPI_Allreduce ( &local_dof_count[0], &dofs_per_component[0], n_target_components, - MPI_UNSIGNED, MPI_SUM, tria->get_communicator()); + MPI_Allreduce ( &local_dof_count[0], &dofs_per_component[0], n_target_components, + MPI_UNSIGNED, MPI_SUM, tria->get_communicator()); } #endif } @@@ -5474,10 -5474,10 +5474,10 @@@ #endif - return n_parameters_on_fine_grid; + return n_parameters_on_fine_grid; } - - + + } } @@@ -5871,8 -5871,8 +5871,8 @@@ template void map_dof_to_boundary_indices ( const DH &dof_handler, - const std::set &boundary_indicators, + const std::set &boundary_indicators, - std::vector &mapping) + std::vector &mapping) { Assert (&dof_handler.get_fe() != 0, ExcNoFESelected()); Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(), @@@ -5963,9 -5963,9 +5963,9 @@@ local_dof_indices.resize(cell->get_fe().dofs_per_cell); cell->get_dof_indices(local_dof_indices); - + const std::vector > & points = - fe_values.get_quadrature_points(); + fe_values.get_quadrature_points(); for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i) // insert the values into the map support_points[local_dof_indices[i]] = points[i]; diff --cc deal.II/source/dofs/dof_tools.inst.in index a629dcc61e,f1ecde11e7..4f3d9e3eea --- a/deal.II/source/dofs/dof_tools.inst.in +++ b/deal.II/source/dofs/dof_tools.inst.in @@@ -604,14 -600,14 +604,14 @@@ templat unsigned int DoFTools::count_dofs_with_subdomain_association > (const hp::DoFHandler &, - const types::subdomain_id_t); + const types::subdomain_id); - + template IndexSet DoFTools::dof_indices_with_subdomain_association > (const hp::DoFHandler &, - const types::subdomain_id_t); + const types::subdomain_id); - + template void DoFTools::get_subdomain_association > diff --cc deal.II/source/hp/dof_handler.cc index 6362c78416,75f2668392..6db6f20fa7 --- a/deal.II/source/hp/dof_handler.cc +++ b/deal.II/source/hp/dof_handler.cc @@@ -2745,7 -2745,7 +2745,7 @@@ namespace h template <> - types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set &boundary_indicators) const - unsigned int DoFHandler<1>::n_boundary_dofs (const std::set &boundary_indicators) const ++ types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set &boundary_indicators) const { Assert (finite_elements != 0, ExcNoFESelected()); @@@ -2797,7 -2797,7 +2797,7 @@@ } template <> - types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set &) const - unsigned int DoFHandler<1,2>::n_boundary_dofs (const std::set &) const ++ types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set &) const { Assert(false,ExcNotImplemented()); return 0; @@@ -2820,7 -2819,7 +2820,7 @@@ types::global_dof_index DoFHandler<1,3> } template <> - types::global_dof_index DoFHandler<1,3>::n_boundary_dofs (const std::set &) const - unsigned int DoFHandler<1,3>::n_boundary_dofs (const std::set &) const ++ types::global_dof_index DoFHandler<1,3>::n_boundary_dofs (const std::set &) const { Assert(false,ExcNotImplemented()); return 0; @@@ -2902,8 -2901,8 +2902,8 @@@ template - unsigned int + types::global_dof_index - DoFHandler::n_boundary_dofs (const std::set &boundary_indicators) const + DoFHandler::n_boundary_dofs (const std::set &boundary_indicators) const { Assert (finite_elements != 0, ExcNoFESelected()); Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(), @@@ -2955,7 -2954,7 +2955,7 @@@ template <> - types::global_dof_index DoFHandler<2,3>::n_boundary_dofs (const std::set &) const - unsigned int DoFHandler<2,3>::n_boundary_dofs (const std::set &) const ++ types::global_dof_index DoFHandler<2,3>::n_boundary_dofs (const std::set &) const { Assert(false,ExcNotImplemented()); return 0; @@@ -4143,31 -4140,6 +4143,32 @@@ template < } + template + template + types::global_dof_index - DoFHandler::get_dof_index (const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index) const ++ DoFHandler::get_dof_index (const unsigned int, ++ const unsigned int, ++ const unsigned int, ++ const unsigned int) const + { - ++ Assert (false, ExcNotImplemented()); ++ return numbers::invalid_unsigned_int; + } + + + template + template + void - DoFHandler::set_dof_index (const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - const types::global_dof_index global_index) const ++ DoFHandler::set_dof_index (const unsigned int, ++ const unsigned int, ++ const unsigned int, ++ const unsigned int, ++ const types::global_dof_index) const + { - ++ Assert (false, ExcNotImplemented()); + } + + template void DoFHandler::clear_space () { diff --cc deal.II/source/multigrid/mg_tools.cc index d60d62e20a,411b60a7e0..c81d632032 --- a/deal.II/source/multigrid/mg_tools.cc +++ b/deal.II/source/multigrid/mg_tools.cc @@@ -867,9 -866,9 +867,9 @@@ namespace MGTool { const FiniteElement& fe = dof.get_fe(); const unsigned int n_comp = fe.n_components(); - + Assert ((level>=1) && (level::cell_iterator - cell = dof.begin(), - endc = dof.end(); + cell = dof.begin(), + endc = dof.end(); for (; cell!=endc; ++cell) - { - const FiniteElement &fe = cell->get_fe(); - const unsigned int level = cell->level(); - local_dofs.resize(fe.dofs_per_face); - - for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; - ++face_no) - if (cell->at_boundary(face_no) == true) - { - const typename MGDoFHandler::face_iterator - face = cell->face(face_no); - const types::boundary_id bi = face->boundary_indicator(); + { + const FiniteElement &fe = cell->get_fe(); + const unsigned int level = cell->level(); + local_dofs.resize(fe.dofs_per_face); + + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; + ++face_no) + if (cell->at_boundary(face_no) == true) + { + const typename MGDoFHandler::face_iterator + face = cell->face(face_no); - const types::boundary_id_t bi = face->boundary_indicator(); ++ const types::boundary_id bi = face->boundary_indicator(); // Face is listed in // boundary map - if (function_map.find(bi) != function_map.end()) - { - face->get_mg_dof_indices(level, local_dofs); - for (unsigned int i=0;iget_mg_dof_indices(level, local_dofs); + for (unsigned int i=0;iget_fe().dofs_per_cell; ++i) - { - const std::vector &nonzero_component_array - = cell->get_fe().get_nonzero_components (i); - for (unsigned int c=0; cget_fe().is_primitive (i), - ExcMessage ("This function can only deal with requested boundary " - "values that correspond to primitive (scalar) base " - "elements")); - } - - typename MGDoFHandler::face_iterator face = cell->face(face_no); - const types::boundary_id boundary_component = face->boundary_indicator(); - if (function_map.find(boundary_component) != function_map.end()) + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + { + const std::vector &nonzero_component_array + = cell->get_fe().get_nonzero_components (i); + for (unsigned int c=0; cget_fe().is_primitive (i), + ExcMessage ("This function can only deal with requested boundary " + "values that correspond to primitive (scalar) base " + "elements")); + } + + typename MGDoFHandler::face_iterator face = cell->face(face_no); - const types::boundary_id_t boundary_component = face->boundary_indicator(); ++ const types::boundary_id boundary_component = face->boundary_indicator(); + if (function_map.find(boundary_component) != function_map.end()) // face is of the right component - { + { // get indices, physical location and // boundary values of dofs on this // face @@@ -1644,23 -1433,23 +1644,23 @@@ std::fill (cell_dofs.begin(), cell_dofs.end(), false); for (unsigned int face_nr=0; face_nr::faces_per_cell; ++face_nr) - { - const typename DoFHandler::face_iterator face = cell->face(face_nr); - if (!face->at_boundary()) - { + { + const typename DoFHandler::face_iterator face = cell->face(face_nr); + if (!face->at_boundary()) + { //interior face - const typename MGDoFHandler::cell_iterator - neighbor = cell->neighbor(face_nr); + const typename MGDoFHandler::cell_iterator + neighbor = cell->neighbor(face_nr); - if (neighbor->level() < cell->level()) - { - for (unsigned int j=0; jlevel() < cell->level()) + { + for (unsigned int j=0; jlevel(); cell->get_mg_dof_indices (local_dof_indices); @@@ -1671,69 -1460,69 +1671,69 @@@ } -template -void -extract_non_interface_dofs (const MGDoFHandler &mg_dof_handler, - std::vector > &non_interface_dofs) -{ - Assert (non_interface_dofs.size() == mg_dof_handler.get_tria().n_levels(), - ExcDimensionMismatch (non_interface_dofs.size(), - mg_dof_handler.get_tria().n_levels())); + template + void + extract_non_interface_dofs (const MGDoFHandler &mg_dof_handler, + std::vector > &non_interface_dofs) + { + Assert (non_interface_dofs.size() == mg_dof_handler.get_tria().n_levels(), + ExcDimensionMismatch (non_interface_dofs.size(), + mg_dof_handler.get_tria().n_levels())); - const FiniteElement &fe = mg_dof_handler.get_fe(); + const FiniteElement &fe = mg_dof_handler.get_fe(); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int dofs_per_face = fe.dofs_per_face; + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int dofs_per_face = fe.dofs_per_face; - std::vector local_dof_indices (dofs_per_cell); - std::vector cell_dofs(dofs_per_cell, false); - std::vector cell_dofs_interface(dofs_per_cell, false); + std::vector local_dof_indices (dofs_per_cell); + std::vector cell_dofs(dofs_per_cell, false); + std::vector cell_dofs_interface(dofs_per_cell, false); - typename MGDoFHandler::cell_iterator cell = mg_dof_handler.begin(), - endc = mg_dof_handler.end(); + typename MGDoFHandler::cell_iterator cell = mg_dof_handler.begin(), + endc = mg_dof_handler.end(); - for (; cell!=endc; ++cell) - { - std::fill (cell_dofs.begin(), cell_dofs.end(), false); - std::fill (cell_dofs_interface.begin(), cell_dofs_interface.end(), false); + for (; cell!=endc; ++cell) + { + std::fill (cell_dofs.begin(), cell_dofs.end(), false); + std::fill (cell_dofs_interface.begin(), cell_dofs_interface.end(), false); - + - for (unsigned int face_nr=0; face_nr::faces_per_cell; ++face_nr) - { - const typename DoFHandler::face_iterator face = cell->face(face_nr); - if (!face->at_boundary()) - { + for (unsigned int face_nr=0; face_nr::faces_per_cell; ++face_nr) + { + const typename DoFHandler::face_iterator face = cell->face(face_nr); + if (!face->at_boundary()) + { //interior face - const typename MGDoFHandler::cell_iterator - neighbor = cell->neighbor(face_nr); - - if ((neighbor->level() < cell->level())) + const typename MGDoFHandler::cell_iterator + neighbor = cell->neighbor(face_nr); + + if ((neighbor->level() < cell->level())) + { + for (unsigned int j=0; jlevel(); - cell->get_mg_dof_indices (local_dof_indices); + const unsigned int level = cell->level(); + cell->get_mg_dof_indices (local_dof_indices); - for(unsigned int i=0; i diff --cc deal.II/source/multigrid/mg_tools.inst.in index 0f78643c97,b4d74691ff..ad42c3ce05 --- a/deal.II/source/multigrid/mg_tools.inst.in +++ b/deal.II/source/multigrid/mg_tools.inst.in @@@ -77,15 -77,15 +77,18 @@@ for (deal_II_dimension : DIMENSIONS const Table<2,DoFTools::Coupling>&, const Table<2,DoFTools::Coupling>&); #endif - template void count_dofs_per_component ( + template void count_dofs_per_component ( const MGDoFHandler&, std::vector >&, bool, std::vector); - template void count_dofs_per_component ( + template void count_dofs_per_component ( const MGDoFHandler&, std::vector >&, std::vector); - template void count_dofs_per_block ( + template void count_dofs_per_block ( const MGDoFHandler&, std::vector >&, std::vector); ++ template void count_dofs_per_block ( ++ const DoFHandler&, std::vector >&, ++ std::vector); #if deal_II_dimension > 1 template void make_boundary_list( @@@ -122,9 -122,9 +125,19 @@@ #endif #if deal_II_dimension == 3 - template void count_dofs_per_block<1,3> ( + template void count_dofs_per_block ( const MGDoFHandler<1,3>&, std::vector >&, std::vector); ++ ++ template void count_dofs_per_block ( ++ const DoFHandler<1,3>&, ++ std::vector >&, std::vector); ++ template void count_dofs_per_block ( ++ const DoFHandler<2,3>&, ++ std::vector >&, std::vector); ++ template void count_dofs_per_block ( ++ const DoFHandler<1,2>&, ++ std::vector >&, std::vector); #endif \} }