diff --git a/include/mesh/boundary_info.h b/include/mesh/boundary_info.h index 8eed8dbdc40..e29a6165209 100644 --- a/include/mesh/boundary_info.h +++ b/include/mesh/boundary_info.h @@ -480,6 +480,13 @@ class BoundaryInfo : public ParallelObject unsigned int n_boundary_ids (const Elem * const elem, const unsigned short int side) const; + /** + * \returns The number of raw (excludes ancestors) boundary ids associated with the \p side + * side of element \p elem. + */ + unsigned int n_raw_boundary_ids (const Elem * const elem, + const unsigned short int side) const; + /** * \returns The list of boundary ids associated with the \p side side of * element \p elem. @@ -550,6 +557,21 @@ class BoundaryInfo : public ParallelObject */ void build_shellface_boundary_ids(std::vector & b_ids) const; +#ifdef LIBMESH_ENABLE_AMR + /** + * Update parent's boundary id list so that this information is consistent with + * its children + * + * This is useful when `_children_on_boundary = true`, and is used when the + * element is about to get coarsened i.e., in MeshRefinement::_coarsen_elements() + * + * Specifically, when we coarsen an element whose children have different boundary ids. + * In such scenarios, the parent will inherit the children's boundaries if at + * least 50% them own a boundary while sharing the side of the parent. + */ + void transfer_boundary_ids_from_children(const Elem * const parent); +#endif + /** * \returns The number of element-side-based boundary conditions. * @@ -877,6 +899,18 @@ class BoundaryInfo : public ParallelObject const std::multimap> & get_sideset_map() const { return _boundary_side_id; } + /** + * \returns Whether or not there may be child elements directly assigned boundary sides + */ + bool is_children_on_boundary_side() const + { return _children_on_boundary; } + + /** + * Whether or not to allow directly setting boundary sides on child elements + */ + void allow_children_on_boundary_side(const bool children_on_boundary) + { _children_on_boundary = children_on_boundary; } + private: /** @@ -927,6 +961,13 @@ class BoundaryInfo : public ParallelObject std::pair> _boundary_side_id; + /* + * Whether or not children elements are associated with any boundary + * It is false by default. The flag will be turned on if `add_side` + * function is called with a child element + */ + bool _children_on_boundary; + /** * A collection of user-specified boundary ids for sides, edges, nodes, * and shell faces. diff --git a/include/mesh/mesh_refinement.h b/include/mesh/mesh_refinement.h index 6050cfa9c65..c754b31e101 100644 --- a/include/mesh/mesh_refinement.h +++ b/include/mesh/mesh_refinement.h @@ -241,6 +241,10 @@ class MeshRefinement : public ParallelObject * * \note This function used to take an argument, \p maintain_level_one, * new code should use face_level_mismatch_limit() instead. + * + * \note When we allow boundaries to be directly associated with child elements, + * i.e., `_children_on_boundary = true`. A child's boundary ID may be + * lost during coarsening if it differs from its siblings on that parent side. */ bool coarsen_elements (); diff --git a/src/mesh/boundary_info.C b/src/mesh/boundary_info.C index 7eade125404..18bbf79da1a 100644 --- a/src/mesh/boundary_info.C +++ b/src/mesh/boundary_info.C @@ -101,7 +101,8 @@ const boundary_id_type BoundaryInfo::invalid_id = -123; // BoundaryInfo functions BoundaryInfo::BoundaryInfo(MeshBase & m) : ParallelObject(m.comm()), - _mesh (&m) + _mesh (&m), + _children_on_boundary(false) { } @@ -961,10 +962,6 @@ void BoundaryInfo::add_side(const Elem * elem, const boundary_id_type id) { libmesh_assert(elem); - - // Only add BCs for level-0 elements. - libmesh_assert_equal_to (elem->level(), 0); - libmesh_error_msg_if(id == invalid_id, "ERROR: You may not set a boundary ID of " << invalid_id << "\n That is reserved for internal use."); @@ -975,6 +972,26 @@ void BoundaryInfo::add_side(const Elem * elem, pr.second.second == id) return; +#ifdef LIBMESH_ENABLE_AMR + // Users try to mark boundary on child elements + // If this happens, we will allow users to remove + // side from child elements as well + if (elem->level()) + { + _children_on_boundary = true; + + // Here we have to stop and check if we already have this boundary defined on the + // parent (if yes, no need to add) + std::vector bd_ids; + this->boundary_ids(elem,side,bd_ids); + + if(std::find(bd_ids.begin(), bd_ids.end(), id) != bd_ids.end()) + libmesh_not_implemented_msg("Trying to add boundary ID " + + std::to_string(id) + + " which already exists on the ancestors."); + } +#endif + _boundary_side_id.emplace(elem, std::make_pair(side, id)); _boundary_ids.insert(id); _side_boundary_ids.insert(id); // Also add this ID to the set of side boundary IDs @@ -991,8 +1008,26 @@ void BoundaryInfo::add_side(const Elem * elem, libmesh_assert(elem); - // Only add BCs for level-0 elements. - libmesh_assert_equal_to (elem->level(), 0); +#ifdef LIBMESH_ENABLE_AMR + // Users try to mark boundary on child elements + // If this happens, we will allow users to remove + // side from child elements as well + if (elem->level()) + { + _children_on_boundary = true; + + // Here we have to stop and check if we already have this boundary defined on the + // parent (if yes, no need to add) + std::vector bd_ids; + this->boundary_ids(elem,side,bd_ids); + + for (const auto id : ids) + if(std::find(bd_ids.begin(), bd_ids.end(), id) != bd_ids.end()) + libmesh_not_implemented_msg("Trying to add boundary ID " + + std::to_string(id) + + " which already exists on the ancestors."); + } +#endif // Don't add the same ID twice auto bounds = _boundary_side_id.equal_range(elem); @@ -1240,25 +1275,61 @@ void BoundaryInfo::boundary_ids (const Elem * const elem, // Clear out any previous contents vec_to_fill.clear(); - // Only level-0 elements store BCs. If this is not a level-0 - // element get its level-0 parent and infer the BCs. + // In most cases only level-0 elements store BCs. + // In certain application (such as time-dependent domains), however, children + // need to store BCs too. This case is covered with the _children_on_boundary + // flag. const Elem * searched_elem = elem; + +#ifdef LIBMESH_ENABLE_AMR + if (elem->level() != 0) + { + // If we have children on the boundaries, we need to search for boundary IDs on the + // child and its ancestors too if they share the side. + if (_children_on_boundary) { - if (elem->neighbor_ptr(side) == nullptr) - searched_elem = elem->top_parent (); -#ifdef LIBMESH_ENABLE_AMR - else - while (searched_elem->parent() != nullptr) - { - const Elem * parent = searched_elem->parent(); - if (parent->is_child_on_side(parent->which_child_am_i(searched_elem), side) == false) - return; - searched_elem = parent; - } -#endif + // Loop over ancestors to check if they have boundary ids on the same side + while (searched_elem) + { + for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem))) + // Here we need to check if the boundary id already exists + if (pr.second.first == side && + std::find(vec_to_fill.begin(), vec_to_fill.end(), pr.second.second) == + vec_to_fill.end()) + vec_to_fill.push_back(pr.second.second); + + + const Elem * parent = searched_elem->parent(); + // If the parent doesn't exist or if the child is not on the correct side of the + // parent we are done checking the ancestors + if (!parent || parent->is_child_on_side(parent->which_child_am_i(searched_elem), side) == false) + return; + + searched_elem = parent; + } + + return; } + // If we don't have children on boundaries and we are on an external boundary, + // we just look for the top parent + if (elem->neighbor_ptr(side) == nullptr) + searched_elem = elem->top_parent(); + // Otherwise we loop over the ancestors and check if they have a different BC for us + else + while (searched_elem->parent() != nullptr) + { + const Elem * parent = searched_elem->parent(); + if (parent->is_child_on_side(parent->which_child_am_i(searched_elem), side) == false) + return; + + searched_elem = parent; + } + } + +#endif + // Check each element in the range to see if its side matches the requested side. for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem))) if (pr.second.first == side) @@ -1277,6 +1348,15 @@ unsigned int BoundaryInfo::n_boundary_ids (const Elem * const elem, } +unsigned int BoundaryInfo::n_raw_boundary_ids (const Elem * const elem, + const unsigned short int side) const +{ + std::vector ids; + this->raw_boundary_ids(elem, side, ids); + return cast_int(ids.size()); +} + + void BoundaryInfo::raw_boundary_ids (const Elem * const elem, const unsigned short int side, @@ -1288,7 +1368,7 @@ void BoundaryInfo::raw_boundary_ids (const Elem * const elem, vec_to_fill.clear(); // Only level-0 elements store BCs. - if (elem->parent()) + if (elem->parent() && !_children_on_boundary) return; // Check each element in the range to see if its side matches the requested side. @@ -1438,9 +1518,6 @@ void BoundaryInfo::remove_side (const Elem * elem, { libmesh_assert(elem); - // Only level 0 elements are stored in BoundaryInfo. - libmesh_assert_equal_to (elem->level(), 0); - // Erase (elem, side, *) entries from map. erase_if(_boundary_side_id, elem, [side](decltype(_boundary_side_id)::mapped_type & pr) @@ -1455,6 +1532,25 @@ void BoundaryInfo::remove_side (const Elem * elem, { libmesh_assert(elem); +#ifdef LIBMESH_ENABLE_AMR + // Here we have to stop and check if somebody tries to remove an ancestor's boundary ID + // through a child + if (elem->level()) + { + std::vector bd_ids; + this->boundary_ids(elem,side,bd_ids); + if(std::find(bd_ids.begin(), bd_ids.end(), id) != bd_ids.end()) + { + std::vector raw_bd_ids; + this->raw_boundary_ids(elem, side, raw_bd_ids); + if(std::find(raw_bd_ids.begin(), raw_bd_ids.end(), id) == raw_bd_ids.end()) + libmesh_not_implemented_msg("We cannot delete boundary ID " + + std::to_string(id) + + " using a child because it is inherited from an ancestor."); + } + } +#endif + // Erase (elem, side, id) entries from map. erase_if(_boundary_side_id, elem, [side, id](decltype(_boundary_side_id)::mapped_type & pr) @@ -1576,18 +1672,25 @@ unsigned int BoundaryInfo::side_with_boundary_id(const Elem * const elem, const boundary_id_type boundary_id_in) const { const Elem * searched_elem = elem; - if (elem->level() != 0) - searched_elem = elem->top_parent(); + + // If we don't have a time-dependent domain, we can just go ahead and use the top parent + // (since only those contain boundary conditions). Otherwise, we keep the element + if (elem->level() != 0 && !_children_on_boundary) + searched_elem = elem->top_parent(); // elem may have zero or multiple occurrences for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem))) - { + { // if this is true we found the requested boundary_id // of the element and want to return the side if (pr.second.second == boundary_id_in) - { - unsigned int side = pr.second.first; + { + unsigned int side = pr.second.first; + // Here we branch out. If we don't allow time-dependent boundary domains, + // we need to check if our parents are consistent. + if (!_children_on_boundary) + { // If we're on this external boundary then we share this // external boundary id if (elem->neighbor_ptr(side) == nullptr) @@ -1600,18 +1703,56 @@ unsigned int BoundaryInfo::side_with_boundary_id(const Elem * const elem, #ifdef LIBMESH_ENABLE_AMR while (p != nullptr) - { - const Elem * parent = p->parent(); - if (parent && !parent->is_child_on_side(parent->which_child_am_i(p), side)) - break; - p = parent; - } + { + const Elem * parent = p->parent(); + if (parent && !parent->is_child_on_side(parent->which_child_am_i(p), side)) + break; + p = parent; + } #endif // We're on that side of our top_parent; return it if (!p) return side; } + // Otherwise we need to check if the child's ancestors have something on + // the side of the child + else + return side; + } + } + +#ifdef LIBMESH_ENABLE_AMR + // We might have instances (especially with moving boundary domains) when we + // query the paren't boundary ID on a child. We only do this till we find the + // the first side, for multiple sides see above. + if (_children_on_boundary && elem->level() != 0) + { + for (auto side : make_range(elem->n_sides())) + { + const Elem * p = elem; + while (p->parent() != nullptr) + { + const Elem * parent = p->parent(); + + // First we make sure the parent shares this side + if (parent->is_child_on_side(parent->which_child_am_i(p), side)) + { + // parent may have multiple boundary ids + for (const auto & pr : as_range(_boundary_side_id.equal_range(parent))) + // if this is true we found the requested boundary_id + // of the element and want to return the side + if (pr.second.first == side && pr.second.second == boundary_id_in) + return side; + + p = parent; + } + // If the parent is not on the same side, other ancestors won't be on the same side either + else + break; + } } + } +#endif // if we get here, we found elem in the data structure but not // the requested boundary id, so return the default value @@ -1626,18 +1767,22 @@ BoundaryInfo::sides_with_boundary_id(const Elem * const elem, std::vector returnval; const Elem * searched_elem = elem; - if (elem->level() != 0) + if (elem->level() != 0 && !_children_on_boundary) searched_elem = elem->top_parent(); // elem may have zero or multiple occurrences for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem))) - { + { // if this is true we found the requested boundary_id // of the element and want to return the side if (pr.second.second == boundary_id_in) - { - unsigned int side = pr.second.first; + { + unsigned int side = pr.second.first; + // Here we branch out. If we don't allow time-dependent boundary domains, + // we need to check if our parents are consistent. + if (!_children_on_boundary) + { // If we're on this external boundary then we share this // external boundary id if (elem->neighbor_ptr(side) == nullptr) @@ -1653,23 +1798,61 @@ BoundaryInfo::sides_with_boundary_id(const Elem * const elem, #ifdef LIBMESH_ENABLE_AMR while (p != nullptr) - { - const Elem * parent = p->parent(); - if (parent && !parent->is_child_on_side(parent->which_child_am_i(p), side)) - break; - p = parent; - } + { + const Elem * parent = p->parent(); + if (parent && !parent->is_child_on_side(parent->which_child_am_i(p), side)) + break; + p = parent; + } #endif // We're on that side of our top_parent; return it if (!p) returnval.push_back(side); } + // Otherwise we trust what we got and return the side + else + returnval.push_back(side); + } + } + +#ifdef LIBMESH_ENABLE_AMR + // We might have instances (especially with moving boundary domains) when we + // query the parent boundary ID on a child. + if (_children_on_boundary && elem->level() != 0) + { + for (auto side : make_range(elem->n_sides())) + { + const Elem * p = elem; + while (p->parent() != nullptr) + { + const Elem * parent = p->parent(); + // First we make sure the parent shares this side + if (parent->is_child_on_side(parent->which_child_am_i(p), side)) + { + // parent may have multiple boundary ids + for (const auto & pr : as_range(_boundary_side_id.equal_range(parent))) + { + // if this is true we found the requested boundary_id + // of the element and want to add the side to the vector. We + // also need to check if the side is already in the vector. This might + // happen if the child inherits the boundary from the parent. + if (pr.second.first == side && pr.second.second == boundary_id_in && + std::find(returnval.begin(), returnval.end(), side) == returnval.end()) + returnval.push_back(side); + } + } + // If the parent is not on the same side, other ancestors won't be on the same side either + else + break; + p = parent; + } } + } +#endif return returnval; } - void BoundaryInfo::build_node_boundary_ids(std::vector & b_ids) const { @@ -1712,6 +1895,57 @@ BoundaryInfo::build_shellface_boundary_ids(std::vector & b_ids } } +#ifdef LIBMESH_ENABLE_AMR +void +BoundaryInfo::transfer_boundary_ids_from_children(const Elem * const parent) +{ + // this is only needed when we allow boundary to be associated with children elements + // also, we only transfer the parent's boundary ids when we are actually coarsen the child element + if (!_children_on_boundary || + !(!parent->active() && parent->refinement_flag() == Elem::COARSEN_INACTIVE)) + return; + + // We assume that edges can be divided ito two pieces, while triangles and + // quads can be divided into four smaller areas. This is double because we'll need + // to convert the ratio of the children with given boundary id to a double. + const double number_of_sides_on_children = std::pow(2, parent->dim()-1); + + // In this case the input argument elem is the parent element. We need to check all of its sides + // to grab any potential boundary ids. + for (unsigned int side_i = 0; side_i < parent->n_sides(); ++side_i) + { + // An temporary storage to count how many times the children's boundaries occur. the general + // consensus is that if the boundary occurs more than once we propagate upon coarsening. Otherwise, + // it will get deleted. + std::map boundary_counts; + + for (const auto & child_i : make_range(parent->n_children())) + { + // We only need to check the children which share the side + if (parent->is_child_on_side(child_i, side_i)) + { + // Fetching the boundary tags on the child's side + for (const auto & pr : as_range(_boundary_side_id.equal_range(parent->child_ptr(child_i)))) + { + // Making sure we are on the same boundary + if (pr.second.first == side_i) + ++boundary_counts[pr.second.second]; + } + } + } + + // This is where the decision is made. If 50% of the children have the tags, + // we propagate them upwards upon coarsening. Otherwise, they are deleted. + for (const auto & boundary : boundary_counts) + if (boundary.second / number_of_sides_on_children > 0.5) + this->add_side(parent, side_i, boundary.first); + } + + for (const auto & child_i : make_range(parent->n_children())) + this->remove(parent->child_ptr(child_i)); +} +#endif + std::size_t BoundaryInfo::n_boundary_conds () const { // in serial we know the number of bcs from the diff --git a/src/mesh/mesh_refinement.C b/src/mesh/mesh_refinement.C index 71d755ef61f..eb276efc478 100644 --- a/src/mesh/mesh_refinement.C +++ b/src/mesh/mesh_refinement.C @@ -705,8 +705,6 @@ bool MeshRefinement::refine_elements () elem->set_refinement_flag(Elem::DO_NOTHING); } - - // Parallel consistency has to come first, or coarsening // along processor boundaries might occasionally be falsely // prevented @@ -735,6 +733,7 @@ bool MeshRefinement::refine_elements () if (_face_level_mismatch_limit) libmesh_assert(test_level_one(true)); libmesh_assert(this->make_refinement_compatible()); + // FIXME: This won't pass unless we add a redundant find_neighbors() // call or replace find_neighbors() with on-the-fly neighbor updating // libmesh_assert(!this->eliminate_unrefined_patches()); @@ -1371,6 +1370,10 @@ bool MeshRefinement::_coarsen_elements () for (auto & elem : _mesh.element_ptr_range()) { + // Make sure we transfer the children's boundary id(s) + // up to its parent when necessary before coarsening. + _mesh.get_boundary_info().transfer_boundary_ids_from_children(elem); + // active elements flagged for coarsening will // no longer be deleted until MeshRefinement::contract() if (elem->refinement_flag() == Elem::COARSEN) @@ -1384,8 +1387,11 @@ bool MeshRefinement::_coarsen_elements () elem->nullify_neighbors(); // Remove any boundary information associated - // with this element - _mesh.get_boundary_info().remove (elem); + // with this element if we do not allow children to have boundary info. + // Otherwise, we will do the removal in `transfer_boundary_ids_from_children` + // to make sure we don't delete the information before it is transferred + if (!_mesh.get_boundary_info().is_children_on_boundary_side()) + _mesh.get_boundary_info().remove (elem); // Add this iterator to the _unused_elements // data structure so we might fill it. diff --git a/src/parallel/parallel_elem.C b/src/parallel/parallel_elem.C index 64706b880fd..3afd0375fc0 100644 --- a/src/parallel/parallel_elem.C +++ b/src/parallel/parallel_elem.C @@ -88,23 +88,37 @@ Packing::packed_size (std::vector::const_iterator const unsigned int indexing_size = DofObject::unpackable_indexing_size(in+pre_indexing_size); - unsigned int total_packed_bc_data = 0; - if (level == 0) + // We communicate if we are on the boundary or not + unsigned int total_packed_bc_data = 1; + largest_id_type on_boundary = *(in + pre_indexing_size + indexing_size); + + if (on_boundary) + { + // Extracting if the children are allowed on the boundary + total_packed_bc_data++; + largest_id_type allow_children_on_boundary = *(in + pre_indexing_size + indexing_size + 1); + + // For now, children are only supported on sides, the nodes and shell faces are + // treated using the top parents only + if (level == 0 || allow_children_on_boundary) { for (unsigned int s = 0; s != n_sides; ++s) - { - const int n_bcs = cast_int - (*(in + pre_indexing_size + indexing_size + - total_packed_bc_data++)); - libmesh_assert_greater_equal (n_bcs, 0); - total_packed_bc_data += n_bcs; - } - + { + const int n_bcs = cast_int + (*(in + pre_indexing_size + indexing_size + + total_packed_bc_data++)); + libmesh_assert_greater_equal (n_bcs, 0); + total_packed_bc_data += n_bcs; + } + } + } + if (level == 0) + { for (unsigned int e = 0; e != n_edges; ++e) { const int n_bcs = cast_int (*(in + pre_indexing_size + indexing_size + - total_packed_bc_data++)); + total_packed_bc_data++)); libmesh_assert_greater_equal (n_bcs, 0); total_packed_bc_data += n_bcs; } @@ -113,7 +127,7 @@ Packing::packed_size (std::vector::const_iterator { const int n_bcs = cast_int (*(in + pre_indexing_size + indexing_size + - total_packed_bc_data++)); + total_packed_bc_data++)); libmesh_assert_greater_equal (n_bcs, 0); total_packed_bc_data += n_bcs; } @@ -142,16 +156,35 @@ unsigned int Packing::packable_size (const Elem * const & elem, const MeshBase * mesh) { - unsigned int total_packed_bcs = 0; + // We always communicate if we are on a boundary or not + unsigned int total_packed_bcs = 1; const unsigned short n_sides = elem->n_sides(); - if (elem->level() == 0) + largest_id_type on_boundary = 0; + for (auto s : elem->side_index_range()) + if (mesh->get_boundary_info().n_raw_boundary_ids(elem,s)) + { + on_boundary = 1; + break; + } + + if (on_boundary) + { + // In this case we need another entry to check if we allow children on the boundary. + // We only allow children on sides, edges and sheel faces are treated normally using + // their top parents. + total_packed_bcs++; + if (elem->level() == 0 || mesh->get_boundary_info().is_children_on_boundary_side()) { total_packed_bcs += n_sides; for (unsigned short s = 0; s != n_sides; ++s) total_packed_bcs += - mesh->get_boundary_info().n_boundary_ids(elem,s); + mesh->get_boundary_info().n_raw_boundary_ids(elem,s); + } + } + if (elem->level() == 0) + { const unsigned short n_edges = elem->n_edges(); total_packed_bcs += n_edges; for (unsigned short e = 0; e != n_edges; ++e) @@ -288,41 +321,64 @@ Packing::pack (const Elem * const & elem, // Add any DofObject indices elem->pack_indexing(data_out); - // If this is a coarse element, - // Add any element side boundary condition ids - if (elem->level() == 0) + // We check if this is a boundary cell. We use the raw + // IDs because we also communicate the parents which + // will bring their associated IDs + largest_id_type on_boundary = 0; + for (auto s : elem->side_index_range()) + if (mesh->get_boundary_info().n_raw_boundary_ids(elem,s)) + { + on_boundary = 1; + break; + } + + *data_out++ = on_boundary; + + if (on_boundary) + { + *data_out++ = mesh->get_boundary_info().is_children_on_boundary_side(); + // Again, only do this if we allow children to hold boundary sides, the edges and + // shell faces are treated normally using their top parents + if (elem->level() == 0 || mesh->get_boundary_info().is_children_on_boundary_side()) { std::vector bcs; for (auto s : elem->side_index_range()) - { - mesh->get_boundary_info().boundary_ids(elem, s, bcs); + { + mesh->get_boundary_info().raw_boundary_ids(elem, s, bcs); - *data_out++ =(bcs.size()); + *data_out++ =(bcs.size()); - for (const auto & bid : bcs) - *data_out++ = bid; - } + for (const auto & bid : bcs) + *data_out++ = bid; + } + } + } - for (auto e : elem->edge_index_range()) - { - mesh->get_boundary_info().edge_boundary_ids(elem, e, bcs); + // If this is a coarse element, + // Add any element side boundary condition ids + if (elem->level() == 0) + { + std::vector bcs; + for (auto e : elem->edge_index_range()) + { + mesh->get_boundary_info().edge_boundary_ids(elem, e, bcs); - *data_out++ =(bcs.size()); + *data_out++ =(bcs.size()); - for (const auto & bid : bcs) - *data_out++ = bid; - } + for (const auto & bid : bcs) + *data_out++ = bid; + } - for (unsigned short sf=0; sf != 2; ++sf) - { - mesh->get_boundary_info().shellface_boundary_ids(elem, sf, bcs); + for (unsigned short sf=0; sf != 2; ++sf) + { + mesh->get_boundary_info().shellface_boundary_ids(elem, sf, bcs); - *data_out++ =(bcs.size()); + *data_out++ =(bcs.size()); - for (const auto & bid : bcs) - *data_out++ = bid; - } - } + for (const auto & bid : bcs) + *data_out++ = bid; + } + } } @@ -611,6 +667,10 @@ Packing::unpack (std::vector::const_iterator in, libmesh_assert(neigh->level() < elem->level() || neigh->neighbor_ptr(neighbor_side) == elem); } + else + // We skip these to go to the boundary information if the element is + // actually subactive + in += 2*elem->n_sides(); // Our p level and refinement flags should be "close to" correct // if we're not an active element - we might have a p level @@ -776,9 +836,15 @@ Packing::unpack (std::vector::const_iterator in, in += elem->packed_indexing_size(); - // If this is a coarse element, - // add any element side or edge boundary condition ids - if (level == 0) + // We check if this is cell holds a boundary ID or not + auto on_boundary = *in++; + if (on_boundary) + { + // Only treat the sides with caution. This is because we might hold boundary IDs + // on the sides of the children. This is not supported for edges and shell faces, thus + // they are treated assuming that only top parents can hold the IDs. + auto children_on_boundary = *in++; + if (elem->level() == 0 || children_on_boundary) { for (auto s : elem->side_index_range()) { @@ -789,7 +855,13 @@ Packing::unpack (std::vector::const_iterator in, mesh->get_boundary_info().add_side (elem, s, cast_int(*in++)); } + } + } + // If this is a coarse element, + // add any element side or edge boundary condition ids + if (level == 0) + { for (auto e : elem->edge_index_range()) { const boundary_id_type num_bcs = diff --git a/src/partitioning/partitioner.C b/src/partitioning/partitioner.C index f136ee398ec..485ff4dae44 100644 --- a/src/partitioning/partitioner.C +++ b/src/partitioning/partitioner.C @@ -28,6 +28,7 @@ #include "libmesh/mesh_communication.h" #include "libmesh/parallel_ghost_sync.h" #include "libmesh/wrapped_petsc.h" +#include "libmesh/boundary_info.h" // Subclasses to build() #include "libmesh/enum_partitioner_type.h" diff --git a/src/solvers/petsc_diff_solver.C b/src/solvers/petsc_diff_solver.C index fa330db8fa7..32728a0fafd 100644 --- a/src/solvers/petsc_diff_solver.C +++ b/src/solvers/petsc_diff_solver.C @@ -23,6 +23,7 @@ #include "libmesh/petsc_matrix.h" #include "libmesh/petsc_vector.h" #include "libmesh/petsc_auto_fieldsplit.h" +#include "libmesh/boundary_info.h" #ifdef LIBMESH_HAVE_PETSC @@ -285,6 +286,14 @@ unsigned int PetscDiffSolver::solve() { LOG_SCOPE("solve()", "PetscDiffSolver"); +#if !PETSC_VERSION_LESS_THAN(3,7,3) +#if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_METAPHYSICL) + // GMG is currently not supported if we enable children to be associated with + // boundary sides + libmesh_assert(!_system.get_mesh().get_boundary_info().is_children_on_boundary_side()); +#endif +#endif + PetscVector & x = *(cast_ptr *>(_system.solution.get())); PetscMatrix & jac = diff --git a/tests/mesh/boundary_info.C b/tests/mesh/boundary_info.C index d0780485bb8..b37d0aa3146 100644 --- a/tests/mesh/boundary_info.C +++ b/tests/mesh/boundary_info.C @@ -8,10 +8,12 @@ #include #include #include +#include #include "test_comm.h" #include "libmesh_cppunit.h" +#include using namespace libMesh; @@ -27,6 +29,14 @@ public: #if LIBMESH_DIM > 1 CPPUNIT_TEST( testMesh ); CPPUNIT_TEST( testRenumber ); +# ifdef LIBMESH_ENABLE_AMR +# ifdef LIBMESH_ENABLE_EXCEPTIONS + CPPUNIT_TEST( testBoundaryOnChildrenErrors ); +# endif + CPPUNIT_TEST( testBoundaryOnChildrenElementsRefineCoarsen ); + CPPUNIT_TEST( testBoundaryOnChildrenBoundaryIDs ); + CPPUNIT_TEST( testBoundaryOnChildrenBoundarySides ); +# endif # ifdef LIBMESH_ENABLE_DIRICHLET CPPUNIT_TEST( testShellFaceConstraints ); # endif @@ -474,6 +484,360 @@ public: } #endif // LIBMESH_ENABLE_DIRICHLET +#if LIBMESH_ENABLE_AMR +# if LIBMESH_ENABLE_EXCEPTIONS + void testBoundaryOnChildrenErrors() + { + LOG_UNIT_TEST; + + // We create one cell only. The default boundaries of the cell are below. + // ___2___ + // 3 | | 1 + // |_____| + // 0 + + auto mesh = std::make_unique(*TestCommWorld); + MeshTools::Generation::build_square(*mesh, + 1, 1, + 0., 1., + 0., 1., + QUAD4); + + BoundaryInfo & bi = mesh->get_boundary_info(); + + // We only have one element, but for easy access we use the iterator + for (auto & elem : mesh->active_element_ptr_range()) + elem->set_refinement_flag(Elem::REFINE); + mesh->prepare_for_use(); + + MeshRefinement(*mesh).refine_elements(); + mesh->prepare_for_use(); + + // Now we try to add boundary id 3 to a child on side 3. This should + // result in a "not implemented" error message + bool threw_desired_exception = false; + try { + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 0.5 && c(1) > 0.5) + bi.add_side(elem, 3, 3); + } + } + catch (libMesh::NotImplemented & e) { + std::regex msg_regex("Trying to add boundary ID 3 which already exists on the ancestors"); + CPPUNIT_ASSERT(std::regex_search(e.what(), msg_regex)); + threw_desired_exception = true; + } + // If we have more than 4 processors, or a poor partitioner, we + // might not get an exception on every processor + mesh->comm().max(threw_desired_exception); + + CPPUNIT_ASSERT(threw_desired_exception); + + threw_desired_exception = false; + try { + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 0.5 && c(1) > 0.5) + bi.add_side(elem, 3, {3,4}); + } + } + catch (libMesh::NotImplemented & e) { + std::regex msg_regex("Trying to add boundary ID 3 which already exists on the ancestors"); + CPPUNIT_ASSERT(std::regex_search(e.what(), msg_regex)); + threw_desired_exception = true; + } + + // If we have more than 4 processors, or a poor partitioner, we + // might not get an exception on every processor + mesh->comm().max(threw_desired_exception); + + CPPUNIT_ASSERT(threw_desired_exception); + + // We tested the side addition errors, now we move to the removal parts. + // We will attempt the removal of boundary 3 through the child + threw_desired_exception = false; + bi.allow_children_on_boundary_side(true); + try { + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 0.5 && c(1) > 0.5) + bi.remove_side(elem, 3, 3); + } + } + catch (libMesh::NotImplemented & e) { + std::regex msg_regex("We cannot delete boundary ID 3 using a child because it is inherited from an ancestor"); + CPPUNIT_ASSERT(std::regex_search(e.what(), msg_regex)); + threw_desired_exception = true; + } + + // If we have more than 4 processors, or a poor partitioner, we + // might not get an exception on every processor + mesh->comm().max(threw_desired_exception); + + CPPUNIT_ASSERT(threw_desired_exception); + } +# endif // LIBMESH_ENABLE_EXCEPTIONS + + void testBoundaryOnChildrenElementsRefineCoarsen() + { + LOG_UNIT_TEST; + + // Set subdomain ids for specific elements, we will refine/coarsen + // the cell on subdomain 1 + // _____________ + // | 1 | 2 | + // |_____|_____| + + auto mesh = std::make_unique(*TestCommWorld); + MeshTools::Generation::build_square(*mesh, + 2, 1, + 0., 2., + 0., 1., + QUAD4); + + BoundaryInfo & bi = mesh->get_boundary_info(); + + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 1) + { + elem->subdomain_id() = 1; + elem->set_refinement_flag(Elem::REFINE); + } + else + elem->subdomain_id() = 2; + } + mesh->prepare_for_use(); + + // Refine the elements once in subdomain 1, and + // add the right side subdomain 1 as boundary 5 + MeshRefinement(*mesh).refine_elements(); + mesh->prepare_for_use(); + + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 1 && c(0) > 0.5) + bi.add_side(elem, 1, 5); + } + mesh->prepare_for_use(); + + // Check the middle boundary, we expect to have two sides in boundary 5 + unsigned int count = 0; + for (auto & elem : mesh->active_element_ptr_range()) + if (bi.has_boundary_id(elem, 1, 5)) + count++; + + if (mesh->n_active_local_elem()) + { + CPPUNIT_ASSERT_EQUAL((unsigned int) 2, count); + CPPUNIT_ASSERT(bi.is_children_on_boundary_side()); + } + + // First, we will coarsen the the elements on subdomain 1. This + // is to check if the boundary information propagates upward upon + // coarsening. + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 1) + elem->set_refinement_flag(Elem::COARSEN); + } + mesh->prepare_for_use(); + + // The coarsened element should have its side on boundary 5 + // This is boundary info transferred from this child element + MeshRefinement(*mesh).coarsen_elements(); + mesh->prepare_for_use(); + + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 1) + { + CPPUNIT_ASSERT(bi.has_boundary_id(elem, 1, 5)); + // We clean up this boundary ID for the next round of tests + bi.remove_side(elem, 1, 5); + // we will refine this element again + elem->set_refinement_flag(Elem::REFINE); + } + } + + MeshRefinement(*mesh).refine_elements(); + mesh->prepare_for_use(); + + // This time we remove boundary 5 from one of the children. We expect + // the boundary not to propagate to the next level. Furthermore we + // expect boundary 5 to be deleted from the parent's boundaries + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 1) + elem->set_refinement_flag(Elem::COARSEN); + if (c(0) > 0.5 && c(0) < 1 && c(1) < 0.5) + bi.add_side(elem, 1, 5); + } + mesh->prepare_for_use(); + + MeshRefinement(*mesh).coarsen_elements(); + + mesh->prepare_for_use(); + + // The parent element should not have any side associated with boundary 5 + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 1) + CPPUNIT_ASSERT(!bi.has_boundary_id(elem, 1, 5)); + } + } + + void testBoundaryOnChildrenBoundaryIDs() + { + LOG_UNIT_TEST; + + // We create one cell only. The default boundaries of the cell are below. + // We will refine the mesh and add a new boundary id to the left side (side 3). + // Then will query the available boundary ids on the added side. It should return + // both the parent's and the child's boundaries. + // ___2___ + // 3 | | 1 + // |_____| + // 0 + + auto mesh = std::make_unique(*TestCommWorld); + MeshTools::Generation::build_square(*mesh, + 1, 1, + 0., 1., + 0., 1., + QUAD4); + + BoundaryInfo & bi = mesh->get_boundary_info(); + + // We only have one element, but for easy access we use the iterator + for (auto & elem : mesh->active_element_ptr_range()) + elem->set_refinement_flag(Elem::REFINE); + mesh->prepare_for_use(); + + MeshRefinement(*mesh).refine_elements(); + + // Now we add the extra boundary ID (5) to the element in the top + // left corner + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 0.5 && c(1) > 0.5) + bi.add_side(elem, 3, 5); + } + mesh->prepare_for_use(); + + // Okay, now we query the boundary ids on side 3 of the child and check if it has + // the right elements + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 0.5 && c(1) > 0.5) + { + std::vector container; + bi.boundary_ids(elem, 3, container); + + CPPUNIT_ASSERT_EQUAL((unsigned long) 2, container.size()); + CPPUNIT_ASSERT_EQUAL((short int) 5, container[0]); + CPPUNIT_ASSERT_EQUAL((short int) 3, container[1]); + } + } + } + + void testBoundaryOnChildrenBoundarySides() + { + LOG_UNIT_TEST; + + // We create one cell only. The default boundaries of the cell are below. + // We will refine mesh and see if we can get back the correct sides + // for a given boundary id on an internal boundary. + // ___2___ + // 3 | | 1 + // |_____| + // 0 + + auto mesh = std::make_unique(*TestCommWorld); + MeshTools::Generation::build_square(*mesh, + 1, 1, + 0., 1., + 0., 1., + QUAD4); + + BoundaryInfo & bi = mesh->get_boundary_info(); + + + // We only have one element, but for easy access we use the iterator + for (auto & elem : mesh->active_element_ptr_range()) + elem->set_refinement_flag(Elem::REFINE); + mesh->prepare_for_use(); + MeshRefinement(*mesh).refine_elements(); + + // Now we add the extra boundary ID (5) to two sides of + // the element in the bottom left corner. then we refine again + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 0.5 && c(1) < 0.5) + { + bi.add_side(elem, 1, 5); + bi.add_side(elem, 2, 5); + elem->set_refinement_flag(Elem::REFINE); + } + } + mesh->prepare_for_use(); + MeshRefinement(*mesh).refine_elements(); + + // Okay, now we add another boundary id (6) to the cell which is in the bottom + // right corner of the refined element + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 0.5 && c(0) > 0.25 && c(1) < 0.25) + bi.add_side(elem, 1, 6); + } + + // Time to test if we can get back the boundary sides, first we + // check if we can get back boundary from the ancestors of (5) on + // the cell which only has boundary (6) registered. We also check + // if we can get boundary (6) back. + + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 0.5 && c(0) > 0.25 && c(1) < 0.25) + { + const auto side_5 = bi.side_with_boundary_id(elem, 5); + const auto side_6 = bi.side_with_boundary_id(elem, 6); + CPPUNIT_ASSERT_EQUAL((unsigned int) 1, side_5); + CPPUNIT_ASSERT_EQUAL((unsigned int) 1, side_6); + } + } + + // Now we go and try to query the sides with boundary id (5) using + // the element which is at the top right corner of the bottom + // right parent. + for (auto & elem : mesh->active_element_ptr_range()) + { + const Point c = elem->vertex_average(); + if (c(0) < 0.5 && c(0) > 0.25 && c(1) > 0.25 && c(1) < 0.5) + { + const auto sides = bi.sides_with_boundary_id(elem, 5); + CPPUNIT_ASSERT_EQUAL((unsigned long) 2, sides.size()); + CPPUNIT_ASSERT_EQUAL((unsigned int) 1, sides[0]); + CPPUNIT_ASSERT_EQUAL((unsigned int) 2, sides[1]); + } + } + } +#endif //LIBMESH_ENABLE_AMR }; CPPUNIT_TEST_SUITE_REGISTRATION( BoundaryInfoTest );