Line data Source code
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 "MortarScalarBase.h"
11 :
12 : // MOOSE includes
13 : #include "Assembly.h"
14 : #include "SystemBase.h"
15 : #include "MooseVariable.h"
16 : #include "MooseVariableScalar.h"
17 :
18 : InputParameters
19 43191 : MortarScalarBase::validParams()
20 : {
21 43191 : InputParameters params = MortarConstraint::validParams();
22 : // This parameter can get renamed in derived class to a more relevant variable name
23 43191 : params.addCoupledVar("scalar_variable", "Primary coupled scalar variable");
24 43191 : params.addParam<bool>("compute_scalar_residuals", true, "Whether to compute scalar residuals");
25 43191 : return params;
26 0 : }
27 :
28 199 : MortarScalarBase::MortarScalarBase(const InputParameters & parameters)
29 : : MortarConstraint(parameters),
30 199 : _use_scalar(isParamValid("scalar_variable") ? true : false),
31 199 : _compute_scalar_residuals(!_use_scalar ? false : getParam<bool>("compute_scalar_residuals")),
32 199 : _kappa_var_ptr(_use_scalar ? getScalarVar("scalar_variable", 0) : nullptr),
33 199 : _kappa_var(_use_scalar ? _kappa_var_ptr->number() : 0),
34 199 : _k_order(_use_scalar ? _kappa_var_ptr->order() : 0),
35 398 : _kappa(_use_scalar ? (_is_implicit ? _kappa_var_ptr->sln() : _kappa_var_ptr->slnOld()) : _zero)
36 : {
37 199 : }
38 :
39 : void
40 6468 : MortarScalarBase::computeResidual()
41 : {
42 6468 : MortarConstraintBase::computeResidual();
43 :
44 6468 : if (!_compute_scalar_residuals)
45 48 : return;
46 :
47 6420 : std::vector<Real> scalar_residuals(_k_order);
48 34976 : for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
49 : {
50 28556 : initScalarQpResidual();
51 107844 : for (_h = 0; _h < _k_order; _h++)
52 79288 : scalar_residuals[_h] += _JxW_msm[_qp] * _coord[_qp] * computeScalarQpResidual();
53 : }
54 :
55 6420 : addResiduals(
56 6420 : _assembly, scalar_residuals, _kappa_var_ptr->dofIndices(), _kappa_var_ptr->scalingFactor());
57 6420 : }
58 :
59 : void
60 3212 : MortarScalarBase::computeJacobian()
61 : {
62 : // d-_var-residual / d-_var and d-_var-residual / d-jvar
63 3212 : MortarConstraintBase::computeJacobian();
64 :
65 3212 : if (!_use_scalar)
66 0 : return;
67 :
68 : // Get the list of coupled scalar vars and compute their off-diag jacobians
69 3212 : const auto & coupled_scalar_vars = getCoupledMooseScalarVars();
70 :
71 : // Handle ALL d-_var-residual / d-scalar columns like computeOffDiagJacobianScalar
72 3212 : if (_compute_primal_residuals)
73 : // Do: dvar / dscalar_var, only want to process only nl-variables (not aux ones)
74 9708 : for (const auto & svariable : coupled_scalar_vars)
75 6496 : if (_sys.hasScalarVariable(svariable->name()))
76 : {
77 : // Compute the jacobian for the secondary interior primal dofs
78 3252 : computeOffDiagJacobianScalar(Moose::MortarType::Secondary, svariable->number());
79 : // Compute the jacobian for the primary interior primal dofs.
80 3252 : computeOffDiagJacobianScalar(Moose::MortarType::Primary, svariable->number());
81 : }
82 :
83 3212 : if (_compute_lm_residuals)
84 : // Do: dvar / dscalar_var, only want to process only nl-variables (not aux ones)
85 6180 : for (const auto & svariable : coupled_scalar_vars)
86 4120 : if (_sys.hasScalarVariable(svariable->name()))
87 : // Compute the jacobian for the lower dimensional LM dofs (if we even have an LM variable)
88 2028 : computeOffDiagJacobianScalar(Moose::MortarType::Lower, svariable->number());
89 :
90 3212 : if (_compute_scalar_residuals)
91 : {
92 : // Handle ALL d-_kappa-residual / d-_var and d-_kappa-residual / d-jvar columns
93 3180 : computeScalarOffDiagJacobian();
94 :
95 : // Do: d-_kappa-residual / d-_kappa and d-_kappa-residual / d-jvar,
96 : // only want to process only nl-variables (not aux ones)
97 9612 : for (const auto & svariable : coupled_scalar_vars)
98 : {
99 6432 : if (_sys.hasScalarVariable(svariable->name()))
100 : {
101 3252 : const unsigned int svar_num = svariable->number();
102 3252 : if (svar_num == _kappa_var)
103 3180 : computeScalarJacobian(); // d-_kappa-residual / d-_kappa
104 : else
105 72 : computeScalarOffDiagJacobianScalar(svar_num); // d-_kappa-residual / d-svar
106 : }
107 : }
108 : }
109 : }
110 :
111 : void
112 3180 : MortarScalarBase::computeScalarJacobian()
113 : {
114 3180 : _local_ke.resize(_k_order, _k_order);
115 19776 : for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
116 : {
117 16596 : initScalarQpJacobian(_kappa_var);
118 64332 : for (_h = 0; _h < _k_order; _h++)
119 187128 : for (_l = 0; _l < _k_order; _l++)
120 139392 : _local_ke(_h, _l) += _JxW_msm[_qp] * _coord[_qp] * computeScalarQpJacobian();
121 : }
122 :
123 12720 : addJacobian(_assembly,
124 3180 : _local_ke,
125 3180 : _kappa_var_ptr->dofIndices(),
126 3180 : _kappa_var_ptr->dofIndices(),
127 3180 : _kappa_var_ptr->scalingFactor());
128 3180 : }
129 :
130 : void
131 3180 : MortarScalarBase::computeScalarOffDiagJacobian()
132 : {
133 : typedef Moose::MortarType MType;
134 3180 : std::array<MType, 3> mortar_types = {{MType::Secondary, MType::Primary, MType::Lower}};
135 :
136 3180 : auto & ce = _assembly.scalarFieldCouplingEntries();
137 12072 : for (const auto & it : ce)
138 : {
139 8892 : MooseVariableScalar & ivariable = *(it.first);
140 8892 : MooseVariableFEBase & jvariable = *(it.second);
141 :
142 8892 : const unsigned int ivar_num = ivariable.number();
143 8892 : const unsigned int jvar_num = jvariable.number();
144 :
145 8892 : if (ivar_num != _kappa_var) // only do the row for _kappa_var in this object
146 72 : continue;
147 :
148 : // Load shape functions of different types for easy access; identical to MortarConstraint.C
149 8820 : std::array<size_t, 3> shape_space_sizes{{jvariable.dofIndices().size(),
150 8820 : jvariable.dofIndicesNeighbor().size(),
151 17640 : jvariable.dofIndicesLower().size()}};
152 : std::array<const VariablePhiValue *, 3> phis;
153 : std::array<const VariablePhiGradient *, 3> grad_phis;
154 : std::array<const VectorVariablePhiValue *, 3> vector_phis;
155 : std::array<const VectorVariablePhiGradient *, 3> vector_grad_phis;
156 8820 : if (jvariable.isVector())
157 : {
158 0 : const auto & temp_var = static_cast<MooseVariableFE<RealVectorValue> &>(jvariable);
159 0 : vector_phis = {{&temp_var.phiFace(), &temp_var.phiFaceNeighbor(), &temp_var.phiLower()}};
160 0 : vector_grad_phis = {
161 0 : {&temp_var.gradPhiFace(), &temp_var.gradPhiFaceNeighbor(), &temp_var.gradPhiLower()}};
162 : }
163 : else
164 : {
165 8820 : const auto & temp_var = static_cast<MooseVariableFE<Real> &>(jvariable);
166 8820 : phis = {{&temp_var.phiFace(), &temp_var.phiFaceNeighbor(), &temp_var.phiLower()}};
167 8820 : grad_phis = {
168 8820 : {&temp_var.gradPhiFace(), &temp_var.gradPhiFaceNeighbor(), &temp_var.gradPhiLower()}};
169 : }
170 :
171 : // Loop over 3 types of spatial variables, find out what jvar_num is
172 35280 : for (MooseIndex(3) type_index = 0; type_index < 3; ++type_index)
173 : {
174 26460 : const auto mortar_type = mortar_types[type_index];
175 26460 : const auto shape_space_size = shape_space_sizes[type_index];
176 26460 : std::vector<dof_id_type> dof_indices;
177 26460 : switch (mortar_type)
178 : {
179 8820 : case MType::Secondary:
180 8820 : dof_indices = jvariable.dofIndices();
181 8820 : break;
182 :
183 8820 : case MType::Primary:
184 8820 : dof_indices = jvariable.dofIndicesNeighbor();
185 8820 : break;
186 :
187 8820 : case MType::Lower:
188 8820 : dof_indices = jvariable.dofIndicesLower();
189 8820 : break;
190 : }
191 :
192 : /// Set the proper phis
193 26460 : if (jvariable.isVector())
194 : {
195 0 : _vector_phi = vector_phis[type_index];
196 0 : _vector_grad_phi = vector_grad_phis[type_index];
197 : }
198 : else
199 : {
200 26460 : _phi = phis[type_index];
201 26460 : _grad_phi = grad_phis[type_index];
202 : }
203 :
204 26460 : _local_ke.resize(_k_order, shape_space_size);
205 :
206 189648 : for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
207 : {
208 163188 : initScalarQpOffDiagJacobian(mortar_type, jvar_num);
209 163188 : const Real dV = _JxW_msm[_qp] * _coord[_qp];
210 642060 : for (_h = 0; _h < _k_order; _h++)
211 : {
212 3120048 : for (_j = 0; _j < shape_space_size; _j++)
213 : {
214 2641176 : _local_ke(_h, _j) += computeScalarQpOffDiagJacobian(mortar_type, jvar_num) * dV;
215 : }
216 : }
217 : }
218 :
219 79380 : addJacobian(_assembly,
220 26460 : _local_ke,
221 26460 : _kappa_var_ptr->dofIndices(),
222 : dof_indices,
223 26460 : _kappa_var_ptr->scalingFactor());
224 26460 : }
225 : }
226 3180 : }
227 :
228 : void
229 8532 : MortarScalarBase::computeOffDiagJacobianScalar(Moose::MortarType mortar_type, unsigned int svar_num)
230 : {
231 8532 : unsigned int test_space_size = 0;
232 8532 : std::vector<dof_id_type> dof_indices;
233 8532 : Real scaling_factor = 1;
234 8532 : switch (mortar_type)
235 : {
236 3252 : case Moose::MortarType::Secondary:
237 3252 : test_space_size = _test_secondary.size();
238 3252 : dof_indices = _secondary_var.dofIndices();
239 3252 : scaling_factor = _secondary_var.scalingFactor();
240 3252 : break;
241 :
242 3252 : case Moose::MortarType::Primary:
243 3252 : test_space_size = _test_primary.size();
244 3252 : dof_indices = _primary_var.dofIndicesNeighbor();
245 3252 : scaling_factor = _primary_var.scalingFactor();
246 3252 : break;
247 :
248 2028 : case Moose::MortarType::Lower:
249 : mooseAssert(_var, "LM variable is null");
250 2028 : test_space_size = _test.size();
251 2028 : dof_indices = _var->dofIndicesLower();
252 2028 : scaling_factor = _var->scalingFactor();
253 2028 : break;
254 : }
255 :
256 : // Get dofs and order of this scalar; at least one will be _kappa_var
257 8532 : const auto & svar = _sys.getScalarVariable(_tid, svar_num);
258 8532 : const unsigned int s_order = svar.order();
259 :
260 8532 : _local_ke.resize(test_space_size, s_order);
261 :
262 55008 : for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
263 : {
264 46476 : initScalarQpOffDiagJacobian(mortar_type, svar_num);
265 46476 : const Real dV = _JxW_msm[_qp] * _coord[_qp];
266 180324 : for (_h = 0; _h < s_order; _h++)
267 2415096 : for (_i = 0; _i < test_space_size; _i++)
268 : { // This assumes Galerkin, i.e. the test and trial functions are the
269 : // same
270 2281248 : _j = _i;
271 2281248 : _local_ke(_i, _h) += computeQpOffDiagJacobianScalar(mortar_type, svar_num) * dV;
272 : }
273 : }
274 :
275 8532 : addJacobian(_assembly, _local_ke, dof_indices, svar.dofIndices(), scaling_factor);
276 8532 : }
277 :
278 : void
279 72 : MortarScalarBase::computeScalarOffDiagJacobianScalar(const unsigned int svar_num)
280 : {
281 : // Get dofs and order of this scalar; will NOT be _kappa_var
282 72 : const auto & svar = _sys.getScalarVariable(_tid, svar_num);
283 72 : const unsigned int s_order = svar.order();
284 :
285 72 : _local_ke.resize(_k_order, s_order);
286 216 : for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
287 : {
288 144 : initScalarQpJacobian(svar_num);
289 288 : for (_h = 0; _h < _k_order; _h++)
290 288 : for (_l = 0; _l < s_order; _l++)
291 144 : _local_ke(_h, _l) +=
292 144 : _JxW_msm[_qp] * _coord[_qp] * computeScalarQpOffDiagJacobianScalar(svar_num);
293 : }
294 :
295 288 : addJacobian(_assembly,
296 72 : _local_ke,
297 72 : _kappa_var_ptr->dofIndices(),
298 72 : svar.dofIndices(),
299 72 : _kappa_var_ptr->scalingFactor());
300 72 : }
|
