https://mooseframework.inl.gov
InitialConditionTempl.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
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 #include "SystemBase.h"
15 
16 #include "libmesh/fe_interface.h"
17 #include "libmesh/quadrature.h"
18 
19 using namespace libMesh;
20 
21 template <typename T>
23  : InitialConditionBase(parameters),
24  _fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
25  _tid(getParam<THREAD_ID>("_tid")),
26  _t(_fe_problem.time()),
27  _var(_sys.getActualFieldVariable<T>(parameters.get<THREAD_ID>("_tid"),
28  parameters.get<VariableName>("variable"))),
29  _fe_var(dynamic_cast<MooseVariableFE<T> *>(&_var)),
30  _assembly(
31  _fe_problem.assembly(_tid, _var.kind() == Moose::VAR_SOLVER ? _var.sys().number() : 0)),
32  _coord_sys(_assembly.coordSystem()),
33  _current_elem(_var.currentElem()),
34  _current_elem_volume(_assembly.elemVolume()),
35  _current_node(nullptr),
36  _qp(0),
37  _fe_type(_var.feType())
38 {
39 }
40 
41 template <typename T>
43 {
44 }
45 
46 template <typename T>
47 void
49 {
50  // -- NOTE ----
51  // The following code is a copy from libMesh project_vector.C plus it adds some features, so we
52  // can couple variable values
53  // and we also do not call any callbacks, but we use our initial condition system directly.
54  // ------------
55 
56  // The dimension of the current element
57  _dim = _current_elem->dim();
58  // The number of nodes on the new element
59  const unsigned int n_nodes = _current_elem->n_nodes();
60 
61  // Get FE objects of the appropriate type
62  // We cannot use the FE object in Assembly, since the following code is messing with the
63  // quadrature rules
64  // for projections and would screw it up. However, if we implement projections from one mesh to
65  // another,
66  // this code should use that implementation.
67  std::unique_ptr<FEBaseType> fe(FEBaseType::build(_dim, _fe_type));
68 
69  // Prepare variables for projection
70  std::unique_ptr<QBase> qrule(_fe_type.default_quadrature_rule(_dim));
71  std::unique_ptr<QBase> qedgerule(_fe_type.default_quadrature_rule(1));
72  std::unique_ptr<QBase> qsiderule(_fe_type.default_quadrature_rule(_dim - 1));
73 
74  // The values of the shape functions at the quadrature points
75  _phi = &fe->get_phi();
76 
77  // The gradients of the shape functions at the quadrature points on the child element.
78  _dphi = nullptr;
79 
80  _cont = fe->get_continuity();
81 
82  if (_cont == C_ONE)
83  {
84  const std::vector<std::vector<GradientShapeType>> & ref_dphi = fe->get_dphi();
85  _dphi = &ref_dphi;
86  }
87 
88  // The Jacobian * quadrature weight at the quadrature points
89  _JxW = &fe->get_JxW();
90  // The XYZ locations of the quadrature points
91  _xyz_values = &fe->get_xyz();
92 
93  // Update the DOF indices for this element based on the current mesh
94  _var.prepareIC();
95  _dof_indices = _var.dofIndices();
96 
97  // The number of DOFs on the element
98  const unsigned int n_dofs = _dof_indices.size();
99  if (n_dofs == 0)
100  return;
101 
102  // Fixed vs. free DoFs on edge/face projections
103  _dof_is_fixed.clear();
104  _dof_is_fixed.resize(n_dofs, false);
105  _free_dof.clear();
106  _free_dof.resize(n_dofs, 0);
107 
108  // Zero the interpolated values
109  _Ue.resize(n_dofs);
110  _Ue.zero();
111 
112  DenseVector<char> mask(n_dofs, true);
113 
114  // In general, we need a series of
115  // projections to ensure a unique and continuous
116  // solution. We start by interpolating nodes, then
117  // hold those fixed and project edges, then
118  // hold those fixed and project faces, then
119  // hold those fixed and project interiors
120 
121  // Interpolate node values first
122  _current_dof = 0;
123 
124  for (_n = 0; _n != n_nodes; ++_n)
125  {
126  _nc = FEInterface::n_dofs_at_node(_fe_type, _current_elem, _n);
127 
128  // for nodes that are in more than one subdomain, only compute the initial
129  // condition once on the lowest numbered block
130  auto curr_node = _current_elem->node_ptr(_n);
131  const auto & block_ids = _sys.mesh().getNodeBlockIds(*curr_node);
132 
133  auto priority_block = *(block_ids.begin());
134  for (auto id : block_ids)
135  if (_var.hasBlocks(id))
136  {
137  priority_block = id;
138  break;
139  }
140 
141  if (!hasBlocks(priority_block) && _var.isNodal())
142  {
143  for (decltype(_nc) i = 0; i < _nc; ++i)
144  {
145  mask(_current_dof) = false;
146  _current_dof++;
147  }
148  continue;
149  }
150 
151  // FIXME: this should go through the DofMap,
152  // not duplicate _dof_indices code badly!
153  if (!_current_elem->is_vertex(_n))
154  {
155  _current_dof += _nc;
156  continue;
157  }
158 
159  if (_cont == DISCONTINUOUS || _cont == H_CURL || _cont == H_DIV)
160  libmesh_assert(_nc == 0);
161  else if (_cont == C_ZERO)
162  setCZeroVertices();
163  else if (_fe_type.family == HERMITE)
164  setHermiteVertices();
165  else if (_cont == C_ONE)
166  setOtherCOneVertices();
167  else if (_cont == SIDE_DISCONTINUOUS)
168  continue;
169  else
170  libmesh_error();
171  } // loop over nodes
172 
173  // From here on out we won't be sampling at nodes anymore
174  _current_node = nullptr;
175 
176  auto & dof_map = _var.dofMap();
177  const bool add_p_level =
178  dof_map.should_p_refine(dof_map.var_group_from_var_number(_var.number()));
179 
180  // In 3D, project any edge values next
181  if (_dim > 2 && _cont != DISCONTINUOUS)
182  for (unsigned int e = 0; e != _current_elem->n_edges(); ++e)
183  {
184  FEInterface::dofs_on_edge(_current_elem, _dim, _fe_type, e, _side_dofs, add_p_level);
185 
186  // Some edge dofs are on nodes 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 edge
198  fe->attach_quadrature_rule(qedgerule.get());
199  fe->edge_reinit(_current_elem, e);
200  _n_qp = qedgerule->n_points();
201 
202  choleskySolve(false);
203  }
204 
205  // Project any side values (edges in 2D, faces in 3D)
206  if (_dim > 1 && _cont != DISCONTINUOUS)
207  for (unsigned int s = 0; s != _current_elem->n_sides(); ++s)
208  {
209  FEInterface::dofs_on_side(_current_elem, _dim, _fe_type, s, _side_dofs, add_p_level);
210 
211  // Some side dofs are on nodes/edges and already
212  // fixed, others are free to calculate
213  _free_dofs = 0;
214  for (unsigned int i = 0; i != _side_dofs.size(); ++i)
215  if (!_dof_is_fixed[_side_dofs[i]])
216  _free_dof[_free_dofs++] = i;
217 
218  // There may be nothing to project
219  if (!_free_dofs)
220  continue;
221 
222  // Initialize FE data on the side
223  fe->attach_quadrature_rule(qsiderule.get());
224  fe->reinit(_current_elem, s);
225  _n_qp = qsiderule->n_points();
226 
227  choleskySolve(false);
228  }
229 
230  // Project the interior values, finally
231 
232  // Some interior dofs are on nodes/edges/sides and
233  // already fixed, others are free to calculate
234  _free_dofs = 0;
235  for (unsigned int i = 0; i != n_dofs; ++i)
236  if (!_dof_is_fixed[i])
237  _free_dof[_free_dofs++] = i;
238 
239  // There may be nothing to project
240  if (_free_dofs)
241  {
242  // Initialize FE data
243  fe->attach_quadrature_rule(qrule.get());
244  fe->reinit(_current_elem);
245  _n_qp = qrule->n_points();
246 
247  choleskySolve(true);
248  } // if there are free interior dofs
249 
250  // Make sure every DoF got reached!
251  for (unsigned int i = 0; i != n_dofs; ++i)
252  libmesh_assert(_dof_is_fixed[i]);
253 
254  // Lock the new_vector since it is shared among threads.
255  {
256  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
257  for (size_t i = 0; i < mask.size(); i++)
258  if (mask(i))
259  _var.setDofValue(_Ue(i), i);
260  }
261 }
262 
263 template <typename T>
264 void
266 {
267  // Assume that C_ZERO elements have a single nodal
268  // value shape function
269  libmesh_assert(_nc == 1);
270  _qp = _n;
271  _current_node = _current_elem->node_ptr(_n);
272  _Ue(_current_dof) = value(*_current_node);
273  _dof_is_fixed[_current_dof] = true;
274  _current_dof++;
275 }
276 
277 template <>
278 void
280 {
281  _qp = _n;
282  _current_node = _current_elem->node_ptr(_n);
283  auto point_value = value(*_current_node);
284  for (decltype(_nc) i = 0; i < _nc; ++i)
285  {
286  _Ue(_current_dof) = point_value(i);
287  _dof_is_fixed[_current_dof] = true;
288  _current_dof++;
289  }
290 }
291 
292 template <typename T>
293 T
295 {
296  return grad(i);
297 }
298 
299 template <>
302 {
303  return grad.row(i);
304 }
305 
306 template <>
309 {
310  return grad.col(i);
311 }
312 
313 template <typename T>
314 void
316 {
317  // The hermite element vertex shape functions are weird
318  _qp = _n;
319  _current_node = _current_elem->node_ptr(_n);
320  _Ue(_current_dof) = value(*_current_node);
321  _dof_is_fixed[_current_dof] = true;
322  _current_dof++;
323  GradientType grad = gradient(*_current_node);
324  // x derivative
325  _Ue(_current_dof) = gradientComponent(grad, 0);
326  _dof_is_fixed[_current_dof] = true;
327  _current_dof++;
328  if (_dim > 1)
329  {
330  // We'll finite difference mixed derivatives
331  Point nxminus = _current_elem->point(_n), nxplus = _current_elem->point(_n);
332  nxminus(0) -= TOLERANCE;
333  nxplus(0) += TOLERANCE;
334  GradientType gxminus = gradient(nxminus);
335  GradientType gxplus = gradient(nxplus);
336  // y derivative
337  _Ue(_current_dof) = gradientComponent(grad, 1);
338  _dof_is_fixed[_current_dof] = true;
339  _current_dof++;
340  // xy derivative
341  _Ue(_current_dof) =
342  (gradientComponent(gxplus, 1) - gradientComponent(gxminus, 1)) / 2. / TOLERANCE;
343  _dof_is_fixed[_current_dof] = true;
344  _current_dof++;
345 
346  if (_dim > 2)
347  {
348  // z derivative
349  _Ue(_current_dof) = gradientComponent(grad, 2);
350  _dof_is_fixed[_current_dof] = true;
351  _current_dof++;
352  // xz derivative
353  _Ue(_current_dof) =
354  (gradientComponent(gxplus, 2) - gradientComponent(gxminus, 2)) / 2. / TOLERANCE;
355  _dof_is_fixed[_current_dof] = true;
356  _current_dof++;
357  // We need new points for yz
358  Point nyminus = _current_elem->point(_n), nyplus = _current_elem->point(_n);
359  nyminus(1) -= TOLERANCE;
360  nyplus(1) += TOLERANCE;
361  GradientType gyminus = gradient(nyminus);
362  GradientType gyplus = gradient(nyplus);
363  // xz derivative
364  _Ue(_current_dof) =
365  (gradientComponent(gyplus, 2) - gradientComponent(gyminus, 2)) / 2. / TOLERANCE;
366  _dof_is_fixed[_current_dof] = true;
367  _current_dof++;
368  // Getting a 2nd order xyz is more tedious
369  Point nxmym = _current_elem->point(_n), nxmyp = _current_elem->point(_n),
370  nxpym = _current_elem->point(_n), nxpyp = _current_elem->point(_n);
371  nxmym(0) -= TOLERANCE;
372  nxmym(1) -= TOLERANCE;
373  nxmyp(0) -= TOLERANCE;
374  nxmyp(1) += TOLERANCE;
375  nxpym(0) += TOLERANCE;
376  nxpym(1) -= TOLERANCE;
377  nxpyp(0) += TOLERANCE;
378  nxpyp(1) += TOLERANCE;
379  GradientType gxmym = gradient(nxmym);
380  GradientType gxmyp = gradient(nxmyp);
381  GradientType gxpym = gradient(nxpym);
382  GradientType gxpyp = gradient(nxpyp);
383  DataType gxzplus =
384  (gradientComponent(gxpyp, 2) - gradientComponent(gxmyp, 2)) / 2. / TOLERANCE;
385  DataType gxzminus =
386  (gradientComponent(gxpym, 2) - gradientComponent(gxmym, 2)) / 2. / TOLERANCE;
387  // xyz derivative
388  _Ue(_current_dof) = (gxzplus - gxzminus) / 2. / TOLERANCE;
389  _dof_is_fixed[_current_dof] = true;
390  _current_dof++;
391  }
392  }
393 }
394 
395 template <>
396 void
398 {
399 }
400 
401 template <typename T>
402 void
404 {
405  // Assume that other C_ONE elements have a single nodal
406  // value shape function and nodal gradient component
407  // shape functions
408  libmesh_assert(_nc == 1 + _dim);
409  _current_node = _current_elem->node_ptr(_n);
410  _Ue(_current_dof) = value(*_current_node);
411  _dof_is_fixed[_current_dof] = true;
412  _current_dof++;
413  GradientType grad = gradient(*_current_node);
414  for (unsigned int i = 0; i != _dim; ++i)
415  {
416  _Ue(_current_dof) = gradientComponent(grad, i);
417  _dof_is_fixed[_current_dof] = true;
418  _current_dof++;
419  }
420 }
421 
422 template <>
423 void
425 {
426 }
427 
428 template <typename T>
429 void
431 {
432  // Loop over the quadrature points
433  for (_qp = 0; _qp < _n_qp; _qp++)
434  {
435  // solution at the quadrature point
436  auto fineval = value((*_xyz_values)[_qp]);
437  // solution grad at the quadrature point
438  GradientType finegrad;
439  if (_cont == C_ONE)
440  finegrad = gradient((*_xyz_values)[_qp]);
441 
442  auto dofs_size = is_volume ? _dof_indices.size() : _side_dofs.size();
443 
444  // Form edge projection matrix
445  for (decltype(dofs_size) geomi = 0, freei = 0; geomi != dofs_size; ++geomi)
446  {
447  auto i = is_volume ? geomi : _side_dofs[geomi];
448 
449  // fixed DoFs aren't test functions
450  if (_dof_is_fixed[i])
451  continue;
452  for (decltype(dofs_size) geomj = 0, freej = 0; geomj != dofs_size; ++geomj)
453  {
454  auto j = is_volume ? geomj : _side_dofs[geomj];
455  if (_dof_is_fixed[j])
456  _Fe(freei) -= (*_phi)[i][_qp] * (*_phi)[j][_qp] * (*_JxW)[_qp] * _Ue(j);
457  else
458  _Ke(freei, freej) += (*_phi)[i][_qp] * (*_phi)[j][_qp] * (*_JxW)[_qp];
459  if (_cont == C_ONE)
460  {
461  if (_dof_is_fixed[j])
462  _Fe(freei) -= dotHelper((*_dphi)[i][_qp], (*_dphi)[j][_qp]) * (*_JxW)[_qp] * _Ue(j);
463  else
464  _Ke(freei, freej) += dotHelper((*_dphi)[i][_qp], (*_dphi)[j][_qp]) * (*_JxW)[_qp];
465  }
466  if (!_dof_is_fixed[j])
467  freej++;
468  }
469  _Fe(freei) += (*_phi)[i][_qp] * fineval * (*_JxW)[_qp];
470  if (_cont == C_ONE)
471  _Fe(freei) += dotHelper(finegrad, (*_dphi)[i][_qp]) * (*_JxW)[_qp];
472  freei++;
473  }
474  }
475 }
476 
477 template <typename T>
478 void
480 {
481  _Ke.resize(_free_dofs, _free_dofs);
482  _Ke.zero();
483  _Fe.resize(_free_dofs);
484  _Fe.zero();
485 
486  choleskyAssembly(is_volume);
487 
488  // The new edge coefficients
489  DenseVector<DataType> U(_free_dofs);
490 
491  _Ke.cholesky_solve(_Fe, U);
492 
493  // Transfer new edge solutions to element
494  for (unsigned int i = 0; i != _free_dofs; ++i)
495  {
496  auto the_dof = is_volume ? _free_dof[i] : _side_dofs[_free_dof[i]];
497  DataType & ui = _Ue(the_dof);
498  libmesh_assert(std::abs(ui) < TOLERANCE || std::abs(ui - U(i)) < TOLERANCE);
499  ui = U(i);
500  _dof_is_fixed[the_dof] = true;
501  }
502 }
503 
504 template <>
505 void
507 {
508  _Ke.resize(_free_dofs, _free_dofs);
509  _Ke.zero();
510  _Fe.resize(_free_dofs);
511  for (unsigned int i = 0; i < _free_dofs; ++i)
512  _Fe(i).setZero(_var.count());
513 
514  choleskyAssembly(is_volume);
515 
516  // The new edge coefficients
517  DenseVector<DataType> U = _Fe;
518 
519  for (unsigned int i = 0; i < _var.count(); ++i)
520  {
521  DenseVector<Real> v(_free_dofs), x(_free_dofs);
522  for (unsigned int j = 0; j < _free_dofs; ++j)
523  v(j) = _Fe(j)(i);
524 
525  _Ke.cholesky_solve(v, x);
526 
527  for (unsigned int j = 0; j < _free_dofs; ++j)
528  U(j)(i) = x(j);
529  }
530 
531  // Transfer new edge solutions to element
532  for (unsigned int i = 0; i != _free_dofs; ++i)
533  {
534  auto the_dof = is_volume ? _free_dof[i] : _side_dofs[_free_dof[i]];
535  DataType & ui = _Ue(the_dof);
536  libmesh_assert(ui.matrix().norm() < TOLERANCE || (ui - U(i)).matrix().norm() < TOLERANCE);
537  ui = U(i);
538  _dof_is_fixed[the_dof] = true;
539  }
540 }
541 
542 template <typename T>
543 void
545 {
546  _var.reinitNode();
547  _var.computeNodalValues(); // has to call this to resize the internal array
548  auto return_value = value(p);
549 
550  _var.setNodalValue(return_value); // update variable data, which is referenced by others, so the
551  // value is up-to-date
552 
553  // We are done, so update the solution vector
554  {
555  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
556  _var.insert(_var.sys().solution());
557  }
558 }
559 
560 template class InitialConditionTempl<Real>;
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: EigenADReal.h:42
OutputTools< T >::OutputData DataType
Class for stuff related to variables.
Definition: Adaptivity.h:31
HERMITE
void setOtherCOneVertices()
set the temporary solution vector for node projections of non-Hermitian C1 variables ...
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:1155
This is a template class that implements the workhorse compute and computeNodal methods.
void choleskyAssembly(bool is_volume)
Assemble a small local system for cholesky solve.
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...
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
DISCONTINUOUS
H_CURL
SIDE_DISCONTINUOUS
const dof_id_type n_nodes
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
void setCZeroVertices()
set the temporary solution vector for node projections of C0 variables
C_ZERO
libmesh_assert(ctx)
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.
H_DIV
void choleskySolve(bool is_volume)
Perform the cholesky solves for edge, side, and interior projections.
C_ONE
virtual unsigned int size() const override final
Eigen::Matrix< Real, Eigen::Dynamic, 1 > RealEigenVector
Definition: MooseTypes.h:146
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
T gradientComponent(GradientType grad, unsigned int i)
unsigned int THREAD_ID
Definition: MooseTypes.h:209