www.mooseframework.org
InitialConditionTempl.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 "InitialConditionTempl.h"
11 #include "FEProblem.h"
12 #include "Assembly.h"
13 #include "MooseVariableFE.h"
14 
15 #include "libmesh/fe_interface.h"
16 #include "libmesh/quadrature.h"
17 
18 template <>
20 validParams<InitialConditionTempl<Real>>()
21 {
23 }
24 template <>
26 validParams<InitialConditionTempl<RealVectorValue>>()
27 {
29 }
30 
31 template <typename T>
33  : InitialConditionBase(parameters),
34  _fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
35  _tid(getParam<THREAD_ID>("_tid")),
36  _assembly(_fe_problem.assembly(_tid)),
37  _t(_fe_problem.time()),
38  _coord_sys(_assembly.coordSystem()),
39  _var(_sys.getFieldVariable<T>(parameters.get<THREAD_ID>("_tid"),
40  parameters.get<VariableName>("variable"))),
41  _current_elem(_var.currentElem()),
42  _current_node(nullptr),
43  _qp(0),
44  _fe_type(_var.feType())
45 {
46 }
47 
48 template <typename T>
50 {
51 }
52 
53 template <typename T>
54 void
56 {
57  // -- NOTE ----
58  // The following code is a copy from libMesh project_vector.C plus it adds some features, so we
59  // can couple variable values
60  // and we also do not call any callbacks, but we use our initial condition system directly.
61  // ------------
62 
63  // The dimension of the current element
64  _dim = _current_elem->dim();
65  // The element type
66  const ElemType elem_type = _current_elem->type();
67  // The number of nodes on the new element
68  const unsigned int n_nodes = _current_elem->n_nodes();
69 
70  // Get FE objects of the appropriate type
71  // We cannot use the FE object in Assembly, since the following code is messing with the
72  // quadrature rules
73  // for projections and would screw it up. However, if we implement projections from one mesh to
74  // another,
75  // this code should use that implementation.
76  std::unique_ptr<FEBaseType> fe(FEBaseType::build(_dim, _fe_type));
77 
78  // Prepare variables for projection
79  std::unique_ptr<QBase> qrule(_fe_type.default_quadrature_rule(_dim));
80  std::unique_ptr<QBase> qedgerule(_fe_type.default_quadrature_rule(1));
81  std::unique_ptr<QBase> qsiderule(_fe_type.default_quadrature_rule(_dim - 1));
82 
83  // The values of the shape functions at the quadrature points
84  _phi = &fe->get_phi();
85 
86  // The gradients of the shape functions at the quadrature points on the child element.
87  _dphi = nullptr;
88 
89  _cont = fe->get_continuity();
90 
91  if (_cont == C_ONE)
92  {
93  const std::vector<std::vector<GradientType>> & ref_dphi = fe->get_dphi();
94  _dphi = &ref_dphi;
95  }
96 
97  // The Jacobian * quadrature weight at the quadrature points
98  _JxW = &fe->get_JxW();
99  // The XYZ locations of the quadrature points
100  _xyz_values = &fe->get_xyz();
101 
102  // Update the DOF indices for this element based on the current mesh
103  _var.prepareIC();
104  _dof_indices = _var.dofIndices();
105 
106  // The number of DOFs on the element
107  const unsigned int n_dofs = _dof_indices.size();
108  if (n_dofs == 0)
109  return;
110 
111  // Fixed vs. free DoFs on edge/face projections
112  _dof_is_fixed.clear();
113  _dof_is_fixed.resize(n_dofs, false);
114  _free_dof.clear();
115  _free_dof.resize(n_dofs, 0);
116 
117  // Zero the interpolated values
118  _Ue.resize(n_dofs);
119  _Ue.zero();
120 
121  // In general, we need a series of
122  // projections to ensure a unique and continuous
123  // solution. We start by interpolating nodes, then
124  // hold those fixed and project edges, then
125  // hold those fixed and project faces, then
126  // hold those fixed and project interiors
127 
128  // Interpolate node values first
129  _current_dof = 0;
130  for (_n = 0; _n != n_nodes; ++_n)
131  {
132  // FIXME: this should go through the DofMap,
133  // not duplicate _dof_indices code badly!
134  _nc = FEInterface::n_dofs_at_node(_dim, _fe_type, elem_type, _n);
135  if (!_current_elem->is_vertex(_n))
136  {
137  _current_dof += _nc;
138  continue;
139  }
140  if (_cont == DISCONTINUOUS)
141  libmesh_assert(_nc == 0);
142  else if (_cont == C_ZERO)
143  setCZeroVertices();
144  else if (_fe_type.family == HERMITE)
145  setHermiteVertices();
146  else if (_cont == C_ONE)
147  setOtherCOneVertices();
148  else
149  libmesh_error();
150  } // loop over nodes
151 
152  // From here on out we won't be sampling at nodes anymore
153  _current_node = nullptr;
154 
155  // In 3D, project any edge values next
156  if (_dim > 2 && _cont != DISCONTINUOUS)
157  for (unsigned int e = 0; e != _current_elem->n_edges(); ++e)
158  {
159  FEInterface::dofs_on_edge(_current_elem, _dim, _fe_type, e, _side_dofs);
160 
161  // Some edge dofs are on nodes and already
162  // fixed, others are free to calculate
163  _free_dofs = 0;
164  for (unsigned int i = 0; i != _side_dofs.size(); ++i)
165  if (!_dof_is_fixed[_side_dofs[i]])
166  _free_dof[_free_dofs++] = i;
167 
168  // There may be nothing to project
169  if (!_free_dofs)
170  continue;
171 
172  // Initialize FE data on the edge
173  fe->attach_quadrature_rule(qedgerule.get());
174  fe->edge_reinit(_current_elem, e);
175  _n_qp = qedgerule->n_points();
176 
177  choleskySolve(false);
178  }
179 
180  // Project any side values (edges in 2D, faces in 3D)
181  if (_dim > 1 && _cont != DISCONTINUOUS)
182  for (unsigned int s = 0; s != _current_elem->n_sides(); ++s)
183  {
184  FEInterface::dofs_on_side(_current_elem, _dim, _fe_type, s, _side_dofs);
185 
186  // Some side dofs are on nodes/edges and already
187  // fixed, others are free to calculate
188  _free_dofs = 0;
189  for (unsigned int i = 0; i != _side_dofs.size(); ++i)
190  if (!_dof_is_fixed[_side_dofs[i]])
191  _free_dof[_free_dofs++] = i;
192 
193  // There may be nothing to project
194  if (!_free_dofs)
195  continue;
196 
197  // Initialize FE data on the side
198  fe->attach_quadrature_rule(qsiderule.get());
199  fe->reinit(_current_elem, s);
200  _n_qp = qsiderule->n_points();
201 
202  choleskySolve(false);
203  }
204 
205  // Project the interior values, finally
206 
207  // Some interior dofs are on nodes/edges/sides and
208  // already fixed, others are free to calculate
209  _free_dofs = 0;
210  for (unsigned int i = 0; i != n_dofs; ++i)
211  if (!_dof_is_fixed[i])
212  _free_dof[_free_dofs++] = i;
213 
214  // There may be nothing to project
215  if (_free_dofs)
216  {
217  // Initialize FE data
218  fe->attach_quadrature_rule(qrule.get());
219  fe->reinit(_current_elem);
220  _n_qp = qrule->n_points();
221 
222  choleskySolve(true);
223  } // if there are free interior dofs
224 
225  // Make sure every DoF got reached!
226  for (unsigned int i = 0; i != n_dofs; ++i)
227  libmesh_assert(_dof_is_fixed[i]);
228 
229  NumericVector<Number> & solution = _var.sys().solution();
230 
231  // 'first' and 'last' are no longer used, see note about subdomain-restricted variables below
232  // const dof_id_type
233  // first = solution.first_local_index(),
234  // last = solution.last_local_index();
235 
236  // Lock the new_vector since it is shared among threads.
237  {
238  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
239 
240  for (unsigned int i = 0; i < n_dofs; i++)
241  // We may be projecting a new zero value onto
242  // an old nonzero approximation - RHS
243  // if (_Ue(i) != 0.)
244 
245  // This is commented out because of subdomain restricted variables.
246  // It can be the case that if a subdomain restricted variable's boundary
247  // aligns perfectly with a processor boundary that the variable will get
248  // no value. To counteract this we're going to let every processor set a
249  // value at every node and then let PETSc figure it out.
250  // Later we can choose to do something different / better.
251  // if ((_dof_indices[i] >= first) && (_dof_indices[i] < last))
252  {
253  solution.set(_dof_indices[i], _Ue(i));
254  }
255  _var.setDofValues(_Ue);
256  }
257 }
258 
259 template <typename T>
260 void
262 {
263  // Assume that C_ZERO elements have a single nodal
264  // value shape function
265  libmesh_assert(_nc == 1);
266  _qp = _n;
267  _current_node = _current_elem->node_ptr(_n);
268  _Ue(_current_dof) = value(*_current_node);
269  _dof_is_fixed[_current_dof] = true;
270  _current_dof++;
271 }
272 
273 template <>
274 void
276 {
277  _qp = _n;
278  _current_node = _current_elem->node_ptr(_n);
279  auto point_value = value(*_current_node);
280  for (decltype(_nc) i = 0; i < _nc; ++i)
281  {
282  _Ue(_current_dof) = point_value(i);
283  _dof_is_fixed[_current_dof] = true;
284  _current_dof++;
285  }
286 }
287 
288 template <typename T>
289 void
291 {
292  // The hermite element vertex shape functions are weird
293  _qp = _n;
294  _current_node = _current_elem->node_ptr(_n);
295  _Ue(_current_dof) = value(*_current_node);
296  _dof_is_fixed[_current_dof] = true;
297  _current_dof++;
298  Gradient grad = gradient(*_current_node);
299  // x derivative
300  _Ue(_current_dof) = grad(0);
301  _dof_is_fixed[_current_dof] = true;
302  _current_dof++;
303  if (_dim > 1)
304  {
305  // We'll finite difference mixed derivatives
306  Point nxminus = _current_elem->point(_n), nxplus = _current_elem->point(_n);
307  nxminus(0) -= TOLERANCE;
308  nxplus(0) += TOLERANCE;
309  Gradient gxminus = gradient(nxminus);
310  Gradient gxplus = gradient(nxplus);
311  // y derivative
312  _Ue(_current_dof) = grad(1);
313  _dof_is_fixed[_current_dof] = true;
314  _current_dof++;
315  // xy derivative
316  _Ue(_current_dof) = (gxplus(1) - gxminus(1)) / 2. / TOLERANCE;
317  _dof_is_fixed[_current_dof] = true;
318  _current_dof++;
319 
320  if (_dim > 2)
321  {
322  // z derivative
323  _Ue(_current_dof) = grad(2);
324  _dof_is_fixed[_current_dof] = true;
325  _current_dof++;
326  // xz derivative
327  _Ue(_current_dof) = (gxplus(2) - gxminus(2)) / 2. / TOLERANCE;
328  _dof_is_fixed[_current_dof] = true;
329  _current_dof++;
330  // We need new points for yz
331  Point nyminus = _current_elem->point(_n), nyplus = _current_elem->point(_n);
332  nyminus(1) -= TOLERANCE;
333  nyplus(1) += TOLERANCE;
334  Gradient gyminus = gradient(nyminus);
335  Gradient gyplus = gradient(nyplus);
336  // xz derivative
337  _Ue(_current_dof) = (gyplus(2) - gyminus(2)) / 2. / TOLERANCE;
338  _dof_is_fixed[_current_dof] = true;
339  _current_dof++;
340  // Getting a 2nd order xyz is more tedious
341  Point nxmym = _current_elem->point(_n), nxmyp = _current_elem->point(_n),
342  nxpym = _current_elem->point(_n), nxpyp = _current_elem->point(_n);
343  nxmym(0) -= TOLERANCE;
344  nxmym(1) -= TOLERANCE;
345  nxmyp(0) -= TOLERANCE;
346  nxmyp(1) += TOLERANCE;
347  nxpym(0) += TOLERANCE;
348  nxpym(1) -= TOLERANCE;
349  nxpyp(0) += TOLERANCE;
350  nxpyp(1) += TOLERANCE;
351  Gradient gxmym = gradient(nxmym);
352  Gradient gxmyp = gradient(nxmyp);
353  Gradient gxpym = gradient(nxpym);
354  Gradient gxpyp = gradient(nxpyp);
355  Number gxzplus = (gxpyp(2) - gxmyp(2)) / 2. / TOLERANCE;
356  Number gxzminus = (gxpym(2) - gxmym(2)) / 2. / TOLERANCE;
357  // xyz derivative
358  _Ue(_current_dof) = (gxzplus - gxzminus) / 2. / TOLERANCE;
359  _dof_is_fixed[_current_dof] = true;
360  _current_dof++;
361  }
362  }
363 }
364 
365 template <>
366 void
368 {
369 }
370 
371 template <typename T>
372 void
374 {
375  // Assume that other C_ONE elements have a single nodal
376  // value shape function and nodal gradient component
377  // shape functions
378  libmesh_assert(_nc == 1 + _dim);
379  _current_node = _current_elem->node_ptr(_n);
380  _Ue(_current_dof) = value(*_current_node);
381  _dof_is_fixed[_current_dof] = true;
382  _current_dof++;
383  Gradient grad = gradient(*_current_node);
384  for (unsigned int i = 0; i != _dim; ++i)
385  {
386  _Ue(_current_dof) = grad(i);
387  _dof_is_fixed[_current_dof] = true;
388  _current_dof++;
389  }
390 }
391 
392 template <>
393 void
395 {
396 }
397 
398 template <typename T>
399 Real
401 {
402  return op1 * op2;
403 }
404 
405 template <>
406 Real
408  const GradientType & op2)
409 {
410  return op1.contract(op2);
411 }
412 
413 template <typename T>
414 void
416 {
417  _Ke.resize(_free_dofs, _free_dofs);
418  _Ke.zero();
419  _Fe.resize(_free_dofs);
420  _Fe.zero();
421  // The new edge coefficients
422  DenseVector<Number> U(_free_dofs);
423 
424  // Loop over the quadrature points
425  for (_qp = 0; _qp < _n_qp; _qp++)
426  {
427  // solution at the quadrature point
428  auto fineval = value((*_xyz_values)[_qp]);
429  // solution grad at the quadrature point
430  GradientType finegrad;
431  if (_cont == C_ONE)
432  finegrad = gradient((*_xyz_values)[_qp]);
433 
434  auto dofs_size = is_volume ? _dof_indices.size() : _side_dofs.size();
435 
436  // Form edge projection matrix
437  for (decltype(dofs_size) geomi = 0, freei = 0; geomi != dofs_size; ++geomi)
438  {
439  auto i = is_volume ? geomi : _side_dofs[geomi];
440 
441  // fixed DoFs aren't test functions
442  if (_dof_is_fixed[i])
443  continue;
444  for (decltype(dofs_size) geomj = 0, freej = 0; geomj != dofs_size; ++geomj)
445  {
446  auto j = is_volume ? geomj : _side_dofs[geomj];
447  if (_dof_is_fixed[j])
448  _Fe(freei) -= (*_phi)[i][_qp] * (*_phi)[j][_qp] * (*_JxW)[_qp] * _Ue(j);
449  else
450  _Ke(freei, freej) += (*_phi)[i][_qp] * (*_phi)[j][_qp] * (*_JxW)[_qp];
451  if (_cont == C_ONE)
452  {
453  if (_dof_is_fixed[j])
454  _Fe(freei) -= dotHelper((*_dphi)[i][_qp], (*_dphi)[j][_qp]) * (*_JxW)[_qp] * _Ue(j);
455  else
456  _Ke(freei, freej) += dotHelper((*_dphi)[i][_qp], (*_dphi)[j][_qp]) * (*_JxW)[_qp];
457  }
458  if (!_dof_is_fixed[j])
459  freej++;
460  }
461  _Fe(freei) += (*_phi)[i][_qp] * fineval * (*_JxW)[_qp];
462  if (_cont == C_ONE)
463  _Fe(freei) += dotHelper(finegrad, (*_dphi)[i][_qp]) * (*_JxW)[_qp];
464  freei++;
465  }
466  }
467 
468  _Ke.cholesky_solve(_Fe, U);
469 
470  // Transfer new edge solutions to element
471  for (unsigned int i = 0; i != _free_dofs; ++i)
472  {
473  auto the_dof = is_volume ? _free_dof[i] : _side_dofs[_free_dof[i]];
474  Number & ui = _Ue(the_dof);
475  libmesh_assert(std::abs(ui) < TOLERANCE || std::abs(ui - U(i)) < TOLERANCE);
476  ui = U(i);
477  _dof_is_fixed[the_dof] = true;
478  }
479 }
480 
481 template <typename T>
482 void
484 {
485  _var.reinitNode();
486  _var.computeNodalValues(); // has to call this to resize the internal array
487  auto return_value = value(p);
488 
489  _var.setNodalValue(return_value); // update variable data, which is referenced by others, so the
490  // value is up-to-date
491 
492  // We are done, so update the solution vector
493  {
494  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
495  _var.insert(_var.sys().solution());
496  }
497 }
498 
499 template class InitialConditionTempl<Real>;
MetaPhysicL::DualNumber< T, D > abs(const MetaPhysicL::DualNumber< T, D > &in)
OutputTools< T >::OutputGradient GradientType
void setOtherCOneVertices()
set the temporary solution vector for node projections of non-Hermitian C1 variables ...
This is a template class that implements the workhorse compute and computeNodal methods.
virtual void computeNodal(const Point &p) override
Workhorse method for projecting the initial conditions for boundary restricted initial conditions...
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
InitialConditionBase serves as the abstract base class for InitialConditions and VectorInitialConditi...
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
void setCZeroVertices()
set the temporary solution vector for node projections of C0 variables
void setHermiteVertices()
set the temporary solution vector for node projections of Hermite variables
InitialConditionTempl(const InputParameters &parameters)
Constructor.
virtual void compute() override
Workhorse method for projecting the initial conditions for block initial conditions.
void choleskySolve(bool is_volume)
Perform the cholesky solves for edge, side, and interior projections.
Real dotHelper(const GradientType &op1, const GradientType &op2)
Helps perform multiplication of GradientTypes: a normal dot product for vectors and a contraction for...
InputParameters validParams< InitialConditionBase >()
unsigned int THREAD_ID
Definition: MooseTypes.h:161