From e3a08cfb4459c4de01ada5fb5fb9c2511e7478a9 Mon Sep 17 00:00:00 2001 From: Yaqi Wang Date: Mon, 25 Nov 2024 21:35:38 -0600 Subject: [PATCH] split lower-d mesh blocks based on side types so that exodus output can work properly #29142 --- .../include/loops/ThreadedElementLoopBase.h | 4 +- framework/include/mesh/MooseMesh.h | 13 ++ framework/include/utils/MooseTypes.h | 2 - .../src/actions/SetupMeshCompleteAction.C | 3 +- framework/src/bcs/ArrayHFEMDirichletBC.C | 10 +- framework/src/bcs/ArrayLowerDIntegratedBC.C | 21 ++-- framework/src/bcs/HFEMDirichletBC.C | 10 +- framework/src/bcs/LowerDIntegratedBC.C | 21 ++-- framework/src/dgkernels/ArrayDGLowerDKernel.C | 21 ++-- framework/src/dgkernels/DGLowerDKernel.C | 21 ++-- framework/src/interfaces/BlockRestrictable.C | 5 +- .../src/loops/ComputeFullJacobianThread.C | 2 +- framework/src/mesh/MooseMesh.C | 86 +++++++++++-- framework/src/problems/DisplacedProblem.C | 4 +- framework/src/problems/FEProblemBase.C | 10 +- framework/src/systems/NonlinearSystemBase.C | 11 +- framework/src/utils/MooseTypes.C | 2 - test/tests/kernels/hfem/3d-lower-d-volumes.i | 12 +- test/tests/kernels/hfem/3d-prism.i | 119 ++++++++++++++++++ test/tests/kernels/hfem/array_dirichlet.i | 6 +- .../kernels/hfem/array_dirichlet_pjfnk.i | 6 +- .../kernels/hfem/array_dirichlet_transform.i | 6 +- .../hfem/array_dirichlet_transform_bc.i | 6 +- test/tests/kernels/hfem/array_neumann.i | 4 +- test/tests/kernels/hfem/array_robin.i | 4 +- test/tests/kernels/hfem/dirichlet.i | 12 +- test/tests/kernels/hfem/gold/3d-prism_out.e | Bin 0 -> 75264 bytes test/tests/kernels/hfem/lower-d-volumes.i | 12 +- test/tests/kernels/hfem/neumann.i | 4 +- test/tests/kernels/hfem/robin.i | 4 +- test/tests/kernels/hfem/robin_adapt.i | 4 +- test/tests/kernels/hfem/robin_displaced.i | 4 +- test/tests/kernels/hfem/robin_dist.i | 4 +- test/tests/kernels/hfem/tests | 8 ++ test/tests/kernels/hfem/variable_dirichlet.i | 8 +- test/tests/kernels/hfem/variable_robin.i | 14 +-- .../side-uo-with-lower-d-use.i | 2 +- .../lower_d_side_boundary.i | 10 +- 38 files changed, 356 insertions(+), 139 deletions(-) create mode 100644 test/tests/kernels/hfem/3d-prism.i create mode 100644 test/tests/kernels/hfem/gold/3d-prism_out.e diff --git a/framework/include/loops/ThreadedElementLoopBase.h b/framework/include/loops/ThreadedElementLoopBase.h index 07d2eaade18e..bfecd634ff30 100644 --- a/framework/include/loops/ThreadedElementLoopBase.h +++ b/framework/include/loops/ThreadedElementLoopBase.h @@ -255,8 +255,8 @@ ThreadedElementLoopBase::operator()(const RangeType & range, bool byp onElement(elem); - if (elem->subdomain_id() == Moose::INTERNAL_SIDE_LOWERD_ID || - elem->subdomain_id() == Moose::BOUNDARY_SIDE_LOWERD_ID) + if (_mesh.interiorLowerDBlocks().count(elem->subdomain_id()) > 0 || + _mesh.boundaryLowerDBlocks().count(elem->subdomain_id()) > 0) { postElement(elem); continue; diff --git a/framework/include/mesh/MooseMesh.h b/framework/include/mesh/MooseMesh.h index 09cc6554bb73..74f806a76c7b 100644 --- a/framework/include/mesh/MooseMesh.h +++ b/framework/include/mesh/MooseMesh.h @@ -1365,6 +1365,15 @@ class MooseMesh : public MooseObject, public Restartable, public PerfGraphInterf */ bool hasLowerD() const { return _has_lower_d; } + /** + * @return The set of lower-dimensional blocks for interior sides + */ + const std::set & interiorLowerDBlocks() const { return _lower_d_interior_blocks; } + /** + * @return The set of lower-dimensional blocks for boundary sides + */ + const std::set & boundaryLowerDBlocks() const { return _lower_d_boundary_blocks; } + protected: /// Deprecated (DO NOT USE) std::vector> _ghosting_functors; @@ -1752,6 +1761,10 @@ class MooseMesh : public MooseObject, public Restartable, public PerfGraphInterf /// Holds a map from neighbor subomdain ids to the boundary ids that are attached to it std::unordered_map> _neighbor_subdomain_boundary_ids; + /// Mesh blocks for interior lower-d elements in different types + std::set _lower_d_interior_blocks; + /// Mesh blocks for boundary lower-d elements in different types + std::set _lower_d_boundary_blocks; /// Holds a map from a high-order element side to its corresponding lower-d element std::unordered_map, const Elem *> _higher_d_elem_side_to_lower_d_elem; diff --git a/framework/include/utils/MooseTypes.h b/framework/include/utils/MooseTypes.h index 939f3e4422cd..76346e4fd4fb 100644 --- a/framework/include/utils/MooseTypes.h +++ b/framework/include/utils/MooseTypes.h @@ -659,8 +659,6 @@ namespace Moose { extern const processor_id_type INVALID_PROCESSOR_ID; extern const SubdomainID ANY_BLOCK_ID; -extern const SubdomainID INTERNAL_SIDE_LOWERD_ID; -extern const SubdomainID BOUNDARY_SIDE_LOWERD_ID; extern const SubdomainID INVALID_BLOCK_ID; extern const BoundaryID ANY_BOUNDARY_ID; extern const BoundaryID INVALID_BOUNDARY_ID; diff --git a/framework/src/actions/SetupMeshCompleteAction.C b/framework/src/actions/SetupMeshCompleteAction.C index b4995badf9f2..f948f9351789 100644 --- a/framework/src/actions/SetupMeshCompleteAction.C +++ b/framework/src/actions/SetupMeshCompleteAction.C @@ -71,8 +71,7 @@ SetupMeshCompleteAction::act() if (_mesh->uniformRefineLevel()) { - if (_mesh->meshSubdomains().count(Moose::INTERNAL_SIDE_LOWERD_ID) || - _mesh->meshSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID)) + if (!_mesh->interiorLowerDBlocks().empty() || !_mesh->boundaryLowerDBlocks().empty()) mooseError("HFEM does not support mesh uniform refinement currently."); Adaptivity::uniformRefine(_mesh.get()); diff --git a/framework/src/bcs/ArrayHFEMDirichletBC.C b/framework/src/bcs/ArrayHFEMDirichletBC.C index cb03dc7ee64f..59b26c0e7fba 100644 --- a/framework/src/bcs/ArrayHFEMDirichletBC.C +++ b/framework/src/bcs/ArrayHFEMDirichletBC.C @@ -34,10 +34,12 @@ ArrayHFEMDirichletBC::ArrayHFEMDirichletBC(const InputParameters & parameters) if (_uhat_var) { - if (!_uhat_var->activeSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID)) - paramError("uhat", - "Must be defined on BOUNDARY_SIDE_LOWERD_SUBDOMAIN subdomain that is added by " - "Mesh/build_all_side_lowerd_mesh=true"); + for (const auto & id : _uhat_var->activeSubdomains()) + if (_mesh.boundaryLowerDBlocks().count(id) == 0) + paramError("uhat", + "Must be defined on lower-dimensional boundary subdomain '" + + _mesh.getSubdomainName(id) + + "' that is added by Mesh/build_all_side_lowerd_mesh=true"); if (_uhat_var->count() != _count) paramError("uhat", diff --git a/framework/src/bcs/ArrayLowerDIntegratedBC.C b/framework/src/bcs/ArrayLowerDIntegratedBC.C index 456a9af6d817..9d8e605f88e1 100644 --- a/framework/src/bcs/ArrayLowerDIntegratedBC.C +++ b/framework/src/bcs/ArrayLowerDIntegratedBC.C @@ -35,15 +35,18 @@ ArrayLowerDIntegratedBC::ArrayLowerDIntegratedBC(const InputParameters & paramet _work_vector(_count) { const auto & lower_domains = _lowerd_var.activeSubdomains(); - if (!lower_domains.count(Moose::BOUNDARY_SIDE_LOWERD_ID) && lower_domains.size() != 1) - paramError( - "lowerd_variable", - "Must be only defined on the subdomain BOUNDARY_SIDE_LOWERD_SUBDOMAIN subdomain that is " - "added by Mesh/build_all_side_lowerd_mesh=true"); - - if (_var.activeSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID)) - paramError("variable", - "Must not be defined on the subdomain BOUNDARY_SIDE_LOWERD_SUBDOMAIN subdomain"); + for (const auto & id : _mesh.boundaryLowerDBlocks()) + if (lower_domains.count(id) == 0) + paramError("lowerd_variable", + "Must be defined on the boundary lower-dimensional subdomain '" + + _mesh.getSubdomainName(id) + + "' that is added by Mesh/build_all_side_lowerd_mesh=true"); + + for (const auto & id : _var.activeSubdomains()) + if (_mesh.boundaryLowerDBlocks().count(id) > 0) + paramError("variable", + "Must not be defined on the boundary lower-dimensional subdomain '" + + _mesh.getSubdomainName(id) + "'"); if (_lowerd_var.count() != _count) paramError("lowerd_variable", diff --git a/framework/src/bcs/HFEMDirichletBC.C b/framework/src/bcs/HFEMDirichletBC.C index 00d3beb37c6a..34e96a75fa8b 100644 --- a/framework/src/bcs/HFEMDirichletBC.C +++ b/framework/src/bcs/HFEMDirichletBC.C @@ -29,10 +29,12 @@ HFEMDirichletBC::HFEMDirichletBC(const InputParameters & parameters) { if (_uhat_var) { - if (!_uhat_var->activeSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID)) - paramError("uhat", - "Must be defined on BOUNDARY_SIDE_LOWERD_SUBDOMAIN subdomain that is added by " - "Mesh/build_all_side_lowerd_mesh=true"); + for (const auto & id : _uhat_var->activeSubdomains()) + if (_mesh.boundaryLowerDBlocks().count(id) == 0) + paramError("uhat", + "Must be defined on lower-dimensional boundary subdomain '" + + _mesh.getSubdomainName(id) + + "' that is added by Mesh/build_all_side_lowerd_mesh=true"); if (isParamValid("value")) paramError("uhat", "'uhat' and 'value' can not be both provided"); diff --git a/framework/src/bcs/LowerDIntegratedBC.C b/framework/src/bcs/LowerDIntegratedBC.C index 0d2fdf4db83e..864ecb4d0c7c 100644 --- a/framework/src/bcs/LowerDIntegratedBC.C +++ b/framework/src/bcs/LowerDIntegratedBC.C @@ -34,15 +34,18 @@ LowerDIntegratedBC::LowerDIntegratedBC(const InputParameters & parameters) _test_lambda(_lowerd_var.phiLower()) { const auto & lower_domains = _lowerd_var.activeSubdomains(); - if (!lower_domains.count(Moose::BOUNDARY_SIDE_LOWERD_ID) && lower_domains.size() != 1) - paramError( - "lowerd_variable", - "Must be only defined on the subdomain BOUNDARY_SIDE_LOWERD_SUBDOMAIN subdomain that is " - "added by Mesh/build_all_side_lowerd_mesh=true"); - - if (_var.activeSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID)) - paramError("variable", - "Must not be defined on the subdomain INTERNAL_SIDE_LOWERD_SUBDOMAIN subdomain"); + for (const auto & id : _mesh.boundaryLowerDBlocks()) + if (lower_domains.count(id) == 0) + paramError("lowerd_variable", + "Must be defined on the boundary lower-dimensional subdomain '" + + _mesh.getSubdomainName(id) + + "' that is added by Mesh/build_all_side_lowerd_mesh=true"); + + for (const auto & id : _var.activeSubdomains()) + if (_mesh.boundaryLowerDBlocks().count(id) > 0) + paramError("variable", + "Must not be defined on the boundary lower-dimensional subdomain '" + + _mesh.getSubdomainName(id) + "'"); // Note: the above two conditions also ensure that the variable and lower-d variable are // different. diff --git a/framework/src/dgkernels/ArrayDGLowerDKernel.C b/framework/src/dgkernels/ArrayDGLowerDKernel.C index 68499733cc54..6fc674af34fd 100644 --- a/framework/src/dgkernels/ArrayDGLowerDKernel.C +++ b/framework/src/dgkernels/ArrayDGLowerDKernel.C @@ -42,15 +42,18 @@ ArrayDGLowerDKernel::ArrayDGLowerDKernel(const InputParameters & parameters) _work_vector(_count) { const auto & lower_domains = _lowerd_var.activeSubdomains(); - if (!lower_domains.count(Moose::INTERNAL_SIDE_LOWERD_ID) && lower_domains.size() != 1) - paramError( - "lowerd_variable", - "Must be only defined on the subdomain INTERNAL_SIDE_LOWERD_SUBDOMAIN subdomain that is " - "added by Mesh/build_all_side_lowerd_mesh=true"); - - if (_var.activeSubdomains().count(Moose::INTERNAL_SIDE_LOWERD_ID)) - paramError("variable", - "Must not be defined on the subdomain INTERNAL_SIDE_LOWERD_SUBDOMAIN subdomain"); + for (const auto & id : _mesh.interiorLowerDBlocks()) + if (lower_domains.count(id) == 0) + paramError("lowerd_variable", + "Must be defined on the interior lower-dimensional subdomain '" + + _mesh.getSubdomainName(id) + + "' that is added by Mesh/build_all_side_lowerd_mesh=true"); + + for (const auto & id : _var.activeSubdomains()) + if (_mesh.interiorLowerDBlocks().count(id) > 0) + paramError("variable", + "Must not be defined on the interior lower-dimensional subdomain '" + + _mesh.getSubdomainName(id) + "'"); if (_lowerd_var.count() != _count) paramError("lowerd_variable", diff --git a/framework/src/dgkernels/DGLowerDKernel.C b/framework/src/dgkernels/DGLowerDKernel.C index ea28a581620e..68c275eb1508 100644 --- a/framework/src/dgkernels/DGLowerDKernel.C +++ b/framework/src/dgkernels/DGLowerDKernel.C @@ -41,15 +41,18 @@ DGLowerDKernel::DGLowerDKernel(const InputParameters & parameters) _test_lambda(_lowerd_var.phiLower()) { const auto & lower_domains = _lowerd_var.activeSubdomains(); - if (!lower_domains.count(Moose::INTERNAL_SIDE_LOWERD_ID) && lower_domains.size() != 1) - paramError( - "lowerd_variable", - "Must be only defined on the subdomain INTERNAL_SIDE_LOWERD_SUBDOMAIN subdomain that is " - "added by Mesh/build_all_side_lowerd_mesh=true"); - - if (_var.activeSubdomains().count(Moose::INTERNAL_SIDE_LOWERD_ID)) - paramError("variable", - "Must not be defined on the subdomain INTERNAL_SIDE_LOWERD_SUBDOMAIN subdomain"); + for (const auto & id : _mesh.interiorLowerDBlocks()) + if (lower_domains.count(id) == 0) + paramError("lowerd_variable", + "Must be defined on the interior lower-dimensional subdomain '" + + _mesh.getSubdomainName(id) + + "'that is added by Mesh/build_all_side_lowerd_mesh=true"); + + for (const auto & id : _var.activeSubdomains()) + if (_mesh.interiorLowerDBlocks().count(id) > 0) + paramError("variable", + "Must not be defined on the interior lower-dimensional subdomain'" + + _mesh.getSubdomainName(id) + "'"); // Note: the above two conditions also ensure that the variable and lower-d variable are // different. diff --git a/framework/src/interfaces/BlockRestrictable.C b/framework/src/interfaces/BlockRestrictable.C index b26bc57e3eeb..4dc510140c1b 100644 --- a/framework/src/interfaces/BlockRestrictable.C +++ b/framework/src/interfaces/BlockRestrictable.C @@ -338,8 +338,9 @@ BlockRestrictable::checkVariable(const MooseVariableFieldBase & variable) const { // a variable defined on all internal sides does not need this check because // it can be coupled with other variables in DG kernels - if (variable.activeSubdomains().count(Moose::INTERNAL_SIDE_LOWERD_ID) > 0) - return; + for (const auto & id : variable.activeSubdomains()) + if (_blk_mesh->interiorLowerDBlocks().count(id) > 0) + return; if (!isBlockSubset(variable.activeSubdomains())) { diff --git a/framework/src/loops/ComputeFullJacobianThread.C b/framework/src/loops/ComputeFullJacobianThread.C index a45dd4eb622b..b1e01332dcc1 100644 --- a/framework/src/loops/ComputeFullJacobianThread.C +++ b/framework/src/loops/ComputeFullJacobianThread.C @@ -300,7 +300,7 @@ ComputeFullJacobianThread::computeOnInternalFace(const Elem * neighbor) if (dg->variable().number() == ivar && dg->isImplicit() && dg->hasBlocks(neighbor->subdomain_id()) && (jvariable.activeOnSubdomain(_subdomain) || - jvariable.activeOnSubdomain(Moose::INTERNAL_SIDE_LOWERD_ID))) + jvariable.activeOnSubdomains(_fe_problem.mesh().interiorLowerDBlocks()))) { dg->prepareShapes(jvar); dg->prepareNeighborShapes(jvar); diff --git a/framework/src/mesh/MooseMesh.C b/framework/src/mesh/MooseMesh.C index bc0ebdff5404..0fb7074e1f63 100644 --- a/framework/src/mesh/MooseMesh.C +++ b/framework/src/mesh/MooseMesh.C @@ -60,6 +60,7 @@ #include "libmesh/default_coupling.h" #include "libmesh/ghost_point_neighbors.h" #include "libmesh/fe_type.h" +#include "libmesh/enum_to_string.h" static const int GRAIN_SIZE = 1; // the grain_size does not have much influence on our execution speed @@ -297,6 +298,8 @@ MooseMesh::MooseMesh(const MooseMesh & other_mesh) _patch_update_strategy(other_mesh._patch_update_strategy), _regular_orthogonal_mesh(false), _is_split(other_mesh._is_split), + _lower_d_interior_blocks(other_mesh._lower_d_interior_blocks), + _lower_d_boundary_blocks(other_mesh._lower_d_boundary_blocks), _has_lower_d(other_mesh._has_lower_d), _allow_recovery(other_mesh._allow_recovery), _construct_node_list_from_side_list(other_mesh._construct_node_list_from_side_list), @@ -632,19 +635,69 @@ MooseMesh::buildLowerDMesh() // remove existing lower-d element first std::set deleteable_elems; for (auto & elem : mesh.element_ptr_range()) - if (elem->subdomain_id() == Moose::INTERNAL_SIDE_LOWERD_ID || - elem->subdomain_id() == Moose::BOUNDARY_SIDE_LOWERD_ID) + if (_lower_d_interior_blocks.count(elem->subdomain_id()) || + _lower_d_boundary_blocks.count(elem->subdomain_id())) deleteable_elems.insert(elem); else if (elem->n_sides() > max_n_sides) max_n_sides = elem->n_sides(); for (auto & elem : deleteable_elems) mesh.delete_elem(elem); + for (const auto & id : _lower_d_interior_blocks) + _mesh_subdomains.erase(id); + for (const auto & id : _lower_d_boundary_blocks) + _mesh_subdomains.erase(id); + _lower_d_interior_blocks.clear(); + _lower_d_boundary_blocks.clear(); mesh.comm().max(max_n_sides); deleteable_elems.clear(); + // get all side types + std::set interior_side_types; + std::set boundary_side_types; + for (const auto & elem : mesh.active_element_ptr_range()) + for (const auto side : elem->side_index_range()) + { + Elem * neig = elem->neighbor_ptr(side); + std::unique_ptr side_elem(elem->build_side_ptr(side)); + if (neig) + interior_side_types.insert(side_elem->type()); + else + boundary_side_types.insert(side_elem->type()); + } + mesh.comm().set_union(interior_side_types); + mesh.comm().set_union(boundary_side_types); + + // assign block ids for different side types + std::map interior_block_ids; + std::map boundary_block_ids; + // we assume this id is not used by the mesh + auto id = libMesh::Elem::invalid_subdomain_id - 2; + for (const auto & tpid : interior_side_types) + { + const auto type = ElemType(tpid); + mesh.subdomain_name(id) = "INTERNAL_SIDE_LOWERD_SUBDOMAIN_" + Utility::enum_to_string(type); + interior_block_ids[type] = id; + _lower_d_interior_blocks.insert(id); + if (_mesh_subdomains.count(id) > 0) + mooseError("Trying to add a mesh block with id ", id, " that has existed in the mesh"); + _mesh_subdomains.insert(id); + --id; + } + for (const auto & tpid : boundary_side_types) + { + const auto type = ElemType(tpid); + mesh.subdomain_name(id) = "BOUNDARY_SIDE_LOWERD_SUBDOMAIN_" + Utility::enum_to_string(type); + boundary_block_ids[type] = id; + _lower_d_boundary_blocks.insert(id); + if (_mesh_subdomains.count(id) > 0) + mooseError("Trying to add a mesh block with id ", id, " that has existed in the mesh"); + _mesh_subdomains.insert(id); + --id; + } + dof_id_type max_elem_id = mesh.max_elem_id(); unique_id_type max_unique_id = mesh.parallel_max_unique_id(); @@ -681,9 +734,9 @@ MooseMesh::buildLowerDMesh() // Add subdomain ID if (neig) - side_elem->subdomain_id() = Moose::INTERNAL_SIDE_LOWERD_ID; + side_elem->subdomain_id() = interior_block_ids.at(side_elem->type()); else - side_elem->subdomain_id() = Moose::BOUNDARY_SIDE_LOWERD_ID; + side_elem->subdomain_id() = boundary_block_ids.at(side_elem->type()); // set ids consistently across processors (these ids will be temporary) side_elem->set_id(max_elem_id + elem->id() * max_n_sides + side); @@ -713,11 +766,6 @@ MooseMesh::buildLowerDMesh() for (auto & elem : side_elems) mesh.add_elem(elem); - _mesh_subdomains.insert(Moose::INTERNAL_SIDE_LOWERD_ID); - mesh.subdomain_name(Moose::INTERNAL_SIDE_LOWERD_ID) = "INTERNAL_SIDE_LOWERD_SUBDOMAIN"; - _mesh_subdomains.insert(Moose::BOUNDARY_SIDE_LOWERD_ID); - mesh.subdomain_name(Moose::BOUNDARY_SIDE_LOWERD_ID) = "BOUNDARY_SIDE_LOWERD_SUBDOMAIN"; - // we do all the stuff in prepare_for_use such as renumber_nodes_and_elements(), // update_parallel_id_counts(), cache_elem_dims(), etc. except partitioning here. const bool skip_partitioning_old = mesh.skip_partitioning(); @@ -1326,9 +1374,13 @@ MooseMesh::cacheInfo() _block_node_list.clear(); _higher_d_elem_side_to_lower_d_elem.clear(); _lower_d_elem_to_higher_d_elem_side.clear(); + _lower_d_interior_blocks.clear(); + _lower_d_boundary_blocks.clear(); + + auto & mesh = getMesh(); // TODO: Thread this! - for (const auto & elem : getMesh().element_ptr_range()) + for (const auto & elem : mesh.element_ptr_range()) { const Elem * ip_elem = elem->interior_parent(); @@ -1346,6 +1398,18 @@ MooseMesh::cacheInfo() std::pair, const Elem *>(pair, elem)); _lower_d_elem_to_higher_d_elem_side.insert( std::pair(elem, ip_side)); + + auto id = elem->subdomain_id(); + if (ip_elem->neighbor_ptr(ip_side)) + { + if (mesh.subdomain_name(id).find("INTERNAL_SIDE_LOWERD_SUBDOMAIN_") != std::string::npos) + _lower_d_interior_blocks.insert(id); + } + else + { + if (mesh.subdomain_name(id).find("BOUNDARY_SIDE_LOWERD_SUBDOMAIN_") != std::string::npos) + _lower_d_boundary_blocks.insert(id); + } } } @@ -1356,7 +1420,7 @@ MooseMesh::cacheInfo() } } - for (const auto & elem : getMesh().active_local_element_ptr_range()) + for (const auto & elem : mesh.active_local_element_ptr_range()) { SubdomainID subdomain_id = elem->subdomain_id(); auto & sub_data = _sub_to_data[subdomain_id]; diff --git a/framework/src/problems/DisplacedProblem.C b/framework/src/problems/DisplacedProblem.C index 900f690defe9..a62b9af62f93 100644 --- a/framework/src/problems/DisplacedProblem.C +++ b/framework/src/problems/DisplacedProblem.C @@ -875,7 +875,7 @@ DisplacedProblem::reinitElemNeighborAndLowerD(const Elem * elem, reinitNeighbor(elem, side, tid); const Elem * lower_d_elem = _mesh.getLowerDElem(elem, side); - if (lower_d_elem && lower_d_elem->subdomain_id() == Moose::INTERNAL_SIDE_LOWERD_ID) + if (lower_d_elem && _mesh.interiorLowerDBlocks().count(lower_d_elem->subdomain_id()) > 0) reinitLowerDElem(lower_d_elem, tid); else { @@ -884,7 +884,7 @@ DisplacedProblem::reinitElemNeighborAndLowerD(const Elem * elem, auto & neighbor_side = _assembly[tid][currentNlSysNum()]->neighborSide(); const Elem * lower_d_elem_neighbor = _mesh.getLowerDElem(neighbor, neighbor_side); if (lower_d_elem_neighbor && - lower_d_elem_neighbor->subdomain_id() == Moose::INTERNAL_SIDE_LOWERD_ID) + _mesh.interiorLowerDBlocks().count(lower_d_elem_neighbor->subdomain_id()) > 0) { auto qps = _assembly[tid][currentNlSysNum()]->qPointsFaceNeighbor().stdVector(); std::vector reference_points; diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C index 31bbac8f9827..30923e011aed 100644 --- a/framework/src/problems/FEProblemBase.C +++ b/framework/src/problems/FEProblemBase.C @@ -2253,7 +2253,7 @@ FEProblemBase::reinitElemNeighborAndLowerD(const Elem * elem, reinitNeighbor(elem, side, tid); const Elem * lower_d_elem = _mesh.getLowerDElem(elem, side); - if (lower_d_elem && lower_d_elem->subdomain_id() == Moose::INTERNAL_SIDE_LOWERD_ID) + if (lower_d_elem && _mesh.interiorLowerDBlocks().count(lower_d_elem->subdomain_id()) > 0) reinitLowerDElem(lower_d_elem, tid); else { @@ -2262,7 +2262,7 @@ FEProblemBase::reinitElemNeighborAndLowerD(const Elem * elem, auto & neighbor_side = _assembly[tid][0]->neighborSide(); const Elem * lower_d_elem_neighbor = _mesh.getLowerDElem(neighbor, neighbor_side); if (lower_d_elem_neighbor && - lower_d_elem_neighbor->subdomain_id() == Moose::INTERNAL_SIDE_LOWERD_ID) + _mesh.interiorLowerDBlocks().count(lower_d_elem_neighbor->subdomain_id()) > 0) { auto qps = _assembly[tid][0]->qPointsFaceNeighbor().stdVector(); std::vector reference_points; @@ -7711,8 +7711,7 @@ FEProblemBase::initialAdaptMesh() _cycles_completed = 0; if (n) { - if (_mesh.meshSubdomains().count(Moose::INTERNAL_SIDE_LOWERD_ID) || - _mesh.meshSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID)) + if (!_mesh.interiorLowerDBlocks().empty() || !_mesh.boundaryLowerDBlocks().empty()) mooseError("HFEM does not support mesh adaptivity currently."); TIME_SECTION("initialAdaptMesh", 2, "Performing Initial Adaptivity"); @@ -7757,8 +7756,7 @@ FEProblemBase::adaptMesh() for (unsigned int i = 0; i < cycles_per_step; ++i) { - if (_mesh.meshSubdomains().count(Moose::INTERNAL_SIDE_LOWERD_ID) || - _mesh.meshSubdomains().count(Moose::BOUNDARY_SIDE_LOWERD_ID)) + if (!_mesh.interiorLowerDBlocks().empty() || !_mesh.boundaryLowerDBlocks().empty()) mooseError("HFEM does not support mesh adaptivity currently."); // Markers were already computed once by Executioner diff --git a/framework/src/systems/NonlinearSystemBase.C b/framework/src/systems/NonlinearSystemBase.C index f380e541d37b..8ca2a1160629 100644 --- a/framework/src/systems/NonlinearSystemBase.C +++ b/framework/src/systems/NonlinearSystemBase.C @@ -3640,8 +3640,10 @@ NonlinearSystemBase::checkKernelCoverage(const std::set & mesh_subd std::inserter(difference, difference.end())); // there supposed to be no kernels on this lower-dimensional subdomain - difference.erase(Moose::INTERNAL_SIDE_LOWERD_ID); - difference.erase(Moose::BOUNDARY_SIDE_LOWERD_ID); + for (const auto & id : _mesh.interiorLowerDBlocks()) + difference.erase(id); + for (const auto & id : _mesh.boundaryLowerDBlocks()) + difference.erase(id); if (!difference.empty()) { @@ -3681,8 +3683,9 @@ NonlinearSystemBase::checkKernelCoverage(const std::set & mesh_subd for (auto & var_name : vars) { auto blks = getSubdomainsForVar(var_name); - if (blks.count(Moose::INTERNAL_SIDE_LOWERD_ID) || blks.count(Moose::BOUNDARY_SIDE_LOWERD_ID)) - difference.erase(var_name); + for (const auto & id : blks) + if (_mesh.interiorLowerDBlocks().count(id) > 0 || _mesh.boundaryLowerDBlocks().count(id) > 0) + difference.erase(var_name); } if (!difference.empty()) diff --git a/framework/src/utils/MooseTypes.C b/framework/src/utils/MooseTypes.C index cf85561965fc..284ac4189ed9 100644 --- a/framework/src/utils/MooseTypes.C +++ b/framework/src/utils/MooseTypes.C @@ -17,8 +17,6 @@ namespace Moose { const processor_id_type INVALID_PROCESSOR_ID = libMesh::DofObject::invalid_processor_id; const SubdomainID ANY_BLOCK_ID = libMesh::Elem::invalid_subdomain_id - 1; -const SubdomainID INTERNAL_SIDE_LOWERD_ID = libMesh::Elem::invalid_subdomain_id - 2; -const SubdomainID BOUNDARY_SIDE_LOWERD_ID = libMesh::Elem::invalid_subdomain_id - 3; const SubdomainID INVALID_BLOCK_ID = libMesh::Elem::invalid_subdomain_id; const BoundaryID ANY_BOUNDARY_ID = static_cast(-1); const BoundaryID INVALID_BOUNDARY_ID = libMesh::BoundaryInfo::invalid_id; diff --git a/test/tests/kernels/hfem/3d-lower-d-volumes.i b/test/tests/kernels/hfem/3d-lower-d-volumes.i index 1d9907212c64..2f4374cb1f22 100644 --- a/test/tests/kernels/hfem/3d-lower-d-volumes.i +++ b/test/tests/kernels/hfem/3d-lower-d-volumes.i @@ -18,17 +18,17 @@ [uhat] order = CONSTANT family = MONOMIAL - block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN + block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN_QUAD4 [] [lambda] order = CONSTANT family = MONOMIAL - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_QUAD4 [] [lambdab] order = CONSTANT family = MONOMIAL - block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN + block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN_QUAD4 [] [] @@ -59,12 +59,12 @@ type = Reaction variable = uhat rate = '1' - block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN + block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN_QUAD4 [] [uhat_coupled] type = CoupledForce variable = uhat - block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN + block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN_QUAD4 v = lambdab coef = '1' [] @@ -99,7 +99,7 @@ [lambdanorm] type = ElementL2Norm variable = lambda - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_QUAD4 [] [] diff --git a/test/tests/kernels/hfem/3d-prism.i b/test/tests/kernels/hfem/3d-prism.i new file mode 100644 index 000000000000..32ff6683f849 --- /dev/null +++ b/test/tests/kernels/hfem/3d-prism.i @@ -0,0 +1,119 @@ +[Mesh] + [square] + type = GeneratedMeshGenerator + nx = 2 + ny = 2 + nz = 2 + dim = 3 + elem_type = PRISM6 + [] + build_all_side_lowerd_mesh = true +[] + +[Variables] + [u] + order = THIRD + family = MONOMIAL + block = 0 + [] + [uhat] + order = CONSTANT + family = MONOMIAL + block = 'BOUNDARY_SIDE_LOWERD_SUBDOMAIN_QUAD4 BOUNDARY_SIDE_LOWERD_SUBDOMAIN_TRI3' + [] + [lambda] + order = CONSTANT + family = MONOMIAL + block = 'INTERNAL_SIDE_LOWERD_SUBDOMAIN_QUAD4 INTERNAL_SIDE_LOWERD_SUBDOMAIN_TRI3' + [] + [lambdab] + order = CONSTANT + family = MONOMIAL + block = 'BOUNDARY_SIDE_LOWERD_SUBDOMAIN_QUAD4 BOUNDARY_SIDE_LOWERD_SUBDOMAIN_TRI3' + [] +[] + +[AuxVariables] + [v] + order = CONSTANT + family = MONOMIAL + block = 0 + initial_condition = '1' + [] +[] + +[Kernels] + [diff] + type = MatDiffusion + variable = u + diffusivity = '1' + block = 0 + [] + [source] + type = CoupledForce + variable = u + v = v + coef = '1' + block = 0 + [] + [reaction] + type = Reaction + variable = uhat + rate = '1' + block = 'BOUNDARY_SIDE_LOWERD_SUBDOMAIN_QUAD4 BOUNDARY_SIDE_LOWERD_SUBDOMAIN_TRI3' + [] + [uhat_coupled] + type = CoupledForce + variable = uhat + block = 'BOUNDARY_SIDE_LOWERD_SUBDOMAIN_QUAD4 BOUNDARY_SIDE_LOWERD_SUBDOMAIN_TRI3' + v = lambdab + coef = '1' + [] +[] + +[DGKernels] + [surface] + type = HFEMDiffusion + variable = u + lowerd_variable = lambda + [] +[] + +[BCs] + [all] + type = HFEMDirichletBC + boundary = 'left right top bottom back front' + variable = u + lowerd_variable = lambdab + uhat = uhat + [] +[] + +[Postprocessors] + [intu] + type = ElementIntegralVariablePostprocessor + variable = u + block = 0 + [] + [lambdanorm] + type = ElementL2Norm + variable = lambda + block = 'INTERNAL_SIDE_LOWERD_SUBDOMAIN_QUAD4 INTERNAL_SIDE_LOWERD_SUBDOMAIN_TRI3' + [] +[] + +[Executioner] + type = Steady + solve_type = 'NEWTON' + petsc_options_iname = '-pc_type -snes_linesearch_type -pc_factor_mat_solver_type' + petsc_options_value = 'lu basic mumps' +[] + +[Outputs] + [out] + # we hide lambda because it may flip sign due to element + # renumbering with distributed mesh + type = Exodus + hide = lambda + [] +[] diff --git a/test/tests/kernels/hfem/array_dirichlet.i b/test/tests/kernels/hfem/array_dirichlet.i index d17f7d7df70e..89c2d9257772 100644 --- a/test/tests/kernels/hfem/array_dirichlet.i +++ b/test/tests/kernels/hfem/array_dirichlet.i @@ -18,13 +18,13 @@ [lambda] order = CONSTANT family = MONOMIAL - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 components = 2 [] [lambdab] order = CONSTANT family = MONOMIAL - block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN + block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN_EDGE2 components = 2 [] [] @@ -88,7 +88,7 @@ [lambdanorm] type = ElementArrayL2Norm variable = lambda - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 [] [] diff --git a/test/tests/kernels/hfem/array_dirichlet_pjfnk.i b/test/tests/kernels/hfem/array_dirichlet_pjfnk.i index 76f204434037..d4f5b1005afd 100644 --- a/test/tests/kernels/hfem/array_dirichlet_pjfnk.i +++ b/test/tests/kernels/hfem/array_dirichlet_pjfnk.i @@ -18,13 +18,13 @@ [lambda] order = CONSTANT family = MONOMIAL - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 components = 2 [] [lambdab] order = CONSTANT family = MONOMIAL - block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN + block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN_EDGE2 components = 2 [] [] @@ -101,7 +101,7 @@ [lambdanorm] type = ElementArrayL2Norm variable = lambda - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 [] [] diff --git a/test/tests/kernels/hfem/array_dirichlet_transform.i b/test/tests/kernels/hfem/array_dirichlet_transform.i index 5695e6534e2a..23e8d132c140 100644 --- a/test/tests/kernels/hfem/array_dirichlet_transform.i +++ b/test/tests/kernels/hfem/array_dirichlet_transform.i @@ -18,13 +18,13 @@ [lambda] order = CONSTANT family = MONOMIAL - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 components = 2 [] [lambdab] order = CONSTANT family = MONOMIAL - block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN + block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN_EDGE2 components = 2 [] [] @@ -88,7 +88,7 @@ [lambdanorm] type = ElementArrayL2Norm variable = lambda - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 [] [] diff --git a/test/tests/kernels/hfem/array_dirichlet_transform_bc.i b/test/tests/kernels/hfem/array_dirichlet_transform_bc.i index bbb71d0bb500..4cd6e1b0f7bb 100644 --- a/test/tests/kernels/hfem/array_dirichlet_transform_bc.i +++ b/test/tests/kernels/hfem/array_dirichlet_transform_bc.i @@ -18,13 +18,13 @@ [lambda] order = CONSTANT family = MONOMIAL - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 components = 2 [] [lambdab] order = CONSTANT family = MONOMIAL - block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN + block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN_EDGE2 components = 2 [] [] @@ -88,7 +88,7 @@ [lambdanorm] type = ElementArrayL2Norm variable = lambda - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 [] [] diff --git a/test/tests/kernels/hfem/array_neumann.i b/test/tests/kernels/hfem/array_neumann.i index b8663bdcb179..89898d779ab9 100644 --- a/test/tests/kernels/hfem/array_neumann.i +++ b/test/tests/kernels/hfem/array_neumann.i @@ -18,7 +18,7 @@ [lambda] order = CONSTANT family = MONOMIAL - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 components = 2 [] [] @@ -92,7 +92,7 @@ [lambdanorm] type = ElementArrayL2Norm variable = lambda - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 [] [] diff --git a/test/tests/kernels/hfem/array_robin.i b/test/tests/kernels/hfem/array_robin.i index 4d396b4b4dab..95fa6b173e9a 100644 --- a/test/tests/kernels/hfem/array_robin.i +++ b/test/tests/kernels/hfem/array_robin.i @@ -18,7 +18,7 @@ [lambda] order = CONSTANT family = MONOMIAL - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 components = 2 [] [] @@ -81,7 +81,7 @@ [lambdanorm] type = ElementArrayL2Norm variable = lambda - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 [] [] diff --git a/test/tests/kernels/hfem/dirichlet.i b/test/tests/kernels/hfem/dirichlet.i index c0d3f947c56d..93399242562b 100644 --- a/test/tests/kernels/hfem/dirichlet.i +++ b/test/tests/kernels/hfem/dirichlet.i @@ -17,17 +17,17 @@ [uhat] order = CONSTANT family = MONOMIAL - block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN + block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN_EDGE2 [] [lambda] order = CONSTANT family = MONOMIAL - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 [] [lambdab] order = CONSTANT family = MONOMIAL - block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN + block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN_EDGE2 [] [] @@ -58,12 +58,12 @@ type = Reaction variable = uhat rate = '1' - block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN + block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN_EDGE2 [] [uhat_coupled] type = CoupledForce variable = uhat - block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN + block = BOUNDARY_SIDE_LOWERD_SUBDOMAIN_EDGE2 v = lambdab coef = '1' [] @@ -96,7 +96,7 @@ [lambdanorm] type = ElementL2Norm variable = lambda - block = INTERNAL_SIDE_LOWERD_SUBDOMAIN + block = INTERNAL_SIDE_LOWERD_SUBDOMAIN_EDGE2 [] [] diff --git a/test/tests/kernels/hfem/gold/3d-prism_out.e b/test/tests/kernels/hfem/gold/3d-prism_out.e new file mode 100644 index 0000000000000000000000000000000000000000..6560b3c51fcf855942d291586b20b2604a33c0cb GIT binary patch literal 75264 zcmeHQZHy$xS>7YDedjnhu>%PoA?=U&Blhjx?Vat|$+FqI-LqHR?b>&H=NN-Jy)!+# zZO=@Pr+fDHEFoDaKa&6<68r%nK>~sx1;h`2fJ8!q@(&4tL;?g0@gbo|L*t+(E)IKenBR@{wyY=av z9Ikg%guk!T4ZS$%@4AC9=*RAraM(t^yiRQLkaUu7M>;_lGAn=3e&g@=_Pm3j$8v(Q z<)V=7#(M}(-*<3-7r>u9g~H~=qe%M=w+pw$LEPmBaj?AJeX@OYFbtxe8;;^d`NKO9 z&WyNIl-VnLVe$gM&r=;N8>APZA8_#hh`*D4l%Y7u7j-a=?~xC8^S#4}KNyAs%A|@* z{__6qvhENxzHMU-Au^5X{U62ik$yDxJ<_BOe(!zq;qK$O9!7lqdyG5jd+Pn?5cfzl z{yoZqI`r>-0)dASU;iHCPW!zB#J#;0LPB^ac|w_02WkIn`Ed8o<0|7`6(QtbMcj8$ zSY+JSC&hjD#}I>g#H^qkQC?Jjs&w4xo(;O)%U*XRWruP2eIB@1(eEK2$$uS7<<*o} zD*Af`7I>>;seCZRdOxx5B%f0>nNFsI$CP-A-477&F1BGVE5Dw{=B}6@l#?z_91n@K zw8DPhZ^bMJ{4UD|>o)c6Bck4luMQ-OMV#ATTwi;39g%PV@6a3KZsDO%;(9p2V+Bh* zw#h8hwRkUFT3r)0t;9og9*>Y#g~vL9c*J>Ua%I}e=(}zE!Up5>AjY5YW?|t`9gm+Q z-uL6Y8>a<*kA4;JDf-wb5RYy7yK&ax73uo~;^X)b@$M4*&gUQLRtN46Ir8DL2XQ?t z;zy#9BR}quzlr~R&ymlMui<)FYlrLNwi!3UWj-?Qy*SlDT&5Gp)ljwxwTpEyQZ-YsZRm?!<60zIt5)j@jR5D$kOUh?mr@OPN`I8<6# zU#mlZ_YDN*-(8ya-TR1_e;41SAmLEoB_GtGefPdM5SV`#-<|H=`-vy(Q?TBRy&XLB zd$-|K{Ku`3osaTnVPT2R_5>-FN8z zH*gO85T~&B`tw?W?KJ~0{)HTIHy?E*L+Sd8ijZOE3o1HHc(!YJPr+;H} zd;P-Z>bYAr|0ksXrqBN!uO+S2H@l!$q;DsT<4~CQDJGXN7FLZ`s|5j$ey#Alp|F^RI&+Gqr{eLUV|4pyp%o8~0Ks^U&uP)lMX=mqn zJ;z8m&aID6b6i=C2fqvV9LK!}=e;=JgOg*59RK53-TQEIyo}>u9N((?7QfH$(Vk9Q zI=@RA7|!qWH-?iAMH9!%8O|_{H!z&@Hd@i#V3aAG%G%me>%;KM}DesMKi<6Z~n$G zMKi<5dxk6ilP;!({9+jIb$*jR#u5HoxDtiD_UMzsYwPm+2=i!<3AvbThqtNAW_5f2YTM?`K{-%p09dd>iMT=CJ&fDif+mWpGlX>ABL&4 z@SeORt&GEG@=B$NJY<+k7x~I(l{P*zKEI*(&ouFw@klS9d9QwF7{950XS#`@o|!J* ztNdepKJ#15FFrHAdS-kjKPvr35uA%SkKsIl^CZp>o}joN%>Rq#W7k3=L0w&#`ze|C7e&-d=lr=I60pC z<2YAwuHjVWiarDP;rt-ZM{z!m^El2){2|ZkPQFjxl8=-vKI;J8-7baC*>dCtviX&?-7$Sq0+%`Q3jZQDh|`4Xraugx`Fvqq=WBK zc2wR`eu}cC;;6cgariCL!f%l`>K%Sd<)!+rs%z9YQk)SY-^xeV>+Q6K+rzENh0&fD z-l1bI7yUjTGkJm^`z>94$hxB+4tu5Sx#thfTuITQAJUcc>}!?Y7;b+|j~`TaO@ z@Y9sE>d~CV#bffX7lx7V#(oqZ4F|DvG;q}1kHvtWy)0+5I2-VKtNGa3A-p;&qN zaoEF*matEi(e{w}kbM0aCp&h7*lB?mxcEyap4d6n+G#)Sgneh}_|mDPCr%t)qNX{1 zoC|i&`_Z0CPzE_3b-Iw}*S9aW#^D_2_PhQ#m`hnjX!T+=C3rtK- z>J0?&7~t_Irvfbd$n_P*LbcH!jA93KWc->U^_G$S@=rJXX~%I#{ou8c?|LoP3u0DH zR-zqe*;#G3FZf+iYeaj&;JnxOcKu-$9S&-Memo4j@EzWbbfn~d+1c2H zg*#5rCvE$F2^fQxok#ki+x9!&s2e|;y;yZez#tC1uG@}FA4;*uL+(uU+x!nUb<|*)glbkSU`Go;{d88@66eVKP$9DNF$2X^mDW15$V9 zJfv$rDlk_dDwymFIlQ#JRh6gX;f_-Qq#(@!#Y~`TRC7Si5L4uoT$i0^S1(<>IG(nO zkj5*AitpG5zq~p+xag5ig4Ou4D(Yb^AUc-o zk`l1CiZw5Kj??jaLH8=5OO9n|-kV$JH&)MSjH^(?G2?;JIw1t z07hFouV50VjkiQdz6J;aglx$KCbgm~H8((+LhFmiTWYhgC9d59niT zq3SSim(kE;__PWSVq-SCN+6**B*m}oLc8A6;H@ZxDuRJJe>ezX)$_5U71o_nb_g3l zQ4MNDsAL7+E+9hcF;&QqJdF1CTm}cFH#3nPtmc25?z1@IM=?&h3r^%r1sV6jMTb(I?KlH z`DXLl7j5DhP5C>U>T7~JXFB|G4b-`7pnbKvBdnoTnsgf;&egf197=U=tVmP@@a?%e zSLC6I)kYd5Zv^dYt}Im)g1094W$8JAAqs!kkSbEEh|t0*y`;1cONB%W(0PGIbQFhJ z>J%$wC8H_%aR&3=Xsu>rV}q@=D-qCm(7(>EL_pIa^;)|U0WZeB*S@MVzV05G(rWnw z41o7x#yG%m$^rUja^c|sy*1&00o1AMup;0fXn9z_*a}M!o`{Xt-uL<~)nU}DSo4OS za!7}?mRma5We&*E6q@D|I50Cusd+9_u)S97VkC~sl6&yka)oyWGDkMaF5ILdWDI9h zOgv;Cpa|HFLp3sBk!0!G0lnDGV6tpMlKrDkOf5AFpilx1IEYtrSW`Ik;S`i6G(#2? z9_J*O=?*HysJdczFpwtL3`*V^6cX`80!AF;$`!hy4~Gv7*D9{gu+^t-KP4|>z5MQF&CR0V0Skh2Jv1m%bI0J@?yF?ruL-+cHk_@rAddq)cibU zUZmg41G0(@@5A5_z_}dYaU5wbPP!q5WtB{dD0%qAB%Jhj0A^KxilW^tJBjaM`bzyx4$vK7 z6bJs48XhW{0k*$r)uU*?vZr_89(1r*!f4RWY`H1svU6(b=~GWX`P97b#ole`W zR&bV`D+7&yiMIfRZ)U*pbc`VAwQNG8d0DOkbkSgA=pHT*X!xfWZpV>Y7fiv-ZY+H} zfJzL<4CQ_Uu`!TBR$qNtRNBFQsX$T=A2Vu~6F|Dxlp&h)FafM8E~2EEQdpH-X0%B1 z36ru63PsINwKc8y6PXc!(}Ph+MWzVldem=~;E_VZ+6EMsG$kpC%Z>p96@{$%8_ngC zrbMFoRI%vr%{bFgt=fRfK}yLB-yAe1IarOriCQ)b)m5-dIAuA<5c4; z0-D-wMAymPkK9+h)Jg!Y36t6^GbQbdVagss!=1qEqd|qA!k$00s2l(od&5{dz(Rm? z1&=&5PcqHUx_J(CWrw*O4*hXJ(QGBU4jWe(Fs4!tQQ%Y8WYpT>h@J4HccExBufC#J zgMsTG&;Z6+NV2!B;K)PM8k==ouY=71nDVkrF>O;sWmv5fFFTv-FK%ybPG29wqbi8G zyI7V?{mxu$675>-vZNy=_vV>Gr-%|D6=a?%aM9{Vs9t-*vKhtj4{G^Q6ecr>%#*AG zj!B_S(UTo0*_lEMsNC7z9e$1*fIY z(r+o2ja_XhyawLRLMZlwxWpyfMh@rh(|Mt^$>qk~r|0f5`O;^qIq+?2H%hp#Q|IJ- zBJ)o+a;VhIJW}#>^EoG{DnwH^LRrG1K4O}46ABYqU+W+vW?0uNONt(mNCM84C2}D< zmYoYpmQ5!Ie6w-gW`Zj6@PzHF#E!WRBhkWHv(Fc|!ku7!Ea3u^d<|5Q`NS(SRvk7q z_ARJ6z+0@+LV%8ssqG5}zlzC(B45nyHUMmV;??!9icFZ+WVG-m8Gn6GYbvQCXA6>Z ze8Letftt`8Xt$YNbBF{zNx}#?*Br=&?0|DJ#t6E8`z&3!CQAU{{9bDgrh^O2Wxtc9 zv2bTo!Ok^@*-Zp;(ojXMIM*B|mOYB`yrFMkLUf@%$UQDQx#hg3%9OdyVIrGqBM8lb zQ}5v&DB7Y~^#|tYtTqRoo&^HVH3xDbJ8lHc;X1H6%vF(0WO!$zecgi7jpxH3MX@pO z1>Eq%ALeTWO;SZ}G9$(f;2y)J^Jmp(bnR(^X2S|rPB>s}emw3^94GCp}R}jt)mN zORZ)@KvNG%HRdQcTF4h<(FBdqCxVCDrX{3U_Tw(*Om0(UHeA;*S!YFeeG#|Z{ zKJ5jVjFskqVtB?4s@Bso?|@6dqAPn!q6HvVXvac(6-1iLp<3^x!Y#-g5G5swha*gG z6|0SRg9v`bTwS*xVtExlqO?rkP&tl}61ZC?lm>o&wt+c?uISZyec6B5CrYK`Bz|#V$ zY8HU>B^D0oY3a~!%yYTcbB>w99WutNdPPz33Qv)KK*wLjE1f#tniG++pJjAdD(;Psr-gE@s=zrP#r(Xd;(hO0RsXa044uQNuQ z;-SH8Xpzf(EWsQFV~tZLHq+~dQfxSruG{6xETx;#EIVuKXIC$sn|2om%*xCv8AX+T zHD4{EJabb_EN)QFLNoM*$+lH^quIA{C^Ivi*CI>r@P;dlpSS+nfC zG!>AFqydV?Wh}SN0a%9%|XVerH}EVqRjJ z*i=}u8&e8Z<}`Md-Yl<;yy5rMn}s_2v#J5hEH=j?K}jxkWLB-^#*~KytTnWl_fxLH zK0f6YA!)~fJ-H}U^rTf$}- zmJnzqHc_=rgqdR60f#v)1uQW>+W;P*BT$xF5oS}t9$>k-8bI@+ceQB>uQpBL>rGSm zdb1SXdb>hU18Q>QF|BEvy152A_B@xH%ACe<=D9pOD_`D{TPQYGoTu?32z$M-FF!1M zq6rbD$-%6E!NA9oO=T-wvMDszZ8ShGGYF5UT$5_v(6_>_ZWYd<2@0n-ijxPT`;GM> za#l?p)Le6@?az~krk)rHpz22G)bZXtQ(95o?w}Tc?Xc7HD%P^LqR>C^vE;XpwZ)66 z70+16Zew#}8`~v2+Z*Sx0q*+q?)v6h<%Cjn9S_3xsEUI&JQQ2;jxTzp$Tb)2-Co0) zvfGCrc9<`3iacD+#bys9w;Dk4)?Bb_4OKQJ<&~&HU?i>PNrw6wV5LN{NEP~1@I=5) zY+1t10c*S8Hkl&(qHfR;E3k2C-WVg-sid3B;lLkuU}*3L`X+d@Bsl|wbItyk1esT` z1wd@r*8_4-4JZN7a2equAeSW55V62)X>(Z>h~W>6%gg-?u*xyAM68zFhGT{&1k@tQ zI-pwPEX`z;Wz$wqQiy@BJuC$1^>^XEq?-%O1zQ4+uUng@+KfiU!wj&HYBM-C^pZl6 z*WGY@ML?S*hr6IE%Z99!@Xj2NR+~IV?#aA@6$@W+vBy1lwS{B5-&#<3wT~5UxdDBS zGWB3iVeGe>ZKOom6)QWBSJxoZTx|B^%AM<&Y;ztSv^~YIva1cC<|Z)?X#v>N+PEeF z-Eij>v6~3jl#Q^*OYH7kb>;drz+c8b6xbPr+Ei7YnhT#WB4v}ywIfmlJOi&$bQ5xP zR6n3w!Ocv>)PkBI3rz%)!enuSd}hxg zHl=_h(u|ATTq$;Yn4}OEo~ zA2^;zusX>39%H5l^H#6#Um=Cu&<0Mm`Q7?tF8fDMOtrP1C@~nJPj~Q1^F_~S^IR?g z)HB*VQ%V5!j5Y^^nWY=Rr(5iagqa|@-2y~?czp50&QhfXGFR8IH>v9&T~rGaout`3mn#h{645+U6kv({m5>&I zc{L86gQ*^72N<;i%>M{F*kd(qNMvfoe{%|j`;Mq;&;=Uz39XlSp zEPAG9X938GK0005y*UpBV9CI8f{A?|yBH3suLQMUlO%IrpdeCxlZ-wl*;KFvg%zq7 zkikP}U(?sxg$2QN+T^`C4`u%(-{`ryTG0rTi)J>ue~RN|wy!Fqu2J=+Wq32U!aWA+ zO@po-aPFuxnwpzTCwPm96}-2Ehf2;-Bhd=5@red?hG&~h(fiw4V&$zB;F#fA^UY`mQ1zW*dofo2TqdqeODGE?IYmphZh~c=7+TuL0dq})O_GZg=DOrIfVHM^n|3)h9cCj8 z8XFCn$V@2#%-Cp?tZRRh7dBWpHu`PGCk;8m$C&dQk;ePdl5D5jh$Ww@`nXm+tQ!qLSvE4J zk@dC4dmF*lnz$ruwTl&Co$=m^hpF*iqfFgeiA@SvG8$l$WMO!&@`Saiv1={a05r5? zPq1X$WC|CwLDR{Va>o@M8y+@l8mDyBTC#15o@nvhy>{mbC4Ab|oQJpG+}r1Kqrr_? z@XLB_c41<00plK{rf8#pOx{@|#}_#ykPC{AktUHO0?x+*$c5~HGahOoxCLx*5g050 zcRAUbF2}@7BKRg8sNts7#V`C zi7V`f`s~h(S5xEzg*!WY11yVe~fiAV=fIHjXL^DRm^468P?EXuf<-O`91bD(Xq*t z+STrD0mX3oAcU(aCICXh)6TOEpc;|+<)04UDoL!+|uSW?B<)_>xLe~aNFw_HdoKN7dO_{-E&(n zu3uPlFJ3ymwsn4WW7Cz}b~6&irY55Ow%v2hArl$i*;KHJDl*p`5?0PNhmtxvA4E>w z@S20^*l5F=!ws-OmFiqZvuh!F`DyA&Q!$z?E0Nb1} zyUPw6{jhNzyHmF$nGJ=F8-m0`?GTjr;Z zr?tyFjWo>U)~xOFj!r%S=k0RdK5kQE&)a26Qg+PSWxbS@-ji9A+U3(*mp0c{FT7N1 zSz;twmTDr}TbAbSawgQgGjF7nLt|MRDgZ6KsyIJ*-Y%CGt7!s7yWH4_KfXxGuy2>; zIDd|C0?ymz8?Rm7es1H!S}rH+9+=DX{xTLPyutd*Cfy!cIn0|;%6fLpo6)*@LHx!> zYVlYj-NuREvWp;_(LJxRnfa8D8Jba0Qb-5^=gsJi*No21&ILg(K2or6e6#bsdt1!c z{4;ZTZgwUYvIDyt-(+Uz>wg41FO}E9)5>=w5pUJ?I$*(^gzJa{sT}^w{mwJr+=y3x zGy2s>?)m&;a(!Xt3l|o?_!pi3T=~7veC6+t{o6;9>(lw2PxT#qI3zL31