void
set_dof_indices(const std::vector<types::global_dof_index> &dof_indices);
+ /**
+ * Set the DoF indices of this cell to the given values. This function
+ * bypasses the DoF cache, if one exists for the given DoF handler class.
+ */
+ void
+ set_non_local_dof_indices(
+ const std::vector<types::global_dof_index> &non_local_dof_indices);
+
/**
* Set the Level DoF indices of this cell to the given values.
*/
{
AssertDimension(dof_indices.size(), n_dof_indices(accessor, fe_index));
- const auto &fe = accessor.get_fe(fe_index);
+ const auto &fe = accessor.get_fe(fe_index);
+ const auto non_local_dofs_per_cell = fe.n_non_local_dofs();
unsigned int index = 0;
dof_indices,
fe_index);
+ // 5) non_local dofs
+ for (unsigned int d = 0; d < non_local_dofs_per_cell; ++d, ++index)
+ {
+ // TODO: Check if we need DoFOperation::process_dof also here!
+ }
+
AssertDimension(dof_indices.size(), index);
}
"been initialized, i.e., it doesn't appear that DoF indices "
"have been distributed on it."));
+ switch (structdim)
+ {
+ case 1:
+ Assert(
+ dof_indices.size() ==
+ (this->n_vertices() *
+ this->dof_handler->get_fe(fe_index).n_dofs_per_vertex() +
+ this->dof_handler->get_fe(fe_index).n_dofs_per_line()) +
+ this->dof_handler->get_fe(fe_index).n_non_local_dofs_per_cell(),
+ ExcVectorDoesNotMatch());
+ break;
+ case 2:
+ Assert(
+ dof_indices.size() ==
+ (this->n_vertices() *
+ this->dof_handler->get_fe(fe_index).n_dofs_per_vertex() +
+ this->n_lines() *
+ this->dof_handler->get_fe(fe_index).n_dofs_per_line() +
+ this->dof_handler->get_fe(fe_index).n_dofs_per_quad()) +
+ this->dof_handler->get_fe(fe_index).n_non_local_dofs_per_cell(),
+ ExcVectorDoesNotMatch());
+ break;
+ case 3:
+ Assert(
+ dof_indices.size() ==
+ (this->n_vertices() *
+ this->dof_handler->get_fe(fe_index).n_dofs_per_vertex() +
+ this->n_lines() *
+ this->dof_handler->get_fe(fe_index).n_dofs_per_line() +
+ this->n_faces() *
+ this->dof_handler->get_fe(fe_index).n_dofs_per_quad() +
+ this->dof_handler->get_fe(fe_index).n_dofs_per_hex()) +
+ this->dof_handler->get_fe(fe_index).n_non_local_dofs_per_cell(),
+ ExcVectorDoesNotMatch());
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+
+
// this function really only makes sense if either a) there are degrees of
// freedom defined on the present object, or b) the object is non-active
- // objects but all degrees of freedom are located on vertices, since
- // otherwise there are degrees of freedom on sub-objects which are not
- // allocated for this non-active thing
+ // objects but all degrees of freedom are located on vertices, since otherwise
+ // there are degrees of freedom on sub-objects which are not allocated for
+ // this non-active thing
Assert(this->fe_index_is_active(fe_index) ||
(this->dof_handler->get_fe(fe_index).n_dofs_per_cell() ==
this->n_vertices() *
+ /**
+ * Implement setting dof indices on a cell. TO CHECK ZHAOWEI + LH
+ */
+ template <int dim, int spacedim, bool level_dof_access>
+ static void
+ set_dof_indices(
+ const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor,
+ const std::vector<types::global_dof_index> &local_dof_indices)
+ {
+ Assert(accessor.has_children() == false, ExcInternalError());
+
+ const unsigned int dofs_per_vertex =
+ accessor.get_fe().n_dofs_per_vertex(),
+ dofs_per_line = accessor.get_fe().n_dofs_per_line(),
+ dofs_per_quad = accessor.get_fe().n_dofs_per_quad(),
+ dofs_per_hex = accessor.get_fe().n_dofs_per_hex();
+
+ Assert(local_dof_indices.size() == accessor.get_fe().dofs_per_cell,
+ ExcInternalError());
+
+ unsigned int index = 0;
+
+ for (unsigned int vertex = 0;
+ vertex < GeometryInfo<dim>::vertices_per_cell;
+ ++vertex)
+ for (unsigned int d = 0; d < dofs_per_vertex; ++d, ++index)
+ accessor.set_vertex_dof_index(vertex,
+ d,
+ local_dof_indices[index],
+ accessor.active_fe_index());
+ // now copy dof numbers into the line. for lines in 3d with the
+ // wrong orientation, we have already made sure that we're ok
+ // by picking the correct vertices (this happens automatically
+ // in the vertex() function). however, if the line is in wrong
+ // orientation, we look at it in flipped orientation and we
+ // will have to adjust the shape function indices that we see
+ // to correspond to the correct (cell-local) ordering.
+ //
+ // of course, if dim<3, then there is nothing to adjust
+ for (unsigned int line = 0; line < GeometryInfo<dim>::lines_per_cell;
+ ++line)
+ for (unsigned int d = 0; d < dofs_per_line; ++d, ++index)
+ accessor.line(line)->set_dof_index(
+ dim < 3 ?
+ d :
+ accessor.get_fe().adjust_line_dof_index_for_line_orientation(
+ d, accessor.line_orientation(line)),
+ local_dof_indices[index],
+ accessor.active_fe_index());
+ // now copy dof numbers into the face. for faces in 3d with the
+ // wrong orientation, we have already made sure that we're ok
+ // by picking the correct lines and vertices (this happens
+ // automatically in the line() and vertex()
+ // functions). however, if the face is in wrong orientation,
+ // we look at it in flipped orientation and we will have to
+ // adjust the shape function indices that we see to correspond
+ // to the correct (cell-local) ordering. The same applies, if
+ // the face_rotation or face_orientation is non-standard
+ //
+ // again, if dim<3, then there is nothing to adjust
+ for (unsigned int quad = 0; quad < GeometryInfo<dim>::quads_per_cell;
+ ++quad)
+ for (unsigned int d = 0; d < dofs_per_quad; ++d, ++index)
+ accessor.quad(quad)->set_dof_index(
+ dim < 3 ?
+ d :
+ accessor.get_fe().adjust_quad_dof_index_for_face_orientation(
+ d,
+ accessor.face_orientation(quad),
+ accessor.face_flip(quad),
+ accessor.face_rotation(quad)),
+ local_dof_indices[index],
+ accessor.active_fe_index());
+ for (unsigned int d = 0; d < dofs_per_hex; ++d, ++index)
+ accessor.set_dof_index(d,
+ local_dof_indices[index],
+ accessor.active_fe_index());
+ Assert(index == accessor.get_fe().dofs_per_cell, ExcInternalError());
+ }
+
+
+
+ /**
+ * Implement setting non-local dof indices on a cell.
+ */
+ template <int dim, int spacedim, bool level_dof_access>
+ static void
+ set_non_local_dof_indices(
+ const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor,
+ const std::vector<types::global_dof_index> &local_non_local_dof_indices)
+ {
+ Assert(accessor.has_children() == false, ExcInternalError());
+
+ const unsigned int dofs_per_cell = accessor.get_fe().dofs_per_cell;
+
+ types::global_dof_index *next_dof_index =
+ const_cast<types::global_dof_index *>(
+ dealii::internal::DoFAccessorImplementation::Implementation::
+ get_cache_ptr(accessor.dof_handler,
+ accessor.present_level,
+ accessor.present_index,
+ dofs_per_cell));
+
+ const unsigned int non_local_dofs =
+ accessor.get_fe().n_non_local_dofs_per_cell(),
+ n_dofs = accessor.get_fe().n_dofs_per_cell();
+
+ unsigned int index = 0;
+
+ for (unsigned int d = 0; d < n_dofs; ++d, ++next_dof_index)
+ if (d >= n_dofs - non_local_dofs)
+ {
+ *next_dof_index = local_non_local_dof_indices[index];
+ ++index;
+ }
+ Assert(index == accessor.get_fe().n_non_local_dofs_per_cell(),
+ ExcInternalError());
+ }
+
+
/**
* Do what the active_fe_index function in the parent class is supposed to
* do.
*/
DoFObjects<dim> dof_object;
+ /**
+ * The offsets for each cell of the cache that holds all DoF indices.
+ * Only used for hp elements.
+ */
+ std::vector<unsigned int> cell_cache_offsets;
+
/**
* Return a pointer to the beginning of the DoF indices cache for a
* given cell.
// ---------------------------------------------------------------------
//
-// Copyright (C) 2000 - 2018 by the deal.II authors
+// Copyright (C) 2000 - 2020 by the deal.II authors
//
// This file is part of the deal.II library.
//
*/
const unsigned int dofs_per_hex;
+ /**
+ * Number of non-zero degrees of freedom (DoFs) on the cell that,
+ * however, can not be associated with the specific geometric object (such as
+ * face or edge). In these cases, DoFs are typically globally distributed over
+ * (patches of) the Triangulation. Typical examples are:
+ * - <a
+ * href="https://www.sciencedirect.com/science/article/pii/S0045782504005171">Isogeometric
+ * Analysis</a>
+ * - <a
+ * href="https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291097-0207%2819990910%2946%3A1%3C131%3A%3AAID-NME726%3E3.0.CO%3B2-J">Extended
+ * finite element method</a>
+ * - <a
+ * href="http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.21.4515&rep=rep1&type=pdf">Catmull-Clark's
+ * subdivision surfaces</a>
+ */
+ const unsigned int non_local_dofs_per_cell;
+
/**
* First index of dof on a line.
*/
* to geometrical objects.
*
* @param[in] dofs_per_object A vector that describes the number of degrees
- * of freedom on geometrical objects for each dimension. This vector must
- * have size dim+1, and entry 0 describes the number of degrees of freedom
- * per vertex, entry 1 the number of degrees of freedom per line, etc. As an
- * example, for the common $Q_1$ Lagrange element in 2d, this vector would
- * have elements <code>(1,0,0)</code>. On the other hand, for a $Q_3$
- * element in 3d, it would have entries <code>(1,2,4,8)</code>.
+ * of freedom on geometrical objects for each dimension.
+ * The first dim+1 elements of the vector describe the number
+ * of degrees of freedom per geometrical object: entry 0 describes the number
+ * of degrees of freedom per vertex, entry 1 the number of degrees of freedom
+ * per line, etc. If the vector is of size dim+2, the last entry specifices
+ * the number of non-local degrees of freedom for this element. In case the
+ * vector size is dim+1, the number of non-local degrees of freedom is set to
+ * zero. As an example, for the common $Q_1$ Lagrange element in 2d, this
+ * vector would have elements <code>(1,0,0)</code>. On the other hand, for a
+ * $Q_3$ element in 3d, it would have entries <code>(1,2,4,8)</code>.
*
* @param[in] n_components Number of vector components of the element.
*
unsigned int
n_dofs_per_cell() const;
+ /**
+ * Return the number of non-zero degrees of freedom (DoFs) on the cell that,
+ * however, can not be associated with the specific geometric object (such as
+ * face or edge). In these cases, DoFs are typically globally distributed over
+ * (patches of) the Triangulation. Typical examples are:
+ * - <a
+ * href="https://www.sciencedirect.com/science/article/pii/S0045782504005171">Isogeometric
+ * Analysis</a>
+ * - <a
+ * href="https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291097-0207%2819990910%2946%3A1%3C131%3A%3AAID-NME726%3E3.0.CO%3B2-J">Extended
+ * finite element method</a>
+ * - <a
+ * href="http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.21.4515&rep=rep1&type=pdf">Catmull-Clark's
+ * subdivision surfaces</a>
+ */
+ unsigned int
+ n_non_local_dofs_per_cell() const;
+
/**
* Return the number of degrees per structdim-dimensional object. For
* structdim==0, the function therefore returns dofs_per_vertex, for
+template <int dim>
+inline unsigned int
+FiniteElementData<dim>::n_non_local_dofs_per_cell() const
+{
+ return non_local_dofs_per_cell;
+}
+
+
+
template <int dim>
template <int structdim>
inline unsigned int
{
AssertDimension(fes.size(), multiplicities.size());
- unsigned int multiplied_dofs_per_vertex = 0;
- unsigned int multiplied_dofs_per_line = 0;
- unsigned int multiplied_dofs_per_quad = 0;
- unsigned int multiplied_dofs_per_hex = 0;
+ unsigned int multiplied_dofs_per_vertex = 0;
+ unsigned int multiplied_dofs_per_line = 0;
+ unsigned int multiplied_dofs_per_quad = 0;
+ unsigned int multiplied_dofs_per_hex = 0;
+ unsigned int multiplied_non_local_dofs_per_cell = 0;
unsigned int multiplied_n_components = 0;
fes[i]->n_dofs_per_quad() * multiplicities[i];
multiplied_dofs_per_hex +=
fes[i]->n_dofs_per_hex() * multiplicities[i];
+ multiplied_non_local_dofs_per_cell +=
+ fes[i]->n_non_local_dofs_per_cell() * multiplicities[i];
multiplied_n_components +=
fes[i]->n_components() * multiplicities[i];
dpo.push_back(multiplied_dofs_per_quad);
if (dim > 2)
dpo.push_back(multiplied_dofs_per_hex);
+ if (multiplied_non_local_dofs_per_cell > 0)
+ dpo.push_back(multiplied_non_local_dofs_per_cell);
BlockIndices block_indices(0, 0);
fes[base]->restriction_is_additive(index_in_base);
}
}
+ unsigned int total_index_offset = total_index;
+ // 5. Non_local
+ unsigned int n_fe = 0;
+ for (unsigned int base = 0; base < fes.size(); ++base)
+ {
+ if (fes[base] != NULL && fes[base]->n_non_local_dofs_per_cell() > 0)
+ {
+ ++n_fe;
+ }
+ }
+ for (unsigned int base = 0; base < n_fe; ++base)
+ for (unsigned int m = 0; m < multiplicities[base]; ++m)
+ for (unsigned int local_index = 0;
+ local_index < fes[base]->n_non_local_dofs_per_cell();
+ ++local_index)
+ {
+ total_index = total_index_offset +
+ local_index * (n_fe * multiplicities[base]) +
+ base * multiplicities[base] + m;
+
+ const unsigned int index_in_base =
+ (local_index + (fes[base]->n_dofs_per_cell() -
+ fes[base]->n_non_local_dofs_per_cell()));
+ Assert(index_in_base < fes[base]->n_dofs_per_cell(),
+ ExcInternalError());
+ retval[total_index] =
+ fes[base]->restriction_is_additive(index_in_base);
+ }
+ // TODO[LH]: Make sure this Assert here is correct
Assert(total_index == n_shape_functions, ExcInternalError());
return retval;
{
AssertDimension(fes.size(), multiplicities.size());
- // first count the number of dofs and components that will emerge from the
- // given FEs
+ // first count the number of dofs and components that will emerge
+ // from the given FEs
unsigned int n_shape_functions = 0;
for (unsigned int i = 0; i < fes.size(); ++i)
if (multiplicities[i] > 0) // needed because fe might be nullptr
std::vector<bool>(n_components,
false));
- // finally go through all the shape functions of the base elements, and
- // copy their flags. this somehow copies the code in build_cell_table,
- // which is not nice as it uses too much implicit knowledge about the
- // layout of the individual bases in the composed FE, but there seems no
- // way around...
+ // finally go through all the shape functions of the base elements,
+ // and copy their flags. this somehow copies the code in
+ // build_cell_table, which is not nice as it uses too much implicit
+ // knowledge about the layout of the individual bases in the
+ // composed FE, but there seems no way around...
//
- // for each shape function, copy the non-zero flags from the base element
- // to this one, taking into account multiplicities, multiple components in
- // base elements, and other complications
+ // for each shape function, copy the non-zero flags from the base
+ // element to this one, taking into account multiplicities, multiple
+ // components in base elements, and other complications
unsigned int total_index = 0;
for (const unsigned int vertex_number :
ReferenceCell::internal::Info::get_cell(
}
}
}
+ unsigned int total_index_offset = total_index;
+ // 5. non_local
+ unsigned int comp_start = 0;
+ unsigned int n_fe = 0;
+ for (unsigned int base = 0; base < fes.size(); ++base)
+ {
+ if (fes[base] != NULL && fes[base]->n_non_local_dofs_per_cell() > 0)
+ {
+ ++n_fe;
+ }
+ }
+ for (unsigned int base = 0; base < n_fe; ++base)
+ {
+ for (unsigned int m = 0; m < multiplicities[base];
+ ++m, comp_start += fes[base]->n_components() * do_tensor_product)
+ {
+ for (unsigned int non_local_index = 0;
+ non_local_index < fes[base]->n_non_local_dofs_per_cell();
+ ++non_local_index)
+ {
+ total_index =
+ total_index_offset +
+ non_local_index * (n_fe * multiplicities[base]) +
+ base * multiplicities[base] + m;
+
+ const unsigned int index_in_base =
+ (non_local_index +
+ (fes[base]->n_dofs_per_cell() -
+ fes[base]->n_non_local_dofs_per_cell()));
+
+ Assert(comp_start + fes[base]->n_components() <=
+ retval[total_index].size(),
+ ExcInternalError());
+
+ for (unsigned int c = 0; c < fes[base]->n_components(); ++c)
+ {
+ Assert(c < fes[base]
+ ->get_nonzero_components(index_in_base)
+ .size(),
+ ExcInternalError());
+ retval[total_index][comp_start + c] =
+ fes[base]->get_nonzero_components(index_in_base)[c];
+ }
+ }
+ }
+ }
+ // TODO[LH]: Make sure this assert here is correct
Assert(total_index == n_shape_functions, ExcInternalError());
// now copy the vector<vector<bool> > into a vector<ComponentMask>.
non_primitive_index;
}
}
+ unsigned int total_index_offset = total_index;
+ // 5. non-local
+ for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int m = 0; m < fe.element_multiplicity(base);
+ ++m,
+ comp_start += fe.base_element(base).n_components() *
+ do_tensor_product)
+ {
+ for (unsigned int non_local_index = 0;
+ non_local_index <
+ fe.base_element(base).n_non_local_dofs_per_cell();
+ ++non_local_index)
+ {
+ total_index =
+ total_index_offset +
+ non_local_index *
+ (fe.n_base_elements() * fe.element_multiplicity(base)) +
+ base * fe.element_multiplicity(base) + m;
+
+ const unsigned int index_in_base =
+ (non_local_index +
+ (fe.base_element(base).n_dofs_per_cell() -
+ fe.base_element(base).n_non_local_dofs_per_cell()));
+
+ system_to_base_table[total_index] =
+ std::make_pair(std::make_pair(base, m), index_in_base);
+
+ if (fe.base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base =
+ fe.base_element(base)
+ .system_to_component_index(index_in_base)
+ .first;
+ const unsigned int comp = comp_start + comp_in_base;
+ const unsigned int index_in_comp =
+ fe.base_element(base)
+ .system_to_component_index(index_in_base)
+ .second;
+ system_to_component_table[total_index] =
+ std::make_pair(comp, index_in_comp);
+ }
+ else
+ {
+ system_to_component_table[total_index] =
+ non_primitive_index;
+ }
+ }
+ }
+ }
}
{
// get (cell) index of this shape function inside the base
// element to see whether the shape function is primitive
- // (assume that all shape functions on vertices share the same
- // primitivity property; assume likewise for all shape
+ // (assume that all shape functions on vertices share the
+ // same primitivity property; assume likewise for all shape
// functions located on lines, quads, etc. this way, we can
- // ask for primitivity of only _one_ shape function, which is
- // taken as representative for all others located on the same
- // type of object):
+ // ask for primitivity of only _one_ shape function, which
+ // is taken as representative for all others located on the
+ // same type of object):
const unsigned int index_in_base =
(fe.base_element(base).n_dofs_per_vertex() * vertex_number +
local_index);
+template <int dim, int spacedim, bool lda>
+void
+DoFCellAccessor<dim, spacedim, lda>::set_non_local_dof_indices(
+ const std::vector<types::global_dof_index> &local_non_local_dof_indices)
+{
+ Assert(static_cast<unsigned int>(this->present_level) <
+ this->dof_handler->object_dof_indices.size(),
+ ExcMessage("DoFHandler not initialized"));
+
+ Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject());
+
+ internal::DoFCellAccessorImplementation::Implementation::
+ set_non_local_dof_indices(*this, local_non_local_dof_indices);
+}
+
+
+
template <int dim, int spacedim, bool lda>
TriaIterator<DoFCellAccessor<dim, spacedim, lda>>
DoFCellAccessor<dim, spacedim, lda>::neighbor_child_on_subface(
// ---------------------------------------------------------------------
//
-// Copyright (C) 2001 - 2018 by the deal.II authors
+// Copyright (C) 2001 - 2020 by the deal.II authors
//
// This file is part of the deal.II library.
//
, dofs_per_line(dofs_per_object[1])
, dofs_per_quad(dim > 1 ? dofs_per_object[2] : 0)
, dofs_per_hex(dim > 2 ? dofs_per_object[3] : 0)
+ , non_local_dofs_per_cell(
+ dofs_per_object.size() == dim + 2 ? dofs_per_object[dim + 1] : 0)
, first_line_index(
ReferenceCell::internal::Info::get_cell(cell_type).n_vertices() *
dofs_per_vertex)
ReferenceCell::internal::Info::get_cell(cell_type).n_faces() :
0)) *
dofs_per_quad +
- (dim == 3 ? 1 : 0) * dofs_per_hex)
+ (dim == 3 ? 1 : 0) * dofs_per_hex + non_local_dofs_per_cell)
, components(n_components)
, degree(degree)
, conforming_space(conformity)
BlockIndices(1, dofs_per_cell) :
block_indices)
{
- Assert(dofs_per_object.size() == dim + 1,
- ExcDimensionMismatch(dofs_per_object.size() - 1, dim));
+ Assert(dofs_per_object.size() == dim + 1 || dofs_per_object.size() == dim + 2,
+ ExcMessage("dofs_per_object should have size of either " +
+ std::to_string(dim + 1) + " or " +
+ std::to_string(dim + 2)));
}
(dofs_per_line == f.dofs_per_line) &&
(dofs_per_quad == f.dofs_per_quad) &&
(dofs_per_hex == f.dofs_per_hex) && (components == f.components) &&
- (degree == f.degree) && (conforming_space == f.conforming_space));
+ (degree == f.degree) && (conforming_space == f.conforming_space) &&
+ (non_local_dofs_per_cell == f.non_local_dofs_per_cell));
}
projects / dealii.git / commitdiff