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 
13 template <>
16 {
18  MooseEnum orders("CONSTANT FIRST SECOND THIRD FOURTH");
19  params.addParam<MooseEnum>("patch_polynomial_order",
20  orders,
21  "Polynomial order used in least squares fitting of material property "
22  "over the local patch of elements connected to a given node");
23  params.addParamNamesToGroup("patch_polynomial_order", "Advanced");
24 
25  // TODO make this class work with relationship manager
26  // params.registerRelationshipManagers("ElementSideNeighborLayers");
27  // params.addPrivateParam<unsigned int short>("element_side_neighbor_layers", 2);
28 
29  return params;
30 }
31 
33  : AuxKernel(parameters),
34  _patch_polynomial_order(
35  isParamValid("patch_polynomial_order")
36  ? static_cast<unsigned int>(getParam<MooseEnum>("patch_polynomial_order"))
37  : static_cast<unsigned int>(_var.order())),
38  _fe_problem(*parameters.get<FEProblemBase *>("_fe_problem_base")),
39  _multi_index(multiIndex(_mesh.dimension(), _patch_polynomial_order))
40 {
41  // if the patch polynomial order is lower than the variable order,
42  // it is very likely that the patch recovery is not used at its max accuracy
43  if (_patch_polynomial_order < static_cast<unsigned int>(_var.order()))
44  mooseWarning("Specified 'patch_polynomial_order' is lower than the AuxVariable's order");
45 
46  // TODO remove the manual ghosting once relationship manager is working correctly
47  // no need to ghost if this aux is elemental
48  if (isNodal())
49  {
50  MeshBase & meshhelper = _mesh.getMesh();
51  meshhelper.allow_renumbering(false);
52  for (const auto & elem :
53  as_range(meshhelper.semilocal_elements_begin(), meshhelper.semilocal_elements_end()))
54  _fe_problem.addGhostedElem(elem->id());
55  }
56 }
57 
58 std::vector<std::vector<unsigned int>>
59 NodalPatchRecovery::nChooseK(unsigned int N, unsigned int K)
60 {
61  std::vector<std::vector<unsigned int>> n_choose_k;
62  std::vector<unsigned int> row;
63  std::string bitmask(K, 1); // K leading 1's
64  bitmask.resize(N, 0); // N-K trailing 0's
65 
66  do
67  {
68  row.clear();
69  row.push_back(0);
70  for (unsigned int i = 0; i < N; ++i) // [0..N-1] integers
71  if (bitmask[i])
72  row.push_back(i + 1);
73  row.push_back(N + 1);
74  n_choose_k.push_back(row);
75  } while (std::prev_permutation(bitmask.begin(), bitmask.end()));
76 
77  return n_choose_k;
78 }
79 
80 std::vector<std::vector<unsigned int>>
81 NodalPatchRecovery::multiIndex(unsigned int dim, unsigned int order)
82 {
83  // first row all zero
84  std::vector<std::vector<unsigned int>> multi_index;
85  std::vector<std::vector<unsigned int>> n_choose_k;
86  std::vector<unsigned int> row(dim, 0);
87  multi_index.push_back(row);
88 
89  if (dim == 1)
90  for (unsigned int q = 1; q <= order; q++)
91  {
92  row[0] = q;
93  multi_index.push_back(row);
94  }
95  else
96  for (unsigned int q = 1; q <= order; q++)
97  {
98  n_choose_k = nChooseK(dim + q - 1, dim - 1);
99  for (unsigned int r = 0; r < n_choose_k.size(); r++)
100  {
101  row.clear();
102  for (unsigned int c = 1; c < n_choose_k[0].size(); c++)
103  row.push_back(n_choose_k[r][c] - n_choose_k[r][c - 1] - 1);
104  multi_index.push_back(row);
105  }
106  }
107 
108  return multi_index;
109 }
110 
111 void
113 {
114  _P.resize(_multi_index.size());
115  Real polynomial;
116  for (unsigned int r = 0; r < _multi_index.size(); r++)
117  {
118  polynomial = 1.0;
119  for (unsigned int c = 0; c < _multi_index[0].size(); c++)
120  for (unsigned int p = 0; p < _multi_index[r][c]; p++)
121  polynomial *= q_point(c);
122  _P(r) = polynomial;
123  }
124 }
125 
126 void
128 {
129  DenseMatrix<Number> PxP;
130  PxP.resize(_multi_index.size(), _multi_index.size());
131  PxP.outer_product(_P, _P);
132  _A.add(1.0, PxP);
133 }
134 
135 void
137 {
138  _B.add(val, _P);
139 }
140 
141 void
143 {
144  // clean and resize P, A, B, coef
145  _P.resize(_multi_index.size());
146  _A.resize(_multi_index.size(), _multi_index.size());
147  _B.resize(_multi_index.size());
148  _coef.resize(_multi_index.size());
149 
150  // activate dependent material properties
151  std::set<unsigned int> needed_mat_props;
152  const std::set<unsigned int> & mp_deps = getMatPropDependencies();
153  needed_mat_props.insert(mp_deps.begin(), mp_deps.end());
154  _fe_problem.setActiveMaterialProperties(needed_mat_props, _tid);
155 }
156 
157 void
159 {
160  // Fall back on standard procedure for element variables
161  if (!isNodal())
162  {
164  return;
165  }
166 
167  // Limit current use of NodalPatchRecovery to a single processor
168  if (_communicator.size() > 1)
169  mooseError("The nodal patch recovery option, which calculates the Zienkiewicz-Zhu patch "
170  "recovery for nodal variables (family = LAGRANGE), is not currently implemented for "
171  "parallel runs. Run in serial if you must use the nodal patch capability");
172 
173  // Use Zienkiewicz-Zhu patch recovery for nodal variables
174  reinitPatch();
175 
176  // get node-to-conneted-elem map
177  const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map = _mesh.nodeToElemMap();
178  auto node_to_elem_pair = node_to_elem_map.find(_current_node->id());
179  mooseAssert(node_to_elem_pair != node_to_elem_map.end(), "Missing entry in node to elem map");
180  std::vector<dof_id_type> elem_ids = node_to_elem_pair->second;
181 
182  // consider the case for corner node
183  if (elem_ids.size() == 1)
184  {
185  dof_id_type elem_id = elem_ids[0];
186  for (auto & n : _mesh.elemPtr(elem_id)->node_ref_range())
187  {
188  node_to_elem_pair = node_to_elem_map.find(n.id());
189  std::vector<dof_id_type> elem_ids_candidate = node_to_elem_pair->second;
190  if (elem_ids_candidate.size() > elem_ids.size())
191  elem_ids = elem_ids_candidate;
192  }
193  }
194 
195  // before we go, check if we have enough sample points for solving the least square fitting
196  if (_q_point.size() * elem_ids.size() < _multi_index.size())
197  mooseError("There are not enough sample points to recover the nodal value, try reducing the "
198  "polynomial order or using a higher-order quadrature scheme.");
199 
200  // general treatment for side nodes and internal nodes
201  for (auto elem_id : elem_ids)
202  {
203  const Elem * elem = _mesh.elemPtr(elem_id);
204 
205  if (!elem->is_semilocal(_mesh.processor_id()))
206  mooseError("skipped non local elem!");
207 
208  _fe_problem.prepare(elem, _tid);
209  _fe_problem.reinitElem(elem, _tid);
210 
211  // Set up Sentinel class so that, even if reinitMaterials() throws, we
212  // still remember to swap back during stack unwinding.
214  _fe_problem.reinitMaterials(elem->subdomain_id(), _tid);
215 
216  for (_qp = 0; _qp < _q_point.size(); _qp++)
217  {
221  }
222  }
223 
225 
226  // solve for the coefficients of least square fitting
227  // as we are on the physical coordinates, A is not always SPD
228  _A.svd_solve(_B, _coef);
229 
230  // compute fitted nodal value
232  Real nodal_value = _P.dot(_coef);
233 
234  // set nodal value
236  dof_id_type dof = _var.nodalDofIndex();
237  {
238  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
239  _solution.set(dof, nodal_value);
240  }
241  _var.setNodalValue(nodal_value);
242 }
NodalPatchRecovery(const InputParameters &parameters)
MooseVariableFE< ComputeValueType > & _var
This is a regular kernel so we cast to a regular MooseVariable.
Definition: AuxKernel.h:167
virtual void prepare(const Elem *elem, THREAD_ID tid) override
virtual void addGhostedElem(dof_id_type elem_id) override
Will make sure that all dofs connected to elem_id are ghosted to this processor.
THREAD_ID _tid
Thread ID.
Definition: AuxKernel.h:164
PetscInt N
virtual Elem * elemPtr(const dof_id_type i)
Definition: MooseMesh.C:2267
void mooseWarning(Args &&... args) const
Definition: MooseObject.h:155
MooseMesh & _mesh
Mesh this kernel is active on.
Definition: AuxKernel.h:190
const Node *const & _current_node
Current node (valid only for nodal kernels)
Definition: AuxKernel.h:213
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
void mooseError(Args &&... args) const
Definition: MooseObject.h:147
virtual void reinitElem(const Elem *elem, THREAD_ID tid) override
virtual void reinitMaterials(SubdomainID blk_id, THREAD_ID tid, bool swap_stateful=true)
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
FEProblemBase & _fe_problem
virtual void setActiveMaterialProperties(const std::set< unsigned int > &mat_prop_ids, THREAD_ID tid) override
Record and set the material properties required by the current computing thread.
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:276
const std::set< unsigned int > & getMatPropDependencies() const
Retrieve the set of material properties that this object depends on.
unsigned int size() const
The number of elements that can currently be stored in the array.
Definition: MooseArray.h:259
virtual ComputeValueType computeValue()=0
Compute and return the value of the aux variable.
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
Definition: MooseMesh.C:2567
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...
virtual void clearActiveMaterialProperties(THREAD_ID tid) override
Clear the active material properties.
InputParameters validParams< AuxKernel >()
Definition: AuxKernel.C:25
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:31
bool isNodal()
Nodal or elemental kernel?
Definition: AuxKernel.h:84
virtual void swapBackMaterials(THREAD_ID tid)
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:216
const dof_id_type & nodalDofIndex() const override
PetscInt n
void setNodalValue(OutputType value, unsigned int idx=0)
InputParameters validParams< NodalPatchRecovery >()
unsigned int _qp
Quadrature point index.
Definition: AuxKernel.h:219
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...
const unsigned int _patch_polynomial_order
polynomial order, default is variable order
Base class for creating new auxiliary kernels and auxiliary boundary conditions.
Definition: AuxKernel.h:33
virtual void accumulateAMatrix()
Calculate the patch stiffness matrix as ^n P^TP where n is the number of quadrature points in the ele...
virtual void reinitNode(const Node *node, THREAD_ID tid) override
virtual void reinitPatch()
reserve space for A, B, P, and prepare required material properties
DenseMatrix< Number > _A
const MooseArray< Point > & _q_point
Dimension of the problem being solved.
Definition: AuxKernel.h:195
std::vector< std::vector< unsigned int > > _multi_index
virtual void computePVector(Point q_point)
compute the P vector at given point i.e.
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...
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:690
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