diff --git a/include/base/dof_map.h b/include/base/dof_map.h index 6a61f304c4..a13fbd7a44 100644 --- a/include/base/dof_map.h +++ b/include/base/dof_map.h @@ -1764,8 +1764,9 @@ class DofMap : public DofMapBase, bool use_condensed_system = false) const; /** - * Describe whether the given variable group should be p-refined. If this API is not called with - * \p false, the default is to p-refine + * Set whether the given variable group should be p-refined on a + * p-refined Elem. This changes the FEType of the variable group to + * enable or disable p-refinement. */ void should_p_refine(unsigned int g, bool p_refine); @@ -2257,11 +2258,6 @@ class DofMap : public DofMapBase, * non-SCALAR variable v */ std::vector _first_old_scalar_df; - - /** - * A container of variable groups that we should not p-refine - */ - std::unordered_set _dont_p_refine; #endif #ifdef LIBMESH_ENABLE_CONSTRAINTS @@ -2568,14 +2564,9 @@ inline void DofMap::should_p_refine(const unsigned int g, const bool p_refine) { #ifdef LIBMESH_ENABLE_AMR - if (p_refine) - { - auto it = _dont_p_refine.find(g); - if (it != _dont_p_refine.end()) - _dont_p_refine.erase(it); - } - else - _dont_p_refine.insert(g); + VariableGroup & var = _variable_groups[g]; + var.type().p_refinement = p_refine; + #else libmesh_ignore(g, p_refine); #endif @@ -2585,7 +2576,8 @@ inline bool DofMap::should_p_refine(const unsigned int g) const { #ifdef LIBMESH_ENABLE_AMR - return !_dont_p_refine.count(g); + const VariableGroup & var = this->variable_group(g); + return var.type().p_refinement; #else libmesh_ignore(g); return false; @@ -2604,7 +2596,7 @@ bool DofMap::should_p_refine_var(const unsigned int var) const { #ifdef LIBMESH_ENABLE_AMR const auto vg = this->var_group_from_var_number(var); - return !_dont_p_refine.count(vg); + return this->should_p_refine(vg); #else libmesh_ignore(var); return false; @@ -2640,12 +2632,7 @@ void DofMap::_dof_indices (const Elem & elem, FEType fe_type = var.type(); - const bool add_p_level = -#ifdef LIBMESH_ENABLE_AMR - !_dont_p_refine.count(vg); -#else - false; -#endif + const bool add_p_level = fe_type.p_refinement; #ifdef DEBUG // The number of dofs per element is non-static for subdivision FE diff --git a/include/base/variable.h b/include/base/variable.h index 84735df415..cd67937c72 100644 --- a/include/base/variable.h +++ b/include/base/variable.h @@ -144,6 +144,13 @@ class Variable const FEType & type() const { return _type; } + /** + * \returns The \p FEType for this variable. Altering this while + * this Variable is already in use may be dangerous! + */ + FEType & type() + { return _type; } + /** * \returns The number of components of this variable if the \p FEFamily is \p SCALAR or if the * associated \p FEFieldType is \p TYPE_SCALAR. Otherwise this will error because determination of diff --git a/include/fe/fe_type.h b/include/fe/fe_type.h index f499f9270b..480ff57f3c 100644 --- a/include/fe/fe_type.h +++ b/include/fe/fe_type.h @@ -200,17 +200,24 @@ class FEType #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS /** - * Constructor. Optionally takes the approximation \p Order - * and the finite element family \p FEFamily + * Constructor. Optionally takes the approximation \p Order, + * the finite element family \p FEFamily. + * + * We can't set p_refinement in this argument list without + * conflicting with the \p ro parameter in the InfFE-compatible + * constructor version below, so we use the with_p_refinement API + * below to potentially convert an FEType with p-refinement (the + * default) to one without. */ FEType(const int o = 1, const FEFamily f = LAGRANGE) : order(o), - family(f) + family(f), + p_refinement(true) {} /** - * The approximation order of the element. + * The approximation order of the element (at 0 p-refinement level). */ OrderWrapper order; @@ -220,6 +227,12 @@ class FEType */ FEFamily family; + /** + * Whether or not the finite elements for this type increase their p + * refinement level on geometric elements of increased p level. + */ + bool p_refinement; + #else /** @@ -230,17 +243,24 @@ class FEType * are the same, as with the \p family and \p base_family. It must * be so, otherwise what we switch on would change when infinite * elements are not compiled in. + * + * We can't set p_refinement in this argument list in a way + * that matches the non-InfFE-enabled constructor version above, so + * we use the with_p_refinement API below to potentially convert an + * FEType with p-refinement (the default) to one without. */ FEType(const int o = 1, const FEFamily f = LAGRANGE, const int ro = THIRD, const FEFamily rf = JACOBI_20_00, - const InfMapType im = CARTESIAN) : + const InfMapType im = CARTESIAN + ) : order(o), radial_order(ro), family(f), radial_family(rf), - inf_map(im) + inf_map(im), + p_refinement(true) {} /** @@ -274,8 +294,25 @@ class FEType */ InfMapType inf_map; + /** + * Whether or not the finite elements for this type increase their p + * refinement level on geometric elements of increased p level. + */ + bool p_refinement; + #endif // ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS + /** + * "Fluent API" for constructing a non-default p_refinement, for + * easier compatibility between non-InfFE and InfFE code + */ + FEType with_p_refinement(bool p) + { + FEType returnval = *this; + returnval.p_refinement = p; + return returnval; + } + /** * Tests equality */ @@ -283,6 +320,7 @@ class FEType { return (order == f2.order && family == f2.family + && p_refinement == f2.p_refinement #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS && radial_order == f2.radial_order && radial_family == f2.radial_family @@ -308,7 +346,8 @@ class FEType return (order < f2.order); if (family != f2.family) return (family < f2.family); - + if (p_refinement != f2.p_refinement) + return (p_refinement < f2.p_refinement); #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS if (radial_order != f2.radial_order) return (radial_order < f2.radial_order); @@ -321,16 +360,19 @@ class FEType } /** - * \returns The default quadrature order for this \p FEType. The - * default quadrature order is calculated assuming a polynomial of - * degree \p order and is based on integrating the mass matrix for - * such an element exactly on affine elements. + * \returns The default quadrature order for this \p FEType, on an + * element without p refinement. The default quadrature order is + * calculated assuming a polynomial of degree \p order and is based + * on integrating the mass matrix for such an element exactly on + * affine elements. */ Order default_quadrature_order () const; /** * \returns A quadrature rule of appropriate type and order for this \p - * FEType. The default quadrature rule is based on integrating the mass + * FEType, on an element without p refinement. + * + * The default quadrature rule is based on integrating the mass * matrix for such an element exactly, with an additional power on * the basis order to help account for nonlinearities and/or * nonuniform coefficients. Higher or lower degree rules can be @@ -341,7 +383,9 @@ class FEType /** * \returns The default quadrature order for integrating unweighted - * basis functions of this \p FEType. + * basis functions of this \p FEType, on an element without p + * refinement. + * * The unweighted quadrature order is calculated assuming a * polynomial of degree \p order and is based on integrating the * shape functions for such an element exactly on affine elements. @@ -350,10 +394,12 @@ class FEType /** * \returns A quadrature rule of appropriate type and order for - * unweighted integration of this \p FEType. The default quadrature - * rule is based on integrating the shape functions on an affine - * element exactly. Higher or lower degree rules can be chosen by - * changing the extraorder parameter. + * unweighted integration of this \p FEType, on an element without p + * refinement + * + * The default quadrature rule is based on integrating the shape + * functions on an affine element exactly. Higher or lower degree + * rules can be chosen by changing the extraorder parameter. */ std::unique_ptr unweighted_quadrature_rule (const unsigned int dim, const int extraorder=0) const; @@ -392,7 +438,13 @@ struct hash // Old compiler versions seem to need the static_cast libMesh::boostcopy::hash_combine(seed, static_cast(fe_type.family)); libMesh::boostcopy::hash_combine(seed, fe_type.order._order); - return seed; + libMesh::boostcopy::hash_combine(seed, static_cast(fe_type.p_refinement)); +#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS + libMesh::boostcopy::hash_combine(seed, static_cast(fe_type.radial_family)); + libMesh::boostcopy::hash_combine(seed, fe_type.radial_order._order); + libMesh::boostcopy::hash_combine(seed, static_cast(fe_type.inf_map)); +#endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS + return seed; } }; } diff --git a/include/systems/generic_projector.h b/include/systems/generic_projector.h index 4721b2381b..4656c3ff30 100644 --- a/include/systems/generic_projector.h +++ b/include/systems/generic_projector.h @@ -1628,9 +1628,7 @@ void GenericProjector::SortAndCopy if (!var.active_on_subdomain(elem->subdomain_id())) continue; const FEType fe_type = var.type(); - const auto & dof_map = this->projector.system.get_dof_map(); - const auto vg = dof_map.var_group_from_var_number(v_num); - const bool add_p_level = dof_map.should_p_refine(vg); + const bool add_p_level = fe_type.p_refinement; // If we're trying to do projections on an isogeometric // analysis mesh, only the finite element nodes constrained @@ -2305,8 +2303,7 @@ void GenericProjector::ProjectVert base_fe_type, &elem, elem.get_node_index(&vertex), - this->projector.system.get_dof_map().should_p_refine( - this->projector.system.get_dof_map().var_group_from_var_number(var))), + base_fe_type.p_refinement), (unsigned int)(1 + dim)); const FValue val = @@ -2561,8 +2558,7 @@ void GenericProjector::ProjectSide continue; FEType fe_type = base_fe_type; - const auto & dof_map = system.get_dof_map(); - const bool add_p_level = dof_map.should_p_refine_var(var); + const bool add_p_level = base_fe_type.p_refinement; // This may be a p refined element fe_type.order = fe_type.order + add_p_level*elem.p_level(); @@ -2721,9 +2717,7 @@ void GenericProjector::ProjectInte context.get_element_fe( var, fe, dim ); FEType fe_type = base_fe_type; - const auto & dof_map = system.get_dof_map(); - const auto vg = dof_map.var_group_from_var_number(var); - const bool add_p_level = dof_map.should_p_refine(vg); + const bool add_p_level = fe_type.p_refinement; // This may be a p refined element fe_type.order = fe_type.order + add_p_level * elem->p_level(); diff --git a/include/systems/system.h b/include/systems/system.h index d4b9c50aaa..41e65a8797 100644 --- a/include/systems/system.h +++ b/include/systems/system.h @@ -1189,12 +1189,15 @@ class System : public ReferenceCountedObject, * for this system. Same as before, but assumes \p LAGRANGE * as default value for \p FEType.family. If \p active_subdomains is either * \p nullptr (the default) or points to an empty set, then it will be assumed - * that \p var has no subdomain restrictions + * that \p var has no subdomain restrictions. If \p p_refinement is + * false, then even when on an \p Elem with non-zero \p p_level() + * this variable will not be p-refined. */ unsigned int add_variable (std::string_view var, const Order order = FIRST, const FEFamily = LAGRANGE, - const std::set * const active_subdomains = nullptr); + const std::set * const active_subdomains = nullptr, + const bool p_refinement = true); /** * Adds the variables \p vars to the list of variables @@ -1235,12 +1238,15 @@ class System : public ReferenceCountedObject, * for this system. Same as before, but assumes \p LAGRANGE * as default value for \p FEType.family. If \p active_subdomains is either * \p nullptr (the default) or points to an empty set, then it will be assumed that - * \p var has no subdomain restrictions + * \p var has no subdomain restrictions. If \p p_refinement is + * false, then even when on an \p Elem with non-zero \p p_level() + * this variable will not be p-refined. */ unsigned int add_variables (const std::vector & vars, const Order order = FIRST, const FEFamily = LAGRANGE, - const std::set * const active_subdomains = nullptr); + const std::set * const active_subdomains = nullptr, + const bool p_refinement = true); /** * Return a constant reference to \p Variable \p var. diff --git a/src/base/dof_map.C b/src/base/dof_map.C index 56d284dbc6..c80ba82ac6 100644 --- a/src/base/dof_map.C +++ b/src/base/dof_map.C @@ -579,12 +579,7 @@ void DofMap::reinit const unsigned int n_var_in_group = vg_description.n_variables(); const FEType & base_fe_type = vg_description.type(); - const bool add_p_level = -#ifdef LIBMESH_ENABLE_AMR - !_dont_p_refine.count(vg); -#else - false; -#endif + const bool add_p_level = base_fe_type.p_refinement; // Don't need to loop over elements for a SCALAR variable // Just increment _n_SCALAR_dofs @@ -2510,12 +2505,7 @@ void DofMap::_node_dof_indices (const Elem & elem, const bool extra_hanging_dofs = FEInterface::extra_hanging_dofs(fe_type); - const bool add_p_level = -#ifdef LIBMESH_ENABLE_AMR - !_dont_p_refine.count(vg); -#else - false; -#endif + const bool add_p_level = fe_type.p_refinement; // There is a potential problem with h refinement. Imagine a // quad9 that has a linear FE on it. Then, on the hanging side, @@ -2765,12 +2755,8 @@ void DofMap::old_dof_indices (const Elem * const elem, // or just for the specified variable FEType fe_type = var.type(); - const bool add_p_level = -#ifdef LIBMESH_ENABLE_AMR - !_dont_p_refine.count(vg); -#else - false; -#endif + const bool add_p_level = fe_type.p_refinement; + // Increase the polynomial order on p refined elements, // but make sure you get the right polynomial order for // the OLD degrees of freedom diff --git a/src/fe/fe_base.C b/src/fe/fe_base.C index 217c0f0b6d..86f9e99ce7 100644 --- a/src/fe/fe_base.C +++ b/src/fe/fe_base.C @@ -1550,7 +1550,7 @@ FEGenericBase::compute_proj_constraints (DofConstraints & constraint const Variable & var = dof_map.variable(variable_number); const FEType & base_fe_type = var.type(); - const bool add_p_level = dof_map.should_p_refine_var(variable_number); + const bool add_p_level = base_fe_type.p_refinement; // Construct FE objects for this element and its neighbors. std::unique_ptr> my_fe diff --git a/src/fe/fe_lagrange.C b/src/fe/fe_lagrange.C index 88cf0fb679..87407502e8 100644 --- a/src/fe/fe_lagrange.C +++ b/src/fe/fe_lagrange.C @@ -907,7 +907,7 @@ void lagrange_compute_constraints (DofConstraints & constraints, const Variable & var = dof_map.variable(variable_number); const FEType fe_type = var.type(); - const bool add_p_level = dof_map.should_p_refine_var(variable_number); + const bool add_p_level = fe_type.p_refinement; // Pull objects out of the loop to reduce heap operations std::vector my_dof_indices, parent_dof_indices; diff --git a/src/systems/equation_systems.C b/src/systems/equation_systems.C index 90b80d5abb..9f8227283c 100644 --- a/src/systems/equation_systems.C +++ b/src/systems/equation_systems.C @@ -833,8 +833,7 @@ EquationSystems::build_parallel_solution_vector(const std::set * sy const FEType & fe_type = system.variable_type(var); const Variable & var_description = system.variable(var); unsigned int n_vec_dim = FEInterface::n_vec_dim( sys_ptr->get_mesh(), fe_type ); - const auto vg = dof_map.var_group_from_var_number(var); - const bool add_p_level = dof_map.should_p_refine(vg); + const bool add_p_level = fe_type.p_refinement; for (const auto & elem : _mesh.active_local_element_ptr_range()) { @@ -1374,8 +1373,7 @@ EquationSystems::build_discontinuous_solution_vector const FEType & fe_type = system->variable_type(var); const Variable & var_description = system->variable(var); - const auto vg = dof_map.var_group_from_var_number(var); - const bool add_p_level = dof_map.should_p_refine(vg); + const bool add_p_level = fe_type.p_refinement; unsigned int nn=0; @@ -1383,7 +1381,7 @@ EquationSystems::build_discontinuous_solution_vector { if (var_description.active_on_subdomain(elem->subdomain_id())) { - system->get_dof_map().dof_indices (elem, dof_indices, var); + dof_map.dof_indices (elem, dof_indices, var); soln_coeffs.resize(dof_indices.size()); @@ -1442,7 +1440,7 @@ EquationSystems::build_discontinuous_solution_vector { if (var_description.active_on_subdomain(elem->subdomain_id())) { - system->get_dof_map().dof_indices (elem, dof_indices, var); + dof_map.dof_indices (elem, dof_indices, var); soln_coeffs.resize(dof_indices.size()); @@ -1486,7 +1484,7 @@ EquationSystems::build_discontinuous_solution_vector var_description.active_on_subdomain(neigh->subdomain_id())) { std::vector neigh_indices; - system->get_dof_map().dof_indices (neigh, neigh_indices, var); + dof_map.dof_indices (neigh, neigh_indices, var); std::vector neigh_coeffs(neigh_indices.size()); for (auto i : index_range(neigh_indices)) diff --git a/src/systems/fem_context.C b/src/systems/fem_context.C index 3009e7b822..7baa6098a6 100644 --- a/src/systems/fem_context.C +++ b/src/systems/fem_context.C @@ -244,8 +244,7 @@ void FEMContext::init_internal_data(const System & sys) unsigned int i) { FEType fe_type = sys.variable_type(i); - const auto & dof_map = sys.get_dof_map(); - const bool add_p_level = dof_map.should_p_refine(dof_map.var_group_from_var_number(i)); + const bool add_p_level = fe_type.p_refinement; auto & element_fe = _element_fe[dim][fe_type]; auto & side_fe = _side_fe[dim][fe_type]; diff --git a/src/systems/system.C b/src/systems/system.C index e2a7f0f30d..c85a0f0b03 100644 --- a/src/systems/system.C +++ b/src/systems/system.C @@ -1349,10 +1349,11 @@ unsigned int System::add_variable (std::string_view var, unsigned int System::add_variable (std::string_view var, const Order order, const FEFamily family, - const std::set * const active_subdomains) + const std::set * const active_subdomains, + const bool p_refinement) { return this->add_variable(var, - FEType(order, family), + FEType(order, family).with_p_refinement(p_refinement), active_subdomains); } @@ -1370,10 +1371,11 @@ unsigned int System::add_variables (const std::vector & vars, unsigned int System::add_variables (const std::vector & vars, const Order order, const FEFamily family, - const std::set * const active_subdomains) + const std::set * const active_subdomains, + const bool p_refinement) { return this->add_variables(vars, - FEType(order, family), + FEType(order, family).with_p_refinement(p_refinement), active_subdomains); }