From 4dfe2df05e97d5c118691e278c9f7e427b071e9f Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Thu, 9 Jul 2020 17:28:27 +0200 Subject: [PATCH] Non local dofs implementation. --- include/deal.II/dofs/dof_accessor.h | 8 + include/deal.II/dofs/dof_accessor.templates.h | 175 +++++++++++++++++- include/deal.II/dofs/dof_levels.h | 6 + include/deal.II/fe/fe_base.h | 62 ++++++- include/deal.II/fe/fe_tools.templates.h | 170 +++++++++++++++-- source/dofs/dof_accessor.cc | 17 ++ source/fe/fe_data.cc | 15 +- 7 files changed, 418 insertions(+), 35 deletions(-) diff --git a/include/deal.II/dofs/dof_accessor.h b/include/deal.II/dofs/dof_accessor.h index 997e5befff..d258c02890 100644 --- a/include/deal.II/dofs/dof_accessor.h +++ b/include/deal.II/dofs/dof_accessor.h @@ -1939,6 +1939,14 @@ public: void set_dof_indices(const std::vector &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 &non_local_dof_indices); + /** * Set the Level DoF indices of this cell to the given values. */ diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h index c97f88f247..9f2d6dd98c 100644 --- a/include/deal.II/dofs/dof_accessor.templates.h +++ b/include/deal.II/dofs/dof_accessor.templates.h @@ -925,7 +925,8 @@ namespace internal { 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; @@ -983,6 +984,12 @@ namespace internal 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); } @@ -1682,11 +1689,51 @@ DoFAccessor::get_dof_indices( "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() * @@ -2205,6 +2252,126 @@ namespace internal + /** + * Implement setting dof indices on a cell. TO CHECK ZHAOWEI + LH + */ + template + static void + set_dof_indices( + const DoFCellAccessor &accessor, + const std::vector &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::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::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::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 + static void + set_non_local_dof_indices( + const DoFCellAccessor &accessor, + const std::vector &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( + 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. diff --git a/include/deal.II/dofs/dof_levels.h b/include/deal.II/dofs/dof_levels.h index 4512e520e1..8170dc23fd 100644 --- a/include/deal.II/dofs/dof_levels.h +++ b/include/deal.II/dofs/dof_levels.h @@ -80,6 +80,12 @@ namespace internal */ DoFObjects dof_object; + /** + * The offsets for each cell of the cache that holds all DoF indices. + * Only used for hp elements. + */ + std::vector cell_cache_offsets; + /** * Return a pointer to the beginning of the DoF indices cache for a * given cell. diff --git a/include/deal.II/fe/fe_base.h b/include/deal.II/fe/fe_base.h index 3b69a65a82..d89a618c89 100644 --- a/include/deal.II/fe/fe_base.h +++ b/include/deal.II/fe/fe_base.h @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// 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. // @@ -249,6 +249,23 @@ public: */ 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: + * - Isogeometric + * Analysis + * - Extended + * finite element method + * - Catmull-Clark's + * subdivision surfaces + */ + const unsigned int non_local_dofs_per_cell; + /** * First index of dof on a line. */ @@ -321,12 +338,16 @@ public: * 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 (1,0,0). On the other hand, for a $Q_3$ - * element in 3d, it would have entries (1,2,4,8). + * 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 (1,0,0). On the other hand, for a + * $Q_3$ element in 3d, it would have entries (1,2,4,8). * * @param[in] n_components Number of vector components of the element. * @@ -415,6 +436,24 @@ public: 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: + * - Isogeometric + * Analysis + * - Extended + * finite element method + * - Catmull-Clark's + * subdivision surfaces + */ + 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 @@ -621,6 +660,15 @@ FiniteElementData::n_dofs_per_cell() const +template +inline unsigned int +FiniteElementData::n_non_local_dofs_per_cell() const +{ + return non_local_dofs_per_cell; +} + + + template template inline unsigned int diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h index 34c7d92789..046cd7d6c5 100644 --- a/include/deal.II/fe/fe_tools.templates.h +++ b/include/deal.II/fe/fe_tools.templates.h @@ -85,10 +85,11 @@ namespace FETools { 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; @@ -114,6 +115,8 @@ namespace FETools 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]; @@ -152,6 +155,8 @@ namespace FETools 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); @@ -319,7 +324,36 @@ namespace FETools 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; @@ -378,8 +412,8 @@ namespace FETools { 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 @@ -413,15 +447,15 @@ namespace FETools std::vector(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( @@ -561,7 +595,54 @@ namespace FETools } } } + 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 > into a vector. @@ -843,6 +924,57 @@ namespace FETools 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; + } + } + } + } } @@ -885,12 +1017,12 @@ namespace FETools { // 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); diff --git a/source/dofs/dof_accessor.cc b/source/dofs/dof_accessor.cc index 579536a127..4f649e9906 100644 --- a/source/dofs/dof_accessor.cc +++ b/source/dofs/dof_accessor.cc @@ -138,6 +138,23 @@ DoFCellAccessor::set_dof_indices( +template +void +DoFCellAccessor::set_non_local_dof_indices( + const std::vector &local_non_local_dof_indices) +{ + Assert(static_cast(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 TriaIterator> DoFCellAccessor::neighbor_child_on_subface( diff --git a/source/fe/fe_data.cc b/source/fe/fe_data.cc index f9fff482ff..81ebd0bb65 100644 --- a/source/fe/fe_data.cc +++ b/source/fe/fe_data.cc @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// 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. // @@ -51,6 +51,8 @@ FiniteElementData::FiniteElementData( , 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) @@ -94,7 +96,7 @@ FiniteElementData::FiniteElementData( 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) @@ -102,8 +104,10 @@ FiniteElementData::FiniteElementData( 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))); } @@ -116,7 +120,8 @@ FiniteElementData::operator==(const FiniteElementData &f) const (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)); } -- 2.39.5