From 9408cfb9f32ccad8b9922549fc783991a60bf282 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 | 9 +- framework/src/bcs/ArrayLowerDIntegratedBC.C | 18 +-- framework/src/bcs/HFEMDirichletBC.C | 9 +- framework/src/bcs/LowerDIntegratedBC.C | 18 +-- framework/src/dgkernels/ArrayDGLowerDKernel.C | 18 +-- framework/src/dgkernels/DGLowerDKernel.C | 18 +-- framework/src/interfaces/BlockRestrictable.C | 5 +- .../src/loops/ComputeFullJacobianThread.C | 2 +- framework/src/mesh/MooseMesh.C | 73 +++++++++-- 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 -> 82600 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, 331 insertions(+), 137 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..15513844a9db 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().size() > 0 || _mesh->boundaryLowerDBlocks().size() > 0) 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..12914f86f4bf 100644 --- a/framework/src/bcs/ArrayHFEMDirichletBC.C +++ b/framework/src/bcs/ArrayHFEMDirichletBC.C @@ -34,10 +34,11 @@ 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 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..eaf1dd1250b5 100644 --- a/framework/src/bcs/ArrayLowerDIntegratedBC.C +++ b/framework/src/bcs/ArrayLowerDIntegratedBC.C @@ -35,15 +35,15 @@ 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 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"); if (_lowerd_var.count() != _count) paramError("lowerd_variable", diff --git a/framework/src/bcs/HFEMDirichletBC.C b/framework/src/bcs/HFEMDirichletBC.C index 00d3beb37c6a..cd7c70009df6 100644 --- a/framework/src/bcs/HFEMDirichletBC.C +++ b/framework/src/bcs/HFEMDirichletBC.C @@ -29,10 +29,11 @@ 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 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..eafaf67d8385 100644 --- a/framework/src/bcs/LowerDIntegratedBC.C +++ b/framework/src/bcs/LowerDIntegratedBC.C @@ -34,15 +34,15 @@ 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 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"); // 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..32ed086e7b56 100644 --- a/framework/src/dgkernels/ArrayDGLowerDKernel.C +++ b/framework/src/dgkernels/ArrayDGLowerDKernel.C @@ -42,15 +42,15 @@ 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 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"); if (_lowerd_var.count() != _count) paramError("lowerd_variable", diff --git a/framework/src/dgkernels/DGLowerDKernel.C b/framework/src/dgkernels/DGLowerDKernel.C index ea28a581620e..dd2a3c21f95e 100644 --- a/framework/src/dgkernels/DGLowerDKernel.C +++ b/framework/src/dgkernels/DGLowerDKernel.C @@ -41,15 +41,15 @@ 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 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"); // 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..92e296c73b2c 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,6 +1374,8 @@ 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(); // TODO: Thread this! for (const auto & elem : getMesh().element_ptr_range()) @@ -1346,6 +1396,11 @@ MooseMesh::cacheInfo() std::pair, const Elem *>(pair, elem)); _lower_d_elem_to_higher_d_elem_side.insert( std::pair(elem, ip_side)); + + if (ip_elem->neighbor_ptr(ip_side)) + _lower_d_interior_blocks.insert(elem->subdomain_id()); + else + _lower_d_boundary_blocks.insert(elem->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..a48b03abc23d 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().size() > 0 || _mesh.boundaryLowerDBlocks().size() > 0) 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().size() > 0 || _mesh.boundaryLowerDBlocks().size() > 0) 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..7a6760b27f84 --- /dev/null +++ b/test/tests/kernels/hfem/3d-prism.i @@ -0,0 +1,119 @@ +[Mesh] + [square] + type = GeneratedMeshGenerator + nx = 3 + ny = 3 + nz = 3 + 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..52d3cb9f15f1c8258eac242623316f0933e9921c GIT binary patch literal 82600 zcmeHwd6*o>b>D&~z(W8@QKU$TZi@GTYjF@Dxm+x;3vwf17X%ihDavln&h+jMITz;u z7PJx{vLri{NJUqo*ivZAwrtt4D4M1$E3q7j)Dg#tEnBf2#fkifBotyF3+lM1N+3==|0uJMTn|d* zR|)4JKAwLYv!~@(2`qO!|8|vs9RCyeDxnj9hsvMAk4hl^MJ!h+hXxCEX(En{uS>*{ zxnvv}NXGFd$vDcA_qb5xSy8tg*$|IY@=C-hnIz&=xf5}!Y>7DDHhzy|xjN%$fJz3m zBR`C!ToiwcrK*ln_N(lzSXrpToc|hmlpU9f71;5}{lo#C54tc{EBeQRdZSpaAgJtc z&BU6|KF#!u_$wDm)m*b!S@dhwVx{RHuht9D%M-+=UdiWL(7J_U2{s!h7->wutejge zmRV0ww|rEx#pV*Cqu+b*Uj*=-+fmta%~r!mgLmZIX0ch~4RNu)uHRPJSgRKsWxv{L zZVD4FL_8I7fvU5|m#X>)ygx}U)(!HD*wHKTe_b%&zA0lt=Yr&78DAX+*7J93k-t{2 z)@YNGm-40moUrXsJf1y0iI}jA^842zyspydxF30vi}xd6xYlpMzqQD3-j8{E-A{gh zAM&nibl#7)AeVW+M-aId`OW(=Z@>3jM&9%CRagjpp-gDAa)EcMkJtKl<6oHfgo>d2 zGV)&1&ih2Kyq7UAbwpLrj(9`Sryygx7aCHv(H3-ULyL|1)c zN%pH07G!H;i5@JmJYroyIY($=IYqdA;)%F#B;G~r!+bV=y^bR`WzH6tT*4-h)$<_f8=LF5sQJA@tsA;<14x9{Xh0X(>x8 z-u*`>CRIzT@Xp$U{J1thgMVu^e{2(oN1Thb%YvlgojZ4Un)!LLyw8T-tzUyHc)T6) zUIkhYvf%IRmysXWHnY6!6o|J0bP1>%UIu@6`~-ie-^>$Y6*pZ(KeVD}cdYgPV+ z>f|U7_I+~kcWgUw-TgfNtu^}LZh61Uyc(B!WZug_auN6XFfcX=S8#jfU1DC1OZ^a6 zx2f=jZLHwF6Jf#amG@cZ)wt9vaW4nSW%B!{h!osjd0$~(jZ6I!_cb874BUIzCkSq@ zynC2e<8t6Y+$%tG5m)Q&o&x-W+bi#HFmG7@BYV(&h=L~W9_1Uk=i6{4ujuWyVVe4H zWsY48+k=$l%rfG58p9ve$s#AGIEfEm~B8`q3_2w$xDP|DUL!ti%L z#LT$*w9 z+qs*bk&#KtH?*Dh;GDnlK)m|Y7k|QQl5NvPN?&n2+KmU23uRV)2Sk47A^0(lUa0!@ zAdgX417QMQ3wi@cF2wyWwyC)G^FG8(5~9BUccA|VLc9L`zfZ{{zd@4vJ?}+>W1XU` zBt7XODK$UiN5;Q$R9z}y{7c45IC8%p{%$UduNZ% zOimnr$UibYIprUieQ@gVq<`e--pSd66Vo&P{b#yddX#_eOfUcR%-q!BnTZ3dRsTK8 zf5z*7K9>(yht70E&nVw~wb`tyWwne9XLxkIxVXfJ$klo^tF_fybXD$Bu&|odKeAA- zR#vkBhVB2b{a?-YKWzVp?f+`F|5ZKs85Vvhf%G*Yu5WRl4A&;Pzr@@Z!u=ew4`T!T z+&^(C=rYi&L6?KL?#?xGu4i+dn(NP8SLS*#*MYgV%QZa@bS;SMWY>YN2i*X=5p)yi zX3#C5TS2#hZU=1yZ31luZ2@sPmCK@BX50?q+BDak?*Q!r?FQWmx(md0!!gi3Anu#o z1L8Vaat)2;XB=g~wKcA*Q5F(UnJ|s4R@w_k77|;7Mjd9c&;~7WY zF`oA$2KC7M5rg_<9CgZg#!+YWJ|-=rmJDjJwf@7acMFCY@&Lh;cXTcvDAW!?ambN95&J zki?6Q7|;7M&D0U&ZZqOVM~vrvnP%#UakuMuQ%6#FSZ7Qf@qU{?5-&PpJh7Q(>WFcy zZxSy$Vm$B5G*d^6+oI!59hr62)X^yNZU=1xnL1+p4xMJ|h;ciOc+nB#S?`%<>WFc> zbiAn}-jDZX8)52*acnarUUbBGVl&Ou5#!jtNWAEX@w_k7OdTjxX zh;eMQBwlpHcw#Wk)Dh#@zDd03i1EBH(@Y&PZjX+a_V8SU&jGPLVjTVFfzAT016=@O z9{Smqv3@X4+BMoIgTJM1N;^Y)VIH<0%)_*+Kv#npM?c#K)_=y)&-yRzf(Jh_ zq>W%2GHQE`e`ft9qozvNc7X5 z=%>qkB>HJ<^wZ^UN%YgkD1XX(n|7tWV%!c8!*p5a#ZUP&Oqa5zpTB21E31d`CEqRQs(sY_smN@(WQRG&vb_AQZM4?ZyBb$NBgNK%9=LLy1_8zy;g_mmvqKa zf0E9+!!Y%55r|>>C7p59tE9gQe(FZTtV_&G9bF0{KEu?_t96)uNoO2?E9t~%m^yn6 zh++CAopE9d^z(kyspP*Be&(f3Tm>RN!&ie|r^C!EVPY^$o45wVu-FL0#G?$Qy?BET zvs~AM7^Yv+8Aq8(I_omSEXR!?hUu4d#!-fneiQsGlZ0=EpLr=u)@kB1%yQkT!}Lo! z-iCqBb0*KHt%>6diI@qUud`?1Utz8!w%r7YQo5T9Ynl5GgX%qwByF-+NQ1~E)L z2@{VpyAHGs#Ja{X%e7sH>6diIQD&09Lx)+8ogjwkmvqKahLX;D$}&lK7yQgiSxS9n znB`!d<5@tPRUTKes&#<(|3^T8UiN`Qy z$$HK(@gz*Vb3kIFbZH;w>9`AYnC?Y}Uu=%9*c*St-%7oup2g4KGEDb1hF@%luGkHK z!@O5%m$o5({+407v}Ul)QX z?~6gy(|Qo)EN!Q>iEO_pPqtf>pR|2f!%z82+r~c61DSFaxv~5rGbyL=&oZ(+EC6d&=a7yfc^>St)QO;y%Y3q(9eUO0euAY ztDxTi{cF(2L7xD<0Qxte&wxG$`cI%g0sUvtS3&&=uuDw^f;&mS_VA_`fbqffIbcSJ-dKc(v(0f4d1-&2i&q4nJ^e;gl0{tt{FM@s<^ij~qK+l4H z3-mnblb}z5ei!s_LB9|BEa)ZB=Rsct{d>?Kfc_BlN1!i*z5@Dl(0>8_CFm~DIA|Z} zUeFssv!ElO2SINF{U|603P43r8B_x`LB~Nq2KsT(Pl4VB`WeuZpr=6p4D>$G2S6VL z{T%4SpkDy}66jYzzXtks&~Jj)fX)J)4LS#OF6cba`JlC+b)XAC7lJMVT?~2^Xgz2H z=n~MSApUj`M1IH%`M3;pIp_+|4WOGqw}6KL+83k`(gQLFSHewlQu&8F!jiCUk!Rah~<~x({{*@ z*qhkdwTPq4S#FkD%ENNp0g`-TkF*u)o7f&mu9H;7nLo}_++Pv*N5bQg#`$L`-k z#A?1ux$=BLy-u?_rO;Z+u}@gF*RfSw@{WSD9TH-q>F>zQ2%L|ws+*6n?Wk1i6;ukZT47Z~laLaCY{{l*f+G{YqhQ#em+vYPp=N6ubjC&cvIjFSg1- zrP=UsF9$}wO`A4{f8}bm5%|rZ(cD9R3J)FG~bg5GRVBA0aBLlv*=lQKl@$pvR=kjb9)M+j48tr*w z-bA5rI4G%Bqp?)19n4j7i$Oh$4-YLsrCG0*@J{+-!$i{lF>iY2fr$h9=%Fl_cD&~m zE9C81kOIc2G4JL|)h`4KxmKykN1`4V0gKIIuH+Y*sgU*@ma4~bve|sKg-#&lvga^Y z0;FFF79sXyf&cjNfNx$z;vDl-8`h(s?PMsW)xeu8IF_xNe=LG&vue+w2E-fu=!M!q z;rD1y>@M4RRe3tWkR*XMnq$C z?2(j!U7J|@s^@tNxpJ{|g3u|~7)IVRvj?Xq4)iBtLUe7j?OU66n~kCu(0XG~V_f#_ zApje*^N-?iuu$Fx7`-ll*)4+t+&4ROWNu<+u74Re@;MxkxR_u~2VJz>wY`Cs`Me?3 zb8oH@bY+!%NJonIqRPHHT#+Xm@-R{)j(GyHI2X2AFD{ENC<=8sai{>}$`w4PKBZW$m5O=Hjiy~L-Kd&% zP-j4@2dX&im}j?;$zNwkF9aU*&>J?@dGV30b@T97JfzuQ44O^$rLBd9;&OAzT1eGiWds^SV?$jQLY6TEjwy-OR2djX%q%#UUbwbC|Kjb4$lh(<9i?ueQXo2$%8tR*=&RAHFQAhe$Ow$K(R=hIatX&^5^8MLoR|hgz~ShNQfOae-W+BN z4siB-@&Syg!sGBwsm1;L<={tOoD*;J6!*C_UpT^b6)ZbMb_PB^?0h&eN0af4Uk6oSa!Y#m9%k;avY_`~%4O zaJaVvizq>TK!c^OLnY3^hWogT8;3(CBf=~=y7^YKijVx@jCj4?6JH6jZx6eUyj+AjKco%tf*k!Fx9oG=C7%@3vI^Fa+u^c76BEMviE8BdSlC%1}t zo>W^zfZF%mcoMKy%;)gsE+;I-xHg+y;aIMcmq%KoO96AXqc7Sax6I}J)a5vPK(6<2 z@=(J;m@!W6wM@g#<(od%3i(049CnrYYX1n76V3to*vQpD+4;^!&ntXNNdlIdRav`m zC~~;@1S5C94bxT&iX3d*(YMu%fWr{5#%0Eo%bCM^fNe`rLEFfR!!3NpqrZzob~Ig0 zzg7z;iQ7=2&!Up3I})(LRk&ECs|VPUf<+~Kx}{$B^GgAqWBNpH1L{h*0xZp8*M|z1 z=F4nfQf$ZGY+`_=#cI9STq?I!Gvh~u$_}1<*ApW3GzO55#28i$R@TJ zU(~d3983<&_%7!_0aL1reN8DYLS|e`(c6N22-0V(EA6daZT7s9pJFVgmTg;`lIrae zl8XUu!4+w-mxmmtH8K%mC_}!C*53udcGDlB7@jdt?~d#bZzKlj*Raln9iIjs8krhC z{F?6*jx~kGyu5&~gtlsh_H1*6Ip%F2-MxMH&O1hT_f>8HE-V!M>;cZ0cf8g&d&met z>>e$0#UE=NW-F_Cw6DuqfId2GY%ahi4Rrkd8@K0$qb^KVQ{DLJQw>96Y#oseLNFT( zD`b%?VO1#?^d0rn3@dsRwNgX>l$pMpmYEn6FKn9`L*ERJ~=0u&C}@()zYDP4)dAzEy!(O0{U zn03tnw0#QTH!BzG#X_1!01Ks1d3|5loV>7rS)1HqwZc#Dp<`ChyVB(W(AI3aZh<0) zwnbJgypI=;FintB>=go{$apI$(uk7bO&e?Zuu0$8s00o5K%<{$R71~XtnXm6#imW0 zM!44R>wQ{Gj%|YITN6e`%pS5+lOlHx+_@8@g>4s&VjjIkxz^}$O8|xsv{~a3`{qjn zri?X}DDBNg17=Bx2ovE~QM zoB-o{QsL9Kz)^;_EjHV@@&Z1C$&#lj#kNmTjbZkfc+8ubdT?%b#%S`odLct>SZX%o z6R+(OJD80gU90bl8$jZan$gfjed!=9#S*3yIwX9|~8%-=I<%33}s^4+4uPG*A zhZVY1z3>4le3(KHC>!qm`CyTo{p^dZE1rMk;4F+WoKZ4lX%C2pn|ZuRU|__&wt%O> z=Hyd+a}3)!`mf|OBfoLLKmuT!8hm-LZ zQnw`<`@xduWzS$h9(W}u*PMjbS`BnLk*G2~l;LE@8@!Cwmx~NBcpdC+=UhGYz2J~V z7Y;?Ch^)hP;qWBv$5uvhtk_K1H11-DvV$Nzxgt?$mm+sv4RGOb748Tc_APJm#+ovP759445eUt{Y_3JK)_1OTv#T z6aj~WB}$=vjd_Q)%KFnge>w26O^>0-O+7v!egDyk$?YC9afY+Ym$v5)PjBln`S$03{?vZ)Dy4~BSycEIQniJhK!X@JaNcI+BZr3KM;DBM!;u4}(7v!a8EXWkps=4? zxO!UvfBAzRIoK{PP|LwWyNn5+O9ML`ImG9!Tsr!25+kEkbOBt2BZpx}KY&-Z6YW7P zIOfF;9P8O2^2ftDhn`~Unn8^mcsKxd2dchkrQ-)GWupFB~OZ9JE&{a^tVa%H)oI`y#V}RmDz1g<*R@ZlPqQYXnWS5AkZE zIJE6kOu&>)+4eO>0k-A~I5L4R8z*^ODPRYW4uI;*8u(OBt2`gnQw-a7@r z6^AHdlV*pi?KzwenqBZJ)Cb9O9h@Wp}PQ{-)1~q*gxLpar}9`TEt-@ z(r%&63Xf2rOzn%@?xdn#DeR=8J`14zYntOj9Pb}~KBYs&PE{`i%xuKQydQpupY84U z!p4pq%%)Js=VCDM3UIPd4ll~Jd9{~ep=g|q&;}Bpjr>>&&+ED+)TJ9mpN4^EN?}7D z?Tfsa!-w1kS(jx548v~XSe8cdiP)wGdk#xMWwE)$iQUvwX6*s_=3gnuEp+fEGPHOb zT4}k0L!ukSj?SrVHvNwcL$l$VbfpqcW=Z{A&X_kjwSVI1fquU;fwwZ-jf`+*y7s3f zq$@YV#Nh^3WhRHT7X79VM%&mcvuJGe#)-87Ub3Nh`no)?{sh()*n;Jj(^NC&J=7PF zWHJC%=hC5VdpnFcjGFxz=n4RF%nRO(ilIR@_JC@P9S)@Z<4b5|=`f$nDEnEhVqz}ho5n;VLgSm`LNt}{2bG9+Nvi5C0kl#}qK zq4b83yyL=G!$RQ&g_tB`f-?gA4c;$q_~zN{jFijXur;93VhNUVodKVS!L$QTfbEGE zABQpEk*bXc8Co^iBm8M>qJ^(U8+2RQOOz}e!-(yA~HvYVcn zp2G*7y}9Xw_(bp2A%ALSGV_H}Jn1c03#}{(y2wy{EV%Q~D`D5(uoue%R?1=p`>?C= z`X=nbwKtn3tlSy^Nw58LO4o@ha%fvq9&Kb8OvtN!kuiKN;-o}%NEM!^uoD5_ifb3# z9$ze$N1>9784QwB!OdX5^HCIuzLi zlrzpEC8MVFzZ6Wag+nZKE#V-*TxAj4m(1w`d&3Tk+<9(oRBAgOg~K*rqSSV9?IJ5( zm#T8+^ok6-UR?xNZ8fyP=(Cq0XKmsbv5XG$G?mq+#^QUdNKq-aTJE|cMZkUGV?vCN>;P8v z_%RdpX`xO5Iw|Eiw<{^@PKCBFGPb&@9hTUEk~ZA7A32agd3PcH9T6o0d&AZbk7#3g z_-}_C?E$;4_qdo1AC}qJbnq5N*vpWBRPx{e^^=rdlY+8EMXCw6Hyg2)Vg+A;>SS|; zfSpf#4l>k5#?HB9amlYchc&zn&^}K$V%6noG<>|pKIdv34;76FO9|8+kcUQN$Qe`d zii#WpCitTa?H|CyBI^^BIoRZEknA{hS4t_9g2cWSpxVPNn|67lnKO{FwuVom0w_&&V!&kI z9_oo8j-*vxOZ>Ap8x5#tR`8mEzkm<=MC~<;Y`f?m0rKNd3HA_l<}gxr0@E2V+4Kvp zvr-!QdJv>cW!nSRIBJi+m?exwE_UiBEUkMG(=43sTf$dV6Bex@?X{da!4e|cD@6cP z=3fcv02m*}Vdr2p4s!!+REg1nS%qR^5E5{~sO-Kwk7cl=49Iw_N^EmFqH>#K0ow*(WV??fdF&ysC ziEymAw#q4u7ptE2Gz$UwqL0Z}+ix#J#bL_Avc|-B9!ppb>D~w$sk+GgE>Mx5ENHV9 z0ooe2;xJ?A1yINkqp$v3ZM^1bo`})GphJ5YKA*%pJwMheTIq#{W)AfGDW2D!`^uWs zwQ9Z$g%8H7w#Q(2Gr+AK*xb==(sT?DgZP*-7ohtp0_~*tXRzmiVT+dYXh)z3(pET ztfHA{0E4f#X@-o52q4-;>fDrf+Gq-(nPTm_CzjDZ&w!N^=<6 zHzUkz1x3b!D30{OqBXyB7T^Rdwa;n=fRXi>46Z~e@orgAUHX*t;lD<7Dl(oYj@OIn zUuAIt>^j%(@**mPr*_c$!V%SC(r&Vh> zPJrDe@10~AP2L-viPy<&BEXc@0GA@GiRY}Hur4ij*O_bpTFztln8|ig3a7@=b@HU# z&H;{#3ca@3nxKT*4)RSX5&(c~V6p;BYN~QfOb;jE7bT+XAL}2n-H@dmx9AJzA!i5Ze~P zUt$9_wrLgc+Fmh_@q^x$&G$F(hd0rhcr@sH&_uwl2WQ$h4xDx#^gOnngUy&fG~e!Q z%PPyqB8T&bN@jKJX~COp_#6!rc)1iAi_3ZRDrhCRTAsL80Em~tI`LM*VLd3OeDvLg z!-0R3D|}wg*dpoDFZ1~r=Vmr}G<2i2V4hS(E!(~~iV^I|_*Q288+(T7xF{vI?2mUI zv3;X`xdN)?^jZ~LP4NOCEIfL6b^(;-m^NF6&;GOB10ITGH5$(Gs$VYU%JYR>1$)_7 z-aAuTgXotPPyGI};>)>wDL*8K!;YR(XkSrBe_(V5W{@Jacocv6gYM|@vc}xh;hBj8 z{*md)DgVIigHwkm{Ub;BPR<^jn4a;&&vr8t&8AC6_s@0@M-FY1q0gm()kBft$U%#8 zIC6+LEh!;N>a-s@gq>_k=~AQEFK^ns^6oeSKMg)mB@dU40Tw$d{~$cZU}-Qc(y)6E zmP9V${Dp)Xn^P=`fCC@fEab|?tPSJn0uIj3%pSxkM*T^c9Ckf2$K}CCPb=`0lJIqK z-bF5#2DaxD;%nZMzhROUN0E=BHocmT9qtu_L& zy}iA5?2-+KIxscl~n2n0H`e4rgBVXEb4vI{}ZiU%_%Ia;{#_WxNRCHn!<8 zRdouCVDr^-ty;m1UfMi{i(v;veEEawPP<&jhXWsDFLXPS%!R{&PXwv=kT3^yB3&t! zlFgDe^`HhzE%T5sLp%0G9`I14-#XuSKBHfrA4tbcX|1$ho;TSi;ILng&yTyb*u#FA zmejt6{jyohQvDQGuYP&&?9rLYiNg^0zln|Dn$RLYRO1D?0p~4u41g{z?7& zYxVsW_5EI_`Fp$O=WjLs4vqimI6vF}kIw(DD1TVqcj)q6s@-dJ`uDW^clv%Cbh*yc z{OpeIqsqG-`C0Dmf1}f{(!YDP#`(WG|BE``w{-j~`o33;3qPU!>zd#9>vG?#%k$e> zE}zony+X%7uicO8^tb5WO=+CB==^u+e9!6lr}TZFj_jcOc*3?n%*rqE(vM15T#3(e zC7-xi`6OQQiJKLlOqYD(N_x_j^rV~2C*fqi{=<44G5+GsU;UXEp1eGAU;6TAe*gY& zERMhU`b#I*U-dxbhVcvIp_|k3)7ssyZOSHRw{PSmj?{iP?`EcZZ;m`laSAXee zpB{hdmFq6N{Y!UicjtH*|90&@spD_d?oB%Wts3X&A~&6XU-JciR5=c5{(S8|tmEIP z-TQR>n&^pej8eV+}HTtrTx9Nw+VJ*M3qy8Lg`?mIK&^HY&r zLc9-llqasgxrBUvRP+C)4EaEA%3fX*xu5;oH?R4`z3&}=X{mhZEsLKYe`)74A9}}M z{_XhZKj59c;Wu9Z&IfDn`=3Sa=5+eEb^3?Kzwnh?zW&JiP1?Od$3Lg@y{vzCRJ(81 z=^xSQU(@&fo$-*5r?vYY9sijKC*yiWh$8uuB^$Bo*(Qpf+Q&i76IJIaUU z5`I1wrH69(UCqb)BDcGI9@Th{Xq=>cZq@u;ulada=lhoa9rZ|l-lWriQKyT1KB@V5 zDj}cW)A{~J|LzkSC(v$|e6G>_yfz`9uj}7El;G#r^?m+Y>*dptT)z0FPkrn!KKPCC PKlsEafBaqFz3~47=6U