// leaves out is faces in 3d --
// let's hope that this never is
// a problem
- Assert ((dim==structdim
- ?
- typename
- dealii::DoFHandler<dim,spacedim>::
- raw_cell_iterator (&dof_handler.get_tria(),
- obj_level,
- obj_index,
- &dof_handler)->used()
- :
- (structdim==1
- ?
- typename
- dealii::DoFHandler<dim,spacedim>::
- raw_line_iterator (&dof_handler.get_tria(),
- obj_level,
- obj_index,
- &dof_handler)->used()
- :
- true))
- == true,
- ExcMessage ("This cell is not active and therefore can't be "
- "queried for its active FE indices"));
- Assert ((dim==structdim
- ?
- typename
- dealii::internal::DoFHandler::
- Iterators<dealii::DoFHandler<dim,spacedim> >::
- raw_cell_iterator (&dof_handler.get_tria(),
- obj_level,
- obj_index,
- &dof_handler)->used()
- :
- (structdim==1
- ?
- typename
- dealii::internal::DoFHandler::
- Iterators<dealii::DoFHandler<dim,spacedim> >::
- raw_line_iterator (&dof_handler.get_tria(),
- obj_level,
- obj_index,
- &dof_handler)->used()
- :
- true))
- == true,
- ExcMessage ("This cell is not active and therefore can't be "
- "queried for its active FE indices"));
- return 1;
- }
-
-
-
- template <int structdim, int dim, int spacedim>
- static
- unsigned int
- nth_active_fe_index (const dealii::DoFHandler<dim,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int n,
- const dealii::internal::int2type<structdim> &)
- {
++ Assert ((dim==structdim
++ ?
++ typename
++ dealii::internal::DoFHandler::
++ Iterators<dealii::DoFHandler<dim,spacedim> >::
++ raw_cell_iterator (&dof_handler.get_tria(),
++ obj_level,
++ obj_index,
++ &dof_handler)->used()
++ :
++ (structdim==1
++ ?
++ typename
++ dealii::internal::DoFHandler::
++ Iterators<dealii::DoFHandler<dim,spacedim> >::
++ raw_line_iterator (&dof_handler.get_tria(),
++ obj_level,
++ obj_index,
++ &dof_handler)->used()
++ :
++ true))
++ == true,
++ ExcMessage ("This cell is not active and therefore can't be "
++ "queried for its active FE indices"));
+ return 1;
+ }
+
+
+
+ template <int structdim, int dim, int spacedim>
+ static
+ unsigned int
+ nth_active_fe_index (const dealii::DoFHandler<dim,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int n,
+ const dealii::internal::int2type<structdim> &)
+ {
// check that the object we look
// at is in fact active. the
// problem is that we have
// leaves out is faces in 3d --
// let's hope that this never is
// a problem
- Assert ((dim==structdim
- ?
- typename
- dealii::internal::DoFHandler::
- Iterators<dealii::DoFHandler<dim,spacedim> >::
- raw_cell_iterator (&dof_handler.get_tria(),
- obj_level,
- obj_index,
- &dof_handler)->used()
- :
- (structdim==1
- ?
- typename
- dealii::internal::DoFHandler::
- Iterators<dealii::DoFHandler<dim,spacedim> >::
- raw_line_iterator (&dof_handler.get_tria(),
- obj_level,
- obj_index,
- &dof_handler)->used()
- :
- true))
- == true,
- ExcMessage ("This cell is not active and therefore can't be "
- "queried for its active FE indices"));
- Assert (n == 0, ExcIndexRange (n, 0, 1));
-
- return dealii::DoFHandler<dim,spacedim>::default_fe_index;
- }
-
-
- template <int spacedim>
- static
- bool
- fe_index_is_active (const dealii::hp::DoFHandler<1,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int fe_index,
- const dealii::internal::int2type<1> &)
- {
- return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
- obj_index,
- fe_index,
- obj_level);
- }
-
-
- template <int spacedim>
- static
- unsigned int
- n_active_fe_indices (const dealii::hp::DoFHandler<1,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const dealii::internal::int2type<1> &)
- {
- return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
- obj_index);
- }
-
-
-
- template <int spacedim>
- static
- unsigned int
- nth_active_fe_index (const dealii::hp::DoFHandler<1,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int n,
- const dealii::internal::int2type<1> &)
- {
- return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
- obj_level,
- obj_index,
- n);
- }
-
-
- template <int spacedim>
- static
- bool
- fe_index_is_active (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int fe_index,
- const dealii::internal::int2type<1> &)
- {
- return dof_handler.faces->lines.fe_index_is_active(dof_handler,
- obj_index,
- fe_index,
- obj_level);
- }
-
-
- template <int spacedim>
- static
- unsigned int
- n_active_fe_indices (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
- const unsigned int ,
- const unsigned int obj_index,
- const dealii::internal::int2type<1> &)
- {
- return dof_handler.faces->lines.n_active_fe_indices (dof_handler,
- obj_index);
- }
-
-
- template <int spacedim>
- static
- unsigned int
- nth_active_fe_index (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int n,
- const dealii::internal::int2type<1> &)
- {
- return dof_handler.faces->lines.nth_active_fe_index (dof_handler,
- obj_level,
- obj_index,
- n);
- }
-
-
-
- template <int spacedim>
- static
- bool
- fe_index_is_active (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int fe_index,
- const dealii::internal::int2type<2> &)
- {
- return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
- obj_index,
- fe_index,
- obj_level);
- }
-
-
- template <int spacedim>
- static
- unsigned int
- n_active_fe_indices (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const dealii::internal::int2type<2> &)
- {
- return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
- obj_index);
- }
-
-
-
- template <int spacedim>
- static
- unsigned int
- nth_active_fe_index (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int n,
- const dealii::internal::int2type<2> &)
- {
- return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
- obj_level,
- obj_index,
- n);
- }
-
-
-
- template <int spacedim>
- static
- bool
- fe_index_is_active (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int fe_index,
- const dealii::internal::int2type<1> &)
- {
- return dof_handler.faces->lines.fe_index_is_active(dof_handler,
- obj_index,
- fe_index,
- obj_level);
- }
-
-
- template <int spacedim>
- static
- unsigned int
- n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
- const unsigned int,
- const unsigned int obj_index,
- const dealii::internal::int2type<1> &)
- {
- return dof_handler.faces->lines.n_active_fe_indices (dof_handler,
- obj_index);
- }
-
-
-
- template <int spacedim>
- static
- unsigned int
- nth_active_fe_index (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int n,
- const dealii::internal::int2type<1> &)
- {
- return dof_handler.faces->lines.nth_active_fe_index (dof_handler,
- obj_level,
- obj_index,
- n);
- }
-
-
-
- template <int spacedim>
- static
- bool
- fe_index_is_active (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int fe_index,
- const dealii::internal::int2type<2> &)
- {
- return dof_handler.faces->quads.fe_index_is_active(dof_handler,
- obj_index,
- fe_index,
- obj_level);
- }
-
- template <int spacedim>
- static
- bool
- fe_index_is_active (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int fe_index,
- const dealii::internal::int2type<3> &)
- {
- return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
- obj_index,
- fe_index,
- obj_level);
- }
-
-
- template <int spacedim>
- static
- unsigned int
- n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
- const unsigned int ,
- const unsigned int obj_index,
- const dealii::internal::int2type<2> &)
- {
- return dof_handler.faces->quads.n_active_fe_indices (dof_handler,
- obj_index);
- }
-
-
-
- template <int spacedim>
- static
- unsigned int
- nth_active_fe_index (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int n,
- const dealii::internal::int2type<2> &)
- {
- return dof_handler.faces->quads.nth_active_fe_index (dof_handler,
- obj_level,
- obj_index,
- n);
- }
-
-
-
- template <int spacedim>
- static
- unsigned int
- n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const dealii::internal::int2type<3> &)
- {
- return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
- obj_index);
- }
-
-
-
- template <int spacedim>
- static
- unsigned int
- nth_active_fe_index (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int n,
- const dealii::internal::int2type<3> &)
- {
- return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
- obj_level,
- obj_index,
- n);
- }
+ Assert ((dim==structdim
- ?
- typename
- dealii::DoFHandler<dim,spacedim>::
- raw_cell_iterator (&dof_handler.get_tria(),
- obj_level,
- obj_index,
- &dof_handler)->used()
- :
- (structdim==1
- ?
- typename
- dealii::DoFHandler<dim,spacedim>::
- raw_line_iterator (&dof_handler.get_tria(),
- obj_level,
- obj_index,
- &dof_handler)->used()
- :
- true))
- == true,
- ExcMessage ("This cell is not active and therefore can't be "
- "queried for its active FE indices"));
++ ?
++ typename
++ dealii::internal::DoFHandler::
++ Iterators<dealii::DoFHandler<dim,spacedim> >::
++ raw_cell_iterator (&dof_handler.get_tria(),
++ obj_level,
++ obj_index,
++ &dof_handler)->used()
++ :
++ (structdim==1
++ ?
++ typename
++ dealii::internal::DoFHandler::
++ Iterators<dealii::DoFHandler<dim,spacedim> >::
++ raw_line_iterator (&dof_handler.get_tria(),
++ obj_level,
++ obj_index,
++ &dof_handler)->used()
++ :
++ true))
++ == true,
++ ExcMessage ("This cell is not active and therefore can't be "
++ "queried for its active FE indices"));
+ Assert (n == 0, ExcIndexRange (n, 0, 1));
+
+ return dealii::DoFHandler<dim,spacedim>::default_fe_index;
+ }
+
+
+ template <int spacedim>
+ static
+ bool
+ fe_index_is_active (const dealii::hp::DoFHandler<1,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int fe_index,
+ const dealii::internal::int2type<1> &)
+ {
+ return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
+ obj_index,
+ fe_index,
+ obj_level);
+ }
+
+
+ template <int spacedim>
+ static
+ unsigned int
+ n_active_fe_indices (const dealii::hp::DoFHandler<1,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const dealii::internal::int2type<1> &)
+ {
+ return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
+ obj_index);
+ }
+
+
+
+ template <int spacedim>
+ static
+ unsigned int
+ nth_active_fe_index (const dealii::hp::DoFHandler<1,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int n,
+ const dealii::internal::int2type<1> &)
+ {
+ return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
+ obj_level,
+ obj_index,
+ n);
+ }
+
+
+ template <int spacedim>
+ static
+ bool
+ fe_index_is_active (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int fe_index,
+ const dealii::internal::int2type<1> &)
+ {
+ return dof_handler.faces->lines.fe_index_is_active(dof_handler,
+ obj_index,
+ fe_index,
+ obj_level);
+ }
+
+
+ template <int spacedim>
+ static
+ unsigned int
+ n_active_fe_indices (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
+ const unsigned int ,
+ const unsigned int obj_index,
+ const dealii::internal::int2type<1> &)
+ {
+ return dof_handler.faces->lines.n_active_fe_indices (dof_handler,
+ obj_index);
+ }
+
+
+ template <int spacedim>
+ static
+ unsigned int
+ nth_active_fe_index (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int n,
+ const dealii::internal::int2type<1> &)
+ {
+ return dof_handler.faces->lines.nth_active_fe_index (dof_handler,
+ obj_level,
+ obj_index,
+ n);
+ }
+
+
+
+ template <int spacedim>
+ static
+ bool
+ fe_index_is_active (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int fe_index,
+ const dealii::internal::int2type<2> &)
+ {
+ return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
+ obj_index,
+ fe_index,
+ obj_level);
+ }
+
+
+ template <int spacedim>
+ static
+ unsigned int
+ n_active_fe_indices (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const dealii::internal::int2type<2> &)
+ {
+ return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
+ obj_index);
+ }
+
+
+
+ template <int spacedim>
+ static
+ unsigned int
+ nth_active_fe_index (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int n,
+ const dealii::internal::int2type<2> &)
+ {
+ return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
+ obj_level,
+ obj_index,
+ n);
+ }
+
+
+
+ template <int spacedim>
+ static
+ bool
+ fe_index_is_active (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int fe_index,
+ const dealii::internal::int2type<1> &)
+ {
+ return dof_handler.faces->lines.fe_index_is_active(dof_handler,
+ obj_index,
+ fe_index,
+ obj_level);
+ }
+
+
+ template <int spacedim>
+ static
+ unsigned int
+ n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
+ const unsigned int,
+ const unsigned int obj_index,
+ const dealii::internal::int2type<1> &)
+ {
+ return dof_handler.faces->lines.n_active_fe_indices (dof_handler,
+ obj_index);
+ }
+
+
+
+ template <int spacedim>
+ static
+ unsigned int
+ nth_active_fe_index (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int n,
+ const dealii::internal::int2type<1> &)
+ {
+ return dof_handler.faces->lines.nth_active_fe_index (dof_handler,
+ obj_level,
+ obj_index,
+ n);
+ }
+
+
+
+ template <int spacedim>
+ static
+ bool
+ fe_index_is_active (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int fe_index,
+ const dealii::internal::int2type<2> &)
+ {
+ return dof_handler.faces->quads.fe_index_is_active(dof_handler,
+ obj_index,
+ fe_index,
+ obj_level);
+ }
+
+ template <int spacedim>
+ static
+ bool
+ fe_index_is_active (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int fe_index,
+ const dealii::internal::int2type<3> &)
+ {
+ return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
+ obj_index,
+ fe_index,
+ obj_level);
+ }
+
+
+ template <int spacedim>
+ static
+ unsigned int
+ n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
+ const unsigned int ,
+ const unsigned int obj_index,
+ const dealii::internal::int2type<2> &)
+ {
+ return dof_handler.faces->quads.n_active_fe_indices (dof_handler,
+ obj_index);
+ }
+
+
+
+ template <int spacedim>
+ static
+ unsigned int
+ nth_active_fe_index (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int n,
+ const dealii::internal::int2type<2> &)
+ {
+ return dof_handler.faces->quads.nth_active_fe_index (dof_handler,
+ obj_level,
+ obj_index,
+ n);
+ }
+
+
+
+ template <int spacedim>
+ static
+ unsigned int
+ n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const dealii::internal::int2type<3> &)
+ {
+ return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
+ obj_index);
+ }
+
+
+
+ template <int spacedim>
+ static
+ unsigned int
+ nth_active_fe_index (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int n,
+ const dealii::internal::int2type<3> &)
+ {
+ return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
+ obj_level,
+ obj_index,
+ n);
+ }
/**
* Set the @p local_index-th
// applies, if the face_rotation or
// face_orientation is non-standard
for (unsigned int quad=0; quad<6; ++quad)
- for (unsigned int d=0; d<dofs_per_quad; ++d)
- *next++ = accessor.quad(quad)->dof_index(accessor.get_fe(fe_index).
- adjust_quad_dof_index_for_face_orientation(d,
- accessor.face_orientation(quad),
- accessor.face_flip(quad),
- accessor.face_rotation(quad)),
- fe_index);
+ for (unsigned int d=0; d<dofs_per_quad; ++d)
+ *next++ = accessor.quad(quad)->dof_index(accessor.get_fe(fe_index).
+ adjust_quad_dof_index_for_face_orientation(d,
+ accessor.face_orientation(quad),
+ accessor.face_flip(quad),
+ accessor.face_rotation(quad)),
+ fe_index);
for (unsigned int d=0; d<dofs_per_hex; ++d)
- *next++ = accessor.dof_index(d,fe_index);
+ *next++ = accessor.dof_index(d,fe_index);
+ }
-
++
+ template<class DH>
+ void get_mg_dof_indices (const dealii::DoFAccessor<1, DH>& accessor, const int level, std::vector<types::global_dof_index>& dof_indices, const unsigned int fe_index) {
+ const FiniteElement<DH::dimension, DH::space_dimension>& fe = accessor.get_dof_handler ().get_fe ()[fe_index];
+ std::vector<types::global_dof_index>::iterator next = dof_indices.begin ();
-
++
+ for (unsigned int vertex = 0; vertex < GeometryInfo<1>::vertices_per_cell; ++vertex)
+ for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
+ *next++ = accessor.mg_vertex_dof_index (level, vertex, dof);
-
++
+ for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof)
+ *next++ = accessor.mg_dof_index (level, dof);
-
++
+ Assert (next == dof_indices.end (), ExcInternalError ());
+ }
-
++
+ template<class DH>
+ void get_mg_dof_indices (const dealii::DoFAccessor<2, DH>& accessor, const int level, std::vector<types::global_dof_index>& dof_indices, const unsigned int fe_index) {
+ const FiniteElement<DH::dimension, DH::space_dimension>& fe = accessor.get_dof_handler ().get_fe ()[fe_index];
+ std::vector<types::global_dof_index>::iterator next = dof_indices.begin ();
-
++
+ for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_cell; ++vertex)
+ for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
+ *next++ = accessor.mg_vertex_dof_index (level, vertex, dof);
-
++
+ for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell; ++line)
+ for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof)
+ *next++ = accessor.line (line)->mg_dof_index (level, dof);
-
++
+ for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof)
+ *next++ = accessor.mg_dof_index (level, dof);
-
++
+ Assert (next == dof_indices.end (), ExcInternalError ());
+ }
-
++
+ template<class DH>
+ void get_mg_dof_indices (const dealii::DoFAccessor<3, DH>& accessor, const int level, std::vector<types::global_dof_index>& dof_indices, const unsigned int fe_index) {
+ const FiniteElement<DH::dimension, DH::space_dimension>& fe = accessor.get_dof_handler ().get_fe ()[fe_index];
+ std::vector<types::global_dof_index>::iterator next = dof_indices.begin ();
-
++
+ for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_cell; ++vertex)
+ for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
+ *next++ = accessor.mg_vertex_dof_index (level, vertex, dof);
-
++
+ for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell; ++line)
+ for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof)
+ *next++ = accessor.line (line)->mg_dof_index (level, dof);
-
++
+ for (unsigned int quad = 0; quad < GeometryInfo<3>::quads_per_cell; ++quad)
+ for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof)
+ *next++ = accessor.quad (quad)->mg_dof_index (level, dof);
-
++
+ for (unsigned int dof = 0; dof < fe.dofs_per_hex; ++dof)
+ *next++ = accessor.mg_dof_index (level, dof);
-
++
+ Assert (next == dof_indices.end (), ExcInternalError ());
}
}
}
dealii::internal::DoFAccessor::get_dof_indices (*this, dof_indices, fe_index);
}
-
+template<int structdim, class DH>
+inline
+void DoFAccessor<structdim, DH>::get_mg_dof_indices (const int level, std::vector<types::global_dof_index>& dof_indices, const unsigned int fe_index) const {
+ Assert (this->dof_handler != 0, ExcInvalidObject ());
+ Assert (&this->dof_handler->get_fe () != 0, ExcInvalidObject ());
-
++
+ switch (structdim) {
+ case 1: {
+ Assert (dof_indices.size () == 2 * this->dof_handler->get_fe ()[fe_index].dofs_per_vertex + this->dof_handler->get_fe ()[fe_index].dofs_per_line, ExcVectorDoesNotMatch ());
+ break;
+ }
-
++
+ case 2: {
+ Assert (dof_indices.size () == 4 * (this->dof_handler->get_fe ()[fe_index].dofs_per_vertex + this->dof_handler->get_fe ()[fe_index].dofs_per_line) + this->dof_handler->get_fe ()[fe_index].dofs_per_quad, ExcVectorDoesNotMatch ());
+ break;
+ }
-
+
+ case 3: {
+ Assert (dof_indices.size () == 8 * this->dof_handler->get_fe ()[fe_index].dofs_per_vertex + 12 * this->dof_handler->get_fe ()[fe_index].dofs_per_line + 6 * this->dof_handler->get_fe ()[fe_index].dofs_per_quad + this->dof_handler->get_fe ()[fe_index].dofs_per_hex, ExcVectorDoesNotMatch ());
+ break;
+ }
-
+
+ default:
+ Assert (false, ExcNotImplemented ());
+ }
+
+ internal::DoFAccessor::get_mg_dof_indices (*this, level, dof_indices, fe_index);
+}
template <int structdim, class DH>
inline
* constraint line <i>i</i>, according
* to the discussion in the general
* class description.
- * Note: the line needs to be added with
+ *
+ * @note the line needs to be added with
* one of the add_line() calls first.
*/
- void set_inhomogeneity (const unsigned int line,
+ void set_inhomogeneity (const types::global_dof_index line,
const double value);
/**
// that is actually wholly on
// the boundary, not only by
// one line or one vertex
- typename dealii::DoFHandler<dim,spacedim>::active_face_iterator
- face = dof_handler.begin_active_face(),
- endf = dof_handler.end_face();
+ typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
- std::vector<unsigned int> face_dof_indices (fe.dofs_per_face);
+ std::vector<types::global_dof_index> face_dof_indices (fe.dofs_per_face);
- for (; face!=endf; ++face)
- if (face->at_boundary())
- {
- face->get_dof_indices (face_dof_indices);
- for (unsigned int i=0; i<fe.dofs_per_face; ++i)
- // enter zero boundary values
- // for all boundary nodes
- //
- // we need not care about
- // vector valued elements here,
- // since we set all components
- boundary_values[face_dof_indices[i]] = 0.;
- }
+ for (; cell!=endc; ++cell)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell->at_boundary(f))
+ {
+ cell->face(f)->get_dof_indices (face_dof_indices);
+ for (unsigned int i=0; i<fe.dofs_per_face; ++i)
+ // enter zero boundary values
+ // for all boundary nodes
+ //
+ // we need not care about
+ // vector valued elements here,
+ // since we set all components
+ boundary_values[face_dof_indices[i]] = 0.;
+ }
}
}
interpolate_boundary_values (const DH &dof,
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 ComponentMask &component_mask)
{
- Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
dof, boundary_component,
boundary_function, boundary_values, component_mask);
void
interpolate_boundary_values (const DH &dof,
const typename FunctionMap<DH::space_dimension>::type &function_map,
- std::map<unsigned int,double> &boundary_values,
+ std::map<types::global_dof_index,double> &boundary_values,
const ComponentMask &component_mask)
{
- Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
dof, function_map,
boundary_values, component_mask);
project_boundary_values (const DoFHandler<dim,spacedim> &dof,
const typename FunctionMap<spacedim>::type &boundary_functions,
const Quadrature<dim-1> &q,
- std::map<unsigned int,double> &boundary_values,
+ std::map<types::global_dof_index,double> &boundary_values,
std::vector<unsigned int> component_mapping)
{
- Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
project_boundary_values(StaticMappingQ1<dim,spacedim>::mapping, dof, boundary_functions, q,
boundary_values, component_mapping);
}
* 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()->dof_object.dofs
- .resize (dof_handler.tria->n_raw_cells(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_cells(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>);
+ 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()->dof_object.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);
- }
++ dof_handler.levels
++ .push_back (new internal::DoFHandler::DoFLevel<1>);
++
++ dof_handler.levels.back()->dof_object.dofs
++ .resize (dof_handler.tria->n_raw_cells(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_cells(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()->dof_object.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()->dof_object.dofs
+ .resize (dof_handler.tria->n_raw_cells(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.levels.back()->cell_dof_indices_cache
+ .resize (dof_handler.tria->n_raw_cells(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);
- }
+ 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);
+ 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>);
+ 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()->dof_object.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()->dof_object.dofs
+ .resize (dof_handler.tria->n_raw_cells(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.levels.back()->cell_dof_indices_cache
+ .resize (dof_handler.tria->n_raw_cells(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);
- }
+ 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 ()->dof_object.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 ()->dof_object.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 ()->dof_object.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.dof_object.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.dof_object.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.dof_object.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.dof_object.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.dof_object.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.dof_object.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
+ }
};
}
}
std::set<int> boundary_dofs;
const unsigned int dofs_per_face = get_fe().dofs_per_face;
- std::vector<unsigned int> dofs_on_face(dofs_per_face);
+ std::vector<types::global_dof_index> dofs_on_face(dofs_per_face);
- // loop over all faces to check
- // whether they are at a
- // boundary. note that we need not
- // take special care of single
- // lines (using
+ // loop over all faces of all cells
+ // and see whether they are at a
+ // boundary. note (i) that we visit
+ // interior faces twice (which we
+ // don't care about) but exterior
+ // faces only once as is
+ // appropriate, and (ii) that we
+ // need not take special care of
+ // single lines (using
// @p{cell->has_boundary_lines}),
// since we do not support
// boundaries of dimension dim-2,
std::set<int> boundary_dofs;
const unsigned int dofs_per_face = get_fe().dofs_per_face;
- std::vector<unsigned int> dofs_on_face(dofs_per_face);
+ std::vector<types::global_dof_index> dofs_on_face(dofs_per_face);
- active_face_iterator face = begin_active_face (),
- endf = end_face();
- for (; face!=endf; ++face)
- if (boundary_indicators.find(face->boundary_indicator()) !=
- boundary_indicators.end())
- {
- face->get_dof_indices (dofs_on_face);
- for (unsigned int i=0; i<dofs_per_face; ++i)
- boundary_dofs.insert(dofs_on_face[i]);
- }
+
+ // same as in the previous
+ // function, but with an additional
+ // check for the boundary indicator
+ active_cell_iterator cell = begin_active (),
+ endc = end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell->at_boundary(f)
+ &&
+ (boundary_indicators.find(cell->face(f)->boundary_indicator()) !=
+ boundary_indicators.end()))
+ {
+ cell->face(f)->get_dof_indices (dofs_on_face);
+ for (unsigned int i=0; i<dofs_per_face; ++i)
+ boundary_dofs.insert(dofs_on_face[i]);
+ }
+
return boundary_dofs.size();
}
std::set<int> boundary_dofs;
const unsigned int dofs_per_face = get_fe().dofs_per_face;
- std::vector<unsigned int> dofs_on_face(dofs_per_face);
+ std::vector<types::global_dof_index> dofs_on_face(dofs_per_face);
- active_face_iterator face = begin_active_face (),
- endf = end_face();
- for (; face!=endf; ++face)
- if (std::find (boundary_indicators.begin(),
- boundary_indicators.end(),
- face->boundary_indicator()) !=
- boundary_indicators.end())
- {
- face->get_dof_indices (dofs_on_face);
- for (unsigned int i=0; i<dofs_per_face; ++i)
- boundary_dofs.insert(dofs_on_face[i]);
- }
+
+ // same as in the previous
+ // function, but with a different
+ // check for the boundary indicator
+ active_cell_iterator cell = begin_active (),
+ endc = end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell->at_boundary(f)
+ &&
+ (std::find (boundary_indicators.begin(),
+ boundary_indicators.end(),
+ cell->face(f)->boundary_indicator()) !=
+ boundary_indicators.end()))
+ {
+ cell->face(f)->get_dof_indices (dofs_on_face);
+ for (unsigned int i=0; i<dofs_per_face; ++i)
+ boundary_dofs.insert(dofs_on_face[i]);
+ }
+
return boundary_dofs.size();
}
delete faces;
faces = 0;
- std::vector<unsigned int> tmp;
+ std::vector<types::global_dof_index> tmp;
std::swap (vertex_dofs, tmp);
+
+ number_cache.clear ();
}
+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"
#include <deal.II/base/thread_management.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/types.h>
+ #include <deal.II/base/template_constraints.h>
#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/lac/sparsity_tools.h>
namespace boost
{
- #ifndef DEAL_II_BOOST_GRAPH_COMPILER_BUG
- namespace types
+ namespace boosttypes
{
using namespace ::boost;
using namespace std;
const bool reversed_numbering,
const bool use_constraints)
{
- #ifdef DEAL_II_BOOST_GRAPH_COMPILER_BUG
- (void)new_dof_indices;
- (void)dof_handler;
- (void)reversed_numbering;
- (void)use_constraints;
- Assert (false,
- ExcMessage ("Due to a bug in your compiler, this function triggers an internal "
- "compiler error and has been disabled. If you need to use the "
- "function, you need to upgrade to a newer version of the compiler."));
- #else
- types::Graph
+ boosttypes::Graph
graph(dof_handler.n_dofs());
- types::property_map<types::Graph,types::vertex_degree_t>::type
+ boosttypes::property_map<boosttypes::Graph,boosttypes::vertex_degree_t>::type
graph_degree;
internal::create_graph (dof_handler, use_constraints, graph, graph_degree);
const bool reversed_numbering,
const bool use_constraints)
{
- #ifdef DEAL_II_BOOST_GRAPH_COMPILER_BUG
- (void)new_dof_indices;
- (void)dof_handler;
- (void)reversed_numbering;
- (void)use_constraints;
- Assert (false,
- ExcMessage ("Due to a bug in your compiler, this function triggers an internal "
- "compiler error and has been disabled. If you need to use the "
- "function, you need to upgrade to a newer version of the compiler."));
- #else
- types::Graph
+ boosttypes::Graph
graph(dof_handler.n_dofs());
- types::property_map<types::Graph,types::vertex_degree_t>::type
+ boosttypes::property_map<boosttypes::Graph,boosttypes::vertex_degree_t>::type
graph_degree;
internal::create_graph (dof_handler, use_constraints, graph, graph_degree);
const BlockMask&,
std::vector<bool>&);
++template
++void
++DoFTools::extract_level_dofs<DoFHandler<deal_II_dimension, deal_II_dimension+1> >
++(const unsigned int level,
++ const DoFHandler<deal_II_dimension, deal_II_dimension+1>&,
++ const BlockMask&,
++ std::vector<bool>&);
++
#endif
const ComponentMask &,
std::vector<bool>&);
-DoFTools::extract_level_dofs<1,3>
+ template
+ void
++DoFTools::extract_level_dofs<DoFHandler<1,3> >
++(const unsigned int level,
++ const DoFHandler<1,3>&,
++ const BlockMask&,
++ std::vector<bool>&);
++
+template
+void
+DoFTools::extract_level_dofs<MGDoFHandler<1, 3> >
(const unsigned int level,
const MGDoFHandler<1,3>&,
const BlockMask &,
for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; v++)
if (cell->vertex_index(v) == vertex)
{
-- // OK, we found a cell that contains
-- // the given vertex. We add it
-- // to the list.
++ // OK, we found a cell that contains
++ // the given vertex. We add it
++ // to the list.
adjacent_cells.insert(cell);
-- // as explained above, in 2+d we need to check whether
-- // this vertex is on a face behind which there is a
-- // (possibly) coarser neighbor. if this is the case,
-- // then we need to also add this neighbor
++ // as explained above, in 2+d we need to check whether
++ // this vertex is on a face behind which there is a
++ // (possibly) coarser neighbor. if this is the case,
++ // then we need to also add this neighbor
if (dim >= 2)
for (unsigned int vface = 0; vface < dim; vface++)
{
const unsigned int face =
-- GeometryInfo<dim>::vertex_to_face[v][vface];
++ GeometryInfo<dim>::vertex_to_face[v][vface];
if (!cell->at_boundary(face)
&&
cell->neighbor(face)->active())
{
-- // there is a (possibly) coarser cell behind a
-- // face to which the vertex belongs. the
-- // vertex we are looking at is then either a
-- // vertex of that coarser neighbor, or it is a
-- // hanging node on one of the faces of that
-- // cell. in either case, it is adjacent to the
-- // vertex, so add it to the list as well (if
-- // the cell was already in the list then the
-- // std::set makes sure that we get it only
-- // once)
++ // there is a (possibly) coarser cell behind a
++ // face to which the vertex belongs. the
++ // vertex we are looking at is then either a
++ // vertex of that coarser neighbor, or it is a
++ // hanging node on one of the faces of that
++ // cell. in either case, it is adjacent to the
++ // vertex, so add it to the list as well (if
++ // the cell was already in the list then the
++ // std::set makes sure that we get it only
++ // once)
adjacent_cells.insert (cell->neighbor(face));
}
}
-- // in any case, we have found a cell, so go to the next cell
++ // in any case, we have found a cell, so go to the next cell
goto next_cell;
}
-- // in 3d also loop over the edges
++ // in 3d also loop over the edges
if (dim >= 3)
{
for (unsigned int e=0; e<GeometryInfo<dim>::lines_per_cell; ++e)
if (cell->line(e)->has_children())
-- // the only place where this vertex could have been
-- // hiding is on the mid-edge point of the edge we
-- // are looking at
++ // the only place where this vertex could have been
++ // hiding is on the mid-edge point of the edge we
++ // are looking at
if (cell->line(e)->child(0)->vertex_index(1) == vertex)
{
adjacent_cells.insert(cell);
-- // jump out of this tangle of nested loops
++ // jump out of this tangle of nested loops
goto next_cell;
}
}
-- // in more than 3d we would probably have to do the same as
-- // above also for even lower-dimensional objects
++ // in more than 3d we would probably have to do the same as
++ // above also for even lower-dimensional objects
Assert (dim <= 3, ExcNotImplemented());
-- // move on to the next cell if we have found the
-- // vertex on the current one
++ // move on to the next cell if we have found the
++ // vertex on the current one
next_cell:
;
}
-- // if this was an active vertex then there needs to have been
-- // at least one cell to which it is adjacent!
++ // if this was an active vertex then there needs to have been
++ // at least one cell to which it is adjacent!
Assert (adjacent_cells.size() > 0, ExcInternalError());
-- // return the result as a vector, rather than the set we built above
++ // return the result as a vector, rather than the set we built above
return
-- std::vector<typename Container<dim,spacedim>::active_cell_iterator>
-- (adjacent_cells.begin(), adjacent_cells.end());
++ std::vector<typename Container<dim,spacedim>::active_cell_iterator>
++ (adjacent_cells.begin(), adjacent_cells.end());
}
namespace
{
-- template <int dim, template<int, int> class Container, int spacedim>
-- void find_active_cell_around_point_internal(const Container<dim,spacedim>& container,
-- std::set<typename Container<dim,spacedim>::active_cell_iterator>& searched_cells,
-- std::set<typename Container<dim,spacedim>::active_cell_iterator>& adjacent_cells)
++ template <int dim, template<int, int> class Container, int spacedim>
++ void find_active_cell_around_point_internal(const Container<dim,spacedim>& container,
++ std::set<typename Container<dim,spacedim>::active_cell_iterator>& searched_cells,
++ std::set<typename Container<dim,spacedim>::active_cell_iterator>& adjacent_cells)
++ {
++ typedef typename Container<dim,spacedim>::active_cell_iterator cell_iterator;
++
++ // update the searched cells
++ searched_cells.insert(adjacent_cells.begin(), adjacent_cells.end());
++ // now we to collect all neighbors
++ // of the cells in adjacent_cells we
++ // have not yet searched.
++ std::set<cell_iterator> adjacent_cells_new;
++
++ typename std::set<cell_iterator>::const_iterator
++ cell = adjacent_cells.begin(),
++ endc = adjacent_cells.end();
++ for(; cell != endc; ++cell)
{
-- typedef typename Container<dim,spacedim>::active_cell_iterator cell_iterator;
--
-- // update the searched cells
-- searched_cells.insert(adjacent_cells.begin(), adjacent_cells.end());
-- // now we to collect all neighbors
-- // of the cells in adjacent_cells we
-- // have not yet searched.
-- std::set<cell_iterator> adjacent_cells_new;
--
-- typename std::set<cell_iterator>::const_iterator
-- cell = adjacent_cells.begin(),
-- endc = adjacent_cells.end();
-- for(; cell != endc; ++cell)
-- {
-- std::vector<cell_iterator> active_neighbors;
-- get_active_neighbors<Container<dim, spacedim> >(*cell, active_neighbors);
-- for (unsigned int i=0; i<active_neighbors.size(); ++i)
-- if(searched_cells.find(active_neighbors[i]) == searched_cells.end())
-- adjacent_cells_new.insert(active_neighbors[i]);
-- }
-- adjacent_cells.clear();
-- adjacent_cells.insert(adjacent_cells_new.begin(), adjacent_cells_new.end());
-- if (adjacent_cells.size() == 0)
-- {
-- // we haven't found any other cell that would be a
-- // neighbor of a previously found cell, but we know
-- // that we haven't checked all cells yet. that means
-- // that the domain is disconnected. in that case,
-- // choose the first previously untouched cell we
-- // can find
-- cell_iterator it = container.begin_active();
-- for ( ; it!=container.end();++it)
-- if(searched_cells.find(it) == searched_cells.end())
-- {
-- adjacent_cells.insert(it);
-- break;
-- }
-- }
++ std::vector<cell_iterator> active_neighbors;
++ get_active_neighbors<Container<dim, spacedim> >(*cell, active_neighbors);
++ for (unsigned int i=0; i<active_neighbors.size(); ++i)
++ if(searched_cells.find(active_neighbors[i]) == searched_cells.end())
++ adjacent_cells_new.insert(active_neighbors[i]);
++ }
++ adjacent_cells.clear();
++ adjacent_cells.insert(adjacent_cells_new.begin(), adjacent_cells_new.end());
++ if (adjacent_cells.size() == 0)
++ {
++ // we haven't found any other cell that would be a
++ // neighbor of a previously found cell, but we know
++ // that we haven't checked all cells yet. that means
++ // that the domain is disconnected. in that case,
++ // choose the first previously untouched cell we
++ // can find
++ cell_iterator it = container.begin_active();
++ for ( ; it!=container.end();++it)
++ if(searched_cells.find(it) == searched_cells.end())
++ {
++ adjacent_cells.insert(it);
++ break;
++ }
}
++ }
}
template <int dim, template<int, int> class Container, int spacedim>
std::vector<cell_iterator> adjacent_cells_tmp =
find_cells_adjacent_to_vertex(container, vertex);
-- // Make sure that we have found
-- // at least one cell adjacent to vertex.
-- Assert(adjacent_cells_tmp.size()>0, ExcInternalError());
--
-- // Copy all the cells into a std::set
-- std::set<cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(), adjacent_cells_tmp.end());
-- std::set<cell_iterator> searched_cells;
--
-- // Determine the maximal number of cells
-- // in the grid.
-- // As long as we have not found
-- // the cell and have not searched
-- // every cell in the triangulation,
-- // we keep on looking.
-- const unsigned int n_cells =get_tria(container).n_cells();
-- bool found = false;
-- unsigned int cells_searched = 0;
-- while (!found && cells_searched < n_cells)
++ // Make sure that we have found
++ // at least one cell adjacent to vertex.
++ Assert(adjacent_cells_tmp.size()>0, ExcInternalError());
++
++ // Copy all the cells into a std::set
++ std::set<cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(), adjacent_cells_tmp.end());
++ std::set<cell_iterator> searched_cells;
++
++ // Determine the maximal number of cells
++ // in the grid.
++ // As long as we have not found
++ // the cell and have not searched
++ // every cell in the triangulation,
++ // we keep on looking.
++ const unsigned int n_cells =get_tria(container).n_cells();
++ bool found = false;
++ unsigned int cells_searched = 0;
++ while (!found && cells_searched < n_cells)
++ {
++ typename std::set<cell_iterator>::const_iterator
++ cell = adjacent_cells.begin(),
++ endc = adjacent_cells.end();
++ for(; cell != endc; ++cell)
{
-- typename std::set<cell_iterator>::const_iterator
-- cell = adjacent_cells.begin(),
-- endc = adjacent_cells.end();
-- for(; cell != endc; ++cell)
-- {
-- try
-- {
-- const Point<dim> p_cell = mapping.transform_real_to_unit_cell(*cell, p);
--
-- // calculate the infinity norm of
-- // the distance vector to the unit cell.
-- const double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);
--
-- // We compare if the point is inside the
-- // unit cell (or at least not too far
-- // outside). If it is, it is also checked
-- // that the cell has a more refined state
-- if (dist < best_distance ||
-- (dist == best_distance && (*cell)->level() > best_level))
-- {
-- found = true;
-- best_distance = dist;
-- best_level = (*cell)->level();
-- best_cell = std::make_pair(*cell, p_cell);
-- }
-- }
-- catch (typename MappingQ1<dim,spacedim>::ExcTransformationFailed &)
-- {
-- // ok, the transformation
-- // failed presumably
-- // because the point we
-- // are looking for lies
-- // outside the current
-- // cell. this means that
-- // the current cell can't
-- // be the cell around the
-- // point, so just ignore
-- // this cell and move on
-- // to the next
-- }
-- }
-- //udpate the number of cells searched
-- cells_searched += adjacent_cells.size();
-- // if we have not found the cell in
-- // question and have not yet searched every
-- // cell, we expand our search to
-- // all the not already searched neighbors of
-- // the cells in adjacent_cells. This is
-- // what find_active_cell_around_poin_internal
-- // is for.
-- if(!found && cells_searched < n_cells)
++ try
++ {
++ const Point<dim> p_cell = mapping.transform_real_to_unit_cell(*cell, p);
++
++ // calculate the infinity norm of
++ // the distance vector to the unit cell.
++ const double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);
++
++ // We compare if the point is inside the
++ // unit cell (or at least not too far
++ // outside). If it is, it is also checked
++ // that the cell has a more refined state
++ if (dist < best_distance ||
++ (dist == best_distance && (*cell)->level() > best_level))
{
-- find_active_cell_around_point_internal(container, searched_cells, adjacent_cells);
++ found = true;
++ best_distance = dist;
++ best_level = (*cell)->level();
++ best_cell = std::make_pair(*cell, p_cell);
}
++ }
++ catch (typename MappingQ1<dim,spacedim>::ExcTransformationFailed &)
++ {
++ // ok, the transformation
++ // failed presumably
++ // because the point we
++ // are looking for lies
++ // outside the current
++ // cell. this means that
++ // the current cell can't
++ // be the cell around the
++ // point, so just ignore
++ // this cell and move on
++ // to the next
++ }
++ }
++ //udpate the number of cells searched
++ cells_searched += adjacent_cells.size();
++ // if we have not found the cell in
++ // question and have not yet searched every
++ // cell, we expand our search to
++ // all the not already searched neighbors of
++ // the cells in adjacent_cells. This is
++ // what find_active_cell_around_poin_internal
++ // is for.
++ if(!found && cells_searched < n_cells)
++ {
++ find_active_cell_around_point_internal(container, searched_cells, adjacent_cells);
}
++ }
Assert (best_cell.first.state() == IteratorState::valid,
ExcPointNotFound<spacedim>(p));
// all adjacent cells
unsigned int vertex = find_closest_vertex(container, p);
-- std::vector<cell_iterator> adjacent_cells_tmp =
-- find_cells_adjacent_to_vertex(container, vertex);
--
-- // Make sure that we have found
-- // at least one cell adjacent to vertex.
-- Assert(adjacent_cells_tmp.size()>0, ExcInternalError());
--
-- // Copy all the cells into a std::set
-- std::set<cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(), adjacent_cells_tmp.end());
-- std::set<cell_iterator> searched_cells;
--
-- // Determine the maximal number of cells
-- // in the grid.
-- // As long as we have not found
-- // the cell and have not searched
-- // every cell in the triangulation,
-- // we keep on looking.
-- const unsigned int n_cells =get_tria(container).n_cells();
-- bool found = false;
-- unsigned int cells_searched = 0;
-- while (!found && cells_searched < n_cells)
++ std::vector<cell_iterator> adjacent_cells_tmp =
++ find_cells_adjacent_to_vertex(container, vertex);
++
++ // Make sure that we have found
++ // at least one cell adjacent to vertex.
++ Assert(adjacent_cells_tmp.size()>0, ExcInternalError());
++
++ // Copy all the cells into a std::set
++ std::set<cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(), adjacent_cells_tmp.end());
++ std::set<cell_iterator> searched_cells;
++
++ // Determine the maximal number of cells
++ // in the grid.
++ // As long as we have not found
++ // the cell and have not searched
++ // every cell in the triangulation,
++ // we keep on looking.
++ const unsigned int n_cells =get_tria(container).n_cells();
++ bool found = false;
++ unsigned int cells_searched = 0;
++ while (!found && cells_searched < n_cells)
{
typename std::set<cell_iterator>::const_iterator
-- cell = adjacent_cells.begin(),
-- endc = adjacent_cells.end();
-- for(; cell != endc; ++cell)
++ cell = adjacent_cells.begin(),
++ endc = adjacent_cells.end();
++ for(; cell != endc; ++cell)
++ {
++ try
{
-- try
-- {
-- const Point<dim> p_cell = mapping[(*cell)->active_fe_index()].transform_real_to_unit_cell(*cell, p);
--
--
-- // calculate the infinity norm of
-- // the distance vector to the unit cell.
-- const double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);
--
-- // We compare if the point is inside the
-- // unit cell (or at least not too far
-- // outside). If it is, it is also checked
-- // that the cell has a more refined state
-- if (dist < best_distance ||
-- (dist == best_distance && (*cell)->level() > best_level))
-- {
-- found = true;
-- best_distance = dist;
-- best_level = (*cell)->level();
-- best_cell = std::make_pair(*cell, p_cell);
-- }
-- }
-- catch (typename MappingQ1<dim,spacedim>::ExcTransformationFailed &)
-- {
-- // ok, the transformation
-- // failed presumably
-- // because the point we
-- // are looking for lies
-- // outside the current
-- // cell. this means that
-- // the current cell can't
-- // be the cell around the
-- // point, so just ignore
-- // this cell and move on
-- // to the next
-- }
++ const Point<dim> p_cell = mapping[(*cell)->active_fe_index()].transform_real_to_unit_cell(*cell, p);
++
++
++ // calculate the infinity norm of
++ // the distance vector to the unit cell.
++ const double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);
++
++ // We compare if the point is inside the
++ // unit cell (or at least not too far
++ // outside). If it is, it is also checked
++ // that the cell has a more refined state
++ if (dist < best_distance ||
++ (dist == best_distance && (*cell)->level() > best_level))
++ {
++ found = true;
++ best_distance = dist;
++ best_level = (*cell)->level();
++ best_cell = std::make_pair(*cell, p_cell);
++ }
}
-- //udpate the number of cells searched
-- cells_searched += adjacent_cells.size();
-- // if we have not found the cell in
-- // question and have not yet searched every
-- // cell, we expand our search to
-- // all the not already searched neighbors of
-- // the cells in adjacent_cells.
-- if(!found && cells_searched < n_cells)
++ catch (typename MappingQ1<dim,spacedim>::ExcTransformationFailed &)
{
-- find_active_cell_around_point_internal(container, searched_cells, adjacent_cells);
++ // ok, the transformation
++ // failed presumably
++ // because the point we
++ // are looking for lies
++ // outside the current
++ // cell. this means that
++ // the current cell can't
++ // be the cell around the
++ // point, so just ignore
++ // this cell and move on
++ // to the next
}
++ }
++ //udpate the number of cells searched
++ cells_searched += adjacent_cells.size();
++ // if we have not found the cell in
++ // question and have not yet searched every
++ // cell, we expand our search to
++ // all the not already searched neighbors of
++ // the cells in adjacent_cells.
++ if(!found && cells_searched < n_cells)
++ {
++ find_active_cell_around_point_internal(container, searched_cells, adjacent_cells);
++ }
}
-- }
++ }
Assert (best_cell.first.state() == IteratorState::valid,
ExcPointNotFound<spacedim>(p));
const CellIterator &cell2,
const int direction,
const dealii::Tensor<1,CellIterator::AccessorType::space_dimension>
-- &offset)
++ &offset)
{
-- // An orthogonal equality test for
-- // cells:
-- // Two cells are equal if the
-- // corresponding vertices of the
-- // boundary faces can be transformed
-- // into each other parallel to the
-- // given direction.
-- // It is assumed that
-- // cell1->face(2*direction) and
-- // cell2->face(2*direction+1) are the
-- // corresponding boundary faces.
-- // The optional argument offset will
-- // be added to each vertex of cell1
++ // An orthogonal equality test for
++ // cells:
++ // Two cells are equal if the
++ // corresponding vertices of the
++ // boundary faces can be transformed
++ // into each other parallel to the
++ // given direction.
++ // It is assumed that
++ // cell1->face(2*direction) and
++ // cell2->face(2*direction+1) are the
++ // corresponding boundary faces.
++ // The optional argument offset will
++ // be added to each vertex of cell1
static const int dim = CellIterator::AccessorType::dimension;
static const int space_dim = CellIterator::AccessorType::space_dimension;
Assert(dim == space_dim,
ExcNotImplemented());
-- // ... otherwise we would need two
-- // directions: One for selecting the
-- // faces, thus living in dim,
-- // the other direction for selecting the
-- // direction for the orthogonal
-- // projection, thus living in
-- // space_dim.
++ // ... otherwise we would need two
++ // directions: One for selecting the
++ // faces, thus living in dim,
++ // the other direction for selecting the
++ // direction for the orthogonal
++ // projection, thus living in
++ // space_dim.
for (int i = 0; i < space_dim; ++i) {
-- // Only compare
-- // coordinate-components != direction:
++ // Only compare
++ // coordinate-components != direction:
if (i == direction)
continue;
for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
if (fabs ( cell1->face(2*direction)->vertex(v)(i)
-- + offset[i]
-- - cell2->face(2*direction+1)->vertex(v)(i))
-- > 1.e-10)
++ + offset[i]
++ - cell2->face(2*direction+1)->vertex(v)(i))
++ > 1.e-10)
return false;
}
return true;
const types::boundary_id boundary_component,
int direction,
const dealii::Tensor<1,CellIterator::AccessorType::space_dimension>
-- &offset)
++ &offset)
{
static const int space_dim = CellIterator::AccessorType::space_dimension;
Assert (0<=direction && direction<space_dim,
std::set<CellIterator> cells1;
std::set<CellIterator> cells2;
-- // collect boundary cells, i.e. cells
-- // with face(2*direction), resp.
-- // face(2*direction+1) being on the
-- // boundary and having the correct
-- // color:
++ // collect boundary cells, i.e. cells
++ // with face(2*direction), resp.
++ // face(2*direction+1) being on the
++ // boundary and having the correct
++ // color:
CellIterator cell = begin;
for (; cell!= end; ++cell) {
if (cell->face(2*direction)->at_boundary() &&
std::map<CellIterator, CellIterator> matched_cells;
-- // Match with a complexity of O(n^2).
-- // This could be improved...
++ // Match with a complexity of O(n^2).
++ // This could be improved...
typedef typename std::set<CellIterator>::const_iterator SetIterator;
for (SetIterator it1 = cells1.begin(); it1 != cells1.end(); ++it1) {
for (SetIterator it2 = cells2.begin(); it2 != cells2.end(); ++it2) {
if (orthogonal_equality(*it1, *it2, direction, offset)) {
-- // We have a match, so insert the
-- // matching pairs and remove the matched
-- // cell in cells2 to speed up the
-- // matching:
++ // We have a match, so insert the
++ // matching pairs and remove the matched
++ // cell in cells2 to speed up the
++ // matching:
matched_cells[*it1] = *it2;
cells2.erase(it2);
break;
const types::boundary_id boundary_component,
int direction,
const dealii::Tensor<1,DH::space_dimension>
-- &offset)
++ &offset)
{
return collect_periodic_cell_pairs<typename DH::cell_iterator>
-- (dof_handler.begin_active(),
-- dof_handler.end(),
-- boundary_component,
-- direction,
-- offset);
++ (dof_handler.begin_active(),
++ dof_handler.end(),
++ boundary_component,
++ direction,
++ offset);
}
template<int dim, int spacedim>
void
DoFHandler<dim,spacedim>::
- compute_line_dof_identities (std::vector<unsigned int> &new_dof_indices) const
+ compute_line_dof_identities (std::vector<types::global_dof_index> &new_dof_indices) const
{
+ // we will mark lines that we
+ // have already treated, so first
+ // save and clear the user flags
+ // on lines and later restore
+ // them
+ std::vector<bool> user_flags;
+ this->get_tria().save_user_flags_line(user_flags);
+ const_cast<Triangulation<dim,spacedim> &>(this->get_tria()).clear_user_flags_line ();
+
// An implementation of the
// algorithm described in the hp
// paper, including the
// at most 2
Table<2,std_cxx1x::shared_ptr<internal::hp::DoFIdentities> >
line_dof_identities (finite_elements->size(),
- finite_elements->size());
+ finite_elements->size());
- for (line_iterator line=begin_line(); line!=end_line(); ++line)
- {
+ for (active_cell_iterator cell=begin_active(); cell!=end(); ++cell)
+ for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
+ if (cell->line(l)->user_flag_set() == false)
+ {
+ const line_iterator line = cell->line(l);
+ line->set_user_flag ();
+
unsigned int unique_sets_of_dofs
= line->n_active_fe_indices();
- // do a first loop over all sets of
- // dofs and do identity
- // uniquification
+ // do a first loop over all sets of
+ // dofs and do identity
+ // uniquification
for (unsigned int f=0; f<line->n_active_fe_indices(); ++f)
for (unsigned int g=f+1; g<line->n_active_fe_indices(); ++g)
- {
- const unsigned int fe_index_1 = line->nth_active_fe_index (f),
- fe_index_2 = line->nth_active_fe_index (g);
-
- if (((*finite_elements)[fe_index_1].dofs_per_line
- ==
- (*finite_elements)[fe_index_2].dofs_per_line)
- &&
- ((*finite_elements)[fe_index_1].dofs_per_line > 0))
- {
- internal::hp::ensure_existence_of_dof_identities<1>
- ((*finite_elements)[fe_index_1],
- (*finite_elements)[fe_index_2],
- line_dof_identities[fe_index_1][fe_index_2]);
- // see if these sets of dofs
- // are identical. the first
- // condition for this is that
- // indeed there are n
- // identities
- if (line_dof_identities[fe_index_1][fe_index_2]->size()
- ==
- (*finite_elements)[fe_index_1].dofs_per_line)
- {
- unsigned int i=0;
- for (; i<(*finite_elements)[fe_index_1].dofs_per_line; ++i)
- if (((*(line_dof_identities[fe_index_1][fe_index_2]))[i].first != i)
- &&
- ((*(line_dof_identities[fe_index_1][fe_index_2]))[i].second != i))
- // not an identity
- break;
-
- if (i == (*finite_elements)[fe_index_1].dofs_per_line)
- {
- // the dofs of these
- // two finite
- // elements are
- // identical. as a
- // safety check,
- // ensure that none
- // of the two FEs is
- // trying to dominate
- // the other, which
- // wouldn't make any
- // sense in this case
- Assert ((*finite_elements)[fe_index_1].compare_for_face_domination
- ((*finite_elements)[fe_index_2])
- ==
- FiniteElementDomination::either_element_can_dominate,
- ExcInternalError());
-
- --unique_sets_of_dofs;
-
- for (unsigned int j=0; j<(*finite_elements)[fe_index_1].dofs_per_line; ++j)
- {
- const types::global_dof_index master_dof_index
- = line->dof_index (j, fe_index_1);
- const types::global_dof_index slave_dof_index
- = line->dof_index (j, fe_index_2);
-
- // if master dof
- // was already
- // constrained,
- // constrain to
- // that one,
- // otherwise
- // constrain
- // slave to
- // master
- if (new_dof_indices[master_dof_index] !=
- numbers::invalid_unsigned_int)
- {
- Assert (new_dof_indices[new_dof_indices[master_dof_index]] ==
- numbers::invalid_unsigned_int,
- ExcInternalError());
-
- new_dof_indices[slave_dof_index]
- = new_dof_indices[master_dof_index];
- }
- else
- {
- Assert ((new_dof_indices[master_dof_index] ==
- numbers::invalid_unsigned_int)
- ||
- (new_dof_indices[slave_dof_index] ==
- master_dof_index),
- ExcInternalError());
-
- new_dof_indices[slave_dof_index] = master_dof_index;
- }
- }
- }
- }
- }
- }
-
- // if at this point, there is only
- // one unique set of dofs left, then
- // we have taken care of everything
- // above. if there are two, then we
- // need to deal with them here. if
- // there are more, then we punt, as
- // described in the paper (and
- // mentioned above)
+ {
+ const unsigned int fe_index_1 = line->nth_active_fe_index (f),
+ fe_index_2 = line->nth_active_fe_index (g);
+
+ if (((*finite_elements)[fe_index_1].dofs_per_line
+ ==
+ (*finite_elements)[fe_index_2].dofs_per_line)
+ &&
+ ((*finite_elements)[fe_index_1].dofs_per_line > 0))
+ {
+ internal::hp::ensure_existence_of_dof_identities<1>
+ ((*finite_elements)[fe_index_1],
+ (*finite_elements)[fe_index_2],
+ line_dof_identities[fe_index_1][fe_index_2]);
+ // see if these sets of dofs
+ // are identical. the first
+ // condition for this is that
+ // indeed there are n
+ // identities
+ if (line_dof_identities[fe_index_1][fe_index_2]->size()
+ ==
+ (*finite_elements)[fe_index_1].dofs_per_line)
+ {
+ unsigned int i=0;
+ for (; i<(*finite_elements)[fe_index_1].dofs_per_line; ++i)
+ if (((*(line_dof_identities[fe_index_1][fe_index_2]))[i].first != i)
+ &&
+ ((*(line_dof_identities[fe_index_1][fe_index_2]))[i].second != i))
+ // not an identity
+ break;
+
+ if (i == (*finite_elements)[fe_index_1].dofs_per_line)
+ {
+ // the dofs of these
+ // two finite
+ // elements are
+ // identical. as a
+ // safety check,
+ // ensure that none
+ // of the two FEs is
+ // trying to dominate
+ // the other, which
+ // wouldn't make any
+ // sense in this case
+ Assert ((*finite_elements)[fe_index_1].compare_for_face_domination
+ ((*finite_elements)[fe_index_2])
+ ==
+ FiniteElementDomination::either_element_can_dominate,
+ ExcInternalError());
+
+ --unique_sets_of_dofs;
+
+ for (unsigned int j=0; j<(*finite_elements)[fe_index_1].dofs_per_line; ++j)
+ {
- const unsigned int master_dof_index
++ const types::global_dof_index master_dof_index
+ = line->dof_index (j, fe_index_1);
- const unsigned int slave_dof_index
++ const types::global_dof_index slave_dof_index
+ = line->dof_index (j, fe_index_2);
+
+ // if master dof
+ // was already
+ // constrained,
+ // constrain to
+ // that one,
+ // otherwise
+ // constrain
+ // slave to
+ // master
+ if (new_dof_indices[master_dof_index] !=
+ numbers::invalid_unsigned_int)
+ {
+ Assert (new_dof_indices[new_dof_indices[master_dof_index]] ==
+ numbers::invalid_unsigned_int,
+ ExcInternalError());
+
+ new_dof_indices[slave_dof_index]
+ = new_dof_indices[master_dof_index];
+ }
+ else
+ {
+ Assert ((new_dof_indices[master_dof_index] ==
+ numbers::invalid_unsigned_int)
+ ||
+ (new_dof_indices[slave_dof_index] ==
+ master_dof_index),
+ ExcInternalError());
+
+ new_dof_indices[slave_dof_index] = master_dof_index;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // if at this point, there is only
+ // one unique set of dofs left, then
+ // we have taken care of everything
+ // above. if there are two, then we
+ // need to deal with them here. if
+ // there are more, then we punt, as
+ // described in the paper (and
+ // mentioned above)
//TODO: The check for 'dim==2' was inserted by intuition. It fixes
// the previous problems with step-27 in 3D. But an explanation
// for this is still required, and what we do here is not what we
// describe in the paper!.
if ((unique_sets_of_dofs == 2) && (dim == 2))
{
- // find out which is the
- // most dominating finite
- // element of the ones that
- // are used on this line
- const unsigned int most_dominating_fe_index
- = internal::hp::get_most_dominating_fe_index<dim,spacedim> (line);
-
- const unsigned int n_active_fe_indices
- = line->n_active_fe_indices ();
-
- // loop over the indices of
- // all the finite elements
- // that are not dominating,
- // and identify their dofs
- // to the most dominating
- // one
- for (unsigned int f=0; f<n_active_fe_indices; ++f)
- if (line->nth_active_fe_index (f) !=
- most_dominating_fe_index)
- {
- const unsigned int
- other_fe_index = line->nth_active_fe_index (f);
-
- internal::hp::ensure_existence_of_dof_identities<1>
- ((*finite_elements)[most_dominating_fe_index],
- (*finite_elements)[other_fe_index],
- line_dof_identities[most_dominating_fe_index][other_fe_index]);
-
- internal::hp::DoFIdentities &identities
- = *line_dof_identities[most_dominating_fe_index][other_fe_index];
- for (unsigned int i=0; i<identities.size(); ++i)
- {
- const types::global_dof_index master_dof_index
- = line->dof_index (identities[i].first, most_dominating_fe_index);
- const types::global_dof_index slave_dof_index
- = line->dof_index (identities[i].second, other_fe_index);
-
- Assert ((new_dof_indices[master_dof_index] ==
- numbers::invalid_unsigned_int)
- ||
- (new_dof_indices[slave_dof_index] ==
- master_dof_index),
- ExcInternalError());
-
- new_dof_indices[slave_dof_index] = master_dof_index;
- }
- }
+ // find out which is the
+ // most dominating finite
+ // element of the ones that
+ // are used on this line
+ const unsigned int most_dominating_fe_index
+ = internal::hp::get_most_dominating_fe_index<dim,spacedim> (line);
+
+ const unsigned int n_active_fe_indices
+ = line->n_active_fe_indices ();
+
+ // loop over the indices of
+ // all the finite elements
+ // that are not dominating,
+ // and identify their dofs
+ // to the most dominating
+ // one
+ for (unsigned int f=0; f<n_active_fe_indices; ++f)
+ if (line->nth_active_fe_index (f) !=
+ most_dominating_fe_index)
+ {
+ const unsigned int
+ other_fe_index = line->nth_active_fe_index (f);
+
+ internal::hp::ensure_existence_of_dof_identities<1>
+ ((*finite_elements)[most_dominating_fe_index],
+ (*finite_elements)[other_fe_index],
+ line_dof_identities[most_dominating_fe_index][other_fe_index]);
+
+ internal::hp::DoFIdentities &identities
+ = *line_dof_identities[most_dominating_fe_index][other_fe_index];
+ for (unsigned int i=0; i<identities.size(); ++i)
+ {
- const unsigned int master_dof_index
++ const types::global_dof_index master_dof_index
+ = line->dof_index (identities[i].first, most_dominating_fe_index);
- const unsigned int slave_dof_index
++ const types::global_dof_index slave_dof_index
+ = line->dof_index (identities[i].second, other_fe_index);
+
+ Assert ((new_dof_indices[master_dof_index] ==
+ numbers::invalid_unsigned_int)
+ ||
+ (new_dof_indices[slave_dof_index] ==
+ master_dof_index),
+ ExcInternalError());
+
+ new_dof_indices[slave_dof_index] = master_dof_index;
+ }
+ }
}
- }
+ }
+
+ // finally restore the user flags
+ const_cast<Triangulation<dim,spacedim> &>(this->get_tria())
+ .load_user_flags_line(user_flags);
}
template<int dim, int spacedim>
void
DoFHandler<dim,spacedim>::
- compute_quad_dof_identities (std::vector<unsigned int> &new_dof_indices) const
+ compute_quad_dof_identities (std::vector<types::global_dof_index> &new_dof_indices) const
{
+ // we will mark quads that we
+ // have already treated, so first
+ // save and clear the user flags
+ // on quads and later restore
+ // them
+ std::vector<bool> user_flags;
+ this->get_tria().save_user_flags_quad(user_flags);
+ const_cast<Triangulation<dim,spacedim> &>(this->get_tria()).clear_user_flags_quad ();
+
// An implementation of the
// algorithm described in the hp
// paper, including the
// particularly frequent case
Table<2,std_cxx1x::shared_ptr<internal::hp::DoFIdentities> >
quad_dof_identities (finite_elements->size(),
- finite_elements->size());
+ finite_elements->size());
- for (quad_iterator quad=begin_quad(); quad!=end_quad(); ++quad)
- if (quad->n_active_fe_indices() == 2)
- {
- // find out which is the
- // most dominating finite
- // element of the ones that
- // are used on this quad
- const unsigned int most_dominating_fe_index
- = internal::hp::get_most_dominating_fe_index<dim,spacedim> (quad);
-
- const unsigned int n_active_fe_indices
- = quad->n_active_fe_indices ();
-
- // loop over the indices of
- // all the finite elements
- // that are not dominating,
- // and identify their dofs
- // to the most dominating
- // one
- for (unsigned int f=0; f<n_active_fe_indices; ++f)
- if (quad->nth_active_fe_index (f) !=
- most_dominating_fe_index)
- {
- const unsigned int
- other_fe_index = quad->nth_active_fe_index (f);
-
- internal::hp::ensure_existence_of_dof_identities<2>
- ((*finite_elements)[most_dominating_fe_index],
- (*finite_elements)[other_fe_index],
- quad_dof_identities[most_dominating_fe_index][other_fe_index]);
-
- internal::hp::DoFIdentities &identities
- = *quad_dof_identities[most_dominating_fe_index][other_fe_index];
- for (unsigned int i=0; i<identities.size(); ++i)
- {
- const types::global_dof_index master_dof_index
- = quad->dof_index (identities[i].first, most_dominating_fe_index);
- const types::global_dof_index slave_dof_index
- = quad->dof_index (identities[i].second, other_fe_index);
-
- Assert ((new_dof_indices[master_dof_index] ==
- numbers::invalid_unsigned_int)
- ||
- (new_dof_indices[slave_dof_index] ==
- master_dof_index),
- ExcInternalError());
-
- new_dof_indices[slave_dof_index] = master_dof_index;
- }
- }
- }
+ for (active_cell_iterator cell=begin_active(); cell!=end(); ++cell)
+ for (unsigned int q=0; q<GeometryInfo<dim>::quads_per_cell; ++q)
+ if ((cell->quad(q)->user_flag_set() == false)
+ &&
+ (cell->quad(q)->n_active_fe_indices() == 2))
+ {
+ const quad_iterator quad = cell->quad(q);
+ quad->set_user_flag ();
+
+ // find out which is the
+ // most dominating finite
+ // element of the ones that
+ // are used on this quad
+ const unsigned int most_dominating_fe_index
+ = internal::hp::get_most_dominating_fe_index<dim,spacedim> (quad);
+
+ const unsigned int n_active_fe_indices
+ = quad->n_active_fe_indices ();
+
+ // loop over the indices of
+ // all the finite elements
+ // that are not dominating,
+ // and identify their dofs
+ // to the most dominating
+ // one
+ for (unsigned int f=0; f<n_active_fe_indices; ++f)
+ if (quad->nth_active_fe_index (f) !=
+ most_dominating_fe_index)
+ {
+ const unsigned int
+ other_fe_index = quad->nth_active_fe_index (f);
+
+ internal::hp::ensure_existence_of_dof_identities<2>
+ ((*finite_elements)[most_dominating_fe_index],
+ (*finite_elements)[other_fe_index],
+ quad_dof_identities[most_dominating_fe_index][other_fe_index]);
+
+ internal::hp::DoFIdentities &identities
+ = *quad_dof_identities[most_dominating_fe_index][other_fe_index];
+ for (unsigned int i=0; i<identities.size(); ++i)
+ {
- const unsigned int master_dof_index
++ const types::global_dof_index master_dof_index
+ = quad->dof_index (identities[i].first, most_dominating_fe_index);
- const unsigned int slave_dof_index
++ const types::global_dof_index slave_dof_index
+ = quad->dof_index (identities[i].second, other_fe_index);
+
+ Assert ((new_dof_indices[master_dof_index] ==
+ numbers::invalid_unsigned_int)
+ ||
+ (new_dof_indices[slave_dof_index] ==
+ master_dof_index),
+ ExcInternalError());
+
+ new_dof_indices[slave_dof_index] = master_dof_index;
+ }
+ }
+ }
+
+ // finally restore the user flags
+ const_cast<Triangulation<dim,spacedim> &>(this->get_tria())
+ .load_user_flags_quad(user_flags);
}
template<>
void
DoFHandler<2,2>::
- renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
- internal::int2type<2>)
+ renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
+ internal::int2type<2>)
{
+ const unsigned int dim = 2;
+ const unsigned int spacedim = 2;
+
Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
renumber_dofs_internal (new_numbers, internal::int2type<1>());
template<>
void
DoFHandler<2,3>::
- renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
- internal::int2type<2>)
+ renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
+ internal::int2type<2>)
{
+ const unsigned int dim = 2;
+ const unsigned int spacedim = 3;
+
Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
renumber_dofs_internal (new_numbers, internal::int2type<1>());
template<>
void
DoFHandler<3,3>::
- renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
- internal::int2type<2>)
+ renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
+ internal::int2type<2>)
{
+ const unsigned int dim = 3;
+ const unsigned int spacedim = 3;
+
Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
renumber_dofs_internal (new_numbers, internal::int2type<1>());
template<>
void
DoFHandler<3,3>::
- renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
- internal::int2type<3>)
+ renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
+ internal::int2type<3>)
{
+ const unsigned int dim = 3;
+ const unsigned int spacedim = 3;
+
Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
renumber_dofs_internal (new_numbers, internal::int2type<2>());
projects / dealii-svn.git / commitdiff