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 "MortarConstraintBase.h"
11 : #include "FEProblemBase.h"
12 : #include "Assembly.h"
13 : #include "MooseVariableFE.h"
14 :
15 : #include "libmesh/string_to_enum.h"
16 :
17 : InputParameters
18 48342 : MortarConstraintBase::validParams()
19 : {
20 48342 : InputParameters params = Constraint::validParams();
21 48342 : params += MortarConsumerInterface::validParams();
22 48342 : params += TwoMaterialPropertyInterface::validParams();
23 :
24 : // Whether on a displaced or undisplaced mesh, coupling ghosting will only happen for
25 : // cross-interface elements
26 145026 : params.addRelationshipManager("AugmentSparsityOnInterface",
27 : Moose::RelationshipManagerType::COUPLING,
28 0 : [](const InputParameters & obj_params, InputParameters & rm_params)
29 : {
30 3666 : rm_params.set<bool>("use_displaced_mesh") =
31 3666 : obj_params.get<bool>("use_displaced_mesh");
32 7332 : rm_params.set<BoundaryName>("secondary_boundary") =
33 7332 : obj_params.get<BoundaryName>("secondary_boundary");
34 7332 : rm_params.set<BoundaryName>("primary_boundary") =
35 7332 : obj_params.get<BoundaryName>("primary_boundary");
36 7332 : rm_params.set<SubdomainName>("secondary_subdomain") =
37 7332 : obj_params.get<SubdomainName>("secondary_subdomain");
38 7332 : rm_params.set<SubdomainName>("primary_subdomain") =
39 7332 : obj_params.get<SubdomainName>("primary_subdomain");
40 3666 : rm_params.set<bool>("ghost_higher_d_neighbors") =
41 3666 : obj_params.get<bool>("ghost_higher_d_neighbors");
42 3666 : });
43 :
44 : // If the LM is ever a Lagrange variable, we will attempt to obtain its dof indices from the
45 : // process that owns the node. However, in order for the dof indices to get set for the Lagrange
46 : // variable, the process that owns the node needs to have local copies of any lower-d elements
47 : // that have the connected node. Note that the geometric ghosting done here is different than that
48 : // done by the AugmentSparsityOnInterface RM, even when ghost_point_neighbors is true. The latter
49 : // ghosts equal-manifold lower-dimensional secondary element point neighbors and their interface
50 : // couplings. This ghosts lower-dimensional point neighbors of higher-dimensional elements.
51 : // Neither is guaranteed to be a superset of the other. For instance ghosting of lower-d point
52 : // neighbors (AugmentSparsityOnInterface with ghost_point_neighbors = true) is only guaranteed to
53 : // ghost those lower-d point neighbors on *processes that own lower-d elements*. And you may have
54 : // a process that only owns higher-dimensional elements
55 : //
56 : // Note that in my experience it is only important for the higher-d lower-d point neighbors to be
57 : // ghosted when forming sparsity patterns and so I'm putting this here instead of at the
58 : // MortarConsumerInterface level
59 145026 : params.addRelationshipManager("GhostHigherDLowerDPointNeighbors",
60 : Moose::RelationshipManagerType::GEOMETRIC);
61 :
62 193368 : params.addParam<VariableName>("secondary_variable", "Primal variable on secondary surface.");
63 193368 : params.addParam<VariableName>(
64 : "primary_variable",
65 : "Primal variable on primary surface. If this parameter is not provided then the primary "
66 : "variable will be initialized to the secondary variable");
67 96684 : params.makeParamNotRequired<NonlinearVariableName>("variable");
68 193368 : params.setDocString(
69 : "variable",
70 : "The name of the lagrange multiplier variable that this constraint is applied to. This "
71 : "parameter may not be supplied in the case of using penalty methods for example");
72 145026 : params.addParam<bool>(
73 96684 : "compute_primal_residuals", true, "Whether to compute residuals for the primal variable.");
74 145026 : params.addParam<bool>(
75 96684 : "compute_lm_residuals", true, "Whether to compute Lagrange Multiplier residuals");
76 241710 : params.addDeprecatedParam<MooseEnum>(
77 : "quadrature",
78 241710 : MooseEnum("DEFAULT FIRST SECOND THIRD FOURTH", "DEFAULT"),
79 : "Polynomial basis order (think Variable order) to assume when building quadrature rules to "
80 : "use on mortar segments. "
81 : "For 2D mortar DEFAULT is recommended. "
82 : "For 3D mortar, QUAD meshes are integrated using triangular mortar segments. "
83 : "While DEFAULT order is typically sufficiently accurate, exact integration of "
84 : "QUAD mortar faces with first order polynomial bases (bilinears) requires building a "
85 : "quadrature rule based off of a *quadratic* (SECOND) polynomial basis on triangles. "
86 : "Similarly, exact integration of QUAD mortar faces with second order polynomial bases "
87 : "(biquadratics) requires building a quadrature rule based off of a *quartic* (FOURTH) "
88 : "polynomial basis on triangles. Note that the actual quadrature order will be double this "
89 : "parameter plus one.",
90 : "This parameter is deprecated in favor of "
91 : "'segment_quadrature' which directly specifies the quadrature order.");
92 145026 : params.addParam<MooseEnum>(
93 : "segment_quadrature",
94 241710 : MooseEnum("DEFAULT FIRST SECOND THIRD FOURTH FIFTH SIXTH SEVENTH EIGHTH NINTH", "DEFAULT"),
95 : "Mortar segment quadrature order. "
96 : "For 2D mortar DEFAULT is recommended. "
97 : "For 3D mortar, quad faces are integrated using triangular mortar segments. "
98 : "A finite element family of order p based on a quadrilateral element actually has polynomial "
99 : "order of 2*p because of the tensor-product nature of the element. Consequently, for exact "
100 : "integraton of something like a mass matrix term, if one is using a first order Lagrange "
101 : "variable (as an example), the 'segment_quadrature' should be set to 'fourth' because we "
102 : "double the polynomial order on the triangle to match the tensor product order on the quad, "
103 : "and then double again since we are multiplying the test and trial (shape) function "
104 : "polynomials for the mass matrix term.");
105 145026 : params.addParam<bool>(
106 : "use_petrov_galerkin",
107 96684 : false,
108 : "Whether to use the Petrov-Galerkin approach for the mortar-based constraints. If set to "
109 : "true, we use the standard basis as the test function and dual basis as "
110 : "the shape function for the interpolation of the Lagrange multiplier variable.");
111 145026 : params.addCoupledVar("aux_lm",
112 : "Auxiliary Lagrange multiplier variable that is utilized together with the "
113 : "Petrov-Galerkin approach.");
114 48342 : return params;
115 0 : }
116 :
117 1224 : MortarConstraintBase::MortarConstraintBase(const InputParameters & parameters)
118 : : Constraint(parameters),
119 : NeighborCoupleableMooseVariableDependencyIntermediateInterface(this, false, false),
120 : MortarConsumerInterface(this),
121 : TwoMaterialPropertyInterface(this, Moose::EMPTY_BLOCK_IDS, getBoundaryIDs()),
122 : MooseVariableInterface<Real>(this,
123 : true,
124 4896 : isParamValid("variable") ? "variable" : "secondary_variable",
125 : Moose::VarKindType::VAR_SOLVER,
126 : Moose::VarFieldType::VAR_FIELD_STANDARD),
127 3672 : _fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
128 2448 : _var(isParamValid("variable")
129 2550 : ? &_subproblem.getStandardVariable(_tid, parameters.getMooseType("variable"))
130 : : nullptr),
131 1224 : _secondary_var(
132 2448 : isParamValid("secondary_variable")
133 4896 : ? _sys.getActualFieldVariable<Real>(_tid, parameters.getMooseType("secondary_variable"))
134 1224 : : _sys.getActualFieldVariable<Real>(_tid, parameters.getMooseType("primary_variable"))),
135 1224 : _primary_var(
136 2448 : isParamValid("primary_variable")
137 1224 : ? _sys.getActualFieldVariable<Real>(_tid, parameters.getMooseType("primary_variable"))
138 : : _secondary_var),
139 :
140 2448 : _compute_primal_residuals(getParam<bool>("compute_primal_residuals")),
141 1887 : _compute_lm_residuals(!_var ? false : getParam<bool>("compute_lm_residuals")),
142 1224 : _test_dummy(),
143 1224 : _use_dual(_var ? _var->useDual() : false),
144 1224 : _tangents(_assembly.tangents()),
145 1224 : _coord(_assembly.mortarCoordTransformation()),
146 1224 : _q_point(_assembly.qPointsMortar()),
147 2448 : _use_petrov_galerkin(getParam<bool>("use_petrov_galerkin")),
148 2472 : _aux_lm_var(isCoupled("aux_lm") ? getVar("aux_lm", 0) : nullptr),
149 2448 : _test(_var
150 1224 : ? ((_use_petrov_galerkin && _aux_lm_var) ? _aux_lm_var->phiLower() : _var->phiLower())
151 : : _test_dummy),
152 1224 : _test_secondary(_secondary_var.phiFace()),
153 1224 : _test_primary(_primary_var.phiFaceNeighbor()),
154 1224 : _grad_test_secondary(_secondary_var.gradPhiFace()),
155 1224 : _grad_test_primary(_primary_var.gradPhiFaceNeighbor()),
156 1224 : _interior_secondary_elem(_assembly.elem()),
157 1224 : _interior_primary_elem(_assembly.neighbor()),
158 6120 : _displaced(getParam<bool>("use_displaced_mesh"))
159 : {
160 1224 : if (_use_dual)
161 124 : _assembly.activateDual();
162 :
163 1224 : if (_use_petrov_galerkin && (!_use_dual))
164 0 : paramError("use_petrov_galerkin",
165 : "We need to set `use_dual = true` while using the Petrov-Galerkin approach");
166 :
167 1248 : if (_use_petrov_galerkin && ((!isParamValid("aux_lm")) || _aux_lm_var == nullptr))
168 0 : paramError("use_petrov_galerkin",
169 : "We need to specify an auxiliary variable `aux_lm` while using the Petrov-Galerkin "
170 : "approach");
171 :
172 1224 : if (_use_petrov_galerkin && _aux_lm_var->useDual())
173 0 : paramError("aux_lm",
174 : "Auxiliary LM variable needs to use standard shape function, i.e., set `use_dual = "
175 : "false`.");
176 3672 : if (isParamSetByUser("quadrature") && isParamSetByUser("segment_quadrature"))
177 0 : paramError("quadrature", "Only one of 'quadrature' and 'segment_quadrature' should be set.");
178 :
179 : // Note parameter is discretization order, we then convert to quadrature order
180 2448 : const auto & p_order = getParam<MooseEnum>("quadrature");
181 : // If quadrature not DEFAULT, set mortar qrule
182 1224 : if (p_order != "DEFAULT")
183 : {
184 0 : const Order q_order = static_cast<Order>(2 * Utility::string_to_enum<Order>(p_order) + 1);
185 0 : _assembly.setMortarQRule(q_order);
186 : }
187 2448 : const auto & q_order_enum = getParam<MooseEnum>("segment_quadrature");
188 1224 : if (q_order_enum != "DEFAULT")
189 : {
190 0 : const Order q_order = Utility::string_to_enum<Order>(q_order_enum);
191 0 : _assembly.setMortarQRule(q_order);
192 : }
193 :
194 1224 : if (_var)
195 663 : addMooseVariableDependency(_var);
196 1224 : addMooseVariableDependency(&_secondary_var);
197 1224 : addMooseVariableDependency(&_primary_var);
198 1224 : }
199 :
200 : void
201 181890 : MortarConstraintBase::computeResidual()
202 : {
203 181890 : precalculateResidual();
204 :
205 181890 : if (_compute_primal_residuals)
206 : {
207 : // Compute the residual for the secondary interior primal dofs
208 181890 : computeResidual(Moose::MortarType::Secondary);
209 :
210 : // Compute the residual for the primary interior primal dofs.
211 181890 : computeResidual(Moose::MortarType::Primary);
212 : }
213 :
214 181890 : if (_compute_lm_residuals)
215 : // Compute the residual for the lower dimensional LM dofs (if we even have an LM variable)
216 86258 : computeResidual(Moose::MortarType::Lower);
217 181890 : }
218 :
219 : void
220 88234 : MortarConstraintBase::computeJacobian()
221 : {
222 88234 : precalculateResidual();
223 :
224 88234 : if (_compute_primal_residuals)
225 : {
226 : // Compute the jacobian for the secondary interior primal dofs
227 88234 : computeJacobian(Moose::MortarType::Secondary);
228 :
229 : // Compute the jacobian for the primary interior primal dofs.
230 88234 : computeJacobian(Moose::MortarType::Primary);
231 : }
232 :
233 88234 : if (_compute_lm_residuals)
234 : // Compute the jacobian for the lower dimensional LM dofs (if we even have an LM variable)
235 47610 : computeJacobian(Moose::MortarType::Lower);
236 88234 : }
237 :
238 : void
239 15070 : MortarConstraintBase::zeroInactiveLMDofs(const std::unordered_set<const Node *> & inactive_lm_nodes,
240 : const std::unordered_set<const Elem *> & inactive_lm_elems)
241 : {
242 : // If no LM variable has been defined, skip
243 15070 : if (!_var)
244 6274 : return;
245 :
246 8796 : const auto sn = _sys.number();
247 8796 : const auto vn = _var->number();
248 :
249 : // If variable is nodal, zero DoFs based on inactive LM nodes
250 8796 : if (_var->isNodal())
251 : {
252 6044 : for (const auto node : inactive_lm_nodes)
253 : {
254 : // Allow mixed Lagrange orders between primal and LM
255 0 : if (!node->n_comp(sn, vn))
256 0 : continue;
257 :
258 0 : const auto dof_index = node->dof_number(sn, vn, 0);
259 : // No scaling; this is not physics
260 0 : if (_assembly.computingJacobian())
261 0 : addJacobianElement(
262 : _assembly, /*element_value=*/1, dof_index, dof_index, /*scaling_factor=*/1);
263 0 : if (_assembly.computingResidual())
264 : {
265 0 : const Real lm_value = _var->getNodalValue(*node);
266 0 : addResiduals(_assembly,
267 0 : std::array<Real, 1>{{lm_value}},
268 0 : std::array<dof_id_type, 1>{{dof_index}},
269 : /*scaling_factor=*/1);
270 : }
271 : }
272 : }
273 : // If variable is elemental, zero based on inactive LM elems
274 : else
275 : {
276 3121 : for (const auto el : inactive_lm_elems)
277 : {
278 369 : const auto n_comp = el->n_comp(sn, vn);
279 :
280 738 : for (const auto comp : make_range(n_comp))
281 : {
282 369 : const auto dof_index = el->dof_number(sn, vn, comp);
283 : // No scaling; this is not physics
284 369 : if (_assembly.computingJacobian())
285 123 : addJacobianElement(
286 : _assembly, /*element_value=*/1, dof_index, dof_index, /*scaling_factor=*/1);
287 369 : if (_assembly.computingResidual())
288 : {
289 246 : const Real lm_value = _var->getElementalValue(el, comp);
290 246 : addResiduals(_assembly,
291 0 : std::array<Real, 1>{{lm_value}},
292 246 : std::array<dof_id_type, 1>{{dof_index}},
293 : /*scaling_factor=*/1);
294 : }
295 : }
296 : }
297 : }
298 : }
|
