From 8a35ffe013217aa0744990828679b6723e1e7c4a Mon Sep 17 00:00:00 2001 From: Yaqi Wang Date: Tue, 25 Apr 2023 23:24:18 -0500 Subject: [PATCH 01/12] enable p refinement #24141 --- framework/include/base/Adaptivity.h | 7 +++++++ framework/src/actions/AdaptivityAction.C | 3 +++ framework/src/actions/SetAdaptivityOptionsAction.C | 3 +++ framework/src/base/Adaptivity.C | 13 +++++++++++++ test/tests/dgkernels/2d_diffusion_dg/tests | 9 +++++++++ 5 files changed, 35 insertions(+) diff --git a/framework/include/base/Adaptivity.h b/framework/include/base/Adaptivity.h index f87ad8a7dcbc..e67802f17dd5 100644 --- a/framework/include/base/Adaptivity.h +++ b/framework/include/base/Adaptivity.h @@ -128,6 +128,11 @@ class Adaptivity : public ConsoleStreamInterface, */ void setRecomputeMarkersFlag(const bool flag) { _recompute_markers_during_cycles = flag; } + /** + * Switch from h-refinement to p-refinement + */ + void switchHToPRefinement(); + /** * Adapts the mesh based on the error estimator used * @@ -302,6 +307,8 @@ class Adaptivity : public ConsoleStreamInterface, /// Stores pointers to ErrorVectors associated with indicator field names std::map> _indicator_field_to_error_vector; + + bool _p_refinement_flag = false; }; template diff --git a/framework/src/actions/AdaptivityAction.C b/framework/src/actions/AdaptivityAction.C index 5cb1c2dfae65..77a92a5af009 100644 --- a/framework/src/actions/AdaptivityAction.C +++ b/framework/src/actions/AdaptivityAction.C @@ -79,6 +79,7 @@ AdaptivityAction::validParams() "show_initial_progress", true, "Show the progress of the initial adaptivity"); params.addParam( "recompute_markers_during_cycles", false, "Recompute markers during adaptivity cycles"); + params.addParam("switch_h_to_p_refinement", false, "True to perform p-refinement"); return params; } @@ -202,6 +203,8 @@ AdaptivityAction::act() adapt.setTimeActive(getParam("start_time"), getParam("stop_time")); adapt.setInterval(getParam("interval")); + if (getParam("switch_h_to_p_refinement")) + adapt.switchHToPRefinement(); } } diff --git a/framework/src/actions/SetAdaptivityOptionsAction.C b/framework/src/actions/SetAdaptivityOptionsAction.C index 8bb4a840d70f..261d641a8fe3 100644 --- a/framework/src/actions/SetAdaptivityOptionsAction.C +++ b/framework/src/actions/SetAdaptivityOptionsAction.C @@ -49,6 +49,7 @@ SetAdaptivityOptionsAction::validParams() "The number of adaptive steps to use when on each timestep during a Transient simulation."); params.addParam( "recompute_markers_during_cycles", false, "Recompute markers during adaptivity cycles"); + params.addParam("switch_h_to_p_refinement", false, "True to perform p-refinement"); return params; } @@ -137,5 +138,7 @@ SetAdaptivityOptionsAction::act() adapt.setInterval(getParam("interval")); adapt.setRecomputeMarkersFlag(getParam("recompute_markers_during_cycles")); + if (getParam("switch_h_to_p_refinement")) + adapt.switchHToPRefinement(); } } diff --git a/framework/src/base/Adaptivity.C b/framework/src/base/Adaptivity.C index d1c30fadf1ca..f9bac0476c1a 100644 --- a/framework/src/base/Adaptivity.C +++ b/framework/src/base/Adaptivity.C @@ -220,6 +220,9 @@ Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/) if (distributed_adaptivity) _mesh_refinement->make_flags_parallel_consistent(); + if (_p_refinement_flag) + _mesh_refinement->switch_h_to_p_refinement(); + // Perform refinement and coarsening mesh_changed = _mesh_refinement->refine_and_coarsen_elements(); @@ -229,6 +232,10 @@ Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/) // we sync them here. if (distributed_adaptivity) _displaced_mesh_refinement->make_flags_parallel_consistent(); + + if (_p_refinement_flag) + _displaced_mesh_refinement->switch_h_to_p_refinement(); + #ifndef NDEBUG bool displaced_mesh_changed = #endif @@ -378,4 +385,10 @@ Adaptivity::isAdaptivityDue() return _mesh_refinement_on && (_start_time <= _t && _t < _stop_time) && _step % _interval == 0; } +void +Adaptivity::switchHToPRefinement() +{ + _p_refinement_flag = true; +} + #endif // LIBMESH_ENABLE_AMR diff --git a/test/tests/dgkernels/2d_diffusion_dg/tests b/test/tests/dgkernels/2d_diffusion_dg/tests index 8c0b8c9d85f7..60ea9696eec5 100644 --- a/test/tests/dgkernels/2d_diffusion_dg/tests +++ b/test/tests/dgkernels/2d_diffusion_dg/tests @@ -9,6 +9,15 @@ requirement = "The system shall support solving 2D diffusion using the discontinuous Galerkin " "method." [] + [p_refinement_test] + type = 'Exodiff' + input = '2d_diffusion_dg_test.i' + cli_args = 'Outputs/file_base=p_refinement Executioner/Adaptivity/switch_h_to_p_refinement=true' + exodiff = 'p_refinement.e-s003' + max_parallel = 1 + requirement = "The system shall support solving 2D diffusion using the discontinuous Galerkin " + "method with p-refinement." + [] [constMonomial] type = 'Exodiff' input = '2d_diffusion_dg_test.i' From 0c2cb9fd1a3a52a6a6676f439c671abb70f2f027 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Wed, 31 May 2023 15:51:13 -0700 Subject: [PATCH 02/12] subproblem -> fe_problem in Adaptivity class --- framework/include/base/Adaptivity.h | 4 +-- framework/src/base/Adaptivity.C | 52 ++++++++++++++--------------- 2 files changed, 28 insertions(+), 28 deletions(-) diff --git a/framework/include/base/Adaptivity.h b/framework/include/base/Adaptivity.h index e67802f17dd5..2e8fc651fe72 100644 --- a/framework/include/base/Adaptivity.h +++ b/framework/include/base/Adaptivity.h @@ -50,7 +50,7 @@ class Adaptivity : public ConsoleStreamInterface, public libMesh::ParallelObject { public: - Adaptivity(FEProblemBase & subproblem); + Adaptivity(FEProblemBase & fe_problem); virtual ~Adaptivity(); /** @@ -250,7 +250,7 @@ class Adaptivity : public ConsoleStreamInterface, bool isAdaptivityDue(); protected: - FEProblemBase & _subproblem; + FEProblemBase & _fe_problem; MooseMesh & _mesh; /// on/off flag reporting if the adaptivity is being used diff --git a/framework/src/base/Adaptivity.C b/framework/src/base/Adaptivity.C index f9bac0476c1a..2c20b782e9f3 100644 --- a/framework/src/base/Adaptivity.C +++ b/framework/src/base/Adaptivity.C @@ -28,19 +28,19 @@ #ifdef LIBMESH_ENABLE_AMR -Adaptivity::Adaptivity(FEProblemBase & subproblem) - : ConsoleStreamInterface(subproblem.getMooseApp()), - PerfGraphInterface(subproblem.getMooseApp().perfGraph(), "Adaptivity"), - ParallelObject(subproblem.getMooseApp()), - _subproblem(subproblem), - _mesh(_subproblem.mesh()), +Adaptivity::Adaptivity(FEProblemBase & fe_problem) + : ConsoleStreamInterface(fe_problem.getMooseApp()), + PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "Adaptivity"), + ParallelObject(fe_problem.getMooseApp()), + _fe_problem(fe_problem), + _mesh(_fe_problem.mesh()), _mesh_refinement_on(false), _initialized(false), _initial_steps(0), _steps(0), _print_mesh_changed(false), - _t(_subproblem.time()), - _step(_subproblem.timeStep()), + _t(_fe_problem.time()), + _step(_fe_problem.timeStep()), _interval(1), _start_time(-std::numeric_limits::max()), _stop_time(std::numeric_limits::max()), @@ -59,12 +59,12 @@ Adaptivity::init(unsigned int steps, unsigned int initial_steps) // Get the pointer to the DisplacedProblem, this cannot be done at construction because // DisplacedProblem // does not exist at that point. - _displaced_problem = _subproblem.getDisplacedProblem(); + _displaced_problem = _fe_problem.getDisplacedProblem(); _mesh_refinement = std::make_unique(_mesh); _error = std::make_unique(); - EquationSystems & es = _subproblem.es(); + EquationSystems & es = _fe_problem.es(); es.parameters.set("adaptivity") = true; _initial_steps = initial_steps; @@ -72,7 +72,7 @@ Adaptivity::init(unsigned int steps, unsigned int initial_steps) _mesh_refinement_on = true; _mesh_refinement->set_periodic_boundaries_ptr( - _subproblem.getNonlinearSystemBase().dofMap().get_periodic_boundaries()); + _fe_problem.getNonlinearSystemBase().dofMap().get_periodic_boundaries()); // displaced problem if (_displaced_problem != nullptr) @@ -88,7 +88,7 @@ Adaptivity::init(unsigned int steps, unsigned int initial_steps) // i.e. neighbors across periodic boundaries, for the purposes of // refinement. _displaced_mesh_refinement->set_periodic_boundaries_ptr( - _subproblem.getNonlinearSystemBase().dofMap().get_periodic_boundaries()); + _fe_problem.getNonlinearSystemBase().dofMap().get_periodic_boundaries()); // TODO: This is currently an empty function on the DisplacedProblem... could it be removed? _displaced_problem->initAdaptivity(); @@ -141,7 +141,7 @@ Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/) std::vector serialized_solution; - auto distributed_mesh = dynamic_cast(&_subproblem.mesh().getMesh()); + auto distributed_mesh = dynamic_cast(&_fe_problem.mesh().getMesh()); // Element range std::unique_ptr all_elems; @@ -152,7 +152,7 @@ Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/) if (distributed_mesh && !distributed_mesh->is_serial_on_zero()) { // We update here to make sure local solution is up-to-date - _subproblem.getAuxiliarySystem().update(); + _fe_problem.getAuxiliarySystem().update(); distributed_adaptivity = true; // We can not assume that geometric and algebraic ghosting functors cover @@ -166,13 +166,13 @@ Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/) // After we set markers for all local elements, we will do a global // communication to sync markers for ghosted elements from their owners. all_elems = std::make_unique( - _subproblem.mesh().getMesh().active_local_elements_begin(), - _subproblem.mesh().getMesh().active_local_elements_end()); + _fe_problem.mesh().getMesh().active_local_elements_begin(), + _fe_problem.mesh().getMesh().active_local_elements_end()); } else // This is not scalable but it might be useful for small-size problems { - _subproblem.getAuxiliarySystem().solution().close(); - _subproblem.getAuxiliarySystem().solution().localize(serialized_solution); + _fe_problem.getAuxiliarySystem().solution().close(); + _fe_problem.getAuxiliarySystem().solution().localize(serialized_solution); distributed_adaptivity = false; // For a replicated mesh or a serialized distributed mesh, the solution @@ -182,20 +182,20 @@ Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/) // We might not care about much since a replicated mesh // or a serialized distributed mesh is not scalable anyway. all_elems = - std::make_unique(_subproblem.mesh().getMesh().active_elements_begin(), - _subproblem.mesh().getMesh().active_elements_end()); + std::make_unique(_fe_problem.mesh().getMesh().active_elements_begin(), + _fe_problem.mesh().getMesh().active_elements_end()); } FlagElementsThread fet( - _subproblem, serialized_solution, _max_h_level, marker_name, !distributed_adaptivity); + _fe_problem, serialized_solution, _max_h_level, marker_name, !distributed_adaptivity); Threads::parallel_reduce(*all_elems, fet); - _subproblem.getAuxiliarySystem().solution().close(); + _fe_problem.getAuxiliarySystem().solution().close(); } } else { // Compute the error for each active element - _error_estimator->estimate_error(_subproblem.getNonlinearSystemBase().system(), *_error); + _error_estimator->estimate_error(_fe_problem.getNonlinearSystemBase().system(), *_error); // Flag elements to be refined and coarsened _mesh_refinement->flag_elements_by_error_fraction(*_error); @@ -310,7 +310,7 @@ Adaptivity::uniformRefineWithProjection() if (_displaced_problem) displaced_mesh_refinement.uniformly_refine(1); - _subproblem.meshChanged(); + _fe_problem.meshChanged(); } } @@ -371,12 +371,12 @@ Adaptivity::updateErrorVectors() } // Fill the vectors with the local contributions - UpdateErrorVectorsThread uevt(_subproblem, _indicator_field_to_error_vector); + UpdateErrorVectorsThread uevt(_fe_problem, _indicator_field_to_error_vector); Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), uevt); // Now sum across all processors for (const auto & it : _indicator_field_to_error_vector) - _subproblem.comm().sum((std::vector &)*(it.second)); + _fe_problem.comm().sum((std::vector &)*(it.second)); } bool From a04dfe70a8ffacd6aa2b8999a9a783f6cb38eb2b Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 1 Jun 2023 10:01:08 -0700 Subject: [PATCH 03/12] Refactor helpers in Assembly to accomodate p-refinement --- framework/include/base/Assembly.h | 42 ++- framework/include/problems/DisplacedProblem.h | 5 + framework/include/problems/FEProblemBase.h | 5 + framework/src/base/Adaptivity.C | 1 + framework/src/base/Assembly.C | 270 ++++++++++++++---- framework/src/problems/DisplacedProblem.C | 8 + framework/src/problems/FEProblemBase.C | 11 + 7 files changed, 280 insertions(+), 62 deletions(-) diff --git a/framework/include/base/Assembly.h b/framework/include/base/Assembly.h index c96b61d8c847..2912f0542514 100644 --- a/framework/include/base/Assembly.h +++ b/framework/include/base/Assembly.h @@ -1817,6 +1817,11 @@ class Assembly */ const Elem * const & msmElem() const { return _msm_elem; } + /** + * Indicate that we have p-refinement + */ + void havePRefinement(); + private: /** * Just an internal helper function to reinit the volume FE objects. @@ -2150,6 +2155,11 @@ class Assembly return _jacobian_block_nonlocal_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar]; } + /** + * request phi, dphi, xyz, JxW, etc. data through the FE helper functions + */ + void helpersRequestData(); + SystemBase & _sys; SubProblem & _subproblem; @@ -2194,6 +2204,28 @@ class Assembly unsigned int _mesh_dimension; + /// The finite element type of the FE helper classes. The helper class gives us data like JxW, the + /// physical quadrature point locations, etc. + const FEType _helper_type; + + /// Whether user code requested a \p FEType the same as our \p _helper_type + mutable bool _user_added_fe_of_helper_type; + mutable bool _user_added_fe_face_of_helper_type; + mutable bool _user_added_fe_face_neighbor_of_helper_type; + mutable bool _user_added_fe_neighbor_of_helper_type; + mutable bool _user_added_fe_lower_of_helper_type; + + /// Containers for holding unique FE helper types if we are doing p-refinement. If we are not + /// doing p-refinement then the helper data is owned by the \p _fe data members + std::vector> _unique_fe_helper; + std::vector> _unique_fe_face_helper; + std::vector> _unique_fe_face_neighbor_helper; + std::vector> _unique_fe_neighbor_helper; + std::vector> _unique_fe_lower_helper; + + /// Whether we are currently building the FE classes for the helpers + bool _building_helpers; + /// The XFEM controller std::shared_ptr _xfem; @@ -2222,7 +2254,7 @@ class Assembly /// Each dimension's actual vector fe objects indexed on type mutable std::map> _vector_fe; /// Each dimension's helper objects - std::map _holder_fe_helper; + std::map _holder_fe_helper; /// The current helper object for transforming coordinates FEBase * _current_fe_helper; /// The current current quadrature rule being used (could be either volumetric or arbitrary - for dirac kernels) @@ -2334,7 +2366,7 @@ class Assembly /// types of vector finite elements mutable std::map> _vector_fe_face; /// Each dimension's helper objects - std::map _holder_fe_face_helper; + std::map _holder_fe_face_helper; /// helper object for transforming coordinates FEBase * _current_fe_face_helper; /// quadrature rule used on faces @@ -2370,15 +2402,15 @@ class Assembly mutable std::map> _vector_fe_face_neighbor; /// Each dimension's helper objects - std::map _holder_fe_neighbor_helper; - std::map _holder_fe_face_neighbor_helper; + std::map _holder_fe_neighbor_helper; + std::map _holder_fe_face_neighbor_helper; /// FE objects for lower dimensional elements mutable std::map> _fe_lower; /// Vector FE objects for lower dimensional elements mutable std::map> _vector_fe_lower; /// helper object for transforming coordinates for lower dimensional element quadrature points - std::map _holder_fe_lower_helper; + std::map _holder_fe_lower_helper; /// quadrature rule used on neighbors QBase * _current_qrule_neighbor; diff --git a/framework/include/problems/DisplacedProblem.h b/framework/include/problems/DisplacedProblem.h index 874640e74512..56dcc638a2fa 100644 --- a/framework/include/problems/DisplacedProblem.h +++ b/framework/include/problems/DisplacedProblem.h @@ -354,6 +354,11 @@ class DisplacedProblem : public SubProblem virtual const std::vector & currentResidualVectorTags() const override; + /** + * Indicate that we have p-refinement + */ + void havePRefinement(); + protected: FEProblemBase & _mproblem; MooseMesh & _mesh; diff --git a/framework/include/problems/FEProblemBase.h b/framework/include/problems/FEProblemBase.h index 25fe5514e39d..b89be1520a59 100644 --- a/framework/include/problems/FEProblemBase.h +++ b/framework/include/problems/FEProblemBase.h @@ -2102,6 +2102,11 @@ class FEProblemBase : public SubProblem, public Restartable */ void clearCurrentResidualVectorTags(); + /** + * Indicate that we have p-refinement + */ + void havePRefinement(); + protected: /// Create extra tagged vectors and matrices void createTagVectors(); diff --git a/framework/src/base/Adaptivity.C b/framework/src/base/Adaptivity.C index 2c20b782e9f3..064a401879cd 100644 --- a/framework/src/base/Adaptivity.C +++ b/framework/src/base/Adaptivity.C @@ -389,6 +389,7 @@ void Adaptivity::switchHToPRefinement() { _p_refinement_flag = true; + _fe_problem.havePRefinement(); } #endif // LIBMESH_ENABLE_AMR diff --git a/framework/src/base/Assembly.C b/framework/src/base/Assembly.C index 6bb20c57e486..d4f6476cb4be 100644 --- a/framework/src/base/Assembly.C +++ b/framework/src/base/Assembly.C @@ -86,6 +86,13 @@ Assembly::Assembly(SystemBase & sys, THREAD_ID tid) _tid(tid), _mesh(sys.mesh()), _mesh_dimension(_mesh.dimension()), + _helper_type(_mesh.hasSecondOrderElements() ? SECOND : FIRST, LAGRANGE), + _user_added_fe_of_helper_type(false), + _user_added_fe_face_of_helper_type(false), + _user_added_fe_face_neighbor_of_helper_type(false), + _user_added_fe_neighbor_of_helper_type(false), + _user_added_fe_lower_of_helper_type(false), + _building_helpers(false), _current_qrule(nullptr), _current_qrule_volume(nullptr), _current_qrule_arbitrary(nullptr), @@ -131,55 +138,39 @@ Assembly::Assembly(SystemBase & sys, THREAD_ID tid) _calculate_curvatures(false), _calculate_ad_coord(false) { - Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST; + const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST; + _building_helpers = true; // Build fe's for the helpers buildFE(FEType(helper_order, LAGRANGE)); buildFaceFE(FEType(helper_order, LAGRANGE)); buildNeighborFE(FEType(helper_order, LAGRANGE)); buildFaceNeighborFE(FEType(helper_order, LAGRANGE)); buildLowerDFE(FEType(helper_order, LAGRANGE)); + _building_helpers = false; // Build an FE helper object for this type for each dimension up to the dimension of the current // mesh for (unsigned int dim = 0; dim <= _mesh_dimension; dim++) { - _holder_fe_helper[dim] = &_fe[dim][FEType(helper_order, LAGRANGE)]; - (*_holder_fe_helper[dim])->get_phi(); - (*_holder_fe_helper[dim])->get_dphi(); - (*_holder_fe_helper[dim])->get_xyz(); - (*_holder_fe_helper[dim])->get_JxW(); - - _holder_fe_face_helper[dim] = &_fe_face[dim][FEType(helper_order, LAGRANGE)]; - (*_holder_fe_face_helper[dim])->get_phi(); - (*_holder_fe_face_helper[dim])->get_dphi(); - (*_holder_fe_face_helper[dim])->get_xyz(); - (*_holder_fe_face_helper[dim])->get_JxW(); - (*_holder_fe_face_helper[dim])->get_normals(); - - _holder_fe_face_neighbor_helper[dim] = &_fe_face_neighbor[dim][FEType(helper_order, LAGRANGE)]; - (*_holder_fe_face_neighbor_helper[dim])->get_xyz(); - (*_holder_fe_face_neighbor_helper[dim])->get_JxW(); - (*_holder_fe_face_neighbor_helper[dim])->get_normals(); - - _holder_fe_neighbor_helper[dim] = &_fe_neighbor[dim][FEType(helper_order, LAGRANGE)]; - (*_holder_fe_neighbor_helper[dim])->get_xyz(); - (*_holder_fe_neighbor_helper[dim])->get_JxW(); + _holder_fe_helper[dim] = _fe[dim][FEType(helper_order, LAGRANGE)]; + _holder_fe_face_helper[dim] = _fe_face[dim][FEType(helper_order, LAGRANGE)]; + _holder_fe_face_neighbor_helper[dim] = _fe_face_neighbor[dim][FEType(helper_order, LAGRANGE)]; + _holder_fe_neighbor_helper[dim] = _fe_neighbor[dim][FEType(helper_order, LAGRANGE)]; } - for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++) - { - _holder_fe_lower_helper[dim] = &_fe_lower[dim][FEType(helper_order, LAGRANGE)]; - // We need these computations in order to compute correct lower-d element volumes in curvilinear - // coordinates - (*_holder_fe_lower_helper[dim])->get_xyz(); - (*_holder_fe_lower_helper[dim])->get_JxW(); - } + for (unsigned int dim = 0; dim < _mesh_dimension; dim++) + _holder_fe_lower_helper[dim] = _fe_lower[dim][FEType(helper_order, LAGRANGE)]; + + // request phi, dphi, xyz, JxW, etc. data + helpersRequestData(); // For 3D mortar, mortar segments are always TRI3 elements so we want FIRST LAGRANGE regardless // of discretization _fe_msm = (_mesh_dimension == 2) ? FEGenericBase::build(_mesh_dimension - 1, FEType(helper_order, LAGRANGE)) : FEGenericBase::build(_mesh_dimension - 1, FEType(FIRST, LAGRANGE)); + // This FE object should not take part in p-refinement + _fe_msm->add_p_level_in_reinit(false); _JxW_msm = &_fe_msm->get_JxW(); // Prerequest xyz so that it is computed for _fe_msm so that it can be used for calculating // _coord_msm @@ -267,6 +258,9 @@ Assembly::JxWNeighbor() const void Assembly::buildFE(FEType type) const { + if (!_building_helpers && type == _helper_type) + _user_added_fe_of_helper_type = true; + if (!_fe_shape_data[type]) _fe_shape_data[type] = std::make_unique(); @@ -290,6 +284,9 @@ Assembly::buildFE(FEType type) const void Assembly::buildFaceFE(FEType type) const { + if (!_building_helpers && type == _helper_type) + _user_added_fe_face_of_helper_type = true; + if (!_fe_shape_data_face[type]) _fe_shape_data_face[type] = std::make_unique(); @@ -309,6 +306,9 @@ Assembly::buildFaceFE(FEType type) const void Assembly::buildNeighborFE(FEType type) const { + if (!_building_helpers && type == _helper_type) + _user_added_fe_neighbor_of_helper_type = true; + if (!_fe_shape_data_neighbor[type]) _fe_shape_data_neighbor[type] = std::make_unique(); @@ -328,6 +328,9 @@ Assembly::buildNeighborFE(FEType type) const void Assembly::buildFaceNeighborFE(FEType type) const { + if (!_building_helpers && type == _helper_type) + _user_added_fe_face_neighbor_of_helper_type = true; + if (!_fe_shape_data_face_neighbor[type]) _fe_shape_data_face_neighbor[type] = std::make_unique(); @@ -347,6 +350,9 @@ Assembly::buildFaceNeighborFE(FEType type) const void Assembly::buildLowerDFE(FEType type) const { + if (!_building_helpers && type == _helper_type) + _user_added_fe_lower_of_helper_type = true; + if (!_fe_shape_data_lower[type]) _fe_shape_data_lower[type] = std::make_unique(); @@ -632,6 +638,11 @@ Assembly::setVolumeQRule(QBase * qrule, unsigned int dim) it.second->attach_quadrature_rule(qrule); for (auto & it : _vector_fe[dim]) it.second->attach_quadrature_rule(qrule); + if (!_unique_fe_helper.empty()) + { + mooseAssert(dim < _unique_fe_helper.size(), "We should not be indexing out of bounds"); + _unique_fe_helper[dim]->attach_quadrature_rule(qrule); + } } } @@ -644,6 +655,11 @@ Assembly::setFaceQRule(QBase * qrule, unsigned int dim) it.second->attach_quadrature_rule(qrule); for (auto & it : _vector_fe_face[dim]) it.second->attach_quadrature_rule(qrule); + if (!_unique_fe_face_helper.empty()) + { + mooseAssert(dim < _unique_fe_face_helper.size(), "We should not be indexing out of bounds"); + _unique_fe_face_helper[dim]->attach_quadrature_rule(qrule); + } } void @@ -658,6 +674,11 @@ Assembly::setLowerQRule(QBase * qrule, unsigned int dim) it.second->attach_quadrature_rule(qrule); for (auto & it : _vector_fe_lower[dim]) it.second->attach_quadrature_rule(qrule); + if (!_unique_fe_lower_helper.empty()) + { + mooseAssert(dim < _unique_fe_lower_helper.size(), "We should not be indexing out of bounds"); + _unique_fe_lower_helper[dim]->attach_quadrature_rule(qrule); + } } void @@ -669,6 +690,12 @@ Assembly::setNeighborQRule(QBase * qrule, unsigned int dim) it.second->attach_quadrature_rule(qrule); for (auto & it : _vector_fe_face_neighbor[dim]) it.second->attach_quadrature_rule(qrule); + if (!_unique_fe_face_neighbor_helper.empty()) + { + mooseAssert(dim < _unique_fe_face_neighbor_helper.size(), + "We should not be indexing out of bounds"); + _unique_fe_face_neighbor_helper[dim]->attach_quadrature_rule(qrule); + } } void @@ -740,12 +767,17 @@ Assembly::reinitFE(const Elem * elem) fesd._curl_phi.shallowCopy( const_cast>> &>(fe.get_curl_phi())); } + if (!_unique_fe_helper.empty()) + { + mooseAssert(dim < _unique_fe_helper.size(), "We should be in bounds here"); + _unique_fe_helper[dim]->reinit(elem); + } // During that last loop the helper objects will have been reinitialized as well // We need to dig out the q_points and JxW from it. _current_q_points.shallowCopy( - const_cast &>((*_holder_fe_helper[dim])->get_xyz())); - _current_JxW.shallowCopy(const_cast &>((*_holder_fe_helper[dim])->get_JxW())); + const_cast &>(_holder_fe_helper[dim]->get_xyz())); + _current_JxW.shallowCopy(const_cast &>(_holder_fe_helper[dim]->get_JxW())); if (_subproblem.haveADObjects()) { @@ -755,7 +787,7 @@ Assembly::reinitFE(const Elem * elem) { const auto & qw = _current_qrule->get_weights(); for (unsigned int qp = 0; qp != n_qp; qp++) - computeSinglePointMapAD(elem, qw, qp, *_holder_fe_helper[dim]); + computeSinglePointMapAD(elem, qw, qp, _holder_fe_helper[dim]); } else for (unsigned qp = 0; qp < n_qp; ++qp) @@ -1235,15 +1267,20 @@ Assembly::reinitFEFace(const Elem * elem, unsigned int side) fesd._curl_phi.shallowCopy( const_cast>> &>(fe_face.get_curl_phi())); } + if (!_unique_fe_face_helper.empty()) + { + mooseAssert(dim < _unique_fe_face_helper.size(), "We should be in bounds here"); + _unique_fe_face_helper[dim]->reinit(elem, side); + } // During that last loop the helper objects will have been reinitialized as well // We need to dig out the q_points and JxW from it. _current_q_points_face.shallowCopy( - const_cast &>((*_holder_fe_face_helper[dim])->get_xyz())); + const_cast &>(_holder_fe_face_helper[dim]->get_xyz())); _current_JxW_face.shallowCopy( - const_cast &>((*_holder_fe_face_helper[dim])->get_JxW())); + const_cast &>(_holder_fe_face_helper[dim]->get_JxW())); _current_normals.shallowCopy( - const_cast &>((*_holder_fe_face_helper[dim])->get_normals())); + const_cast &>(_holder_fe_face_helper[dim]->get_normals())); _mapped_normals.resize(_current_normals.size(), Eigen::Map(nullptr)); for (unsigned int i = 0; i < _current_normals.size(); i++) @@ -1252,7 +1289,7 @@ Assembly::reinitFEFace(const Elem * elem, unsigned int side) if (_calculate_curvatures) _curvatures.shallowCopy( - const_cast &>((*_holder_fe_face_helper[dim])->get_curvatures())); + const_cast &>(_holder_fe_face_helper[dim]->get_curvatures())); computeADFace(*elem, side); @@ -1277,9 +1314,9 @@ Assembly::computeFaceMap(const Elem & elem, const unsigned int side, const std:: const Elem & side_elem = _compute_face_map_side_elem_builder(elem, side); const auto dim = elem.dim(); const auto n_qp = qw.size(); - const auto & dpsidxi_map = (*_holder_fe_face_helper[dim])->get_fe_map().get_dpsidxi(); - const auto & dpsideta_map = (*_holder_fe_face_helper[dim])->get_fe_map().get_dpsideta(); - const auto & psi_map = (*_holder_fe_face_helper[dim])->get_fe_map().get_psi(); + const auto & dpsidxi_map = _holder_fe_face_helper[dim]->get_fe_map().get_dpsidxi(); + const auto & dpsideta_map = _holder_fe_face_helper[dim]->get_fe_map().get_dpsideta(); + const auto & psi_map = _holder_fe_face_helper[dim]->get_fe_map().get_psi(); std::vector> const * d2psidxi2_map = nullptr; std::vector> const * d2psidxideta_map = nullptr; std::vector> const * d2psideta2_map = nullptr; @@ -1288,9 +1325,9 @@ Assembly::computeFaceMap(const Elem & elem, const unsigned int side, const std:: if (_calculate_curvatures) { - d2psidxi2_map = &(*_holder_fe_face_helper[dim])->get_fe_map().get_d2psidxi2(); - d2psidxideta_map = &(*_holder_fe_face_helper[dim])->get_fe_map().get_d2psidxideta(); - d2psideta2_map = &(*_holder_fe_face_helper[dim])->get_fe_map().get_d2psideta2(); + d2psidxi2_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psidxi2(); + d2psidxideta_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psidxideta(); + d2psideta2_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psideta2(); } switch (dim) @@ -1530,6 +1567,12 @@ Assembly::reinitFEFaceNeighbor(const Elem * neighbor, const std::vector & fesd._curl_phi.shallowCopy(const_cast>> &>( fe_face_neighbor.get_curl_phi())); } + if (!_unique_fe_face_neighbor_helper.empty()) + { + mooseAssert(neighbor_dim < _unique_fe_face_neighbor_helper.size(), + "We should be in bounds here"); + _unique_fe_face_neighbor_helper[neighbor_dim]->reinit(neighbor, &reference_points); + } } void @@ -1577,6 +1620,11 @@ Assembly::reinitFENeighbor(const Elem * neighbor, const std::vector & ref fesd._curl_phi.shallowCopy( const_cast>> &>(fe_neighbor.get_curl_phi())); } + if (!_unique_fe_neighbor_helper.empty()) + { + mooseAssert(neighbor_dim < _unique_fe_neighbor_helper.size(), "We should be in bounds here"); + _unique_fe_neighbor_helper[neighbor_dim]->reinit(neighbor, &reference_points); + } } void @@ -1596,7 +1644,7 @@ Assembly::reinitNeighbor(const Elem * neighbor, const std::vector & refer if (_need_neighbor_elem_volume) { unsigned int dim = neighbor->dim(); - FEBase & fe = **_holder_fe_neighbor_helper[dim]; + FEBase & fe = *_holder_fe_neighbor_helper[dim]; QBase * qrule = qrules(dim).vol.get(); fe.attach_quadrature_rule(qrule); @@ -1691,7 +1739,7 @@ Assembly::reinitAtPhysical(const Elem * elem, const std::vector & physica _current_elem_volume_computed = false; FEInterface::inverse_map(elem->dim(), - (*_holder_fe_helper[elem->dim()])->get_fe_type(), + _holder_fe_helper[elem->dim()]->get_fe_type(), elem, physical_points, _temp_reference_points); @@ -1899,7 +1947,7 @@ Assembly::reinitElemAndNeighbor(const Elem * elem, // compute JxW on the neighbor's face _current_JxW_neighbor.shallowCopy(const_cast &>( - (*_holder_fe_face_neighbor_helper[_current_neighbor_side_elem->dim()])->get_JxW())); + _holder_fe_face_neighbor_helper[_current_neighbor_side_elem->dim()]->get_JxW())); } reinitFEFaceNeighbor(neighbor, *reference_points_ptr); @@ -1971,20 +2019,26 @@ Assembly::reinitElemFaceRef(const Elem * elem, fesd._curl_phi.shallowCopy( const_cast>> &>(fe_face.get_curl_phi())); } + if (!_unique_fe_face_helper.empty()) + { + mooseAssert(elem_dim < _unique_fe_face_helper.size(), "We should be in bounds here"); + _unique_fe_face_helper[elem_dim]->reinit(elem, elem_side, tolerance, pts, weights); + } + // During that last loop the helper objects will have been reinitialized _current_q_points_face.shallowCopy( - const_cast &>((*_holder_fe_face_helper[elem_dim])->get_xyz())); + const_cast &>(_holder_fe_face_helper[elem_dim]->get_xyz())); _current_normals.shallowCopy( - const_cast &>((*_holder_fe_face_helper[elem_dim])->get_normals())); + const_cast &>(_holder_fe_face_helper[elem_dim]->get_normals())); _current_tangents.shallowCopy(const_cast> &>( - (*_holder_fe_face_helper[elem_dim])->get_tangents())); + _holder_fe_face_helper[elem_dim]->get_tangents())); // Note that if the user did pass in points and not weights to this method, JxW will be garbage // and should not be used _current_JxW_face.shallowCopy( - const_cast &>((*_holder_fe_face_helper[elem_dim])->get_JxW())); + const_cast &>(_holder_fe_face_helper[elem_dim]->get_JxW())); if (_calculate_curvatures) _curvatures.shallowCopy( - const_cast &>((*_holder_fe_face_helper[elem_dim])->get_curvatures())); + const_cast &>(_holder_fe_face_helper[elem_dim]->get_curvatures())); computeADFace(*elem, elem_side); } @@ -2012,7 +2066,7 @@ Assembly::computeADFace(const Elem & elem, const unsigned int side) const std::vector dummy_qw(n_qp, 1.); for (unsigned int qp = 0; qp != n_qp; qp++) - computeSinglePointMapAD(&elem, dummy_qw, qp, *_holder_fe_face_helper[dim]); + computeSinglePointMapAD(&elem, dummy_qw, qp, _holder_fe_face_helper[dim]); } else for (unsigned qp = 0; qp < n_qp; ++qp) @@ -2129,12 +2183,19 @@ Assembly::reinitNeighborFaceRef(const Elem * neighbor, fesd._curl_phi.shallowCopy(const_cast>> &>( fe_face_neighbor.get_curl_phi())); } + if (!_unique_fe_face_neighbor_helper.empty()) + { + mooseAssert(neighbor_dim < _unique_fe_face_neighbor_helper.size(), + "We should be in bounds here"); + _unique_fe_face_neighbor_helper[neighbor_dim]->reinit( + neighbor, neighbor_side, tolerance, pts, weights); + } // During that last loop the helper objects will have been reinitialized as well // We need to dig out the q_points from it - _current_q_points_face_neighbor.shallowCopy(const_cast &>( - (*_holder_fe_face_neighbor_helper[neighbor_dim])->get_xyz())); + _current_q_points_face_neighbor.shallowCopy( + const_cast &>(_holder_fe_face_neighbor_helper[neighbor_dim]->get_xyz())); _current_neighbor_normals.shallowCopy(const_cast &>( - (*_holder_fe_face_neighbor_helper[neighbor_dim])->get_normals())); + _holder_fe_face_neighbor_helper[neighbor_dim]->get_normals())); } void @@ -2210,6 +2271,11 @@ Assembly::reinitLowerDElem(const Elem * elem, const_cast>> &>(fe_lower.get_dual_d2phi())); } } + if (!_unique_fe_lower_helper.empty()) + { + mooseAssert(elem_dim < _unique_fe_lower_helper.size(), "We should be in bounds here"); + _unique_fe_lower_helper[elem_dim]->reinit(elem); + } if (!_need_lower_d_elem_volume) return; @@ -2231,7 +2297,7 @@ Assembly::reinitLowerDElem(const Elem * elem, else { // During that last loop the helper objects will have been reinitialized as well - FEBase & helper_fe = **_holder_fe_lower_helper[elem_dim]; + FEBase & helper_fe = *_holder_fe_lower_helper[elem_dim]; const auto & physical_q_points = helper_fe.get_xyz(); const auto & JxW = helper_fe.get_JxW(); MooseArray coord; @@ -2296,7 +2362,7 @@ Assembly::reinitNeighborAtPhysical(const Elem * neighbor, // compute JxW on the neighbor's face unsigned int neighbor_side_dim = _current_neighbor_side_elem->dim(); _current_JxW_neighbor.shallowCopy(const_cast &>( - (*_holder_fe_face_neighbor_helper[neighbor_side_dim])->get_JxW())); + _holder_fe_face_neighbor_helper[neighbor_side_dim]->get_JxW())); reinitFEFaceNeighbor(neighbor, reference_points); reinitNeighbor(neighbor, reference_points); @@ -4605,6 +4671,96 @@ Assembly::adCurvatures() const return _ad_curvatures; } +void +Assembly::helpersRequestData() +{ + for (unsigned int dim = 0; dim <= _mesh_dimension; dim++) + { + _holder_fe_helper[dim]->get_phi(); + _holder_fe_helper[dim]->get_dphi(); + _holder_fe_helper[dim]->get_xyz(); + _holder_fe_helper[dim]->get_JxW(); + + _holder_fe_face_helper[dim]->get_phi(); + _holder_fe_face_helper[dim]->get_dphi(); + _holder_fe_face_helper[dim]->get_xyz(); + _holder_fe_face_helper[dim]->get_JxW(); + _holder_fe_face_helper[dim]->get_normals(); + + _holder_fe_face_neighbor_helper[dim]->get_xyz(); + _holder_fe_face_neighbor_helper[dim]->get_JxW(); + _holder_fe_face_neighbor_helper[dim]->get_normals(); + + _holder_fe_neighbor_helper[dim]->get_xyz(); + _holder_fe_neighbor_helper[dim]->get_JxW(); + } + + for (unsigned int dim = 0; dim < _mesh_dimension; dim++) + { + // We need these computations in order to compute correct lower-d element volumes in curvilinear + // coordinates + _holder_fe_lower_helper[dim]->get_xyz(); + _holder_fe_lower_helper[dim]->get_JxW(); + } +} + +void +Assembly::havePRefinement() +{ + const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST; + const FEType helper_type(helper_order, LAGRANGE); + auto setup_helpers = [&helper_type](auto & unique_helper_container, + auto & helper_container, + const unsigned int num_dimensionalities, + const bool user_added_helper_type, + auto & fe_container) + { + unique_helper_container.resize(num_dimensionalities); + for (const auto dim : make_range(num_dimensionalities)) + { + auto & unique_helper = unique_helper_container[dim]; + unique_helper = FEGenericBase::build(dim, helper_type); + // don't participate in p-refinement + unique_helper->add_p_level_in_reinit(false); + helper_container[dim] = unique_helper.get(); + + // If the user did not request the helper type then we should erase it from our FE container + // so that they're not penalized (in the "we should be able to do p-refinement sense") for our + // perhaps silly helpers + if (!user_added_helper_type) + fe_container[dim].erase(helper_type); + } + }; + + setup_helpers(_unique_fe_helper, + _holder_fe_helper, + _mesh_dimension + 1, + _user_added_fe_of_helper_type, + _fe); + setup_helpers(_unique_fe_face_helper, + _holder_fe_face_helper, + _mesh_dimension + 1, + _user_added_fe_face_of_helper_type, + _fe_face); + setup_helpers(_unique_fe_face_neighbor_helper, + _holder_fe_face_neighbor_helper, + _mesh_dimension + 1, + _user_added_fe_face_neighbor_of_helper_type, + _fe_face_neighbor); + setup_helpers(_unique_fe_neighbor_helper, + _holder_fe_neighbor_helper, + _mesh_dimension + 1, + _user_added_fe_neighbor_of_helper_type, + _fe_neighbor); + setup_helpers(_unique_fe_lower_helper, + _holder_fe_lower_helper, + _mesh_dimension, + _user_added_fe_lower_of_helper_type, + _fe_lower); + + helpersRequestData(); +} + template void coordTransformFactor(const SubProblem & s, SubdomainID sub_id, const Point & point, diff --git a/framework/src/problems/DisplacedProblem.C b/framework/src/problems/DisplacedProblem.C index 2c7a64ce91d1..bc9d3a4c9722 100644 --- a/framework/src/problems/DisplacedProblem.C +++ b/framework/src/problems/DisplacedProblem.C @@ -1303,3 +1303,11 @@ DisplacedProblem::currentResidualVectorTags() const { return _mproblem.currentResidualVectorTags(); } + +void +DisplacedProblem::havePRefinement() +{ + for (auto & assembly_vecs : _assembly) + for (auto & assembly : assembly_vecs) + assembly->havePRefinement(); +} diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C index f91c69f13a02..a561b91809d2 100644 --- a/framework/src/problems/FEProblemBase.C +++ b/framework/src/problems/FEProblemBase.C @@ -8095,3 +8095,14 @@ FEProblemBase::reinitMortarUserObjects(const BoundaryID primary_boundary_id, mortar_uo->reinit(); } } + +void +FEProblemBase::havePRefinement() +{ + for (auto & assembly_vecs : _assembly) + for (auto & assembly : assembly_vecs) + assembly->havePRefinement(); + + if (_displaced_problem) + _displaced_problem->havePRefinement(); +} From 53904f11787a321878c524a04a1c4d54b2054de9 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Mon, 26 Jun 2023 09:25:54 -0700 Subject: [PATCH 04/12] Add gold file for test --- .../2d_diffusion_dg/gold/p_refinement.e-s003 | Bin 0 -> 6800 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 test/tests/dgkernels/2d_diffusion_dg/gold/p_refinement.e-s003 diff --git a/test/tests/dgkernels/2d_diffusion_dg/gold/p_refinement.e-s003 b/test/tests/dgkernels/2d_diffusion_dg/gold/p_refinement.e-s003 new file mode 100644 index 0000000000000000000000000000000000000000..f9fb8f769d700b391c970dd4aacc111dc8759d67 GIT binary patch literal 6800 zcmeHLO=ule6rR{NnR$6>YSWjNBGQFOT~uOQ-3XZ!pG6TD1r-r-nB>jN8#;f^Ow!aY z+!Qx1gz836$-;$OMS=)|bz!pWDqXs8(eJx=&g;!<(mdCb19$Gd=YHoq z=bn4!zDZtapMAu#tP!*Z+95PM3KB2z<3Ot|rwrr3E6O0lVkKkAQtWlYm{Zz2oF{1~ zC^W!9XHyqM0YQ7w8hY-{N9PrYJ}1P^D}u7JC06N+umiuF^!qE0Fd@w|iih~gipR@> z96r?o^-bzEPeCPPHFYPq(4E>sce>Jzi=RZ%S$m*2bfg>9k#0~&y3M*}Kc|j4QGAP0 z>Qy@WXrhgA&Cfz_ImnAJO)$0gp^cmw`ROC{W83JDbfeTS!(`FR(l9B#3u)d#PTs@$ zs&k=F8_ns45%R4Uj5V%L#r{ed^Ndj!9tvc!>_OePKZbD=j7fyp1A5|= z=xyLlf5h(OPwX!P>l6?AV?BtN{yu@y2Hx~X?6Cg61a@>GMU|{5dgg>RD8krZ(i7)f z^cs6nXC4+!R+>1BPe4&;m^wuW~$jG{OEXh|0fG32CREW<_{bLQ0KlKcO_E@t zl>3hF3GZX^?eI?a%8QxtK8~Jy{Y+ckO1jqPfQPhtcr+D$vw!*&;kc^R)>jzY3SX!T zlCko4MTA(!18IGShi%-=55xXQY{knf} z-kq=;#8$jaD|si-L`>e>Pzvu(*t^8)dg$kOzmTFQZxw&*Ujr+7nA{pWt>US8WpCA) zDhq>mntF*r;(c@*-p8PC^2OrELG^>i;yzg6r6z}b@28K$ZerPw@V06`7CD5R%<(q& z)X^EN={WBHu+nAwSouNDeTQAI>tmN=9+0E z>%?2`hU)sgoR$i&o!|1mK6n3k>%{}>&%WmV^6BimSFS$gUca< Date: Mon, 26 Jun 2023 09:29:11 -0700 Subject: [PATCH 05/12] Remove serial restriction on some tests --- test/tests/dgkernels/2d_diffusion_dg/tests | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/tests/dgkernels/2d_diffusion_dg/tests b/test/tests/dgkernels/2d_diffusion_dg/tests index 60ea9696eec5..4db589f741d7 100644 --- a/test/tests/dgkernels/2d_diffusion_dg/tests +++ b/test/tests/dgkernels/2d_diffusion_dg/tests @@ -5,7 +5,6 @@ type = 'Exodiff' input = '2d_diffusion_dg_test.i' exodiff = 'out.e-s003' - max_parallel = 1 requirement = "The system shall support solving 2D diffusion using the discontinuous Galerkin " "method." [] @@ -14,16 +13,15 @@ input = '2d_diffusion_dg_test.i' cli_args = 'Outputs/file_base=p_refinement Executioner/Adaptivity/switch_h_to_p_refinement=true' exodiff = 'p_refinement.e-s003' - max_parallel = 1 requirement = "The system shall support solving 2D diffusion using the discontinuous Galerkin " "method with p-refinement." + issues = '#869 #24141' [] [constMonomial] type = 'Exodiff' input = '2d_diffusion_dg_test.i' exodiff = 'constMonomial.e-s003' cli_args = 'Outputs/file_base=constMonomial Variables/u/order=CONSTANT' - max_parallel = 1 issues = '#21314' requirement = "The system shall support solving 2D diffusion using the discontinuous Galerkin " "method with constant monomial shape functions." From a0307ed2d1fc9453a7b7a7c6193d2c94e24a4a87 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Mon, 26 Jun 2023 09:47:20 -0700 Subject: [PATCH 06/12] Add 3D test --- .../3d_diffusion_p_refinement.i | 107 ++++++++++++++++++ .../gold/3d_diffusion_p_refinement_out.e-s003 | Bin 0 -> 8544 bytes test/tests/dgkernels/3d_diffusion_dg/tests | 15 ++- 3 files changed, 118 insertions(+), 4 deletions(-) create mode 100644 test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_p_refinement.i create mode 100644 test/tests/dgkernels/3d_diffusion_dg/gold/3d_diffusion_p_refinement_out.e-s003 diff --git a/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_p_refinement.i b/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_p_refinement.i new file mode 100644 index 000000000000..f3809bf0171a --- /dev/null +++ b/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_p_refinement.i @@ -0,0 +1,107 @@ +[Mesh] + type = GeneratedMesh + dim = 3 + nx = 1 + ny = 1 + nz = 1 + xmin = 0 + xmax = 1 + ymin = 0 + ymax = 1 + zmin = 0 + zmax = 1 + elem_type = HEX8 +[] + +[Variables] + [u] + order = FIRST + family = MONOMIAL + [InitialCondition] + type = ConstantIC + value = 0.5 + [] + [] +[] + +[Functions] + [forcing_fn] + type = ParsedFunction + expression = 2*pow(e,-x-(y*y))*(1-2*y*y) + [] + [exact_fn] + type = ParsedGradFunction + value = pow(e,-x-(y*y)) + grad_x = -pow(e,-x-(y*y)) + grad_y = -2*y*pow(e,-x-(y*y)) + [] +[] + +[Kernels] + [diff] + type = Diffusion + variable = u + [] + [abs] + type = Reaction + variable = u + [] + [forcing] + type = BodyForce + variable = u + function = forcing_fn + [] +[] + +[DGKernels] + [dg_diff] + type = DGDiffusion + variable = u + epsilon = -1 + sigma = 6 + [] +[] + +[BCs] + [all] + type = DGFunctionDiffusionDirichletBC + variable = u + boundary = '0 1 2 3 4 5' + function = exact_fn + epsilon = -1 + sigma = 6 + [] +[] + +[Executioner] + type = Steady + solve_type = 'PJFNK' + [Adaptivity] + switch_h_to_p_refinement = true + steps = 2 + refine_fraction = 1.0 + coarsen_fraction = 0 + max_h_level = 8 + [] +[] + +[Postprocessors] + [h] + type = AverageElementSize + execute_on = 'initial timestep_end' + [] + [dofs] + type = NumDOFs + execute_on = 'initial timestep_end' + [] + [l2_err] + type = ElementL2Error + variable = u + function = exact_fn + execute_on = 'initial timestep_end' + [] +[] + +[Outputs] + exodus = true +[] diff --git a/test/tests/dgkernels/3d_diffusion_dg/gold/3d_diffusion_p_refinement_out.e-s003 b/test/tests/dgkernels/3d_diffusion_dg/gold/3d_diffusion_p_refinement_out.e-s003 new file mode 100644 index 0000000000000000000000000000000000000000..cfd26d76da41cc032870bfdc498492ff91f0354b GIT binary patch literal 8544 zcmeHL&x;&I6z(-Ywr6KHn>Z$}u*wD^UJ{9$HOV13E5@*aAP5rnR`l-F>@@atH{CUx zWXU0i96SigA%`3kauz}HAR>r~e}TCeym=J#kh>S(S5xP0X$&+}%e%u!jS;zu$TaWIfZ?KM<5kRn$yrA{le@o*rz;Q-A8R2pcJ7~Q*Lbmy(E$^OLJtap&~jt*vvI+!i$V791(*`f|+_v?0o z3_8qL+ue>5VRcxq7L{4_1!*XDWR{0XOik|)m6`Wvez`~WY>9rc7bSrT<86^9VXVZB zBt{NgKJ>Q-!Od`hoS{4ulyKXsPwLV5WqNKA`T12!yg&_+ z8i*mrLJp}=MKDN${C)FWm(sD<8)70N6`AaXu^dQpFBz)0<(d5KSsh@)HmSTtWlrby zjecT@z|H*h=phd=50aDkX~@kzXQ&M<;>800P)M4neZr?%l5*Ywi@0OlLGylT-w#~9 zOd2zm$LsjT`3v^o6AQ3-UrCza^$G9V1GdT@k3U7&M*gBboz&hZ zggv~KQ0Cm!^oR-ak_Eo6nMO0B+F);K&3|FP9$|MYv40Z+YlNjhej$%|J$W8;;z}4r z*Mn$il8yIDkkd0oeTTOS`@rGl zGj*Ya_o2gEh5Z6plfR>UI9lMf?WIQ!me=iA3HDnD%WLs+3HGsr<@mNsu)jK3j@c_E z*gqUB$KiMh_BrWW*kJI{WI4j0s zU9KXYhIn2F(RXUo+WW2QXJi$ws-JndF5-)A#?`gn47SWamzoA(Gf^thfBEVefGphZ zt9njd&srsEjc4z9k?hrZ`MqA2#B~n9bpN03|8>s)>Ha_6|LdIpwf^{-PJS8zbKg-`@( pzx6*_dHiId(Err-e`os}Ypa&$pYM17+_?5tC;ai}yQjbT@?RNdx`_Y) literal 0 HcmV?d00001 diff --git a/test/tests/dgkernels/3d_diffusion_dg/tests b/test/tests/dgkernels/3d_diffusion_dg/tests index 7a1206169c62..1c7a1b9c80b9 100644 --- a/test/tests/dgkernels/3d_diffusion_dg/tests +++ b/test/tests/dgkernels/3d_diffusion_dg/tests @@ -1,13 +1,20 @@ [Tests] - issues = '#869' design = 'source/bcs/DGFunctionDiffusionDirichletBC.md source/dgkernels/DGDiffusion.md' - [./test] + [test] type = 'Exodiff' input = '3d_diffusion_dg_test.i' exodiff = 'out.e' group = 'adaptive' - max_parallel = 1 valgrind = 'HEAVY' requirement = "The system shall support solving 3D diffusion using the discontinuous Galerkin method." - [../] + issues = '#869' + [] + [p_refinement_test] + type = 'Exodiff' + input = '3d_diffusion_p_refinement.i' + exodiff = '3d_diffusion_p_refinement_out.e-s003' + requirement = "The system shall support solving 3D diffusion using the discontinuous Galerkin " + "method with p-refinement." + issues = '#869 #24141' + [] [] From cc9821c10feec4511ebc5eb66cc1424e78f2aa97 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Mon, 26 Jun 2023 10:03:38 -0700 Subject: [PATCH 07/12] Fix valgrind leak of helper objects --- framework/src/base/Assembly.C | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/framework/src/base/Assembly.C b/framework/src/base/Assembly.C index d4f6476cb4be..b62e8f57b142 100644 --- a/framework/src/base/Assembly.C +++ b/framework/src/base/Assembly.C @@ -4728,7 +4728,13 @@ Assembly::havePRefinement() // so that they're not penalized (in the "we should be able to do p-refinement sense") for our // perhaps silly helpers if (!user_added_helper_type) - fe_container[dim].erase(helper_type); + { + auto & fe_container_dim = fe_container[dim]; + auto fe_it = fe_container_dim.find(helper_type); + mooseAssert(fe_it != fe_container_dim.end(), "We should have the helper type"); + delete fe_it->second; + fe_container_dim.erase(fe_it); + } } }; From 8e3152507711f1475193e046dfeb7a6d7a8518f0 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Mon, 26 Jun 2023 13:00:16 -0700 Subject: [PATCH 08/12] Input file fixes - formatting - 'value' to 'expression' for ParsedGradFunction --- .../2d_diffusion_dg/2d_diffusion_dg_test.i | 54 +++++++++---------- .../3d_diffusion_dg/3d_diffusion_dg_test.i | 50 ++++++++--------- .../3d_diffusion_p_refinement.i | 2 +- 3 files changed, 53 insertions(+), 53 deletions(-) diff --git a/test/tests/dgkernels/2d_diffusion_dg/2d_diffusion_dg_test.i b/test/tests/dgkernels/2d_diffusion_dg/2d_diffusion_dg_test.i index ba5eda6628e8..7c2238c39e42 100644 --- a/test/tests/dgkernels/2d_diffusion_dg/2d_diffusion_dg_test.i +++ b/test/tests/dgkernels/2d_diffusion_dg/2d_diffusion_dg_test.i @@ -28,30 +28,30 @@ [Variables] active = 'u' - [./u] + [u] order = FIRST family = MONOMIAL - [./InitialCondition] + [InitialCondition] type = ConstantIC value = 1 - [../] - [../] + [] + [] [] [Functions] active = 'forcing_fn exact_fn' - [./forcing_fn] + [forcing_fn] type = ParsedFunction # function = -4.0+(x*x)+(y*y) # function = x # function = (x*x)-2.0 expression = 2*pow(e,-x-(y*y))*(1-2*y*y) # function = (x*x*x)-6.0*x - [../] + [] - [./exact_fn] + [exact_fn] type = ParsedGradFunction # function = x # grad_x = 1 @@ -65,58 +65,58 @@ # grad_x = 2*x # grad_y = 0 - value = pow(e,-x-(y*y)) + expression = pow(e,-x-(y*y)) grad_x = -pow(e,-x-(y*y)) grad_y = -2*y*pow(e,-x-(y*y)) # function = (x*x*x) # grad_x = 3*x*x # grad_y = 0 - [../] + [] [] [Kernels] active = 'diff abs forcing' - [./diff] + [diff] type = Diffusion variable = u - [../] + [] - [./abs] # u * v + [abs] # u * v type = Reaction variable = u - [../] + [] - [./forcing] + [forcing] type = BodyForce variable = u function = forcing_fn - [../] + [] [] [DGKernels] active = 'dg_diff' - [./dg_diff] + [dg_diff] type = DGDiffusion variable = u epsilon = -1 sigma = 6 - [../] + [] [] [BCs] active = 'all' - [./all] + [all] type = DGFunctionDiffusionDirichletBC variable = u boundary = '0 1 2 3' function = exact_fn epsilon = -1 sigma = 6 - [../] + [] [] [Executioner] @@ -129,12 +129,12 @@ # petsc_options = '-snes_mf' # max_r_steps = 2 - [./Adaptivity] + [Adaptivity] steps = 2 refine_fraction = 1.0 coarsen_fraction = 0 max_h_level = 8 - [../] + [] nl_rel_tol = 1e-10 @@ -144,19 +144,19 @@ [Postprocessors] active = 'h dofs l2_err' - [./h] + [h] type = AverageElementSize - [../] + [] - [./dofs] + [dofs] type = NumDOFs - [../] + [] - [./l2_err] + [l2_err] type = ElementL2Error variable = u function = exact_fn - [../] + [] [] [Outputs] diff --git a/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_dg_test.i b/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_dg_test.i index a0bdff3c19d8..bc34bf5f90e3 100644 --- a/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_dg_test.i +++ b/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_dg_test.i @@ -16,76 +16,76 @@ [Variables] active = 'u' - [./u] + [u] order = FIRST family = MONOMIAL - [./InitialCondition] + [InitialCondition] type = ConstantIC value = 0.5 - [../] - [../] + [] + [] [] [Functions] active = 'forcing_fn exact_fn' - [./forcing_fn] + [forcing_fn] type = ParsedFunction expression = 2*pow(e,-x-(y*y))*(1-2*y*y) - [../] + [] - [./exact_fn] + [exact_fn] type = ParsedGradFunction - value = pow(e,-x-(y*y)) + expression = pow(e,-x-(y*y)) grad_x = -pow(e,-x-(y*y)) grad_y = -2*y*pow(e,-x-(y*y)) - [../] + [] [] [Kernels] active = 'diff abs forcing' - [./diff] + [diff] type = Diffusion variable = u - [../] + [] - [./abs] # u * v + [abs] # u * v type = Reaction variable = u - [../] + [] - [./forcing] + [forcing] type = BodyForce variable = u function = forcing_fn - [../] + [] [] [DGKernels] active = 'dg_diff' - [./dg_diff] + [dg_diff] type = DGDiffusion variable = u epsilon = -1 sigma = 6 - [../] + [] [] [BCs] active = 'all' - [./all] + [all] type = DGFunctionDiffusionDirichletBC variable = u boundary = '0 1 2 3 4 5' function = exact_fn epsilon = -1 sigma = 6 - [../] + [] [] [Executioner] @@ -97,22 +97,22 @@ [Postprocessors] active = 'h dofs l2_err' - [./h] + [h] type = AverageElementSize execute_on = 'initial timestep_end' - [../] + [] - [./dofs] + [dofs] type = NumDOFs execute_on = 'initial timestep_end' - [../] + [] - [./l2_err] + [l2_err] type = ElementL2Error variable = u function = exact_fn execute_on = 'initial timestep_end' - [../] + [] [] [Outputs] diff --git a/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_p_refinement.i b/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_p_refinement.i index f3809bf0171a..83e3bd21ea5f 100644 --- a/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_p_refinement.i +++ b/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_p_refinement.i @@ -31,7 +31,7 @@ [] [exact_fn] type = ParsedGradFunction - value = pow(e,-x-(y*y)) + expression = pow(e,-x-(y*y)) grad_x = -pow(e,-x-(y*y)) grad_y = -2*y*pow(e,-x-(y*y)) [] From 9fa77820801f77ae96912b8bdea459b93b3cca68 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Mon, 26 Jun 2023 13:42:36 -0700 Subject: [PATCH 09/12] Don't add_p_level_in_reinit when computing max qps This is consistent with previous behavior of only using first order Lagrange --- framework/src/loops/MaxQpsThread.C | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/framework/src/loops/MaxQpsThread.C b/framework/src/loops/MaxQpsThread.C index 841523baacd9..47dda65485fa 100644 --- a/framework/src/loops/MaxQpsThread.C +++ b/framework/src/loops/MaxQpsThread.C @@ -64,11 +64,17 @@ MaxQpsThread::operator()(const ConstElemRange & range) // We'll use one for element interiors, which calculates nothing std::unique_ptr fe(FEBase::build(dim, fe_type)); fe->get_nothing(); + // Alex: This seems wrong, e.g. like we'd want the quadrature order to be higher if the + // polynomial basis is higher, but our fe_type up above is first Lagrange, so it seems like the + // person who wrote this code didn't want to consider higher order bases even if if the user + // asked for higher order bases + fe->add_p_level_in_reinit(false); // And another for element sides, which calculates the minimum // libMesh currently allows for that std::unique_ptr side_fe(FEBase::build(dim, fe_type)); side_fe->get_xyz(); + side_fe->add_p_level_in_reinit(false); // figure out the number of qps for volume auto qrule = assem.attachQRuleElem(dim, *fe); From 3c779f04d36121d93305ee759c237d4c5e622af9 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Mon, 26 Jun 2023 13:43:37 -0700 Subject: [PATCH 10/12] Don't add p level when querying the number of shape functions Once libmesh/libmesh#3597 is in, then we can use that. But until then we can do this clumsy thing --- framework/src/loops/MaxQpsThread.C | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/framework/src/loops/MaxQpsThread.C b/framework/src/loops/MaxQpsThread.C index 47dda65485fa..2d871ad4a59b 100644 --- a/framework/src/loops/MaxQpsThread.C +++ b/framework/src/loops/MaxQpsThread.C @@ -11,7 +11,7 @@ #include "FEProblem.h" #include "Assembly.h" -#include "libmesh/fe_base.h" +#include "libmesh/fe.h" #include "libmesh/threads.h" #include LIBMESH_INCLUDE_UNORDERED_SET LIBMESH_DEFINE_HASH_POINTERS @@ -82,7 +82,23 @@ MaxQpsThread::operator()(const ConstElemRange & range) if (qrule->n_points() > _max) _max = qrule->n_points(); - unsigned int n_shape_funcs = fe->n_shape_functions(); + const auto n_shape_funcs = [dim, &fe]() + { + const auto elem_type = fe->get_type(); + switch (dim) + { + case 0: + return static_cast *>(fe.get())->n_shape_functions(elem_type, FIRST); + case 1: + return static_cast *>(fe.get())->n_shape_functions(elem_type, FIRST); + case 2: + return static_cast *>(fe.get())->n_shape_functions(elem_type, FIRST); + case 3: + return static_cast *>(fe.get())->n_shape_functions(elem_type, FIRST); + default: + mooseError("Unhandled dimension ", dim); + } + }(); if (n_shape_funcs > _max_shape_funcs) _max_shape_funcs = n_shape_funcs; From b327bb2cfbb8f06013450e37c8e679c326fa6db1 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Wed, 5 Jul 2023 10:36:45 -0700 Subject: [PATCH 11/12] Remove _max_shape_funcs --- framework/include/loops/MaxQpsThread.h | 5 ---- framework/include/problems/FEProblemBase.h | 8 ------ framework/src/loops/MaxQpsThread.C | 30 ++-------------------- framework/src/problems/FEProblemBase.C | 11 -------- 4 files changed, 2 insertions(+), 52 deletions(-) diff --git a/framework/include/loops/MaxQpsThread.h b/framework/include/loops/MaxQpsThread.h index 11a0c1f146b6..9b3c003c5633 100644 --- a/framework/include/loops/MaxQpsThread.h +++ b/framework/include/loops/MaxQpsThread.h @@ -38,8 +38,6 @@ class MaxQpsThread unsigned int max() const { return _max; } - unsigned int max_shape_funcs() const { return _max_shape_funcs; } - protected: FEProblemBase & _fe_problem; @@ -47,7 +45,4 @@ class MaxQpsThread /// Maximum number of qps encountered unsigned int _max; - - /// Maximum number of shape functions encountered - unsigned int _max_shape_funcs; }; diff --git a/framework/include/problems/FEProblemBase.h b/framework/include/problems/FEProblemBase.h index b89be1520a59..4e71b3065f47 100644 --- a/framework/include/problems/FEProblemBase.h +++ b/framework/include/problems/FEProblemBase.h @@ -339,11 +339,6 @@ class FEProblemBase : public SubProblem, public Restartable */ unsigned int getMaxQps() const; - /** - * @return The maximum number of quadrature points in use on any element in this problem. - */ - unsigned int getMaxShapeFunctions() const; - /** * @return The maximum order for all scalar variables in this problem's systems. */ @@ -2389,9 +2384,6 @@ class FEProblemBase : public SubProblem, public Restartable /// Maximum number of quadrature points used in the problem unsigned int _max_qps; - /// Maximum number of shape functions on any element in the problem - unsigned int _max_shape_funcs; - /// Maximum scalar variable order Order _max_scalar_order; diff --git a/framework/src/loops/MaxQpsThread.C b/framework/src/loops/MaxQpsThread.C index 2d871ad4a59b..5e542dbf3da5 100644 --- a/framework/src/loops/MaxQpsThread.C +++ b/framework/src/loops/MaxQpsThread.C @@ -17,14 +17,11 @@ LIBMESH_DEFINE_HASH_POINTERS #include "libmesh/quadrature.h" -MaxQpsThread::MaxQpsThread(FEProblemBase & fe_problem) - : _fe_problem(fe_problem), _max(0), _max_shape_funcs(0) -{ -} +MaxQpsThread::MaxQpsThread(FEProblemBase & fe_problem) : _fe_problem(fe_problem), _max(0) {} // Splitting Constructor MaxQpsThread::MaxQpsThread(MaxQpsThread & x, Threads::split /*split*/) - : _fe_problem(x._fe_problem), _max(x._max), _max_shape_funcs(x._max_shape_funcs) + : _fe_problem(x._fe_problem), _max(x._max) { } @@ -82,26 +79,6 @@ MaxQpsThread::operator()(const ConstElemRange & range) if (qrule->n_points() > _max) _max = qrule->n_points(); - const auto n_shape_funcs = [dim, &fe]() - { - const auto elem_type = fe->get_type(); - switch (dim) - { - case 0: - return static_cast *>(fe.get())->n_shape_functions(elem_type, FIRST); - case 1: - return static_cast *>(fe.get())->n_shape_functions(elem_type, FIRST); - case 2: - return static_cast *>(fe.get())->n_shape_functions(elem_type, FIRST); - case 3: - return static_cast *>(fe.get())->n_shape_functions(elem_type, FIRST); - default: - mooseError("Unhandled dimension ", dim); - } - }(); - if (n_shape_funcs > _max_shape_funcs) - _max_shape_funcs = n_shape_funcs; - // figure out the number of qps for the face // NOTE: user might specify higher order rule for faces, thus possibly ending up with more qps // than in the volume @@ -126,7 +103,4 @@ MaxQpsThread::join(const MaxQpsThread & y) { if (y._max > _max) _max = y._max; - - if (y._max_shape_funcs > _max_shape_funcs) - _max_shape_funcs = y._max_shape_funcs; } diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C index a561b91809d2..110b10344489 100644 --- a/framework/src/problems/FEProblemBase.C +++ b/framework/src/problems/FEProblemBase.C @@ -362,7 +362,6 @@ FEProblemBase::FEProblemBase(const InputParameters & parameters) _material_dependency_check(getParam("material_dependency_check")), _uo_aux_state_check(getParam("check_uo_aux_state")), _max_qps(std::numeric_limits::max()), - _max_shape_funcs(std::numeric_limits::max()), _max_scalar_order(INVALID_ORDER), _has_time_integrator(false), _has_exception(false), @@ -1362,14 +1361,6 @@ FEProblemBase::getMaxQps() const return _max_qps; } -unsigned int -FEProblemBase::getMaxShapeFunctions() const -{ - if (_max_shape_funcs == std::numeric_limits::max()) - mooseError("Max shape functions uninitialized"); - return _max_shape_funcs; -} - Order FEProblemBase::getMaxScalarOrder() const { @@ -5241,13 +5232,11 @@ FEProblemBase::updateMaxQps() MaxQpsThread mqt(*this); Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), mqt); _max_qps = mqt.max(); - _max_shape_funcs = mqt.max_shape_funcs(); // If we have more shape functions or more quadrature points on // another processor, then we may need to handle those elements // ourselves later after repartitioning. _communicator.max(_max_qps); - _communicator.max(_max_shape_funcs); } unsigned int max_qpts = getMaxQps(); From 0383427ade8643a7c50702f353061ddbb05e8296 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Wed, 5 Jul 2023 16:47:12 -0700 Subject: [PATCH 12/12] Just init quadrature rules in MaxQpsThread --- framework/include/base/Assembly.h | 15 +++++++++ framework/src/base/Assembly.C | 48 ++++++++++++++++++---------- framework/src/loops/MaxQpsThread.C | 50 +++++++++--------------------- 3 files changed, 62 insertions(+), 51 deletions(-) diff --git a/framework/include/base/Assembly.h b/framework/include/base/Assembly.h index 2912f0542514..0e610a2283f4 100644 --- a/framework/include/base/Assembly.h +++ b/framework/include/base/Assembly.h @@ -568,6 +568,11 @@ class Assembly */ bool needDual() const { return _need_dual; } + /** + * Set the cached quadrature rules to nullptr + */ + void clearCachedQRules(); + private: /** * Set the qrule to be used for lower dimensional integration. @@ -595,6 +600,11 @@ class Assembly */ void reinit(const Elem * elem); + /** + * Set the volumetric quadrature rule based on the provided element + */ + void setVolumeQRule(const Elem * elem); + /** * Reinitialize FE data for the given element on the given side, optionally * with a given set of reference points @@ -664,6 +674,11 @@ class Assembly */ void reinit(const Elem * elem, const std::vector & reference_points); + /** + * Set the face quadrature rule based on the provided element and side + */ + void setFaceQRule(const Elem * const elem, const unsigned int side); + /** * Reinitialize the assembly data on an side of an element */ diff --git a/framework/src/base/Assembly.C b/framework/src/base/Assembly.C index b62e8f57b142..e142eda8649b 100644 --- a/framework/src/base/Assembly.C +++ b/framework/src/base/Assembly.C @@ -698,6 +698,15 @@ Assembly::setNeighborQRule(QBase * qrule, unsigned int dim) } } +void +Assembly::clearCachedQRules() +{ + _current_qrule = nullptr; + _current_qrule_face = nullptr; + _current_qrule_lower = nullptr; + _current_qrule_neighbor = nullptr; +} + void Assembly::setMortarQRule(Order order) { @@ -1750,6 +1759,16 @@ Assembly::reinitAtPhysical(const Elem * elem, const std::vector & physica _current_physical_points = physical_points; } +void +Assembly::setVolumeQRule(const Elem * const elem) +{ + unsigned int elem_dimension = elem->dim(); + _current_qrule_volume = qrules(elem_dimension).vol.get(); + // Make sure the qrule is the right one + if (_current_qrule != _current_qrule_volume) + setVolumeQRule(_current_qrule_volume, elem_dimension); +} + void Assembly::reinit(const Elem * elem) { @@ -1759,14 +1778,7 @@ Assembly::reinit(const Elem * elem) "current subdomain has been set incorrectly"); _current_elem_volume_computed = false; - unsigned int elem_dimension = elem->dim(); - - _current_qrule_volume = qrules(elem_dimension).vol.get(); - - // Make sure the qrule is the right one - if (_current_qrule != _current_qrule_volume) - setVolumeQRule(_current_qrule_volume, elem_dimension); - + setVolumeQRule(elem); reinitFE(elem); computeCurrentElemVolume(); @@ -1853,7 +1865,17 @@ Assembly::qruleArbitraryFace(const Elem * elem, unsigned int side) } void -Assembly::reinit(const Elem * elem, unsigned int side) +Assembly::setFaceQRule(const Elem * const elem, const unsigned int side) +{ + const auto elem_dimension = elem->dim(); + //// Make sure the qrule is the right one + auto rule = qruleFace(elem, side); + if (_current_qrule_face != rule) + setFaceQRule(rule, elem_dimension); +} + +void +Assembly::reinit(const Elem * const elem, const unsigned int side) { _current_elem = elem; _current_neighbor_elem = nullptr; @@ -1863,15 +1885,9 @@ Assembly::reinit(const Elem * elem, unsigned int side) _current_elem_volume_computed = false; _current_side_volume_computed = false; - unsigned int elem_dimension = elem->dim(); - _current_side_elem = &_current_side_elem_builder(*elem, side); - //// Make sure the qrule is the right one - auto rule = qruleFace(elem, side); - if (_current_qrule_face != rule) - setFaceQRule(rule, elem_dimension); - + setFaceQRule(elem, side); reinitFEFace(elem, side); computeCurrentFaceVolume(); diff --git a/framework/src/loops/MaxQpsThread.C b/framework/src/loops/MaxQpsThread.C index 5e542dbf3da5..2517a54dc10a 100644 --- a/framework/src/loops/MaxQpsThread.C +++ b/framework/src/loops/MaxQpsThread.C @@ -47,47 +47,22 @@ MaxQpsThread::operator()(const ConstElemRange & range) // This ensures we can access the correct qrules if any block-specific // qrules have been created. assem.setCurrentSubdomainID(elem->subdomain_id()); + assem.setVolumeQRule(elem); - FEType fe_type(FIRST, LAGRANGE); - unsigned int dim = elem->dim(); - unsigned int side = 0; // we assume that any element will have at least one side ;) - - // We cannot mess with the FE objects in Assembly, because we might need to request second - // derivatives - // later on. If we used them, we'd call reinit on them, thus making the call to request second - // derivatives harmful (i.e. leading to segfaults/asserts). Thus, we have to use a locally - // allocated object here. - // - // We'll use one for element interiors, which calculates nothing - std::unique_ptr fe(FEBase::build(dim, fe_type)); - fe->get_nothing(); - // Alex: This seems wrong, e.g. like we'd want the quadrature order to be higher if the - // polynomial basis is higher, but our fe_type up above is first Lagrange, so it seems like the - // person who wrote this code didn't want to consider higher order bases even if if the user - // asked for higher order bases - fe->add_p_level_in_reinit(false); - - // And another for element sides, which calculates the minimum - // libMesh currently allows for that - std::unique_ptr side_fe(FEBase::build(dim, fe_type)); - side_fe->get_xyz(); - side_fe->add_p_level_in_reinit(false); - - // figure out the number of qps for volume - auto qrule = assem.attachQRuleElem(dim, *fe); - fe->reinit(elem); + auto & qrule = assem.writeableQRule(); + qrule->init(elem->type(), elem->p_level()); if (qrule->n_points() > _max) _max = qrule->n_points(); - // figure out the number of qps for the face - // NOTE: user might specify higher order rule for faces, thus possibly ending up with more qps - // than in the volume - auto qrule_face = assem.attachQRuleFace(dim, *side_fe); - if (dim > 0) // side reinit in 0D makes no sense, but we may have NodeElems + if (elem->dim()) { - side_fe->reinit(elem, side); + // We are assuming here that side 0 ends up with at least as many quadrature points as any + // other side + assem.setFaceQRule(elem, /*side=*/0); + auto & qrule_face = assem.writeableQRuleFace(); + qrule_face->init(elem->side_type(0), elem->p_level()); if (qrule_face->n_points() > _max) - _max = qrule_face->n_points(); + _max = qrule->n_points(); } // In initial conditions nodes are enumerated as pretend quadrature points @@ -96,6 +71,11 @@ MaxQpsThread::operator()(const ConstElemRange & range) if (elem->n_nodes() > _max) _max = elem->n_nodes(); } + + // Clear the cached quadrature rules because we may add FE objects in between now and simulation + // start and we only ensure we set all the FE quadrature rules if a quadrature rule is different + // from the cached quadrature rule + assem.clearCachedQRules(); } void