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 43167 : MortarScalarBase::validParams()
20 : {
21 43167 : InputParameters params = MortarConstraint::validParams();
22 : // This parameter can get renamed in derived class to a more relevant variable name
23 43167 : params.addCoupledVar("scalar_variable", "Primary coupled scalar variable");
24 43167 : params.addParam<bool>("compute_scalar_residuals", true, "Whether to compute scalar residuals");
25 43167 : return params;
26 0 : }
27 :
28 187 : MortarScalarBase::MortarScalarBase(const InputParameters & parameters)
29 : : MortarConstraint(parameters),
30 187 : _use_scalar(isParamValid("scalar_variable") ? true : false),
31 187 : _compute_scalar_residuals(!_use_scalar ? false : getParam<bool>("compute_scalar_residuals")),
32 187 : _kappa_var_ptr(_use_scalar ? getScalarVar("scalar_variable", 0) : nullptr),
33 187 : _kappa_var(_use_scalar ? _kappa_var_ptr->number() : 0),
34 187 : _k_order(_use_scalar ? _kappa_var_ptr->order() : 0),
35 374 : _kappa(_use_scalar ? (_is_implicit ? _kappa_var_ptr->sln() : _kappa_var_ptr->slnOld()) : _zero)
36 : {
37 187 : }
38 :
39 : void
40 6172 : MortarScalarBase::computeResidual()
41 : {
42 6172 : MortarConstraintBase::computeResidual();
43 :
44 6172 : if (!_compute_scalar_residuals)
45 48 : return;
46 :
47 6124 : std::vector<Real> scalar_residuals(_k_order);
48 34024 : for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
49 : {
50 27900 : initScalarQpResidual();
51 105908 : for (_h = 0; _h < _k_order; _h++)
52 78008 : scalar_residuals[_h] += _JxW_msm[_qp] * _coord[_qp] * computeScalarQpResidual();
53 : }
54 :
55 6124 : addResiduals(
56 6124 : _assembly, scalar_residuals, _kappa_var_ptr->dofIndices(), _kappa_var_ptr->scalingFactor());
57 6124 : }
58 :
59 : void
60 3124 : MortarScalarBase::computeJacobian()
61 : {
62 : // d-_var-residual / d-_var and d-_var-residual / d-jvar
63 3124 : MortarConstraintBase::computeJacobian();
64 :
65 3124 : if (!_use_scalar)
66 0 : return;
67 :
68 : // Get the list of coupled scalar vars and compute their off-diag jacobians
69 3124 : const auto & coupled_scalar_vars = getCoupledMooseScalarVars();
70 :
71 : // Handle ALL d-_var-residual / d-scalar columns like computeOffDiagJacobianScalar
72 3124 : if (_compute_primal_residuals)
73 : // Do: dvar / dscalar_var, only want to process only nl-variables (not aux ones)
74 9436 : for (const auto & svariable : coupled_scalar_vars)
75 6312 : if (_sys.hasScalarVariable(svariable->name()))
76 : {
77 : // Compute the jacobian for the secondary interior primal dofs
78 3156 : computeOffDiagJacobianScalar(Moose::MortarType::Secondary, svariable->number());
79 : // Compute the jacobian for the primary interior primal dofs.
80 3156 : computeOffDiagJacobianScalar(Moose::MortarType::Primary, svariable->number());
81 : }
82 :
83 3124 : if (_compute_lm_residuals)
84 : // Do: dvar / dscalar_var, only want to process only nl-variables (not aux ones)
85 6084 : for (const auto & svariable : coupled_scalar_vars)
86 4056 : if (_sys.hasScalarVariable(svariable->name()))
87 : // Compute the jacobian for the lower dimensional LM dofs (if we even have an LM variable)
88 1996 : computeOffDiagJacobianScalar(Moose::MortarType::Lower, svariable->number());
89 :
90 3124 : if (_compute_scalar_residuals)
91 : {
92 : // Handle ALL d-_kappa-residual / d-_var and d-_kappa-residual / d-jvar columns
93 3092 : 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 9340 : for (const auto & svariable : coupled_scalar_vars)
98 : {
99 6248 : if (_sys.hasScalarVariable(svariable->name()))
100 : {
101 3156 : const unsigned int svar_num = svariable->number();
102 3156 : if (svar_num == _kappa_var)
103 3092 : computeScalarJacobian(); // d-_kappa-residual / d-_kappa
104 : else
105 64 : computeScalarOffDiagJacobianScalar(svar_num); // d-_kappa-residual / d-svar
106 : }
107 : }
108 : }
109 : }
110 :
111 : void
112 3092 : MortarScalarBase::computeScalarJacobian()
113 : {
114 3092 : _local_ke.resize(_k_order, _k_order);
115 19480 : for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
116 : {
117 16388 : initScalarQpJacobian(_kappa_var);
118 63724 : for (_h = 0; _h < _k_order; _h++)
119 185944 : for (_l = 0; _l < _k_order; _l++)
120 138608 : _local_ke(_h, _l) += _JxW_msm[_qp] * _coord[_qp] * computeScalarQpJacobian();
121 : }
122 :
123 12368 : addJacobian(_assembly,
124 3092 : _local_ke,
125 3092 : _kappa_var_ptr->dofIndices(),
126 3092 : _kappa_var_ptr->dofIndices(),
127 3092 : _kappa_var_ptr->scalingFactor());
128 3092 : }
129 :
130 : void
131 3092 : MortarScalarBase::computeScalarOffDiagJacobian()
132 : {
133 : typedef Moose::MortarType MType;
134 3092 : std::array<MType, 3> mortar_types = {{MType::Secondary, MType::Primary, MType::Lower}};
135 :
136 3092 : auto & ce = _assembly.scalarFieldCouplingEntries();
137 11840 : for (const auto & it : ce)
138 : {
139 8748 : MooseVariableScalar & ivariable = *(it.first);
140 8748 : MooseVariableFEBase & jvariable = *(it.second);
141 :
142 8748 : const unsigned int ivar_num = ivariable.number();
143 8748 : const unsigned int jvar_num = jvariable.number();
144 :
145 8748 : if (ivar_num != _kappa_var) // only do the row for _kappa_var in this object
146 64 : continue;
147 :
148 : // Load shape functions of different types for easy access; identical to MortarConstraint.C
149 8684 : std::array<size_t, 3> shape_space_sizes{{jvariable.dofIndices().size(),
150 8684 : jvariable.dofIndicesNeighbor().size(),
151 17368 : 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 8684 : 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 8684 : const auto & temp_var = static_cast<MooseVariableFE<Real> &>(jvariable);
166 8684 : phis = {{&temp_var.phiFace(), &temp_var.phiFaceNeighbor(), &temp_var.phiLower()}};
167 8684 : grad_phis = {
168 8684 : {&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 34736 : for (MooseIndex(3) type_index = 0; type_index < 3; ++type_index)
173 : {
174 26052 : const auto mortar_type = mortar_types[type_index];
175 26052 : const auto shape_space_size = shape_space_sizes[type_index];
176 26052 : std::vector<dof_id_type> dof_indices;
177 26052 : switch (mortar_type)
178 : {
179 8684 : case MType::Secondary:
180 8684 : dof_indices = jvariable.dofIndices();
181 8684 : break;
182 :
183 8684 : case MType::Primary:
184 8684 : dof_indices = jvariable.dofIndicesNeighbor();
185 8684 : break;
186 :
187 8684 : case MType::Lower:
188 8684 : dof_indices = jvariable.dofIndicesLower();
189 8684 : break;
190 : }
191 :
192 : /// Set the proper phis
193 26052 : 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 26052 : _phi = phis[type_index];
201 26052 : _grad_phi = grad_phis[type_index];
202 : }
203 :
204 26052 : _local_ke.resize(_k_order, shape_space_size);
205 :
206 188136 : for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
207 : {
208 162084 : initScalarQpOffDiagJacobian(mortar_type, jvar_num);
209 162084 : const Real dV = _JxW_msm[_qp] * _coord[_qp];
210 638796 : for (_h = 0; _h < _k_order; _h++)
211 : {
212 3112032 : for (_j = 0; _j < shape_space_size; _j++)
213 : {
214 2635320 : _local_ke(_h, _j) += computeScalarQpOffDiagJacobian(mortar_type, jvar_num) * dV;
215 : }
216 : }
217 : }
218 :
219 78156 : addJacobian(_assembly,
220 26052 : _local_ke,
221 26052 : _kappa_var_ptr->dofIndices(),
222 : dof_indices,
223 26052 : _kappa_var_ptr->scalingFactor());
224 26052 : }
225 : }
226 3092 : }
227 :
228 : void
229 8308 : MortarScalarBase::computeOffDiagJacobianScalar(Moose::MortarType mortar_type, unsigned int svar_num)
230 : {
231 8308 : unsigned int test_space_size = 0;
232 8308 : std::vector<dof_id_type> dof_indices;
233 8308 : Real scaling_factor = 1;
234 8308 : switch (mortar_type)
235 : {
236 3156 : case Moose::MortarType::Secondary:
237 3156 : test_space_size = _test_secondary.size();
238 3156 : dof_indices = _secondary_var.dofIndices();
239 3156 : scaling_factor = _secondary_var.scalingFactor();
240 3156 : break;
241 :
242 3156 : case Moose::MortarType::Primary:
243 3156 : test_space_size = _test_primary.size();
244 3156 : dof_indices = _primary_var.dofIndicesNeighbor();
245 3156 : scaling_factor = _primary_var.scalingFactor();
246 3156 : break;
247 :
248 1996 : case Moose::MortarType::Lower:
249 : mooseAssert(_var, "LM variable is null");
250 1996 : test_space_size = _test.size();
251 1996 : dof_indices = _var->dofIndicesLower();
252 1996 : scaling_factor = _var->scalingFactor();
253 1996 : break;
254 : }
255 :
256 : // Get dofs and order of this scalar; at least one will be _kappa_var
257 8308 : const auto & svar = _sys.getScalarVariable(_tid, svar_num);
258 8308 : const unsigned int s_order = svar.order();
259 :
260 8308 : _local_ke.resize(test_space_size, s_order);
261 :
262 54240 : for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
263 : {
264 45932 : initScalarQpOffDiagJacobian(mortar_type, svar_num);
265 45932 : const Real dV = _JxW_msm[_qp] * _coord[_qp];
266 178756 : for (_h = 0; _h < s_order; _h++)
267 2408440 : for (_i = 0; _i < test_space_size; _i++)
268 : { // This assumes Galerkin, i.e. the test and trial functions are the
269 : // same
270 2275616 : _j = _i;
271 2275616 : _local_ke(_i, _h) += computeQpOffDiagJacobianScalar(mortar_type, svar_num) * dV;
272 : }
273 : }
274 :
275 8308 : addJacobian(_assembly, _local_ke, dof_indices, svar.dofIndices(), scaling_factor);
276 8308 : }
277 :
278 : void
279 64 : MortarScalarBase::computeScalarOffDiagJacobianScalar(const unsigned int svar_num)
280 : {
281 : // Get dofs and order of this scalar; will NOT be _kappa_var
282 64 : const auto & svar = _sys.getScalarVariable(_tid, svar_num);
283 64 : const unsigned int s_order = svar.order();
284 :
285 64 : _local_ke.resize(_k_order, s_order);
286 192 : for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
287 : {
288 128 : initScalarQpJacobian(svar_num);
289 256 : for (_h = 0; _h < _k_order; _h++)
290 256 : for (_l = 0; _l < s_order; _l++)
291 128 : _local_ke(_h, _l) +=
292 128 : _JxW_msm[_qp] * _coord[_qp] * computeScalarQpOffDiagJacobianScalar(svar_num);
293 : }
294 :
295 256 : addJacobian(_assembly,
296 64 : _local_ke,
297 64 : _kappa_var_ptr->dofIndices(),
298 64 : svar.dofIndices(),
299 64 : _kappa_var_ptr->scalingFactor());
300 64 : }
|
