16 #include "libmesh/fe_interface.h" 17 #include "libmesh/quadrature.h" 24 _fe_problem(*getCheckedPointerParam<
FEProblemBase *>(
"_fe_problem_base")),
26 _t(_fe_problem.time()),
27 _var(_sys.getActualFieldVariable<T>(parameters.
get<
THREAD_ID>(
"_tid"),
28 parameters.
get<VariableName>(
"variable"))),
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),
37 _fe_type(_var.feType())
57 _dim = _current_elem->dim();
59 const unsigned int n_nodes = _current_elem->n_nodes();
67 std::unique_ptr<FEBaseType> fe(FEBaseType::build(_dim, _fe_type));
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));
75 _phi = &fe->get_phi();
80 _cont = fe->get_continuity();
84 const std::vector<std::vector<GradientShapeType>> & ref_dphi = fe->get_dphi();
89 _JxW = &fe->get_JxW();
91 _xyz_values = &fe->get_xyz();
95 _dof_indices = _var.dofIndices();
98 const unsigned int n_dofs = _dof_indices.size();
103 _dof_is_fixed.clear();
104 _dof_is_fixed.resize(n_dofs,
false);
106 _free_dof.resize(n_dofs, 0);
124 for (_n = 0; _n !=
n_nodes; ++_n)
126 _nc = FEInterface::n_dofs_at_node(_fe_type, _current_elem, _n);
130 auto curr_node = _current_elem->node_ptr(_n);
131 const auto & block_ids = _sys.mesh().getNodeBlockIds(*curr_node);
133 auto priority_block = *(block_ids.begin());
134 for (
auto id : block_ids)
135 if (_var.hasBlocks(
id))
141 if (!hasBlocks(priority_block) && _var.isNodal())
143 for (decltype(_nc) i = 0; i < _nc; ++i)
145 mask(_current_dof) =
false;
153 if (!_current_elem->is_vertex(_n))
163 else if (_fe_type.family ==
HERMITE)
164 setHermiteVertices();
165 else if (_cont ==
C_ONE)
166 setOtherCOneVertices();
174 _current_node =
nullptr;
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()));
182 for (
unsigned int e = 0; e != _current_elem->n_edges(); ++e)
184 FEInterface::dofs_on_edge(_current_elem, _dim, _fe_type, e, _side_dofs, add_p_level);
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;
198 fe->attach_quadrature_rule(qedgerule.get());
199 fe->edge_reinit(_current_elem, e);
200 _n_qp = qedgerule->n_points();
202 choleskySolve(
false);
207 for (
unsigned int s = 0; s != _current_elem->n_sides(); ++s)
209 FEInterface::dofs_on_side(_current_elem, _dim, _fe_type, s, _side_dofs, add_p_level);
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;
223 fe->attach_quadrature_rule(qsiderule.get());
224 fe->reinit(_current_elem, s);
225 _n_qp = qsiderule->n_points();
227 choleskySolve(
false);
235 for (
unsigned int i = 0; i != n_dofs; ++i)
236 if (!_dof_is_fixed[i])
237 _free_dof[_free_dofs++] = i;
243 fe->attach_quadrature_rule(qrule.get());
244 fe->reinit(_current_elem);
245 _n_qp = qrule->n_points();
251 for (
unsigned int i = 0; i != n_dofs; ++i)
256 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
257 for (
size_t i = 0; i < mask.
size(); i++)
259 _var.setDofValue(_Ue(i), i);
263 template <
typename T>
271 _current_node = _current_elem->node_ptr(_n);
272 _Ue(_current_dof) =
value(*_current_node);
273 _dof_is_fixed[_current_dof] =
true;
282 _current_node = _current_elem->node_ptr(_n);
283 auto point_value =
value(*_current_node);
284 for (decltype(_nc) i = 0; i < _nc; ++i)
286 _Ue(_current_dof) = point_value(i);
287 _dof_is_fixed[_current_dof] =
true;
292 template <
typename T>
313 template <
typename T>
319 _current_node = _current_elem->node_ptr(_n);
320 _Ue(_current_dof) =
value(*_current_node);
321 _dof_is_fixed[_current_dof] =
true;
325 _Ue(_current_dof) = gradientComponent(grad, 0);
326 _dof_is_fixed[_current_dof] =
true;
331 Point nxminus = _current_elem->point(_n), nxplus = _current_elem->point(_n);
332 nxminus(0) -= TOLERANCE;
333 nxplus(0) += TOLERANCE;
337 _Ue(_current_dof) = gradientComponent(grad, 1);
338 _dof_is_fixed[_current_dof] =
true;
342 (gradientComponent(gxplus, 1) - gradientComponent(gxminus, 1)) / 2. / TOLERANCE;
343 _dof_is_fixed[_current_dof] =
true;
349 _Ue(_current_dof) = gradientComponent(grad, 2);
350 _dof_is_fixed[_current_dof] =
true;
354 (gradientComponent(gxplus, 2) - gradientComponent(gxminus, 2)) / 2. / TOLERANCE;
355 _dof_is_fixed[_current_dof] =
true;
358 Point nyminus = _current_elem->point(_n), nyplus = _current_elem->point(_n);
359 nyminus(1) -= TOLERANCE;
360 nyplus(1) += TOLERANCE;
365 (gradientComponent(gyplus, 2) - gradientComponent(gyminus, 2)) / 2. / TOLERANCE;
366 _dof_is_fixed[_current_dof] =
true;
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;
384 (gradientComponent(gxpyp, 2) - gradientComponent(gxmyp, 2)) / 2. / TOLERANCE;
386 (gradientComponent(gxpym, 2) - gradientComponent(gxmym, 2)) / 2. / TOLERANCE;
388 _Ue(_current_dof) = (gxzplus - gxzminus) / 2. / TOLERANCE;
389 _dof_is_fixed[_current_dof] =
true;
401 template <
typename T>
409 _current_node = _current_elem->node_ptr(_n);
410 _Ue(_current_dof) =
value(*_current_node);
411 _dof_is_fixed[_current_dof] =
true;
414 for (
unsigned int i = 0; i != _dim; ++i)
416 _Ue(_current_dof) = gradientComponent(grad, i);
417 _dof_is_fixed[_current_dof] =
true;
428 template <
typename T>
433 for (_qp = 0; _qp < _n_qp; _qp++)
436 auto fineval =
value((*_xyz_values)[_qp]);
440 finegrad = gradient((*_xyz_values)[_qp]);
442 auto dofs_size = is_volume ? _dof_indices.size() : _side_dofs.size();
445 for (decltype(dofs_size) geomi = 0, freei = 0; geomi != dofs_size; ++geomi)
447 auto i = is_volume ? geomi : _side_dofs[geomi];
450 if (_dof_is_fixed[i])
452 for (decltype(dofs_size) geomj = 0, freej = 0; geomj != dofs_size; ++geomj)
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);
458 _Ke(freei, freej) += (*_phi)[i][_qp] * (*_phi)[j][_qp] * (*_JxW)[_qp];
461 if (_dof_is_fixed[j])
462 _Fe(freei) -= dotHelper((*_dphi)[i][_qp], (*_dphi)[j][_qp]) * (*_JxW)[_qp] * _Ue(j);
464 _Ke(freei, freej) += dotHelper((*_dphi)[i][_qp], (*_dphi)[j][_qp]) * (*_JxW)[_qp];
466 if (!_dof_is_fixed[j])
469 _Fe(freei) += (*_phi)[i][_qp] * fineval * (*_JxW)[_qp];
471 _Fe(freei) += dotHelper(finegrad, (*_dphi)[i][_qp]) * (*_JxW)[_qp];
477 template <
typename T>
481 _Ke.resize(_free_dofs, _free_dofs);
483 _Fe.resize(_free_dofs);
486 choleskyAssembly(is_volume);
491 _Ke.cholesky_solve(_Fe, U);
494 for (
unsigned int i = 0; i != _free_dofs; ++i)
496 auto the_dof = is_volume ? _free_dof[i] : _side_dofs[_free_dof[i]];
500 _dof_is_fixed[the_dof] =
true;
508 _Ke.resize(_free_dofs, _free_dofs);
510 _Fe.resize(_free_dofs);
511 for (
unsigned int i = 0; i < _free_dofs; ++i)
512 _Fe(i).setZero(_var.count());
514 choleskyAssembly(is_volume);
519 for (
unsigned int i = 0; i < _var.count(); ++i)
522 for (
unsigned int j = 0; j < _free_dofs; ++j)
525 _Ke.cholesky_solve(v, x);
527 for (
unsigned int j = 0; j < _free_dofs; ++j)
532 for (
unsigned int i = 0; i != _free_dofs; ++i)
534 auto the_dof = is_volume ? _free_dof[i] : _side_dofs[_free_dof[i]];
536 libmesh_assert(ui.matrix().norm() < TOLERANCE || (ui - U(i)).matrix().norm() < TOLERANCE);
538 _dof_is_fixed[the_dof] =
true;
542 template <
typename T>
547 _var.computeNodalValues();
548 auto return_value =
value(p);
550 _var.setNodalValue(return_value);
555 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
556 _var.insert(_var.sys().solution());
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
OutputTools< T >::OutputData DataType
Class for stuff related to variables.
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...
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...
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.
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
void setHermiteVertices()
set the temporary solution vector for node projections of Hermite variables
InitialConditionTempl(const InputParameters ¶meters)
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.
virtual unsigned int size() const override final
Eigen::Matrix< Real, Eigen::Dynamic, 1 > RealEigenVector
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)
virtual ~InitialConditionTempl()