www.mooseframework.org
NodalPatchRecovery.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "NodalPatchRecovery.h"
11 #include "SwapBackSentinel.h"
12 
15 {
17  MooseEnum orders("CONSTANT FIRST SECOND THIRD FOURTH");
18  params.addParam<MooseEnum>("patch_polynomial_order",
19  orders,
20  "Polynomial order used in least squares fitting of material property "
21  "over the local patch of elements connected to a given node");
22  params.addParamNamesToGroup("patch_polynomial_order", "Advanced");
23 
24  // TODO make this class work with relationship manager
25  // params.registerRelationshipManagers("ElementSideNeighborLayers");
26  // params.addPrivateParam<unsigned int short>("element_side_neighbor_layers", 2);
27 
28  return params;
29 }
30 
32  : AuxKernel(parameters),
33  _patch_polynomial_order(
34  isParamValid("patch_polynomial_order")
35  ? static_cast<unsigned int>(getParam<MooseEnum>("patch_polynomial_order"))
36  : static_cast<unsigned int>(_var.order())),
37  _fe_problem(*parameters.get<FEProblemBase *>("_fe_problem_base")),
38  _multi_index(multiIndex(_mesh.dimension(), _patch_polynomial_order))
39 {
40  // if the patch polynomial order is lower than the variable order,
41  // it is very likely that the patch recovery is not used at its max accuracy
42  if (_patch_polynomial_order < static_cast<unsigned int>(_var.order()))
43  mooseWarning("Specified 'patch_polynomial_order' is lower than the AuxVariable's order");
44 
45  // TODO remove the manual ghosting once relationship manager is working correctly
46  // no need to ghost if this aux is elemental
47  if (isNodal())
48  {
49  MeshBase & meshhelper = _mesh.getMesh();
50  meshhelper.allow_renumbering(false);
51  for (const auto & elem :
52  as_range(meshhelper.semilocal_elements_begin(), meshhelper.semilocal_elements_end()))
53  _fe_problem.addGhostedElem(elem->id());
54  }
55 }
56 
57 std::vector<std::vector<unsigned int>>
58 NodalPatchRecovery::nChooseK(unsigned int N, unsigned int K)
59 {
60  std::vector<std::vector<unsigned int>> n_choose_k;
61  std::vector<unsigned int> row;
62  std::string bitmask(K, 1); // K leading 1's
63  bitmask.resize(N, 0); // N-K trailing 0's
64 
65  do
66  {
67  row.clear();
68  row.push_back(0);
69  for (unsigned int i = 0; i < N; ++i) // [0..N-1] integers
70  if (bitmask[i])
71  row.push_back(i + 1);
72  row.push_back(N + 1);
73  n_choose_k.push_back(row);
74  } while (std::prev_permutation(bitmask.begin(), bitmask.end()));
75 
76  return n_choose_k;
77 }
78 
79 std::vector<std::vector<unsigned int>>
80 NodalPatchRecovery::multiIndex(unsigned int dim, unsigned int order)
81 {
82  // first row all zero
83  std::vector<std::vector<unsigned int>> multi_index;
84  std::vector<std::vector<unsigned int>> n_choose_k;
85  std::vector<unsigned int> row(dim, 0);
86  multi_index.push_back(row);
87 
88  if (dim == 1)
89  for (unsigned int q = 1; q <= order; q++)
90  {
91  row[0] = q;
92  multi_index.push_back(row);
93  }
94  else
95  for (unsigned int q = 1; q <= order; q++)
96  {
97  n_choose_k = nChooseK(dim + q - 1, dim - 1);
98  for (unsigned int r = 0; r < n_choose_k.size(); r++)
99  {
100  row.clear();
101  for (unsigned int c = 1; c < n_choose_k[0].size(); c++)
102  row.push_back(n_choose_k[r][c] - n_choose_k[r][c - 1] - 1);
103  multi_index.push_back(row);
104  }
105  }
106 
107  return multi_index;
108 }
109 
110 void
112 {
113  _P.resize(_multi_index.size());
115  for (unsigned int r = 0; r < _multi_index.size(); r++)
116  {
117  polynomial = 1.0;
118  for (unsigned int c = 0; c < _multi_index[0].size(); c++)
119  for (unsigned int p = 0; p < _multi_index[r][c]; p++)
120  polynomial *= q_point(c);
121  _P(r) = polynomial;
122  }
123 }
124 
125 void
127 {
129  PxP.resize(_multi_index.size(), _multi_index.size());
130  PxP.outer_product(_P, _P);
131  _A.add(1.0, PxP);
132 }
133 
134 void
136 {
137  _B.add(val, _P);
138 }
139 
140 void
142 {
143  // clean and resize P, A, B, coef
144  _P.resize(_multi_index.size());
145  _A.resize(_multi_index.size(), _multi_index.size());
146  _B.resize(_multi_index.size());
147  _coef.resize(_multi_index.size());
148 
149  // activate dependent material properties
150  std::unordered_set<unsigned int> needed_mat_props;
151  const auto & mp_deps = getMatPropDependencies();
152  needed_mat_props.insert(mp_deps.begin(), mp_deps.end());
153  _fe_problem.setActiveMaterialProperties(needed_mat_props, _tid);
154 }
155 
156 void
158 {
159  // Fall back on standard procedure for element variables
160  if (!isNodal())
161  {
163  return;
164  }
165 
166  // Limit current use of NodalPatchRecovery to a single processor
167  if (_communicator.size() > 1)
168  mooseError("The nodal patch recovery option, which calculates the Zienkiewicz-Zhu patch "
169  "recovery for nodal variables (family = LAGRANGE), is not currently implemented for "
170  "parallel runs. Run in serial if you must use the nodal patch capability");
171 
172  // Use Zienkiewicz-Zhu patch recovery for nodal variables
173  reinitPatch();
174 
175  // get node-to-conneted-elem map
176  const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map = _mesh.nodeToElemMap();
177  auto node_to_elem_pair = node_to_elem_map.find(_current_node->id());
178  mooseAssert(node_to_elem_pair != node_to_elem_map.end(), "Missing entry in node to elem map");
179  std::vector<dof_id_type> elem_ids = node_to_elem_pair->second;
180 
181  // consider the case for corner node
182  if (elem_ids.size() == 1)
183  {
184  dof_id_type elem_id = elem_ids[0];
185  for (auto & n : _mesh.elemPtr(elem_id)->node_ref_range())
186  {
187  node_to_elem_pair = node_to_elem_map.find(n.id());
188  std::vector<dof_id_type> elem_ids_candidate = node_to_elem_pair->second;
189  if (elem_ids_candidate.size() > elem_ids.size())
190  elem_ids = elem_ids_candidate;
191  }
192  }
193 
194  // before we go, check if we have enough sample points for solving the least square fitting
195  if (_q_point.size() * elem_ids.size() < _multi_index.size())
196  mooseError("There are not enough sample points to recover the nodal value, try reducing the "
197  "polynomial order or using a higher-order quadrature scheme.");
198 
199  // general treatment for side nodes and internal nodes
200  for (auto elem_id : elem_ids)
201  {
202  const Elem * elem = _mesh.elemPtr(elem_id);
203 
204  if (!elem->is_semilocal(_mesh.processor_id()))
205  mooseError("skipped non local elem!");
206 
207  _fe_problem.prepare(elem, _tid);
208  _fe_problem.reinitElem(elem, _tid);
209 
210  // Set up Sentinel class so that, even if reinitMaterials() throws, we
211  // still remember to swap back during stack unwinding.
213  _fe_problem.reinitMaterials(elem->subdomain_id(), _tid);
214 
215  for (_qp = 0; _qp < _q_point.size(); _qp++)
216  {
220  }
221  }
222 
224 
225  // solve for the coefficients of least square fitting
226  // as we are on the physical coordinates, A is not always SPD
227  _A.svd_solve(_B, _coef);
228 
229  // compute fitted nodal value
231  Real nodal_value = _P.dot(_coef);
232 
233  // set nodal value
236  {
237  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
238  _solution.set(dof, nodal_value);
239  }
240  _var.setNodalValue(nodal_value);
241 }
NodalPatchRecovery(const InputParameters &parameters)
void setActiveMaterialProperties(const std::unordered_set< unsigned int > &mat_prop_ids, const THREAD_ID tid)
Record and set the material properties required by the current computing thread.
virtual void addGhostedElem(dof_id_type elem_id) override
Will make sure that all dofs connected to elem_id are ghosted to this processor.
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseVector< T3 > &vec)
THREAD_ID _tid
Thread ID.
Definition: AuxKernel.h:171
virtual void reinitNode(const Node *node, const THREAD_ID tid) override
virtual Elem * elemPtr(const dof_id_type i)
Definition: MooseMesh.C:2863
MooseMesh & _mesh
Mesh this kernel is active on.
Definition: AuxKernel.h:188
virtual void setNodalValue(const OutputType &value, unsigned int idx=0)=0
const Node *const & _current_node
Current node (valid only for nodal kernels)
Definition: AuxKernel.h:214
void resize(const unsigned int n)
T * get(const std::unique_ptr< T > &u)
The MooseUtils::get() specializations are used to support making forwards-compatible code changes fro...
Definition: MooseUtils.h:1147
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:148
const Parallel::Communicator & _communicator
void clearActiveMaterialProperties(const THREAD_ID tid)
Clear the active material properties.
virtual const dof_id_type & nodalDofIndex() const =0
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
FEProblemBase & _fe_problem
void mooseWarning(Args &&... args) const
Emits a warning prefixed with object name and type.
void outer_product(const DenseVector< Number > &a, const DenseVector< Number > &b)
virtual std::vector< std::vector< unsigned int > > multiIndex(unsigned int dim, unsigned int order)
generate a complete multi index table for given dimension and order i.e.
virtual void compute()
Computes the value and stores it in the solution vector.
Definition: AuxKernel.C:296
processor_id_type size() const
unsigned int size() const
The number of elements that can currently be stored in the array.
Definition: MooseArray.h:256
virtual Real computeValue()=0
Compute and return the value of the aux variable.
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
Definition: MooseMesh.C:3198
virtual void reinitElem(const Elem *elem, const THREAD_ID tid) override
virtual void accumulateBVector(Real val)
calculate the patch load vector as ^n P^Tval where n is the number of quadrature points in the elemen...
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseMatrix< T3 > &mat)
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:31
R polynomial(const C &c, const T x)
Evaluate a polynomial with the coefficients c at x.
Definition: MathUtils.h:268
DenseVector< Number > _P
virtual void compute() override
solve the equation Ac = B where c is the coefficients vector from least square fitting nodal value is...
Order order() const
Get the order of this variable Note: Order enum can be implicitly converted to unsigned int...
DenseVector< Number > _B
NumericVector< Number > & _solution
reference to the solution vector of auxiliary system
Definition: AuxKernel.h:220
CompareTypes< Number, T2 >::supertype dot(const DenseVector< T2 > &vec) const
virtual void swapBackMaterials(const THREAD_ID tid)
static InputParameters validParams()
MooseVariableField< Real > & _var
This is a regular kernel so we cast to a regular MooseVariable.
Definition: AuxKernel.h:174
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void reinitMaterials(SubdomainID blk_id, const THREAD_ID tid, bool swap_stateful=true)
virtual void prepare(const Elem *elem, const THREAD_ID tid) override
std::vector< std::vector< unsigned int > > multiIndex(unsigned int dim, unsigned int order)
generate a complete multi index table for given dimension and order i.e.
Definition: MathUtils.C:175
void resize(const unsigned int new_m, const unsigned int new_n)
void mooseError(Args &&... args) const
Emits an error prefixed with object name and type.
unsigned int _qp
Quadrature point index.
Definition: AuxKernel.h:230
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an option parameter and a documentation string to the InputParameters object...
static InputParameters validParams()
Definition: AuxKernel.C:27
const unsigned int _patch_polynomial_order
polynomial order, default is variable order
virtual void set(const numeric_index_type i, const Number value)=0
virtual void accumulateAMatrix()
Calculate the patch stiffness matrix as ^n P^TP where n is the number of quadrature points in the ele...
const std::unordered_set< unsigned int > & getMatPropDependencies() const
Retrieve the set of material properties that this object depends on.
virtual void reinitPatch()
reserve space for A, B, P, and prepare required material properties
DenseMatrix< Number > _A
processor_id_type processor_id() const
const MooseArray< Point > & _q_point
Active quadrature points.
Definition: AuxKernel.h:196
std::vector< std::vector< unsigned int > > _multi_index
virtual void computePVector(Point q_point)
compute the P vector at given point i.e.
void svd_solve(const DenseVector< Number > &rhs, DenseVector< Number > &x, Real rcond=std::numeric_limits< Real >::epsilon()) const
void ErrorVector unsigned int
bool isNodal() const
Nodal or elemental kernel?
Definition: AuxKernel.h:86
virtual std::vector< std::vector< unsigned int > > nChooseK(unsigned int N, unsigned int K)
Find out how many different ways you can choose K items from N items set without repetition and witho...
uint8_t dof_id_type
const std::map< dof_id_type, std::vector< dof_id_type > > & nodeToElemMap()
If not already created, creates a map from every node to all elements to which they are connected...
Definition: MooseMesh.C:980
void addParamNamesToGroup(const std::string &space_delim_names, const std::string group_name)
This method takes a space delimited list of parameter names and adds them to the specified group name...
The "SwapBackSentinel" class&#39;s destructor guarantees that FEProblemBase::swapBackMaterials{Face,Neighbor}() is called even when an exception is thrown from FEProblemBase::reinitMaterials{Face,Neighbor}.
DenseVector< Number > _coef