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