#include <boost/fusion/tuple/detail/preprocessed/tuple50.hpp>
#else
#error "FUSION_MAX_VECTOR_SIZE out of bounds for preprocessed headers"
--#endif
++#endif
}
template <int dim>
- void Obstacle<dim>::vector_value (const Point<dim> &p,
- Vector<double> &values) const
+ void PlasticityContactProblem<dim>::make_grid ()
{
- for (unsigned int c=0; c<this->n_components; ++c)
- values(c) = Obstacle<dim>::value (p, c);
+ std::vector<unsigned int> repet(3);
+ repet[0] = 1;//20;
+ repet[1] = 1;
+ repet[2] = 1;
+
+ Point<dim> p1 (0,0,0);
+ Point<dim> 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<GeometryInfo<dim>::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<n_refinements_local; ++step)
+ {
+ cell = triangulation.begin_active(); // Iterator ueber alle Zellen
+
+ for (; cell!=endc; ++cell)
+ for (unsigned int face=0; face<GeometryInfo<dim>::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 <int dim>
+ void PlasticityContactProblem<dim>::setup_system ()
+ {
+ // setup dofs
+ {
+ dof_handler.distribute_dofs (fe);
- // @sect3{Implementation of the <code>Step4</code> class}
-
- // Next for the implementation of the class
- // template that makes use of the functions
- // above. As before, we will write everything
-
- template <int dim>
- Step4<dim>::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<dim>(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<dim> (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 <int dim>
- void Step4<dim>::make_grid ()
- {
- std::vector<unsigned int> repet(3);
- repet[0] = 1;//20;
- repet[1] = 1;
- repet[2] = 1;
-
- Point<dim> p1 (0,0,0);
- Point<dim> 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<GeometryInfo<dim>::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<n_refinements_local; ++step)
+ // setup sparsity pattern
{
- cell = triangulation.begin_active(); // Iterator ueber alle Zellen
-
- for (; cell!=endc; ++cell)
- for (unsigned int face=0; face<GeometryInfo<dim>::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 <int dim>
- void Step4<dim>::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<end; j++)
+ diag_mass_matrix_vector (j) = mass_matrix.diag_element (j);
+ number_iterations = 0;
+
+ diag_mass_matrix_vector.compress ();
+ }
}
- // setup hanging nodes and dirichlet constraints
+ template <int dim>
+ void PlasticityContactProblem<dim>::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<dim> quadrature_formula(2);
+ QGauss<dim-1> face_quadrature_formula(2);
+
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ UpdateFlags(update_values |
+ update_gradients |
+ update_q_points |
+ update_JxW_values));
+
+ FEFaceValues<dim> 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<dim> right_hand_side;
+ std::vector<Vector<double> > right_hand_side_values (n_q_points,
+ Vector<double>(dim));
+ std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
+ Vector<double>(dim));
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+
- std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
++ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+ typename DoFHandler<dim>::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<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
+ fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ {
+ SymmetricTensor<4,dim> 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 <<std::endl;
+ // std::cout<< stress_strain_tensor <<std::endl;
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ stress_tensor = stress_strain_tensor_linearized * plast_lin_hard->get_strain(fe_values, i, q_point);
+
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ cell_matrix(i,j) += (stress_tensor *
+ plast_lin_hard->get_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<GeometryInfo<dim>::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<n_face_q_points; ++q_point)
+ {
+ Tensor<1,dim> rhs_values;
+ rhs_values = 0;
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
+ rhs_values *
+ fe_values_face.JxW (q_point));
+ }
+ }
+ }
+
+ cell->get_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 <int dim>
+ void PlasticityContactProblem<dim>::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<dim> quadrature_formula(2);
+ QGauss<dim-1> face_quadrature_formula(2);
+
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ UpdateFlags(update_values |
+ update_gradients |
+ update_q_points |
+ update_JxW_values));
+
+ FEFaceValues<dim> 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<dim> right_hand_side;
+ std::vector<Vector<double> > right_hand_side_values (n_q_points,
+ Vector<double>(dim));
+ std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
+ Vector<double>(dim));
+
+ Vector<double> cell_rhs (dofs_per_cell);
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ const FEValuesExtractors::Vector displacement (0);
+
+ typename DoFHandler<dim>::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<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
+ fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ {
+ SymmetricTensor<4,dim> 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; i<dofs_per_cell; ++i)
+ {
+ cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor * //(stress_tensor) *
+ plast_lin_hard->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<GeometryInfo<dim>::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<n_face_q_points; ++q_point)
+ {
+ Tensor<1,dim> rhs_values;
+ rhs_values = 0;
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
+ rhs_values *
+ fe_values_face.JxW (q_point));
+ }
+ }
+ }
+
+ cell->get_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 <<std::endl;
+ pcout<< "Plast-Points = " << sum_plast_points <<std::endl;
}
- // setup sparsity pattern
+ template <int dim>
+ void PlasticityContactProblem<dim>::assemble_mass_matrix_diagonal (TrilinosWrappers::SparseMatrix &mass_matrix)
{
- TrilinosWrappers::SparsityPattern sp (locally_owned_dofs,
- mpi_communicator);
+ QTrapez<dim-1> face_quadrature_formula;
- DoFTools::make_sparsity_pattern (dof_handler, sp, constraints_dirichlet_hanging_nodes, false,
- Utilities::MPI::this_mpi_process(mpi_communicator));
+ FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
+ update_values |
+ update_quadrature_points |
+ update_JxW_values);
- sp.compress();
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_face_q_points = face_quadrature_formula.size();
- system_matrix_newton.reinit (sp);
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- mass_matrix.reinit (sp);
- }
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
- assemble_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<end; j++)
- diag_mass_matrix_vector (j) = mass_matrix.diag_element (j);
- number_iterations = 0;
-
- diag_mass_matrix_vector.compress ();
- }
+ const FEValuesExtractors::Vector displacement (0);
- template <int dim>
- void Step4<dim>::assemble_mass_matrix ()
- {
- QTrapez<dim-1> face_quadrature_formula;
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
- FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
- update_values | update_quadrature_points | update_JxW_values);
+ for (; cell!=endc; ++cell)
+ if (cell->is_locally_owned())
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face (face)->at_boundary()
+ && cell->face (face)->boundary_indicator () == 9)
+ {
+ fe_values_face.reinit (cell, face);
+ cell_matrix = 0;
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
- const unsigned int dofs_per_face = fe.dofs_per_face;
- const unsigned int n_face_q_points = face_quadrature_formula.size();
+ for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ cell_matrix(i,i) += (fe_values_face[displacement].value (i, q_point) *
+ fe_values_face[displacement].value (i, q_point) *
+ fe_values_face.JxW (q_point));
- FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ cell->get_dof_indices (local_dof_indices);
- std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+ constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_matrix,
+ local_dof_indices,
+ mass_matrix);
+ }
- const FEValuesExtractors::Vector displacement (0);
+ mass_matrix.compress ();
+ }
- typename DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ // @sect4{PlasticityContactProblem::update_solution_and_constraints}
- for (; cell!=endc; ++cell)
- if (cell->is_locally_owned())
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- if (cell->face (face)->at_boundary()
- && cell->face (face)->boundary_indicator () == 9)
- {
- fe_values_face.reinit (cell, face);
- cell_matrix = 0;
-
- for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- cell_matrix(i,i) += (fe_values_face[displacement].value (i, q_point) *
- fe_values_face[displacement].value (i, q_point) *
- fe_values_face.JxW (q_point));
-
- cell->get_dof_indices (local_dof_indices);
-
- constraints_dirichlet_hanging_nodes.distribute_local_to_global (cell_matrix,
- local_dof_indices,
- mass_matrix);
- }
-
- mass_matrix.compress ();
- }
+ // Projection and updating of the active set
+ // for the dofs which penetrates the obstacle.
+ template <int dim>
+ void PlasticityContactProblem<dim>::update_solution_and_constraints ()
+ {
+ clock_t start_proj, end_proj;
- template <int dim>
- void Step4<dim>::assemble_nl_system (TrilinosWrappers::MPI::Vector &u)
- {
- QGauss<dim> quadrature_formula(2);
- QGauss<dim-1> face_quadrature_formula(2);
+ const EquationData::Obstacle<dim> obstacle;
+ std::vector<bool> vertex_touched (dof_handler.n_dofs (), false);
- FEValues<dim> fe_values (fe, quadrature_formula,
- UpdateFlags(update_values |
- update_gradients |
- update_q_points |
- update_JxW_values));
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
- FEFaceValues<dim> fe_values_face (fe, face_quadrature_formula,
- update_values | update_quadrature_points |
- update_JxW_values);
+ 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<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face (face)->at_boundary()
+ && cell->face (face)->boundary_indicator () == 9)
+ for (unsigned int v=0; v<GeometryInfo<dim-1>::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<dim> 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
+ // <<std::endl;
+ }
+ }
+ distributed_solution.compress(Insert);
+
+ unsigned int sum_contact_constraints = Utilities::MPI::sum(active_set_locally_owned.n_elements (),
+ mpi_communicator);
+ pcout << "Number of Contact-Constaints: " << sum_contact_constraints <<std::endl;
+
+ solution = distributed_solution;
+
+ constraints.close ();
+
+ const ConstraintMatrix::MergeConflictBehavior
+ merge_conflict_behavior = ConstraintMatrix::left_object_wins;
+ constraints.merge (constraints_dirichlet_hanging_nodes, merge_conflict_behavior);
+ }
- 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();
+ template <int dim>
+ void PlasticityContactProblem<dim>::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<bool> 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<dim>(),
+ 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<dim>(),
+ constraints_dirichlet_hanging_nodes,
+ component_mask);
+ constraints_dirichlet_hanging_nodes.close ();
+ }
- const EquationData::RightHandSide<dim> right_hand_side;
- std::vector<Vector<double> > right_hand_side_values (n_q_points,
- Vector<double>(dim));
- std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
- Vector<double>(dim));
+ template <int dim>
+ void PlasticityContactProblem<dim>::solve ()
+ {
+ pcout << "Solving ..." << std::endl;
+ Timer t;
- FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- Vector<double> cell_rhs (dofs_per_cell);
+ TrilinosWrappers::MPI::Vector distributed_solution (system_rhs_newton);
+ distributed_solution = solution;
- std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+ constraints_hanging_nodes.set_zero (distributed_solution);
- typename DoFHandler<dim>::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<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
- fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);
-
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- {
- SymmetricTensor<4,dim> 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 <<std::endl;
- // std::cout<< stress_strain_tensor <<std::endl;
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- stress_tensor = stress_strain_tensor_linearized * plast_lin_hard->get_strain(fe_values, i, q_point);
-
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- {
- cell_matrix(i,j) += (stress_tensor *
- plast_lin_hard->get_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<GeometryInfo<dim>::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<n_face_q_points; ++q_point)
- {
- Tensor<1,dim> rhs_values;
- rhs_values = 0;
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
- rhs_values *
- fe_values_face.JxW (q_point));
- }
- }
- }
-
- cell->get_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 <int dim>
- void Step4<dim>::residual_nl_system (TrilinosWrappers::MPI::Vector &u,
- Vector<double> &sigma_eff_vector)
- {
- QGauss<dim> quadrature_formula(2);
- QGauss<dim-1> face_quadrature_formula(2);
-
- FEValues<dim> fe_values (fe, quadrature_formula,
- UpdateFlags(update_values |
- update_gradients |
- update_q_points |
- update_JxW_values));
-
- FEFaceValues<dim> 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<dim> right_hand_side;
- std::vector<Vector<double> > right_hand_side_values (n_q_points,
- Vector<double>(dim));
- std::vector<Vector<double> > right_hand_side_values_face (n_face_q_points,
- Vector<double>(dim));
-
- Vector<double> cell_rhs (dofs_per_cell);
- Vector<double> cell_sigma_eff (dofs_per_cell);
-
- std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
- const FEValuesExtractors::Vector displacement (0);
-
- typename DoFHandler<dim>::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<SymmetricTensor<2,dim> > strain_tensor (n_q_points);
- fe_values[displacement].get_function_symmetric_gradients (u, strain_tensor);
-
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- {
- SymmetricTensor<4,dim> 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 <<std::endl;*/
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- cell_rhs(i) -= (strain_tensor[q_point] * stress_strain_tensor * //(stress_tensor) *
- plast_lin_hard->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<GeometryInfo<dim>::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<n_face_q_points; ++q_point)
- {
- Tensor<1,dim> rhs_values;
- rhs_values = 0;
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- cell_rhs(i) += (fe_values_face[displacement].value (i, q_point) *
- rhs_values *
- fe_values_face.JxW (q_point));
- }
- }
- }
-
- cell->get_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<TrilinosWrappers::MPI::Vector>
+ // 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 <<std::endl;
- pcout<< "Plast-Points = " << sum_plast_points <<std::endl;
- }
+ PrimitiveVectorMemory<TrilinosWrappers::MPI::Vector> 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<TrilinosWrappers::MPI::Vector>
+ solver(solver_control, mem,
+ SolverFGMRES<TrilinosWrappers::MPI::Vector>::
+ AdditionalData(30, true));
+ solver.solve(system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u);
- // @sect4{Step4::projection_active_set}
+ pcout << "Initial error: " << solver_control.initial_value() <<std::endl;
+ pcout << " " << solver_control.last_step()
+ << " CG iterations needed to obtain convergence with an error: "
+ << solver_control.last_value()
+ << std::endl;
- // Projection and updating of the active set
- // for the dofs which penetrates the obstacle.
- template <int dim>
- void Step4<dim>::projection_active_set ()
- {
- clock_t start_proj, end_proj;
-
- const EquationData::Obstacle<dim> obstacle;
- std::vector<bool> vertex_touched (dof_handler.n_dofs (), false);
-
- typename DoFHandler<dim>::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<GeometryInfo<dim>::faces_per_cell; ++face)
- if (cell->face (face)->at_boundary()
- && cell->face (face)->boundary_indicator () == 9)
- for (unsigned int v=0; v<GeometryInfo<dim-1>::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<dim> 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
- // <<std::endl;
- }
- }
- distributed_solution.compress(Insert);
-
- unsigned int sum_contact_constraints = Utilities::MPI::sum(active_set_locally_owned.n_elements (),
- mpi_communicator);
- pcout << "Number of Contact-Constaints: " << sum_contact_constraints <<std::endl;
-
- solution = distributed_solution;
-
- constraints.close ();
-
- const ConstraintMatrix::MergeConflictBehavior
- merge_conflict_behavior = ConstraintMatrix::left_object_wins;
- constraints.merge (constraints_dirichlet_hanging_nodes, merge_conflict_behavior);
- }
+ MPI_Barrier (mpi_communicator);
+ t.stop();
+ if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+ run_time[7] += t.wall_time();
- template <int dim>
- void Step4<dim>::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<bool> 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<dim>(),
- 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<dim>(),
- constraints_dirichlet_hanging_nodes,
- component_mask);
- constraints_dirichlet_hanging_nodes.close ();
- }
+ number_iterations += solver_control.last_step();
- template <int dim>
- void Step4<dim>::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<TrilinosWrappers::MPI::Vector>
- // solver (reduction_control, mpi_communicator);
- // solver.solve (system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u);
-
- PrimitiveVectorMemory<TrilinosWrappers::MPI::Vector> 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<TrilinosWrappers::MPI::Vector>
- solver(solver_control, mem,
- SolverFGMRES<TrilinosWrappers::MPI::Vector>::
- AdditionalData(30, true));
- solver.solve(system_matrix_newton, distributed_solution, system_rhs_newton, preconditioner_u);
-
- pcout << "Initial error: " << solver_control.initial_value() <<std::endl;
- pcout << " " << solver_control.last_step()
- << " CG iterations needed to obtain convergence with an error: "
- << solver_control.last_value()
- << std::endl;
-
- MPI_Barrier (mpi_communicator);
- t.stop();
- if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
- run_time[7] += t.wall_time();
-
- number_iterations += solver_control.last_step();
-
- constraints.distribute (distributed_solution);
-
- solution = distributed_solution;
- }
-
- template <int dim>
- void Step4<dim>::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<std::vector<bool> > constant_modes;
- std::vector<bool> 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<double> 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<< " " <<std::endl;
- pcout<< j << ". Iteration of the inexact Newton-method." <<std::endl;
- pcout<< "Update of active set" <<std::endl;
-
- MPI_Barrier (mpi_communicator);
- t.restart();
-
- projection_active_set ();
-
- MPI_Barrier (mpi_communicator);
- t.stop();
- if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
- run_time[5] += t.wall_time();
-
- pcout<< "Assembling ... " <<std::endl;
- MPI_Barrier (mpi_communicator);
- t.restart();
- system_matrix_newton = 0;
- system_rhs_newton = 0;
- assemble_nl_system (solution); //compute Newton-Matrix
- MPI_Barrier (mpi_communicator);
- t.stop();
- if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
- run_time[1] += t.wall_time();
-
- number_assemble_system += 1;
-
- MPI_Barrier (mpi_communicator);
- t.restart();
- solve ();
- MPI_Barrier (mpi_communicator);
- t.stop();
- if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
- run_time[2] += t.wall_time();
-
- TrilinosWrappers::MPI::Vector distributed_solution (system_rhs_newton);
- distributed_solution = solution;
-
- int damped = 0;
- tmp_vector = old_solution;
- double a = 0;
- for (unsigned int i=0; (i<10)&&(!damped); i++)
- {
- a=pow(0.5, static_cast<double>(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<end_res; ++n)
- if (constraints.is_inhomogeneously_constrained (n))
- {
- // pcout<< i << ". " << constraints.get_inhomogeneity (n)
- // << ". " << res (n)
- // << ", start = " << start_res
- // << ", end = " << end_res
- // <<std::endl;
- res(n) = 0;
- }
-
- resid = res.l2_norm ();
- pcout<< "Residual: " << resid <<std::endl;
-
- if (resid<resid_old)
- {
- pcout<< "Newton-damping parameter alpha = " << a <<std::endl;
- damped=1;
- }
- MPI_Barrier (mpi_communicator);
- t.stop();
- if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
- run_time[3] += t.wall_time();
- }
-
- if (resid<1e-8)
- {
- pcout<< "Inexact Newton-method stopped with residual = " << resid <<std::endl;
- pcout<< "Number of Assembling systems = " << number_assemble_system <<std::endl;
- break;
- }
- resid_old=resid;
-
- resid_vector = system_rhs_newton;
-
- if (active_set == active_set_old && resid < 1e-10)
- break;
- active_set_old = active_set;
- } // End of active-set-loop
-
-
- pcout<< "Creating output." <<std::endl;
- MPI_Barrier (mpi_communicator);
- t.restart();
- std::ostringstream filename_solution;
- filename_solution << "solution";
- // filename_solution << "solution_";
- // filename_solution << k;
- output_results (filename_solution.str ());
- // output_results (sigma_eff_vector, "sigma_eff");
- MPI_Barrier (mpi_communicator);
- t.stop();
- if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
- run_time[4] += t.wall_time();
-
- pcout<< "Number of Solver-Iterations = " << number_iterations <<std::endl;
-
- pcout<< "%%%%%% Rechenzeit make grid and setup = " << run_time[0] <<std::endl;
- pcout<< "%%%%%% Rechenzeit projection active set = " << run_time[5] <<std::endl;
- pcout<< "%%%%%% Rechenzeit assemble system = " << run_time[1] <<std::endl;
- pcout<< "%%%%%% Rechenzeit solve system = " << run_time[2] <<std::endl;
- pcout<< "%%%%%% Rechenzeit preconditioner = " << run_time[6] <<std::endl;
- pcout<< "%%%%%% Rechenzeit solve with CG = " << run_time[7] <<std::endl;
- pcout<< "%%%%%% Rechenzeit error and lambda = " << run_time[3] <<std::endl;
- pcout<< "%%%%%% Rechenzeit output = " << run_time[4] <<std::endl;
- }
+ constraints.distribute (distributed_solution);
- template <int dim>
- void Step4<dim>::output_results (const std::string& title) const
- {
- move_mesh (solution);
+ solution = distributed_solution;
+ }
- TrilinosWrappers::MPI::Vector lambda (solution);
- lambda = resid_vector;
+ template <int dim>
+ void PlasticityContactProblem<dim>::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<std::vector<bool> > constant_modes;
+ std::vector<bool> 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);
+ unsigned int j = 0;
+ unsigned int number_assemble_system = 0;
+ for (; j<=100;j++)
+ {
+ pcout<< " " <<std::endl;
+ pcout<< j << ". Iteration of the inexact Newton-method." <<std::endl;
+ pcout<< "Update of active set" <<std::endl;
+
+ MPI_Barrier (mpi_communicator);
+ t.restart();
+
+ update_solution_and_constraints ();
+
+ MPI_Barrier (mpi_communicator);
+ t.stop();
+ if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+ run_time[5] += t.wall_time();
+
+ pcout<< "Assembling ... " <<std::endl;
+ MPI_Barrier (mpi_communicator);
+ t.restart();
+ system_matrix_newton = 0;
+ system_rhs_newton = 0;
+ assemble_nl_system (solution); //compute Newton-Matrix
+ MPI_Barrier (mpi_communicator);
+ t.stop();
+ if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+ run_time[1] += t.wall_time();
+
+ number_assemble_system += 1;
+
+ MPI_Barrier (mpi_communicator);
+ t.restart();
+ solve ();
+ MPI_Barrier (mpi_communicator);
+ t.stop();
+ if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+ run_time[2] += t.wall_time();
+
+ TrilinosWrappers::MPI::Vector distributed_solution (system_rhs_newton);
+ distributed_solution = solution;
+
+ int damped = 0;
+ tmp_vector = old_solution;
+ double a = 0;
+ for (unsigned int i=0; (i<10)&&(!damped); i++)
+ {
+ a=pow(0.5, static_cast<double>(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<end_res; ++n)
+ if (constraints.is_inhomogeneously_constrained (n))
+ {
+ // pcout<< i << ". " << constraints.get_inhomogeneity (n)
+ // << ". " << res (n)
+ // << ", start = " << start_res
+ // << ", end = " << end_res
+ // <<std::endl;
+ res(n) = 0;
+ }
+
+ resid = res.l2_norm ();
+ pcout<< "Residual: " << resid <<std::endl;
+
+ if (resid<resid_old)
+ {
+ pcout<< "Newton-damping parameter alpha = " << a <<std::endl;
+ damped=1;
+ }
+ MPI_Barrier (mpi_communicator);
+ t.stop();
+ if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+ run_time[3] += t.wall_time();
+ }
+
+ if (resid<1e-8)
+ {
+ pcout<< "Inexact Newton-method stopped with residual = " << resid <<std::endl;
+ pcout<< "Number of Assembling systems = " << number_assemble_system <<std::endl;
+ break;
+ }
+ resid_old=resid;
+
+ resid_vector = system_rhs_newton;
+
+ if (active_set == active_set_old && resid < 1e-10)
+ break;
+ active_set_old = active_set;
+ } // End of active-set-loop
+
+
+ pcout<< "Creating output." <<std::endl;
+ MPI_Barrier (mpi_communicator);
+ t.restart();
+ std::ostringstream filename_solution;
+ filename_solution << "solution";
+ // filename_solution << "solution_";
+ // filename_solution << k;
+ output_results (filename_solution.str ());
+ MPI_Barrier (mpi_communicator);
+ t.stop();
+ if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+ run_time[4] += t.wall_time();
+
+ pcout<< "Number of Solver-Iterations = " << number_iterations <<std::endl;
+
+ pcout<< "%%%%%% Rechenzeit make grid and setup = " << run_time[0] <<std::endl;
+ pcout<< "%%%%%% Rechenzeit projection active set = " << run_time[5] <<std::endl;
+ pcout<< "%%%%%% Rechenzeit assemble system = " << run_time[1] <<std::endl;
+ pcout<< "%%%%%% Rechenzeit solve system = " << run_time[2] <<std::endl;
+ pcout<< "%%%%%% Rechenzeit preconditioner = " << run_time[6] <<std::endl;
+ pcout<< "%%%%%% Rechenzeit solve with CG = " << run_time[7] <<std::endl;
+ pcout<< "%%%%%% Rechenzeit error and lambda = " << run_time[3] <<std::endl;
+ pcout<< "%%%%%% Rechenzeit output = " << run_time[4] <<std::endl;
+ }
- DataOut<dim> data_out;
+ template <int dim>
+ void PlasticityContactProblem<dim>::output_results (const std::string& title) const
+ {
+ move_mesh (solution);
- data_out.attach_dof_handler (dof_handler);
+ TrilinosWrappers::MPI::Vector lambda (solution);
+ lambda = resid_vector;
- const std::vector<DataComponentInterpretation::DataComponentInterpretation>
- data_component_interpretation
- (dim, DataComponentInterpretation::component_is_part_of_vector);
- data_out.add_data_vector (solution, std::vector<std::string>(dim, "Displacement"),
- DataOut<dim>::type_dof_data,
- data_component_interpretation);
- data_out.add_data_vector (lambda, std::vector<std::string>(dim, "Residual"),
- DataOut<dim>::type_dof_data,
- data_component_interpretation);
- data_out.add_data_vector (active_set, std::vector<std::string>(dim, "ActiveSet"),
- DataOut<dim>::type_dof_data,
- data_component_interpretation);
+ DataOut<dim> data_out;
- Vector<float> subdomain (triangulation.n_active_cells());
- for (unsigned int i=0; i<subdomain.size(); ++i)
- subdomain(i) = triangulation.locally_owned_subdomain();
- data_out.add_data_vector (subdomain, "subdomain");
+ data_out.attach_dof_handler (dof_handler);
- data_out.build_patches ();
+ const std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ data_component_interpretation
+ (dim, DataComponentInterpretation::component_is_part_of_vector);
+ data_out.add_data_vector (solution, std::vector<std::string>(dim, "Displacement"),
+ DataOut<dim>::type_dof_data,
+ data_component_interpretation);
+ data_out.add_data_vector (lambda, std::vector<std::string>(dim, "Residual"),
+ DataOut<dim>::type_dof_data,
+ data_component_interpretation);
+ data_out.add_data_vector (active_set, std::vector<std::string>(dim, "ActiveSet"),
+ DataOut<dim>::type_dof_data,
+ data_component_interpretation);
- const std::string filename = (title + "-" +
- Utilities::int_to_string
- (triangulation.locally_owned_subdomain(), 4));
+ Vector<float> subdomain (triangulation.n_active_cells());
+ for (unsigned int i=0; i<subdomain.size(); ++i)
+ subdomain(i) = triangulation.locally_owned_subdomain();
+ data_out.add_data_vector (subdomain, "subdomain");
- std::ofstream output_vtu ((filename + ".vtu").c_str ());
- data_out.write_vtu (output_vtu);
+ data_out.build_patches ();
- if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
- {
- std::vector<std::string> filenames;
- for (unsigned int i=0;
- i<Utilities::MPI::n_mpi_processes(mpi_communicator);
- ++i)
- filenames.push_back ("solution-" +
- Utilities::int_to_string (i, 4) +
- ".vtu");
-
- std::ofstream master_output ((filename + ".pvtu").c_str());
- data_out.write_pvtu_record (master_output, filenames);
- }
+ const std::string filename = (title + "-" +
+ Utilities::int_to_string
+ (triangulation.locally_owned_subdomain(), 4));
- TrilinosWrappers::MPI::Vector tmp (solution);
- tmp *= -1;
- move_mesh (tmp);
- }
+ std::ofstream output_vtu ((filename + ".vtu").c_str ());
+ data_out.write_vtu (output_vtu);
- template <int dim>
- void Step4<dim>::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const
- {
- pcout<< "Moving mesh." <<std::endl;
-
- std::vector<bool> vertex_touched (triangulation.n_vertices(),
- false);
-
- for (typename DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active ();
- cell != dof_handler.end(); ++cell)
- if (cell->is_locally_owned())
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
- {
- if (vertex_touched[cell->vertex_index(v)] == false)
- {
- vertex_touched[cell->vertex_index(v)] = true;
-
- Point<dim> vertex_displacement;
- for (unsigned int d=0; d<dim; ++d)
- {
- if (_complete_displacement(cell->vertex_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<std::string> filenames;
+ for (unsigned int i=0;
+ i<Utilities::MPI::n_mpi_processes(mpi_communicator);
+ ++i)
+ filenames.push_back ("solution-" +
+ Utilities::int_to_string (i, 4) +
+ ".vtu");
+
+ std::ofstream master_output ((filename + ".pvtu").c_str());
+ data_out.write_pvtu_record (master_output, filenames);
+ }
- template <int dim>
- void Step4<dim>::output_results (TrilinosWrappers::MPI::Vector vector,
- const std::string& title) const
- {
- DataOut<dim> 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 <int dim>
+ void PlasticityContactProblem<dim>::move_mesh (const TrilinosWrappers::MPI::Vector &_complete_displacement) const
+ {
+ pcout<< "Moving mesh." <<std::endl;
+
+ std::vector<bool> vertex_touched (triangulation.n_vertices(),
+ false);
+
+ for (typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active ();
+ cell != dof_handler.end(); ++cell)
+ if (cell->is_locally_owned())
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+ {
+ if (vertex_touched[cell->vertex_index(v)] == false)
+ {
+ vertex_touched[cell->vertex_index(v)] = true;
+
+ Point<dim> vertex_displacement;
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ if (_complete_displacement(cell->vertex_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 <int dim>
+ void PlasticityContactProblem<dim>::output_results (TrilinosWrappers::MPI::Vector vector,
+ const std::string& title) const
+ {
+ DataOut<dim> 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 <int dim>
- void Step4<dim>::output_results (Vector<double> vector, const std::string& title) const
- {
- DataOut<dim> 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 <int dim>
+ void PlasticityContactProblem<dim>::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 <int dim>
- void Step4<dim>::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 <code>main</code> 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
- // <code>Step4@<2@></code> (forcing
- // the compiler to compile the class template
- // with <code>dim</code> replaced by
- // <code>2</code>) 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
- // <code>laplace_problem_2d</code> 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
- // <code>laplace_problem_2d</code> 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 <code>deallog.depth_console(0)</code>
- // 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;
* the @ref distributed module for
* more information.
*/
- const unsigned int artificial_subdomain_id = static_cast<subdomain_id_t>(-2);
+ const unsigned int artificial_subdomain_id = static_cast<subdomain_id>(-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<global_dof_index>(-1);
+
/**
* The type used to denote boundary indicators associated with every
* piece of the boundary and, in the case of meshes that describe
* 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<material_id_t>(-1);
+ const material_id invalid_material_id = static_cast<material_id>(-1);
+
+
}
* object to the finite element.
*/
virtual void distribute_dofs (const FiniteElement<dim,spacedim> &fe,
- const unsigned int offset = 0);
+ const types::global_dof_index offset = 0);
-
++
+ virtual void distribute_mg_dofs (const FiniteElement<dim, spacedim>& fe, const types::global_dof_index offset = 0);
/**
* After distribute_dofs() with
* 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
* indicators under
* consideration.
*/
- unsigned int
+ types::global_dof_index
- n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const;
+ n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const;
/**
* Access to an object informing
*/
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 <int structdim>
+ 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<int structdim>
+ 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
* Triangulation objects.
*/
std::vector<dealii::internal::DoFHandler::DoFLevel<dim>*> levels;
-
+
+ std::vector<dealii::internal::DoFHandler::DoFLevel<dim>*> mg_levels;
+
/**
* Space to store DoF numbers of
* faces. They are not stored in
* array.
*/
dealii::internal::DoFHandler::DoFFaces<dim> *faces;
-
+
+ dealii::internal::DoFHandler::DoFFaces<dim>* mg_faces;
-
++
+ std::vector<MGVertexDoFs> mg_vertex_dofs;
-
++
+ std::vector<types::global_dof_index> mg_used_dofs;
+
/**
* Make accessor objects friends.
*/
#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<types::boundary_id> &) 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<types::boundary_id_t> &) const;
++template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &) const;
/* ----------------------- Inline functions ---------------------------------- */
template <class DH, class SparsityPattern>
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
template <class DH, class SparsityPattern>
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
template<class DH, class SparsityPattern>
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
extract_boundary_dofs (const DH &dof_handler,
const std::vector<bool> &component_mask,
std::vector<bool> &selected_dofs,
- const std::set<types::boundary_id_t> &boundary_indicators = std::set<types::boundary_id_t>());
+ const std::set<types::boundary_id> &boundary_indicators = std::set<types::boundary_id>());
+ /**
+ * This function does the same as the previous one but it
+ * returns its result as an IndexSet rather than a std::vector@<bool@>.
+ * 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 <class DH>
+ void
+ extract_boundary_dofs (const DH &dof_handler,
+ const std::vector<bool> &component_mask,
+ IndexSet &selected_dofs,
+ const std::set<types::boundary_id_t> &boundary_indicators = std::set<types::boundary_id_t>());
+
/**
* This function does the same as the previous one but it
* returns its result as an IndexSet rather than a std::vector@<bool@>.
template <class DH>
void
extract_dofs_with_support_on_boundary (const DH &dof_handler,
- const std::vector<bool> &component_mask,
- std::vector<bool> &selected_dofs,
- const std::set<types::boundary_id> &boundary_indicators = std::set<types::boundary_id>());
+ const std::vector<bool> &component_mask,
+ std::vector<bool> &selected_dofs,
- const std::set<types::boundary_id_t> &boundary_indicators = std::set<types::boundary_id_t>());
++ const std::set<types::boundary_id> &boundary_indicators = std::set<types::boundary_id>());
/**
* @name Hanging Nodes
template <class DH>
void
extract_subdomain_dofs (const DH &dof_handler,
- const types::subdomain_id_t subdomain_id,
- const types::subdomain_id subdomain_id,
- std::vector<bool> &selected_dofs);
++ const types::subdomain_id subdomain_id,
+ std::vector<bool> &selected_dofs);
/**
template <class DH>
void
get_subdomain_association (const DH &dof_handler,
- std::vector<types::subdomain_id_t> &subdomain);
- std::vector<types::subdomain_id> &subdomain);
++ std::vector<types::subdomain_id> &subdomain);
/**
* Count how many degrees of freedom are
template <class DH>
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
template <class DH>
void
count_dofs_with_subdomain_association (const DH &dof_handler,
- const types::subdomain_id_t subdomain,
- const types::subdomain_id subdomain,
- std::vector<unsigned int> &n_dofs_on_subdomain);
++ const types::subdomain_id subdomain,
+ std::vector<unsigned int> &n_dofs_on_subdomain);
/**
* Return a set of indices that denotes
template <class DH>
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
template <class DH>
void
map_dof_to_boundary_indices (const DH &dof_handler,
- const std::set<types::boundary_id_t> &boundary_indicators,
- const std::set<types::boundary_id> &boundary_indicators,
- std::vector<unsigned int> &mapping);
++ const std::set<types::boundary_id> &boundary_indicators,
+ std::vector<types::global_dof_index> &mapping);
/**
* Return a list of support
* indicators under
* consideration.
*/
- unsigned int
+ types::global_dof_index
- n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const;
+ n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const;
/**
* Return the number of
* 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<int structdim>
+ 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<int structdim>
+ 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
* the present class.
*/
std::vector<unsigned int> vertex_dofs_offsets;
-
+
+ std::vector<MGVertexDoFs> mg_vertex_dofs;
+
/**
* Array to store the
* information, if a cell on
return *tria;
}
- types::global_dof_index DoFHandler<dim, spacedim>::MGVertexDoFs::get_index (const unsigned int level, const unsigned int dof_number) const {
+ template<int dim, int spacedim>
+ inline
-
++ DoFHandler<dim, spacedim>::MGVertexDoFs::MGVertexDoFs()
++ {
++ Assert (false, ExcNotImplemented ());
++ }
++
++ template<int dim, int spacedim>
++ inline
++ DoFHandler<dim, spacedim>::MGVertexDoFs::~MGVertexDoFs()
++ {
++ Assert (false, ExcNotImplemented ());
++ }
++
++ template<int dim, int spacedim>
++ inline
++ types::global_dof_index DoFHandler<dim, spacedim>::MGVertexDoFs::get_index (const unsigned int,
++ const unsigned int) const
++ {
+ Assert (false, ExcNotImplemented ());
+ return invalid_dof_index;
+ }
- void DoFHandler<dim, spacedim>::MGVertexDoFs::set_index (const unsigned int level, const unsigned int dof_number, const types::global_dof_index index) {
+
+ template<int dim, int spacedim>
+ inline
++ void DoFHandler<dim, spacedim>::MGVertexDoFs::set_index (const unsigned int,
++ const unsigned int,
++ types::global_dof_index)
++ {
+ Assert (false, ExcNotImplemented ());
+ }
#endif
{
std::size_t result = sizeof(*this);
result += sizeof(unsigned int) * sizes.size();
- #ifdef DEAL_PREFER_MATRIX_EZ
- std::vector<std_cxx1x::shared_ptr<SparseMatrixEZ<double> > >::const_iterator m;
- const std::vector<std_cxx1x::shared_ptr<SparseMatrixEZ<double> > >::const_iterator end = prolongation_matrices.end();
- for (m = prolongation_matrices.begin(); m != end ; ++m)
- result += *m->memory_consumption();
- #else
+
for (unsigned int i=0;i<prolongation_matrices.size();++i)
result += prolongation_matrices[i]->memory_consumption()
- + prolongation_sparsities[i]->memory_consumption();
+ + prolongation_sparsities[i]->memory_consumption();
- #endif
+
return result;
}
void
interpolate_boundary_values (const Mapping<DH::dimension,DH::space_dimension> &mapping,
const DH &dof,
- const types::boundary_id_t boundary_component,
+ const types::boundary_id boundary_component,
const Function<DH::space_dimension> &boundary_function,
- std::map<unsigned int,double> &boundary_values,
+ std::map<types::global_dof_index,double> &boundary_values,
const std::vector<bool> &component_mask = std::vector<bool>());
/**
template <class DH>
void
interpolate_boundary_values (const DH &dof,
- const types::boundary_id_t boundary_component,
+ const types::boundary_id boundary_component,
const Function<DH::space_dimension> &boundary_function,
- std::map<unsigned int,double> &boundary_values,
+ std::map<types::global_dof_index,double> &boundary_values,
const std::vector<bool> &component_mask = std::vector<bool>());
void
interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension> &mapping,
const DH &dof,
- const types::boundary_id_t boundary_component,
+ const types::boundary_id boundary_component,
const Function<DH::space_dimension> &boundary_function,
- std::map<unsigned int,double> &boundary_values,
+ std::map<types::global_dof_index,double> &boundary_values,
const std::vector<bool> &component_mask)
{
typename FunctionMap<DH::space_dimension>::type function_map;
template <class DH>
void
interpolate_boundary_values (const DH &dof,
- const types::boundary_id_t boundary_component,
+ const types::boundary_id boundary_component,
const Function<DH::space_dimension> &boundary_function,
- std::map<unsigned int,double> &boundary_values,
+ std::map<types::global_dof_index,double> &boundary_values,
const std::vector<bool> &component_mask)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
else
AssertDimension (dof.get_fe().n_components(), component_mapping.size());
- std::vector<unsigned int> dof_to_boundary_mapping;
+ std::vector<types::global_dof_index> dof_to_boundary_mapping;
- std::set<types::boundary_id_t> selected_boundary_components;
+ std::set<types::boundary_id> selected_boundary_components;
for (typename FunctionMap<spacedim>::type::const_iterator i=boundary_functions.begin();
i!=boundary_functions.end(); ++i)
selected_boundary_components.insert (i->first);
* was selected when calling
* @p distribute_dofs the last time.
*/
- template <int spacedim>
- 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; i<dof_handler.tria->n_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 <int spacedim>
- 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; i<dof_handler.tria->n_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 <int spacedim>
- 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; i<dof_handler.tria->n_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 <int spacedim>
+ 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; i<dof_handler.tria->n_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 <int spacedim>
+ 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; i<dof_handler.tria->n_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 <int spacedim>
+ 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; i<dof_handler.tria->n_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<int spacedim>
+ 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<unsigned int> (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<unsigned int> max_level (n_vertices, 0);
+ std::vector<unsigned int> 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<int spacedim>
+ 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<unsigned int> (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<unsigned int> (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<unsigned int> max_level (n_vertices, 0);
+ std::vector<unsigned int> 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<int spacedim>
+ 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<unsigned int> (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<unsigned int> (tria.n_raw_lines () * fe.dofs_per_line, DoFHandler<3>::invalid_dof_index);
+ dof_handler.mg_faces->quads.dofs = std::vector<unsigned int> (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<unsigned int> max_level (n_vertices, 0);
+ std::vector<unsigned int> 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<int spacedim>
+ 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<int spacedim>
+ 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<int spacedim>
+ 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<int spacedim>
+ 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<int spacedim>
+ 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<int spacedim>
+ 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<int spacedim>
+ 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<int spacedim>
+ 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<int spacedim>
+ 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<int spacedim>
+ 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<int spacedim>
+ 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<int spacedim>
+ 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<int spacedim>
+ 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<int spacedim>
+ 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<int spacedim>
+ 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);
+ }
};
}
}
template <>
- types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const
-unsigned int DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
++types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
{
// check that only boundary
// indicators 0 and 1 are allowed
// in 1d
- for (std::set<types::boundary_id_t>::const_iterator i=boundary_indicators.begin();
+ for (std::set<types::boundary_id>::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;
}
template <>
- types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const
-unsigned int DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
++types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
{
// check that only boundary
// indicators 0 and 1 are allowed
// in 1d
- for (std::set<types::boundary_id_t>::const_iterator i=boundary_indicators.begin();
+ for (std::set<types::boundary_id>::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;
}
template<int dim, int spacedim>
-unsigned int
+types::global_dof_index
- DoFHandler<dim,spacedim>::n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const
+ DoFHandler<dim,spacedim>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
{
Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(),
- ExcInvalidBoundaryIndicator());
+ ExcInvalidBoundaryIndicator());
std::set<int> boundary_dofs;
DoFHandler<dim,spacedim>::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; i<levels.size(); ++i)
mem += MemoryConsumption::memory_consumption (*levels[i]);
-
+
+ for (unsigned int level = 0; level < mg_levels.size (); ++level)
+ mem += mg_levels[level]->memory_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;
}
}
-
+template<int dim, int spacedim>
+void DoFHandler<dim, spacedim>::distribute_mg_dofs (const FiniteElement<dim, spacedim>& fe, const types::global_dof_index offset) {
+ Assert ((dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim>*> (&*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<bool> user_flags;
-
++
+ (*tria).save_user_flags (user_flags);
+ const_cast<Triangulation<dim, spacedim>&>(*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<spacedim> (cell, next_free_dof);
-
++
+ mg_used_dofs[level] = next_free_dof;
+ }
-
++
+ const_cast<Triangulation<dim, spacedim>&>(*tria).load_user_flags (user_flags);
+ block_info_object.initialize (*this, true, true);
+}
+
+template<int dim, int spacedim>
+void DoFHandler<dim, spacedim>::reserve_space () {
+ internal::DoFHandler::Implementation::reserve_space_mg (*this);
+}
+
+template<int dim, int spacedim>
+void DoFHandler<dim, spacedim>::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<MGVertexDoFs> tmp;
++
+ std::swap (mg_vertex_dofs, tmp);
+}
+
template<int dim, int spacedim>
void DoFHandler<dim,spacedim>::initialize_local_block_info ()
std::swap (vertex_dofs, tmp);
}
-
+template<int dim, int spacedim>
+template<int structdim>
+types::global_dof_index DoFHandler<dim, spacedim>::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<structdim> ());
+}
+
+template<int dim, int spacedim>
+template<int structdim>
+void DoFHandler<dim, spacedim>::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<structdim> ());
+}
+
+template<int dim, int spacedim>
+DoFHandler<dim, spacedim>::MGVertexDoFs::MGVertexDoFs (): coarsest_level (numbers::invalid_unsigned_int), finest_level (0), indices (0), indices_offset (0) {
+}
+
+template<int dim, int spacedim>
+DoFHandler<dim, spacedim>::MGVertexDoFs::~MGVertexDoFs () {
+ delete[] indices;
+ delete[] indices_offset;
+}
+
+template<int dim, int spacedim>
+void DoFHandler<dim, spacedim>::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<dim, spacedim>::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<int dim, int spacedim>
+unsigned int DoFHandler<dim, spacedim>::MGVertexDoFs::get_coarsest_level () const {
+ return coarsest_level;
+}
+
+template<int dim, int spacedim>
+unsigned int DoFHandler<dim, spacedim>::MGVertexDoFs::get_finest_level () const {
+ return finest_level;
+}
+
/*-------------- Explicit Instantiations -------------------------------*/
#include "dof_handler.inst"
*/
template <int dim, int spacedim>
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<dim,spacedim> &dof_handler)
{
const dealii::Triangulation<dim,spacedim> & tria
// mark all own cells, that miss some
// dof_data and collect the neighbors
// that are going to send stuff to us
- std::set<dealii::types::subdomain_id_t> senders;
+ std::set<dealii::types::subdomain_id> senders;
{
- std::vector<unsigned int> local_dof_indices;
+ std::vector<types::global_dof_index> local_dof_indices;
typename DoFHandler<dim,spacedim>::active_cell_iterator
cell, endc = dof_handler.end();
+ template <class DH>
+ void
+ sort_selected_dofs_back (DH& dof_handler,
+ const std::vector<bool>& selected_dofs,
+ const unsigned int level)
+ {
+ std::vector<unsigned int> 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 <class DH>
void
- compute_sort_selected_dofs_back (std::vector<unsigned int>& new_indices,
+ compute_sort_selected_dofs_back (std::vector<types::global_dof_index>& new_indices,
const DH& dof_handler,
const std::vector<bool>& selected_dofs)
{
template <class DH, class SparsityPattern>
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();
template <class DH, class SparsityPattern>
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();
template <class DH, class SparsityPattern>
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();
template <class DH>
void
extract_dofs_with_support_on_boundary (const DH &dof_handler,
- const std::vector<bool> &component_select,
- std::vector<bool> &selected_dofs,
- const std::set<types::boundary_id> &boundary_indicators)
+ const std::vector<bool> &component_select,
+ std::vector<bool> &selected_dofs,
- const std::set<types::boundary_id_t> &boundary_indicators)
++ const std::set<types::boundary_id> &boundary_indicators)
{
AssertDimension (component_select.size(), n_components(dof_handler));
Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(),
template <class DH>
void
extract_subdomain_dofs (const DH &dof_handler,
- const types::subdomain_id_t subdomain_id,
- const types::subdomain_id subdomain_id,
- std::vector<bool> &selected_dofs)
++ const types::subdomain_id subdomain_id,
+ std::vector<bool> &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);
template <class DH>
void
get_subdomain_association (const DH &dof_handler,
- std::vector<types::subdomain_id_t> &subdomain_association)
- std::vector<types::subdomain_id> &subdomain_association)
++ std::vector<types::subdomain_id> &subdomain_association)
{
// if the Triangulation is distributed, the
// only thing we can usefully ask is for
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
template <class DH>
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<types::subdomain_id_t> subdomain_association (dof_handler.n_dofs());
+ std::vector<types::subdomain_id> subdomain_association (dof_handler.n_dofs());
get_subdomain_association (dof_handler, subdomain_association);
return std::count (subdomain_association.begin(),
template <class DH>
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
template <class DH>
void
count_dofs_with_subdomain_association (const DH &dof_handler,
- const types::subdomain_id_t subdomain,
- const types::subdomain_id subdomain,
- std::vector<unsigned int> &n_dofs_on_subdomain)
++ const types::subdomain_id subdomain,
+ std::vector<unsigned int> &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
#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<dim,spacedim> * tria
- = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
- (&dof_handler.get_tria())))
+ = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
+ (&dof_handler.get_tria())))
{
- std::vector<unsigned int> local_dof_count = dofs_per_component;
+ std::vector<unsigned int> 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
}
#endif
- return n_parameters_on_fine_grid;
+ return n_parameters_on_fine_grid;
}
-
-
+
+
}
}
template <class DH>
void map_dof_to_boundary_indices (
const DH &dof_handler,
- const std::set<types::boundary_id_t> &boundary_indicators,
+ const std::set<types::boundary_id> &boundary_indicators,
- std::vector<unsigned int> &mapping)
+ std::vector<types::global_dof_index> &mapping)
{
Assert (&dof_handler.get_fe() != 0, ExcNoFESelected());
Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(),
local_dof_indices.resize(cell->get_fe().dofs_per_cell);
cell->get_dof_indices(local_dof_indices);
-
+
const std::vector<Point<spacedim> > & 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];
unsigned int
DoFTools::count_dofs_with_subdomain_association<hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> >
(const hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
- const types::subdomain_id_t);
+ const types::subdomain_id);
-
+
template
IndexSet
DoFTools::dof_indices_with_subdomain_association<hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> >
(const hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
- const types::subdomain_id_t);
+ const types::subdomain_id);
-
+
template
void
DoFTools::get_subdomain_association<DoFHandler<deal_II_dimension,deal_II_dimension+1> >
template <>
- types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const
- unsigned int DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
++ types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
{
Assert (finite_elements != 0, ExcNoFESelected());
}
template <>
- types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id_t> &) const
- unsigned int DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id> &) const
++ types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id> &) const
{
Assert(false,ExcNotImplemented());
return 0;
}
template <>
- types::global_dof_index DoFHandler<1,3>::n_boundary_dofs (const std::set<types::boundary_id_t> &) const
- unsigned int DoFHandler<1,3>::n_boundary_dofs (const std::set<types::boundary_id> &) const
++ types::global_dof_index DoFHandler<1,3>::n_boundary_dofs (const std::set<types::boundary_id> &) const
{
Assert(false,ExcNotImplemented());
return 0;
template<int dim, int spacedim>
- unsigned int
+ types::global_dof_index
- DoFHandler<dim,spacedim>::n_boundary_dofs (const std::set<types::boundary_id_t> &boundary_indicators) const
+ DoFHandler<dim,spacedim>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_indicators) const
{
Assert (finite_elements != 0, ExcNoFESelected());
Assert (boundary_indicators.find (types::internal_face_boundary_id) == boundary_indicators.end(),
template <>
- types::global_dof_index DoFHandler<2,3>::n_boundary_dofs (const std::set<types::boundary_id_t> &) const
- unsigned int DoFHandler<2,3>::n_boundary_dofs (const std::set<types::boundary_id> &) const
++ types::global_dof_index DoFHandler<2,3>::n_boundary_dofs (const std::set<types::boundary_id> &) const
{
Assert(false,ExcNotImplemented());
return 0;
}
- DoFHandler<dim,spacedim>::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 <int dim, int spacedim>
+ template <int structdim>
+ types::global_dof_index
-
++ DoFHandler<dim,spacedim>::get_dof_index (const unsigned int,
++ const unsigned int,
++ const unsigned int,
++ const unsigned int) const
+ {
- DoFHandler<dim,spacedim>::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
++ Assert (false, ExcNotImplemented());
++ return numbers::invalid_unsigned_int;
+ }
+
+
+ template <int dim, int spacedim>
+ template <int structdim>
+ void
-
++ DoFHandler<dim,spacedim>::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<int dim, int spacedim>
void DoFHandler<dim,spacedim>::clear_space ()
{
{
const FiniteElement<dim>& fe = dof.get_fe();
const unsigned int n_comp = fe.n_components();
-
+
Assert ((level>=1) && (level<dof.get_tria().n_levels()),
- ExcIndexRange(level, 1, dof.get_tria().n_levels()));
+ ExcIndexRange(level, 1, dof.get_tria().n_levels()));
const unsigned int fine_dofs = dof.n_dofs(level);
const unsigned int coarse_dofs = dof.n_dofs(level-1);
if (component_mask.size() == 0)
{
typename MGDoFHandler<dim,spacedim>::cell_iterator
- cell = dof.begin(),
- endc = dof.end();
+ cell = dof.begin(),
+ endc = dof.end();
for (; cell!=endc; ++cell)
- {
- const FiniteElement<dim> &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<dim>::faces_per_cell;
- ++face_no)
- if (cell->at_boundary(face_no) == true)
- {
- const typename MGDoFHandler<dim,spacedim>::face_iterator
- face = cell->face(face_no);
- const types::boundary_id bi = face->boundary_indicator();
+ {
+ const FiniteElement<dim> &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<dim>::faces_per_cell;
+ ++face_no)
+ if (cell->at_boundary(face_no) == true)
+ {
+ const typename MGDoFHandler<dim,spacedim>::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;i<fe.dofs_per_face;++i)
- boundary_indices[level].insert(local_dofs[i]);
- }
- }
- }
+ if (function_map.find(bi) != function_map.end())
+ {
+ face->get_mg_dof_indices(level, local_dofs);
+ for (unsigned int i=0;i<fe.dofs_per_face;++i)
+ boundary_indices[level].insert(local_dofs[i]);
+ }
+ }
+ }
}
else
{
// that are non-zero for the
// components we are interested in,
// are in fact primitive
- for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
- {
- const std::vector<bool> &nonzero_component_array
- = cell->get_fe().get_nonzero_components (i);
- for (unsigned int c=0; c<n_components; ++c)
- if ((nonzero_component_array[c] == true)
- &&
- (component_mask[c] == true))
- Assert (cell->get_fe().is_primitive (i),
- ExcMessage ("This function can only deal with requested boundary "
- "values that correspond to primitive (scalar) base "
- "elements"));
- }
-
- typename MGDoFHandler<dim,spacedim>::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; i<cell->get_fe().dofs_per_cell; ++i)
+ {
+ const std::vector<bool> &nonzero_component_array
+ = cell->get_fe().get_nonzero_components (i);
+ for (unsigned int c=0; c<n_components; ++c)
+ if ((nonzero_component_array[c] == true)
+ &&
+ (component_mask[c] == true))
+ Assert (cell->get_fe().is_primitive (i),
+ ExcMessage ("This function can only deal with requested boundary "
+ "values that correspond to primitive (scalar) base "
+ "elements"));
+ }
+
+ typename MGDoFHandler<dim,spacedim>::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
std::fill (cell_dofs.begin(), cell_dofs.end(), false);
for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
- {
- const typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_nr);
- if (!face->at_boundary())
- {
+ {
+ const typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_nr);
+ if (!face->at_boundary())
+ {
//interior face
- const typename MGDoFHandler<dim>::cell_iterator
- neighbor = cell->neighbor(face_nr);
+ const typename MGDoFHandler<dim>::cell_iterator
+ neighbor = cell->neighbor(face_nr);
- if (neighbor->level() < cell->level())
- {
- for (unsigned int j=0; j<dofs_per_face; ++j)
- cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
+ if (neighbor->level() < cell->level())
+ {
+ for (unsigned int j=0; j<dofs_per_face; ++j)
+ cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
- }
- }
- }
+ }
+ }
+ }
-
+
const unsigned int level = cell->level();
cell->get_mg_dof_indices (local_dof_indices);
}
-template <int dim, int spacedim>
-void
-extract_non_interface_dofs (const MGDoFHandler<dim,spacedim> &mg_dof_handler,
- std::vector<std::set<unsigned int> > &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 <int dim, int spacedim>
+ void
+ extract_non_interface_dofs (const MGDoFHandler<dim,spacedim> &mg_dof_handler,
+ std::vector<std::set<unsigned int> > &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<dim,spacedim> &fe = mg_dof_handler.get_fe();
+ const FiniteElement<dim,spacedim> &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<unsigned int> local_dof_indices (dofs_per_cell);
- std::vector<bool> cell_dofs(dofs_per_cell, false);
- std::vector<bool> cell_dofs_interface(dofs_per_cell, false);
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+ std::vector<bool> cell_dofs(dofs_per_cell, false);
+ std::vector<bool> cell_dofs_interface(dofs_per_cell, false);
- typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
- endc = mg_dof_handler.end();
+ typename MGDoFHandler<dim>::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<GeometryInfo<dim>::faces_per_cell; ++face_nr)
- {
- const typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_nr);
- if (!face->at_boundary())
- {
+ for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
+ {
+ const typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_nr);
+ if (!face->at_boundary())
+ {
//interior face
- const typename MGDoFHandler<dim>::cell_iterator
- neighbor = cell->neighbor(face_nr);
-
- if ((neighbor->level() < cell->level()))
+ const typename MGDoFHandler<dim>::cell_iterator
+ neighbor = cell->neighbor(face_nr);
+
+ if ((neighbor->level() < cell->level()))
+ {
+ for (unsigned int j=0; j<dofs_per_face; ++j)
+ cell_dofs_interface[fe.face_to_cell_index(j,face_nr)] = true;
+ }
+ else
{
for (unsigned int j=0; j<dofs_per_face; ++j)
- cell_dofs_interface[fe.face_to_cell_index(j,face_nr)] = true;
+ cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
}
- else
+ }
+ else
{
- for (unsigned int j=0; j<dofs_per_face; ++j)
- cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
- }
- }
- else
- {
//boundary face
for (unsigned int j=0; j<dofs_per_face; ++j)
- cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
- }
- }
+ cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
+ }
+ }
- const unsigned int level = cell->level();
- 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<dofs_per_cell; ++i)
- if (cell_dofs[i] && !cell_dofs_interface[i])
- non_interface_dofs[level].insert(local_dof_indices[i]);
- }
-}
+ for(unsigned int i=0; i<dofs_per_cell; ++i)
+ if (cell_dofs[i] && !cell_dofs_interface[i])
+ non_interface_dofs[level].insert(local_dof_indices[i]);
+ }
+ }
template <int dim, int spacedim>
const Table<2,DoFTools::Coupling>&, const Table<2,DoFTools::Coupling>&);
#endif
- template void count_dofs_per_component<deal_II_dimension> (
+ template void count_dofs_per_component (
const MGDoFHandler<deal_II_dimension>&, std::vector<std::vector<unsigned int> >&,
bool, std::vector<unsigned int>);
- template void count_dofs_per_component<deal_II_dimension> (
+ template void count_dofs_per_component (
const MGDoFHandler<deal_II_dimension>&, std::vector<std::vector<unsigned int> >&,
std::vector<unsigned int>);
- template void count_dofs_per_block<deal_II_dimension> (
+ template void count_dofs_per_block (
const MGDoFHandler<deal_II_dimension>&, std::vector<std::vector<unsigned int> >&,
std::vector<unsigned int>);
++ template void count_dofs_per_block (
++ const DoFHandler<deal_II_dimension>&, std::vector<std::vector<unsigned int> >&,
++ std::vector<unsigned int>);
#if deal_II_dimension > 1
template void make_boundary_list(
#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<unsigned int> >&, std::vector<unsigned int>);
++
++ template void count_dofs_per_block (
++ const DoFHandler<1,3>&,
++ std::vector<std::vector<unsigned int> >&, std::vector<unsigned int>);
++ template void count_dofs_per_block (
++ const DoFHandler<2,3>&,
++ std::vector<std::vector<unsigned int> >&, std::vector<unsigned int>);
++ template void count_dofs_per_block (
++ const DoFHandler<1,2>&,
++ std::vector<std::vector<unsigned int> >&, std::vector<unsigned int>);
#endif
\}
}