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 "PetscDMMoose.h"
11 :
12 : // PETSc includes
13 : #include <petscerror.h>
14 : #include <petsc/private/dmimpl.h>
15 :
16 : // MOOSE includes
17 : #include "FEProblem.h"
18 : #include "DisplacedProblem.h"
19 : #include "MooseMesh.h"
20 : #include "NonlinearSystem.h"
21 : #include "PenetrationLocator.h"
22 : #include "NearestNodeLocator.h"
23 : #include "GeometricSearchData.h"
24 : #include "MooseVariableScalar.h"
25 :
26 : #include "libmesh/nonlinear_implicit_system.h"
27 : #include "libmesh/nonlinear_solver.h"
28 : #include "libmesh/petsc_macro.h"
29 : #include "libmesh/petsc_vector.h"
30 : #include "libmesh/petsc_matrix.h"
31 : #include "libmesh/dof_map.h"
32 : #include "libmesh/preconditioner.h"
33 : #include "libmesh/elem_side_builder.h"
34 :
35 : using namespace libMesh;
36 :
37 : template <typename I1, typename I2>
38 : void
39 157 : checkSize(const std::string & split_name, const I1 split_size, const I2 size_expected_by_parent)
40 : {
41 314 : if (libMesh::cast_int<libMesh::numeric_index_type>(split_size) !=
42 157 : libMesh::cast_int<libMesh::numeric_index_type>(size_expected_by_parent))
43 8 : mooseError("Split '",
44 : split_name,
45 : "' has size ",
46 8 : libMesh::cast_int<libMesh::numeric_index_type>(split_size),
47 : " but the parent split expected size ",
48 8 : libMesh::cast_int<libMesh::numeric_index_type>(size_expected_by_parent),
49 : ". Make sure that you have non-overlapping complete sets for variables and "
50 : "blocks as well as consistency in sides/unsides, contacts/uncontacts, etc.");
51 149 : }
52 :
53 : struct DM_Moose
54 : {
55 : NonlinearSystemBase * _nl; // nonlinear system context
56 : const DofMapBase * _dof_map;
57 : const System * _system;
58 : DM_Moose * _parent = nullptr;
59 : std::set<std::string> * _vars; // variables
60 : std::map<std::string, unsigned int> * _var_ids;
61 : std::map<unsigned int, std::string> * _var_names;
62 : bool _all_vars; // whether all system variables are included
63 : std::set<std::string> * _blocks; // mesh blocks
64 : std::map<std::string, subdomain_id_type> * _block_ids;
65 : std::map<unsigned int, std::string> * _block_names;
66 : bool _all_blocks; // all blocks are included
67 : std::set<std::string> * _sides; // mesh surfaces (edges in 2D)
68 : std::map<BoundaryID, std::string> * _side_names;
69 : std::map<std::string, BoundaryID> * _side_ids;
70 : std::set<std::string> * _unsides; // excluded sides
71 : std::map<std::string, BoundaryID> * _unside_ids;
72 : std::map<BoundaryID, std::string> * _unside_names;
73 : std::set<std::string> * _unside_by_var; // excluded sides by variable
74 : std::set<std::pair<BoundaryID, unsigned int>> * _unside_by_var_set;
75 : bool _nosides; // whether to include any sides
76 : bool _nounsides; // whether to exclude any sides
77 : bool _nounside_by_var;
78 : typedef std::pair<std::string, std::string> ContactName;
79 : typedef std::pair<BoundaryID, BoundaryID> ContactID;
80 : std::set<ContactName> * _contacts;
81 : std::map<ContactID, ContactName> * _contact_names;
82 : std::set<ContactName> * _uncontacts;
83 : std::map<ContactID, ContactName> * _uncontact_names;
84 : std::map<ContactName, PetscBool> * _contact_displaced;
85 : std::map<ContactName, PetscBool> * _uncontact_displaced;
86 : bool _nocontacts;
87 : bool _nouncontacts;
88 : bool _include_all_contact_nodes;
89 : // to locate splits without having to search, however,
90 : // maintain a multimap from names to split locations (to enable
91 : // the same split to appear in multiple spots (this might
92 : // break the current implementation of PCFieldSplit, though).
93 : std::multimap<std::string, unsigned int> * _splitlocs;
94 : struct SplitInfo
95 : {
96 : DM _dm;
97 : IS _rembedding; // relative embedding
98 : };
99 : std::map<std::string, SplitInfo> * _splits;
100 :
101 : IS _embedding;
102 : PetscBool _print_embedding;
103 :
104 : /// The name of this DM
105 : std::string * _name;
106 :
107 : /**
108 : * Check whether the size of the child matches the size we expect
109 : */
110 : void checkChildSize(DM child, PetscInt child_size, const std::string & child_name);
111 : };
112 :
113 : void
114 28 : DM_Moose::checkChildSize(DM child, const PetscInt child_size, const std::string & child_name)
115 : {
116 40 : for (const auto & split : *_splits)
117 40 : if (split.second._dm == child)
118 : {
119 : mooseAssert(split.first == child_name, "These should match");
120 : PetscInt parent_expected_size;
121 28 : auto ierr = ISGetLocalSize(split.second._rembedding, &parent_expected_size);
122 28 : if (ierr)
123 0 : mooseError("Unable to get size");
124 28 : checkSize(child_name, child_size, parent_expected_size);
125 24 : return;
126 : }
127 :
128 0 : mooseError("No child DM match");
129 : }
130 :
131 : PetscErrorCode
132 4896 : DMMooseValidityCheck(DM dm)
133 : {
134 : PetscBool ismoose;
135 :
136 : PetscFunctionBegin;
137 : PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
138 4896 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose));
139 4896 : if (!ismoose)
140 0 : LIBMESH_SETERRQ2(((PetscObject)dm)->comm,
141 : PETSC_ERR_ARG_WRONG,
142 : "Got DM of type %s, not of type %s",
143 : ((PetscObject)dm)->type_name,
144 : DMMOOSE);
145 4896 : PetscFunctionReturn(PETSC_SUCCESS);
146 : }
147 :
148 : PetscErrorCode
149 0 : DMMooseGetContacts(DM dm,
150 : std::vector<std::pair<std::string, std::string>> & contact_names,
151 : std::vector<PetscBool> & displaced)
152 : {
153 : PetscFunctionBegin;
154 0 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
155 0 : DM_Moose * dmm = (DM_Moose *)dm->data;
156 0 : for (const auto & it : *(dmm->_contact_names))
157 : {
158 0 : contact_names.push_back(it.second);
159 0 : displaced.push_back((*dmm->_contact_displaced)[it.second]);
160 : }
161 0 : PetscFunctionReturn(PETSC_SUCCESS);
162 : }
163 :
164 : PetscErrorCode
165 0 : DMMooseGetUnContacts(DM dm,
166 : std::vector<std::pair<std::string, std::string>> & uncontact_names,
167 : std::vector<PetscBool> & displaced)
168 : {
169 : PetscFunctionBegin;
170 0 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
171 0 : DM_Moose * dmm = (DM_Moose *)dm->data;
172 0 : for (const auto & it : *(dmm->_uncontact_names))
173 : {
174 0 : uncontact_names.push_back(it.second);
175 0 : displaced.push_back((*dmm->_uncontact_displaced)[it.second]);
176 : }
177 0 : PetscFunctionReturn(PETSC_SUCCESS);
178 : }
179 :
180 : PetscErrorCode
181 0 : DMMooseGetSides(DM dm, std::vector<std::string> & side_names)
182 : {
183 : PetscFunctionBegin;
184 0 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
185 0 : DM_Moose * dmm = (DM_Moose *)dm->data;
186 0 : for (const auto & it : *(dmm->_side_ids))
187 0 : side_names.push_back(it.first);
188 0 : PetscFunctionReturn(PETSC_SUCCESS);
189 : }
190 :
191 : PetscErrorCode
192 0 : DMMooseGetUnSides(DM dm, std::vector<std::string> & side_names)
193 : {
194 : PetscFunctionBegin;
195 0 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
196 0 : DM_Moose * dmm = (DM_Moose *)dm->data;
197 0 : for (const auto & it : *(dmm->_unside_ids))
198 0 : side_names.push_back(it.first);
199 0 : PetscFunctionReturn(PETSC_SUCCESS);
200 : }
201 :
202 : PetscErrorCode
203 0 : DMMooseGetBlocks(DM dm, std::vector<std::string> & block_names)
204 : {
205 : PetscFunctionBegin;
206 0 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
207 0 : DM_Moose * dmm = (DM_Moose *)dm->data;
208 0 : for (const auto & it : *(dmm->_block_ids))
209 0 : block_names.push_back(it.first);
210 0 : PetscFunctionReturn(PETSC_SUCCESS);
211 : }
212 :
213 : PetscErrorCode
214 0 : DMMooseGetVariables(DM dm, std::vector<std::string> & var_names)
215 : {
216 : PetscFunctionBegin;
217 0 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
218 0 : DM_Moose * dmm = (DM_Moose *)(dm->data);
219 0 : for (const auto & it : *(dmm->_var_ids))
220 0 : var_names.push_back(it.first);
221 0 : PetscFunctionReturn(PETSC_SUCCESS);
222 : }
223 :
224 : PetscErrorCode
225 443 : DMMooseSetNonlinearSystem(DM dm, NonlinearSystemBase & nl)
226 : {
227 : PetscFunctionBegin;
228 443 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
229 443 : if (dm->setupcalled)
230 0 : SETERRQ(((PetscObject)dm)->comm,
231 : PETSC_ERR_ARG_WRONGSTATE,
232 : "Cannot reset the NonlinearSystem after DM has been set up.");
233 443 : DM_Moose * dmm = (DM_Moose *)(dm->data);
234 443 : dmm->_nl = &nl;
235 443 : PetscFunctionReturn(PETSC_SUCCESS);
236 : }
237 :
238 : PetscErrorCode
239 443 : DMMooseSetDofMap(DM dm, const DofMapBase & dof_map)
240 : {
241 : PetscFunctionBegin;
242 443 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
243 443 : if (dm->setupcalled)
244 0 : SETERRQ(((PetscObject)dm)->comm,
245 : PETSC_ERR_ARG_WRONGSTATE,
246 : "Cannot reset the degree of freedom map after DM has been set up.");
247 443 : DM_Moose * dmm = (DM_Moose *)(dm->data);
248 443 : dmm->_dof_map = &dof_map;
249 443 : PetscFunctionReturn(PETSC_SUCCESS);
250 : }
251 :
252 : PetscErrorCode
253 443 : DMMooseSetSystem(DM dm, const System & system)
254 : {
255 : PetscFunctionBegin;
256 443 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
257 443 : if (dm->setupcalled)
258 0 : SETERRQ(((PetscObject)dm)->comm,
259 : PETSC_ERR_ARG_WRONGSTATE,
260 : "Cannot reset the degree of freedom map after DM has been set up.");
261 443 : DM_Moose * dmm = (DM_Moose *)(dm->data);
262 443 : dmm->_system = &system;
263 443 : PetscFunctionReturn(PETSC_SUCCESS);
264 : }
265 :
266 : PetscErrorCode
267 443 : DMMooseSetName(DM dm, const std::string & dm_name)
268 : {
269 : PetscFunctionBegin;
270 443 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
271 443 : if (dm->setupcalled)
272 0 : SETERRQ(((PetscObject)dm)->comm,
273 : PETSC_ERR_ARG_WRONGSTATE,
274 : "Cannot reset the MOOSE DM name after DM has been set up.");
275 443 : DM_Moose * dmm = (DM_Moose *)(dm->data);
276 443 : *dmm->_name = dm_name;
277 443 : PetscFunctionReturn(PETSC_SUCCESS);
278 : }
279 :
280 : PetscErrorCode
281 314 : DMMooseSetParentDM(DM dm, DM_Moose * parent)
282 : {
283 : PetscFunctionBegin;
284 314 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
285 314 : if (dm->setupcalled)
286 0 : SETERRQ(((PetscObject)dm)->comm,
287 : PETSC_ERR_ARG_WRONGSTATE,
288 : "Cannot reset the parent DM after the child DM has been set up.");
289 :
290 314 : DM_Moose * dmm = (DM_Moose *)(dm->data);
291 314 : dmm->_parent = parent;
292 314 : PetscFunctionReturn(PETSC_SUCCESS);
293 : }
294 :
295 : PetscErrorCode
296 302 : DMMooseSetVariables(DM dm, const std::set<std::string> & vars)
297 : {
298 302 : DM_Moose * dmm = (DM_Moose *)dm->data;
299 :
300 : PetscFunctionBegin;
301 302 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
302 302 : if (dm->setupcalled)
303 0 : SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
304 302 : if (dmm->_vars)
305 0 : delete dmm->_vars;
306 302 : std::set<std::string> processed_vars;
307 604 : for (const auto & var_name : vars)
308 : {
309 : const auto * const var =
310 302 : dmm->_nl->hasVariable(var_name)
311 302 : ? static_cast<MooseVariableBase *>(&dmm->_nl->getVariable(0, var_name))
312 3 : : static_cast<MooseVariableBase *>(&dmm->_nl->getScalarVariable(0, var_name));
313 302 : if (var->isArray())
314 36 : for (const auto i : make_range(var->count()))
315 24 : processed_vars.insert(var->arrayVariableComponent(i));
316 : else
317 290 : processed_vars.insert(var_name);
318 : }
319 :
320 302 : dmm->_vars = new std::set<std::string>(std::move(processed_vars));
321 302 : PetscFunctionReturn(PETSC_SUCCESS);
322 302 : }
323 :
324 : PetscErrorCode
325 0 : DMMooseSetBlocks(DM dm, const std::set<std::string> & blocks)
326 : {
327 0 : DM_Moose * dmm = (DM_Moose *)dm->data;
328 :
329 : PetscFunctionBegin;
330 0 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
331 0 : if (dm->setupcalled)
332 0 : SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
333 0 : if (dmm->_blocks)
334 0 : delete dmm->_blocks;
335 0 : dmm->_blocks = new std::set<std::string>(blocks);
336 0 : PetscFunctionReturn(PETSC_SUCCESS);
337 : }
338 :
339 : PetscErrorCode
340 40 : DMMooseSetSides(DM dm, const std::set<std::string> & sides)
341 : {
342 40 : DM_Moose * dmm = (DM_Moose *)dm->data;
343 :
344 : PetscFunctionBegin;
345 40 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
346 40 : if (dm->setupcalled)
347 0 : SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
348 40 : if (dmm->_sides)
349 0 : delete dmm->_sides;
350 40 : dmm->_sides = new std::set<std::string>(sides);
351 40 : PetscFunctionReturn(PETSC_SUCCESS);
352 : }
353 :
354 : PetscErrorCode
355 28 : DMMooseSetUnSides(DM dm, const std::set<std::string> & unsides)
356 : {
357 28 : DM_Moose * dmm = (DM_Moose *)dm->data;
358 :
359 : PetscFunctionBegin;
360 28 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
361 28 : if (dm->setupcalled)
362 0 : SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
363 28 : if (dmm->_unsides)
364 0 : delete dmm->_unsides;
365 28 : dmm->_unsides = new std::set<std::string>(unsides);
366 28 : PetscFunctionReturn(PETSC_SUCCESS);
367 : }
368 :
369 : PetscErrorCode
370 12 : DMMooseSetUnSideByVar(DM dm, const std::set<std::string> & unside_by_var)
371 : {
372 12 : DM_Moose * dmm = (DM_Moose *)dm->data;
373 :
374 : PetscFunctionBegin;
375 12 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
376 12 : if (dm->setupcalled)
377 0 : SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
378 12 : if (dmm->_unside_by_var)
379 0 : delete dmm->_unside_by_var;
380 12 : dmm->_unside_by_var = new std::set<std::string>(unside_by_var);
381 12 : PetscFunctionReturn(PETSC_SUCCESS);
382 : }
383 :
384 : PetscErrorCode
385 0 : DMMooseSetContacts(DM dm,
386 : const std::vector<std::pair<std::string, std::string>> & contacts,
387 : const std::vector<PetscBool> & displaced)
388 : {
389 0 : DM_Moose * dmm = (DM_Moose *)dm->data;
390 :
391 : PetscFunctionBegin;
392 0 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
393 0 : if (dm->setupcalled)
394 0 : SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
395 0 : if (contacts.size() != displaced.size())
396 0 : LIBMESH_SETERRQ2(PETSC_COMM_SELF,
397 : PETSC_ERR_ARG_SIZ,
398 : "Nonmatching sizes of the contact and displaced arrays: %" LIBMESH_PETSCINT_FMT
399 : " != %" LIBMESH_PETSCINT_FMT,
400 : static_cast<PetscInt>(contacts.size()),
401 : static_cast<PetscInt>(displaced.size()));
402 0 : if (dmm->_contacts)
403 0 : delete dmm->_contacts;
404 0 : dmm->_contact_displaced->clear();
405 0 : dmm->_contacts = new std::set<DM_Moose::ContactName>();
406 0 : for (unsigned int i = 0; i < contacts.size(); ++i)
407 : {
408 0 : dmm->_contacts->insert(contacts[i]);
409 0 : dmm->_contact_displaced->insert(std::make_pair(contacts[i], displaced[i]));
410 : }
411 0 : PetscFunctionReturn(PETSC_SUCCESS);
412 : }
413 :
414 : PetscErrorCode
415 0 : DMMooseSetUnContacts(DM dm,
416 : const std::vector<std::pair<std::string, std::string>> & uncontacts,
417 : const std::vector<PetscBool> & displaced)
418 : {
419 0 : DM_Moose * dmm = (DM_Moose *)dm->data;
420 :
421 : PetscFunctionBegin;
422 0 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
423 0 : if (dm->setupcalled)
424 0 : SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
425 0 : if (uncontacts.size() != displaced.size())
426 0 : LIBMESH_SETERRQ2(
427 : PETSC_COMM_SELF,
428 : PETSC_ERR_ARG_SIZ,
429 : "Nonmatching sizes of the uncontact and displaced arrays: %" LIBMESH_PETSCINT_FMT
430 : " != %" LIBMESH_PETSCINT_FMT,
431 : static_cast<PetscInt>(uncontacts.size()),
432 : static_cast<PetscInt>(displaced.size()));
433 0 : if (dmm->_uncontacts)
434 0 : delete dmm->_uncontacts;
435 0 : dmm->_uncontact_displaced->clear();
436 0 : dmm->_uncontacts = new std::set<DM_Moose::ContactName>();
437 0 : for (unsigned int i = 0; i < uncontacts.size(); ++i)
438 : {
439 0 : dmm->_uncontacts->insert(uncontacts[i]);
440 0 : dmm->_uncontact_displaced->insert(std::make_pair(uncontacts[i], displaced[i]));
441 : }
442 0 : PetscFunctionReturn(PETSC_SUCCESS);
443 : }
444 :
445 : PetscErrorCode
446 0 : DMMooseGetNonlinearSystem(DM dm, NonlinearSystemBase *& nl)
447 : {
448 : PetscFunctionBegin;
449 0 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
450 0 : DM_Moose * dmm = (DM_Moose *)(dm->data);
451 0 : nl = dmm->_nl;
452 0 : PetscFunctionReturn(PETSC_SUCCESS);
453 : }
454 :
455 : PetscErrorCode
456 161 : DMMooseSetSplitNames(DM dm, const std::vector<std::string> & split_names)
457 : {
458 : PetscFunctionBegin;
459 161 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
460 161 : DM_Moose * dmm = (DM_Moose *)(dm->data);
461 :
462 161 : if (dmm->_splits)
463 : {
464 161 : for (auto & it : *(dmm->_splits))
465 : {
466 0 : LibmeshPetscCallQ(DMDestroy(&(it.second._dm)));
467 0 : LibmeshPetscCallQ(ISDestroy(&(it.second._rembedding)));
468 : }
469 161 : delete dmm->_splits;
470 161 : dmm->_splits = LIBMESH_PETSC_NULLPTR;
471 : }
472 161 : if (dmm->_splitlocs)
473 : {
474 0 : delete dmm->_splitlocs;
475 0 : dmm->_splitlocs = LIBMESH_PETSC_NULLPTR;
476 : }
477 161 : dmm->_splits = new std::map<std::string, DM_Moose::SplitInfo>();
478 161 : dmm->_splitlocs = new std::multimap<std::string, unsigned int>();
479 483 : for (unsigned int i = 0; i < split_names.size(); ++i)
480 : {
481 : DM_Moose::SplitInfo info;
482 322 : info._dm = LIBMESH_PETSC_NULLPTR;
483 322 : info._rembedding = LIBMESH_PETSC_NULLPTR;
484 322 : std::string name = split_names[i];
485 322 : (*dmm->_splits)[name] = info;
486 322 : dmm->_splitlocs->insert(std::make_pair(name, i));
487 322 : }
488 161 : PetscFunctionReturn(PETSC_SUCCESS);
489 : }
490 :
491 : PetscErrorCode
492 0 : DMMooseGetSplitNames(DM dm, std::vector<std::string> & split_names)
493 : {
494 : PetscFunctionBegin;
495 0 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
496 0 : DM_Moose * dmm = (DM_Moose *)(dm->data);
497 0 : if (!dm->setupcalled)
498 0 : SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "DM not set up");
499 0 : split_names.clear();
500 0 : split_names.reserve(dmm->_splitlocs->size());
501 0 : if (dmm->_splitlocs && dmm->_splitlocs->size())
502 0 : for (const auto & lit : *(dmm->_splitlocs))
503 : {
504 0 : std::string sname = lit.first;
505 0 : unsigned int sloc = lit.second;
506 0 : split_names[sloc] = sname;
507 0 : }
508 0 : PetscFunctionReturn(PETSC_SUCCESS);
509 : }
510 :
511 : static PetscErrorCode
512 314 : DMMooseGetEmbedding_Private(DM dm, IS * embedding)
513 : {
514 314 : DM_Moose * dmm = (DM_Moose *)dm->data;
515 :
516 : PetscFunctionBegin;
517 314 : if (!embedding)
518 0 : PetscFunctionReturn(PETSC_SUCCESS);
519 314 : if (!dmm->_embedding)
520 : {
521 : // The rules interpreting the coexistence of blocks (un)sides/(un)contacts are these
522 : // [sides and contacts behave similarly, so 'sides' means 'sides/contacts']
523 : // ['ANY' means 'not NONE' and covers 'ALL' as well, unless there is a specific 'ALL' clause,
524 : // which overrides 'ANY'; 'NOT ALL' means not ALL and not NONE]
525 : // [there are always some blocks, since by default 'ALL' is assumed, unless it is overridden by
526 : // a specific list, which implies ANY]
527 : // In general,
528 : // (1) ALL blocks and ANY sides are interpreted as the INTERSECTION of blocks and sides,
529 : // equivalent to just the sides (since ALL blocks are assumed to be a cover).
530 : // (2) NOT ALL blocks and ANY or NO sides are interpreted as the UNION of blocks and sides.
531 : // (3a) ANY unsides and ANY blocks are interpreted as the DIFFERENCE of blocks and unsides.
532 : // (3b) ANY unsides and ANY sides are interpreted as the DIFFERENCE of sides and unsides.
533 : // (4) NO unsides means NO DIFFERENCE is needed.
534 : // The result is easily computed by first computing the result of (1 & 2) followed by difference
535 : // with the result of (3 & 4).
536 : // To simply (1 & 2) observe the following:
537 : // - The intersection is computed only if ALL blocks and ANY sides, and the result is the sides,
538 : // so block dofs do not need to be computed.
539 : // - Otherwise the union is computed, and initially consists of the blocks' dofs, to which the
540 : // sides' dofs are added, if ANY.
541 : // - The result is called 'indices'
542 : // To satisfy (3 & 4) simply cmpute subtrahend set 'unindices' as all of the unsides' dofs:
543 : // Then take the set difference of 'indices' and 'unindices', putting the result in 'dindices'.
544 314 : if (!dmm->_all_vars || !dmm->_all_blocks || !dmm->_nosides || !dmm->_nounsides ||
545 12 : !dmm->_nounside_by_var || !dmm->_nocontacts || !dmm->_nouncontacts)
546 : {
547 314 : auto & dofmap = *dmm->_dof_map;
548 : // Put this outside the lambda scope to avoid constant memory reallocation
549 314 : std::vector<dof_id_type> node_indices;
550 : auto process_nodal_dof_indices =
551 40325 : [&dofmap, &node_indices](const Node & node,
552 : const unsigned int var_num,
553 : std::set<dof_id_type> & local_indices,
554 : std::set<dof_id_type> * const nonlocal_indices = nullptr)
555 : {
556 40325 : dofmap.dof_indices(&node, node_indices, var_num);
557 81370 : for (const auto index : node_indices)
558 : {
559 41045 : if (index >= dofmap.first_dof() && index < dofmap.end_dof())
560 40937 : local_indices.insert(index);
561 108 : else if (nonlocal_indices)
562 0 : nonlocal_indices->insert(index);
563 : }
564 40325 : };
565 :
566 : auto process_elem_dof_indices =
567 35482 : [&dofmap](const std::vector<dof_id_type> & elem_indices,
568 : std::set<dof_id_type> & local_indices,
569 : std::set<dof_id_type> * const nonlocal_indices = nullptr)
570 : {
571 133528 : for (const auto index : elem_indices)
572 : {
573 98046 : if (index >= dofmap.first_dof() && index < dofmap.end_dof())
574 96023 : local_indices.insert(index);
575 2023 : else if (nonlocal_indices)
576 0 : nonlocal_indices->insert(index);
577 : }
578 35482 : };
579 :
580 314 : std::set<dof_id_type> indices;
581 314 : std::set<dof_id_type> unindices;
582 314 : std::set<dof_id_type> cached_indices;
583 314 : std::set<dof_id_type> cached_unindices;
584 314 : auto & lm_mesh = dmm->_system->get_mesh();
585 314 : const auto & node_to_elem_map = dmm->_nl->feProblem().mesh().nodeToElemMap();
586 652 : for (const auto & vit : *(dmm->_var_ids))
587 : {
588 338 : unsigned int v = vit.second;
589 : // Iterate only over this DM's blocks.
590 338 : if (!dmm->_all_blocks || (dmm->_nosides && dmm->_nocontacts))
591 744 : for (const auto & bit : *(dmm->_block_ids))
592 : {
593 406 : subdomain_id_type b = bit.second;
594 812 : for (const auto & elem : as_range(lm_mesh.active_local_subdomain_elements_begin(b),
595 36700 : lm_mesh.active_local_subdomain_elements_end(b)))
596 : {
597 : // Get the degree of freedom indices for the given variable off the current element.
598 35482 : std::vector<dof_id_type> evindices;
599 35482 : dofmap.dof_indices(elem, evindices, v);
600 35482 : process_elem_dof_indices(evindices, indices);
601 35888 : }
602 :
603 : // Sometime, we own nodes but do not own the elements the nodes are connected to
604 : {
605 406 : bool is_on_current_block = false;
606 83801 : for (auto & node : lm_mesh.local_node_ptr_range())
607 : {
608 83395 : const unsigned int n_comp = node->n_comp(dmm->_system->number(), v);
609 :
610 : // skip it if no dof
611 83395 : if (!n_comp)
612 43538 : continue;
613 :
614 52493 : auto node_to_elem_pair = node_to_elem_map.find(node->id());
615 52493 : is_on_current_block = false;
616 114219 : for (const auto & elem_num : node_to_elem_pair->second)
617 : {
618 : // if one of incident elements belongs to a block, we consider
619 : // the node lives in the block
620 101583 : const Elem & neighbor_elem = lm_mesh.elem_ref(elem_num);
621 101583 : if (neighbor_elem.subdomain_id() == b)
622 : {
623 39857 : is_on_current_block = true;
624 39857 : break;
625 : }
626 : }
627 : // we add indices for the current block only
628 52493 : if (!is_on_current_block)
629 12636 : continue;
630 :
631 39857 : process_nodal_dof_indices(*node, v, indices);
632 406 : }
633 : }
634 : }
635 :
636 : // Iterate over the sides from this split.
637 338 : if (dmm->_side_ids->size())
638 : {
639 : // For some reason the following may return an empty node list
640 : // std::vector<dof_id_type> snodes;
641 : // std::vector<boundary_id_type> sides;
642 : // dmm->nl->system().get_mesh().get_boundary_info().build_node_list(snodes, sides);
643 : // // FIXME: make an array of (snode,side) pairs, sort on side and use std::lower_bound
644 : // from <algorithm>
645 : // for (dof_id_type i = 0; i < sides.size(); ++i) {
646 : // boundary_id_type s = sides[i];
647 : // if (!dmm->sidenames->count(s)) continue;
648 : // const Node& node = dmm->nl->system().get_mesh().node_ref(snodes[i]);
649 : // // determine v's dof on node and insert into indices
650 : // }
651 40 : ConstBndNodeRange & bnodes = *dmm->_nl->mesh().getBoundaryNodeRange();
652 520 : for (const auto & bnode : bnodes)
653 : {
654 480 : BoundaryID boundary_id = bnode->_bnd_id;
655 480 : if (dmm->_side_names->find(boundary_id) == dmm->_side_names->end())
656 252 : continue;
657 :
658 228 : const Node * node = bnode->_node;
659 228 : process_nodal_dof_indices(*node, v, indices);
660 : }
661 : }
662 :
663 : // Iterate over the sides excluded from this split.
664 338 : if (dmm->_unside_ids->size())
665 : {
666 28 : ConstBndNodeRange & bnodes = *dmm->_nl->mesh().getBoundaryNodeRange();
667 364 : for (const auto & bnode : bnodes)
668 : {
669 336 : BoundaryID boundary_id = bnode->_bnd_id;
670 336 : if (dmm->_unside_names->find(boundary_id) == dmm->_unside_names->end())
671 168 : continue;
672 168 : const Node * node = bnode->_node;
673 168 : process_nodal_dof_indices(*node, v, unindices);
674 : }
675 : }
676 338 : if (dmm->_unside_by_var_set->size())
677 : {
678 24 : std::set<BoundaryID> eligible_bids;
679 72 : for (const auto & [bid, var] : *(dmm->_unside_by_var_set))
680 48 : if (var == v)
681 24 : eligible_bids.insert(bid);
682 :
683 24 : ConstBndNodeRange & bnodes = *dmm->_nl->mesh().getBoundaryNodeRange();
684 312 : for (const auto & bnode : bnodes)
685 : {
686 288 : BoundaryID boundary_id = bnode->_bnd_id;
687 288 : if (eligible_bids.count(boundary_id))
688 : {
689 72 : const Node * node = bnode->_node;
690 72 : process_nodal_dof_indices(*node, v, unindices);
691 : }
692 : }
693 24 : }
694 :
695 : auto process_contact_all_nodes =
696 0 : [dmm, process_nodal_dof_indices, v](const auto & contact_names,
697 : auto & indices_to_insert_to)
698 : {
699 0 : std::set<boundary_id_type> bc_id_set;
700 : // loop over contacts
701 0 : for (const auto & [contact_bid_pair, contact_bname_pair] : contact_names)
702 : {
703 0 : libmesh_ignore(contact_bname_pair);
704 0 : bc_id_set.insert(contact_bid_pair.first); // primary
705 0 : bc_id_set.insert(contact_bid_pair.second); // secondary
706 : }
707 : // loop over boundary elements
708 0 : ConstBndElemRange & range = *dmm->_nl->feProblem().mesh().getBoundaryElementRange();
709 0 : for (const auto & belem : range)
710 : {
711 0 : const Elem * elem_bdry = belem->_elem;
712 0 : const auto side = belem->_side;
713 0 : BoundaryID boundary_id = belem->_bnd_id;
714 :
715 0 : if (bc_id_set.find(boundary_id) == bc_id_set.end())
716 0 : continue;
717 :
718 0 : for (const auto node_idx : elem_bdry->node_index_range())
719 0 : if (elem_bdry->is_node_on_side(node_idx, side))
720 0 : process_nodal_dof_indices(elem_bdry->node_ref(node_idx), v, indices_to_insert_to);
721 : }
722 0 : };
723 :
724 : auto process_contact_some_nodes =
725 0 : [dmm, process_nodal_dof_indices, v, &dofmap, &lm_mesh, process_elem_dof_indices](
726 : const auto & contact_names,
727 : auto & indices_to_insert_to,
728 : auto & nonlocal_indices_to_insert_to)
729 : {
730 0 : std::vector<dof_id_type> evindices;
731 0 : for (const auto & it : contact_names)
732 : {
733 0 : PetscBool displaced = (*dmm->_uncontact_displaced)[it.second];
734 : PenetrationLocator * locator;
735 0 : if (displaced)
736 : {
737 0 : std::shared_ptr<DisplacedProblem> displaced_problem =
738 0 : dmm->_nl->feProblem().getDisplacedProblem();
739 0 : if (!displaced_problem)
740 : {
741 0 : std::ostringstream err;
742 0 : err << "Cannot use a displaced uncontact (" << it.second.first << ","
743 0 : << it.second.second << ") with an undisplaced problem";
744 0 : mooseError(err.str());
745 0 : }
746 0 : locator = displaced_problem->geomSearchData()._penetration_locators[it.first];
747 0 : }
748 : else
749 0 : locator = dmm->_nl->feProblem().geomSearchData()._penetration_locators[it.first];
750 :
751 0 : evindices.clear();
752 : // penetration locator
753 0 : auto lend = locator->_penetration_info.end();
754 0 : for (auto lit = locator->_penetration_info.begin(); lit != lend; ++lit)
755 : {
756 0 : const dof_id_type secondary_node_num = lit->first;
757 0 : PenetrationInfo * pinfo = lit->second;
758 0 : if (pinfo && pinfo->isCaptured())
759 : {
760 0 : const Node & secondary_node = lm_mesh.node_ref(secondary_node_num);
761 0 : process_nodal_dof_indices(
762 : secondary_node, v, indices_to_insert_to, &nonlocal_indices_to_insert_to);
763 :
764 : // indices for primary element
765 0 : evindices.clear();
766 0 : const Elem * primary_side = pinfo->_side;
767 0 : dofmap.dof_indices(primary_side, evindices, v);
768 0 : process_elem_dof_indices(
769 : evindices, indices_to_insert_to, &nonlocal_indices_to_insert_to);
770 : } // if pinfo
771 : } // for penetration
772 : } // for contact name
773 0 : };
774 :
775 : // Include all nodes on the contact surfaces
776 338 : if (dmm->_contact_names->size() && dmm->_include_all_contact_nodes)
777 0 : process_contact_all_nodes(*dmm->_contact_names, indices);
778 :
779 : // Iterate over the contacts included in this split.
780 338 : if (dmm->_contact_names->size() && !(dmm->_include_all_contact_nodes))
781 0 : process_contact_some_nodes(*dmm->_contact_names, indices, cached_indices);
782 :
783 : // Exclude all nodes on the contact surfaces
784 338 : if (dmm->_uncontact_names->size() && dmm->_include_all_contact_nodes)
785 0 : process_contact_all_nodes(*dmm->_uncontact_names, unindices);
786 :
787 : // Iterate over the contacts excluded from this split.
788 338 : if (dmm->_uncontact_names->size() && !(dmm->_include_all_contact_nodes))
789 0 : process_contact_some_nodes(*dmm->_uncontact_names, unindices, cached_unindices);
790 : } // variables
791 :
792 314 : std::vector<dof_id_type> local_vec_indices(cached_indices.size());
793 314 : std::copy(cached_indices.begin(), cached_indices.end(), local_vec_indices.begin());
794 314 : if (dmm->_contact_names->size() && !(dmm->_include_all_contact_nodes))
795 0 : dmm->_nl->feProblem().mesh().comm().allgather(local_vec_indices, false);
796 : // insert indices
797 314 : for (const auto & dof : local_vec_indices)
798 0 : if (dof >= dofmap.first_dof() && dof < dofmap.end_dof())
799 0 : indices.insert(dof);
800 :
801 314 : local_vec_indices.clear();
802 314 : local_vec_indices.resize(cached_unindices.size());
803 314 : std::copy(cached_unindices.begin(), cached_unindices.end(), local_vec_indices.begin());
804 314 : if (dmm->_uncontact_names->size() && !(dmm->_include_all_contact_nodes))
805 0 : dmm->_nl->feProblem().mesh().comm().allgather(local_vec_indices, false);
806 : // insert unindices
807 314 : for (const auto & dof : local_vec_indices)
808 0 : if (dof >= dofmap.first_dof() && dof < dofmap.end_dof())
809 0 : unindices.insert(dof);
810 :
811 314 : std::set<dof_id_type> dindices;
812 314 : std::set_difference(indices.begin(),
813 : indices.end(),
814 : unindices.begin(),
815 : unindices.end(),
816 : std::inserter(dindices, dindices.end()));
817 : PetscInt * darray;
818 314 : LibmeshPetscCallQ(PetscMalloc(sizeof(PetscInt) * dindices.size(), &darray));
819 314 : dof_id_type i = 0;
820 38173 : for (const auto & dof : dindices)
821 : {
822 37859 : darray[i] = dof;
823 37859 : ++i;
824 : }
825 314 : LibmeshPetscCallQ(ISCreateGeneral(
826 : ((PetscObject)dm)->comm, dindices.size(), darray, PETSC_OWN_POINTER, &dmm->_embedding));
827 314 : }
828 : else
829 : {
830 : // if (dmm->allblocks && dmm->allvars && dmm->nosides && dmm->nounsides && dmm->nocontacts &&
831 : // dmm->nouncontacts)
832 : // DMCreateGlobalVector is defined()
833 : Vec v;
834 : PetscInt low, high;
835 :
836 0 : LibmeshPetscCallQ(DMCreateGlobalVector(dm, &v));
837 0 : LibmeshPetscCallQ(VecGetOwnershipRange(v, &low, &high));
838 0 : LibmeshPetscCallQ(
839 : ISCreateStride(((PetscObject)dm)->comm, (high - low), low, 1, &dmm->_embedding));
840 : }
841 : }
842 314 : LibmeshPetscCallQ(PetscObjectReference((PetscObject)(dmm->_embedding)));
843 314 : *embedding = dmm->_embedding;
844 :
845 314 : PetscFunctionReturn(PETSC_SUCCESS);
846 : }
847 :
848 : static PetscErrorCode
849 157 : DMCreateFieldDecomposition_Moose(
850 : DM dm, PetscInt * len, char *** namelist, IS ** islist, DM ** dmlist)
851 : {
852 157 : DM_Moose * dmm = (DM_Moose *)(dm->data);
853 :
854 : PetscFunctionBegin;
855 :
856 157 : PetscInt split_size_sum = 0;
857 :
858 : /* Only called after DMSetUp(). */
859 157 : if (!dmm->_splitlocs)
860 0 : PetscFunctionReturn(PETSC_SUCCESS);
861 157 : *len = dmm->_splitlocs->size();
862 157 : if (namelist)
863 157 : LibmeshPetscCallQ(PetscMalloc(*len * sizeof(char *), namelist));
864 157 : if (islist)
865 157 : LibmeshPetscCallQ(PetscMalloc(*len * sizeof(IS), islist));
866 157 : if (dmlist)
867 157 : LibmeshPetscCallQ(PetscMalloc(*len * sizeof(DM), dmlist));
868 471 : for (const auto & dit : *(dmm->_splitlocs))
869 : {
870 314 : unsigned int d = dit.second;
871 314 : std::string dname = dit.first;
872 314 : DM_Moose::SplitInfo & dinfo = (*dmm->_splits)[dname];
873 314 : if (!dinfo._dm)
874 : {
875 314 : LibmeshPetscCallQ(DMCreateMoose(
876 : ((PetscObject)dm)->comm, *dmm->_nl, *dmm->_dof_map, *dmm->_system, dname, &dinfo._dm));
877 314 : LibmeshPetscCallQ(
878 : PetscObjectSetOptionsPrefix((PetscObject)dinfo._dm, ((PetscObject)dm)->prefix));
879 314 : std::string suffix = std::string("fieldsplit_") + dname + "_";
880 314 : LibmeshPetscCallQ(PetscObjectAppendOptionsPrefix((PetscObject)dinfo._dm, suffix.c_str()));
881 314 : LibmeshPetscCallQ(DMMooseSetParentDM(dinfo._dm, dmm));
882 314 : }
883 314 : LibmeshPetscCallQ(DMSetFromOptions(dinfo._dm));
884 314 : LibmeshPetscCallQ(DMSetUp(dinfo._dm));
885 314 : if (namelist)
886 314 : LibmeshPetscCallQ(PetscStrallocpy(dname.c_str(), (*namelist) + d));
887 314 : if (islist)
888 : {
889 314 : if (!dinfo._rembedding)
890 : {
891 : IS dembedding, lembedding;
892 314 : LibmeshPetscCallQ(DMMooseGetEmbedding_Private(dinfo._dm, &dembedding));
893 314 : if (dmm->_embedding)
894 : {
895 : // Create a relative embedding into the parent's index space.
896 56 : LibmeshPetscCallQ(ISEmbed(dembedding, dmm->_embedding, PETSC_TRUE, &lembedding));
897 : const PetscInt * lindices;
898 : PetscInt len, dlen, llen, *rindices, off, i;
899 56 : LibmeshPetscCallQ(ISGetLocalSize(dembedding, &dlen));
900 56 : LibmeshPetscCallQ(ISGetLocalSize(lembedding, &llen));
901 56 : if (llen != dlen)
902 0 : LIBMESH_SETERRQ1(
903 : ((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Failed to embed split %u", d);
904 56 : LibmeshPetscCallQ(ISDestroy(&dembedding));
905 : // Convert local embedding to global (but still relative) embedding
906 56 : LibmeshPetscCallQ(PetscMalloc(llen * sizeof(PetscInt), &rindices));
907 56 : LibmeshPetscCallQ(ISGetIndices(lembedding, &lindices));
908 56 : LibmeshPetscCallQ(PetscMemcpy(rindices, lindices, llen * sizeof(PetscInt)));
909 56 : LibmeshPetscCallQ(ISDestroy(&lembedding));
910 : // We could get the index offset from a corresponding global vector, but subDMs don't yet
911 : // have global vectors
912 56 : LibmeshPetscCallQ(ISGetLocalSize(dmm->_embedding, &len));
913 :
914 56 : MPI_Scan(&len,
915 : &off,
916 : 1,
917 : #ifdef PETSC_USE_64BIT_INDICES
918 : MPI_LONG_LONG_INT,
919 : #else
920 : MPI_INT,
921 : #endif
922 : MPI_SUM,
923 : ((PetscObject)dm)->comm);
924 :
925 56 : off -= len;
926 242 : for (i = 0; i < llen; ++i)
927 186 : rindices[i] += off;
928 56 : LibmeshPetscCallQ(ISCreateGeneral(
929 : ((PetscObject)dm)->comm, llen, rindices, PETSC_OWN_POINTER, &(dinfo._rembedding)));
930 : }
931 : else
932 : {
933 258 : dinfo._rembedding = dembedding;
934 : }
935 : }
936 314 : LibmeshPetscCallQ(PetscObjectReference((PetscObject)(dinfo._rembedding)));
937 314 : (*islist)[d] = dinfo._rembedding;
938 : PetscInt is_size;
939 314 : LibmeshPetscCallQ(ISGetLocalSize(dinfo._rembedding, &is_size));
940 314 : split_size_sum += is_size;
941 : }
942 314 : if (dmlist)
943 : {
944 314 : LibmeshPetscCallQ(PetscObjectReference((PetscObject)dinfo._dm));
945 314 : (*dmlist)[d] = dinfo._dm;
946 : }
947 314 : }
948 :
949 : mooseAssert(islist, "What does it even mean if this is NULL?");
950 :
951 157 : if (dmm->_parent)
952 28 : dmm->_parent->checkChildSize(dm, split_size_sum, *dmm->_name);
953 : else
954 129 : checkSize(*dmm->_name, split_size_sum, dmm->_dof_map->n_local_dofs());
955 :
956 149 : PetscFunctionReturn(PETSC_SUCCESS);
957 : }
958 :
959 : static PetscErrorCode
960 0 : DMCreateDomainDecomposition_Moose(
961 : DM dm, PetscInt * len, char *** namelist, IS ** innerislist, IS ** outerislist, DM ** dmlist)
962 : {
963 : PetscFunctionBegin;
964 : /* Use DMCreateFieldDecomposition_Moose() to obtain everything but outerislist, which is currently
965 : * LIBMESH_PETSC_NULLPTR. */
966 0 : if (outerislist)
967 0 : *outerislist = LIBMESH_PETSC_NULLPTR; /* FIX: allow mesh-based overlap. */
968 0 : LibmeshPetscCallQ(DMCreateFieldDecomposition_Moose(dm, len, namelist, innerislist, dmlist));
969 0 : PetscFunctionReturn(PETSC_SUCCESS);
970 : }
971 :
972 : static PetscErrorCode
973 0 : DMMooseFunction(DM dm, Vec x, Vec r)
974 : {
975 : PetscFunctionBegin;
976 : libmesh_assert(x);
977 : libmesh_assert(r);
978 :
979 0 : NonlinearSystemBase * nl = NULL;
980 0 : LibmeshPetscCallQ(DMMooseGetNonlinearSystem(dm, nl));
981 0 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(nl->system().solution.get());
982 0 : PetscVector<Number> X_global(x, nl->comm()), R(r, nl->comm());
983 :
984 : // Use the system's update() to get a good local version of the
985 : // parallel solution. system.update() does change the residual vector,
986 : // so there's no reason to swap PETSc's residual into the system for
987 : // this step.
988 0 : X_global.swap(X_sys);
989 0 : nl->system().update();
990 0 : X_global.swap(X_sys);
991 :
992 : // Enforce constraints (if any) exactly on the
993 : // current_local_solution. This is the solution vector that is
994 : // actually used in the computation of the residual below, and is
995 : // not locked by debug-enabled PETSc the way that "x" is.
996 0 : nl->system().get_dof_map().enforce_constraints_exactly(nl->system(),
997 0 : nl->system().current_local_solution.get());
998 :
999 : // Zero the residual vector before assembling
1000 0 : R.zero();
1001 :
1002 : // if the user has provided both function pointers and objects only the pointer
1003 : // will be used, so catch that as an error
1004 0 : if (nl->nonlinearSolver()->residual && nl->nonlinearSolver()->residual_object)
1005 : {
1006 0 : std::ostringstream err;
1007 0 : err << "ERROR: cannot specifiy both a function and object to compute the Residual!"
1008 0 : << std::endl;
1009 0 : mooseError(err.str());
1010 0 : }
1011 0 : if (nl->nonlinearSolver()->matvec && nl->nonlinearSolver()->residual_and_jacobian_object)
1012 : {
1013 0 : std::ostringstream err;
1014 : err << "ERROR: cannot specifiy both a function and object to compute the combined Residual & "
1015 0 : "Jacobian!"
1016 0 : << std::endl;
1017 0 : mooseError(err.str());
1018 0 : }
1019 0 : if (nl->nonlinearSolver()->residual != NULL)
1020 0 : nl->nonlinearSolver()->residual(
1021 0 : *(nl->system().current_local_solution.get()), R, nl->nonlinearSolver()->system());
1022 0 : else if (nl->nonlinearSolver()->residual_object != NULL)
1023 0 : nl->nonlinearSolver()->residual_object->residual(
1024 0 : *(nl->system().current_local_solution.get()), R, nl->nonlinearSolver()->system());
1025 0 : else if (nl->nonlinearSolver()->matvec != NULL)
1026 0 : nl->nonlinearSolver()->matvec(
1027 0 : *(nl->system().current_local_solution.get()), &R, NULL, nl->nonlinearSolver()->system());
1028 0 : else if (nl->nonlinearSolver()->residual_and_jacobian_object != NULL)
1029 0 : nl->nonlinearSolver()->residual_and_jacobian_object->residual_and_jacobian(
1030 0 : *(nl->system().current_local_solution.get()), &R, NULL, nl->nonlinearSolver()->system());
1031 : else
1032 : {
1033 0 : std::ostringstream err;
1034 0 : err << "No suitable residual computation routine found";
1035 0 : mooseError(err.str());
1036 0 : }
1037 0 : R.close();
1038 0 : PetscFunctionReturn(PETSC_SUCCESS);
1039 0 : }
1040 :
1041 : static PetscErrorCode
1042 0 : SNESFunction_DMMoose(SNES, Vec x, Vec r, void * ctx)
1043 : {
1044 0 : DM dm = (DM)ctx;
1045 :
1046 : PetscFunctionBegin;
1047 0 : LibmeshPetscCallQ(DMMooseFunction(dm, x, r));
1048 0 : PetscFunctionReturn(PETSC_SUCCESS);
1049 : }
1050 :
1051 : static PetscErrorCode
1052 0 : DMMooseJacobian(DM dm, Vec x, Mat jac, Mat pc)
1053 : {
1054 0 : NonlinearSystemBase * nl = NULL;
1055 :
1056 : PetscFunctionBegin;
1057 0 : LibmeshPetscCallQ(DMMooseGetNonlinearSystem(dm, nl));
1058 :
1059 0 : PetscMatrix<Number> the_pc(pc, nl->comm());
1060 0 : PetscMatrix<Number> Jac(jac, nl->comm());
1061 0 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(nl->system().solution.get());
1062 0 : PetscVector<Number> X_global(x, nl->comm());
1063 :
1064 : // Set the dof maps
1065 0 : the_pc.attach_dof_map(nl->system().get_dof_map());
1066 0 : Jac.attach_dof_map(nl->system().get_dof_map());
1067 :
1068 : // Use the system's update() to get a good local version of the
1069 : // parallel solution. system.update() does change the Jacobian, so
1070 : // there's no reason to swap PETSc's Jacobian into the system for
1071 : // this step.
1072 0 : X_global.swap(X_sys);
1073 0 : nl->system().update();
1074 0 : X_global.swap(X_sys);
1075 :
1076 : // Enforce constraints (if any) exactly on the
1077 : // current_local_solution. This is the solution vector that is
1078 : // actually used in the computation of the Jacobian below, and is
1079 : // not locked by debug-enabled PETSc the way that "x" is.
1080 0 : nl->system().get_dof_map().enforce_constraints_exactly(nl->system(),
1081 0 : nl->system().current_local_solution.get());
1082 :
1083 : // Zero out the preconditioner before computing the Jacobian.
1084 0 : the_pc.zero();
1085 :
1086 : // if the user has provided both function pointers and objects only the pointer
1087 : // will be used, so catch that as an error
1088 0 : if (nl->nonlinearSolver()->jacobian && nl->nonlinearSolver()->jacobian_object)
1089 : {
1090 0 : std::ostringstream err;
1091 0 : err << "ERROR: cannot specifiy both a function and object to compute the Jacobian!"
1092 0 : << std::endl;
1093 0 : mooseError(err.str());
1094 0 : }
1095 0 : if (nl->nonlinearSolver()->matvec && nl->nonlinearSolver()->residual_and_jacobian_object)
1096 : {
1097 0 : std::ostringstream err;
1098 : err << "ERROR: cannot specifiy both a function and object to compute the combined Residual & "
1099 0 : "Jacobian!"
1100 0 : << std::endl;
1101 0 : mooseError(err.str());
1102 0 : }
1103 0 : if (nl->nonlinearSolver()->jacobian != NULL)
1104 0 : nl->nonlinearSolver()->jacobian(
1105 0 : *(nl->system().current_local_solution.get()), the_pc, nl->nonlinearSolver()->system());
1106 0 : else if (nl->nonlinearSolver()->jacobian_object != NULL)
1107 0 : nl->nonlinearSolver()->jacobian_object->jacobian(
1108 0 : *(nl->system().current_local_solution.get()), the_pc, nl->nonlinearSolver()->system());
1109 0 : else if (nl->nonlinearSolver()->matvec != NULL)
1110 0 : nl->nonlinearSolver()->matvec(*(nl->system().current_local_solution.get()),
1111 : NULL,
1112 : &the_pc,
1113 0 : nl->nonlinearSolver()->system());
1114 0 : else if (nl->nonlinearSolver()->residual_and_jacobian_object != NULL)
1115 0 : nl->nonlinearSolver()->residual_and_jacobian_object->residual_and_jacobian(
1116 0 : *(nl->system().current_local_solution.get()),
1117 : NULL,
1118 : &the_pc,
1119 0 : nl->nonlinearSolver()->system());
1120 : else
1121 : {
1122 0 : std::ostringstream err;
1123 0 : err << "No suitable Jacobian routine or object";
1124 0 : mooseError(err.str());
1125 0 : }
1126 0 : the_pc.close();
1127 0 : Jac.close();
1128 0 : PetscFunctionReturn(PETSC_SUCCESS);
1129 0 : }
1130 :
1131 : static PetscErrorCode
1132 0 : SNESJacobian_DMMoose(SNES, Vec x, Mat jac, Mat pc, void * ctx)
1133 : {
1134 0 : DM dm = (DM)ctx;
1135 :
1136 : PetscFunctionBegin;
1137 0 : LibmeshPetscCallQ(DMMooseJacobian(dm, x, jac, pc));
1138 0 : PetscFunctionReturn(PETSC_SUCCESS);
1139 : }
1140 :
1141 : static PetscErrorCode
1142 0 : DMVariableBounds_Moose(DM dm, Vec xl, Vec xu)
1143 : {
1144 0 : NonlinearSystemBase * nl = NULL;
1145 :
1146 : PetscFunctionBegin;
1147 0 : LibmeshPetscCallQ(DMMooseGetNonlinearSystem(dm, nl));
1148 :
1149 0 : PetscVector<Number> XL(xl, nl->comm());
1150 0 : PetscVector<Number> XU(xu, nl->comm());
1151 :
1152 0 : LibmeshPetscCallQ(VecSet(xl, PETSC_NINFINITY));
1153 0 : LibmeshPetscCallQ(VecSet(xu, PETSC_INFINITY));
1154 0 : if (nl->nonlinearSolver()->bounds != NULL)
1155 0 : nl->nonlinearSolver()->bounds(XL, XU, nl->nonlinearSolver()->system());
1156 0 : else if (nl->nonlinearSolver()->bounds_object != NULL)
1157 0 : nl->nonlinearSolver()->bounds_object->bounds(XL, XU, nl->nonlinearSolver()->system());
1158 : else
1159 0 : SETERRQ(
1160 : ((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "No bounds calculation in this Moose object");
1161 0 : PetscFunctionReturn(PETSC_SUCCESS);
1162 0 : }
1163 :
1164 : static PetscErrorCode
1165 52 : DMCreateGlobalVector_Moose(DM dm, Vec * x)
1166 : {
1167 52 : DM_Moose * dmm = (DM_Moose *)(dm->data);
1168 :
1169 : PetscFunctionBegin;
1170 52 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
1171 52 : if (!dmm->_nl)
1172 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1173 :
1174 52 : NumericVector<Number> * nv = (dmm->_system->solution).get();
1175 52 : PetscVector<Number> * pv = dynamic_cast<PetscVector<Number> *>(nv);
1176 52 : Vec v = pv->vec();
1177 : /* Unfortunately, currently this does not produce a ghosted vector, so nonlinear subproblem solves
1178 : aren't going to be easily available.
1179 : Should work fine for getting vectors out for linear subproblem solvers. */
1180 52 : if (dmm->_embedding)
1181 : {
1182 : PetscInt n;
1183 52 : LibmeshPetscCallQ(VecCreate(((PetscObject)v)->comm, x));
1184 52 : LibmeshPetscCallQ(ISGetLocalSize(dmm->_embedding, &n));
1185 52 : LibmeshPetscCallQ(VecSetSizes(*x, n, PETSC_DETERMINE));
1186 52 : LibmeshPetscCallQ(VecSetType(*x, ((PetscObject)v)->type_name));
1187 52 : LibmeshPetscCallQ(VecSetFromOptions(*x));
1188 52 : LibmeshPetscCallQ(VecSetUp(*x));
1189 : }
1190 : else
1191 0 : LibmeshPetscCallQ(VecDuplicate(v, x));
1192 :
1193 : #if PETSC_RELEASE_LESS_THAN(3, 13, 0)
1194 : LibmeshPetscCallQ(PetscObjectCompose((PetscObject)*x, "DM", (PetscObject)dm));
1195 : #else
1196 52 : LibmeshPetscCallQ(VecSetDM(*x, dm));
1197 : #endif
1198 52 : PetscFunctionReturn(PETSC_SUCCESS);
1199 : }
1200 :
1201 : static PetscErrorCode
1202 0 : DMCreateMatrix_Moose(DM dm, Mat * A)
1203 : {
1204 0 : DM_Moose * dmm = (DM_Moose *)(dm->data);
1205 : MatType type;
1206 :
1207 : PetscFunctionBegin;
1208 0 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
1209 0 : if (!dmm->_nl)
1210 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1211 0 : LibmeshPetscCallQ(DMGetMatType(dm, &type));
1212 :
1213 : /*
1214 : The simplest thing for now: compute the sparsity_pattern using dof_map and init the matrix using
1215 : that info.
1216 : TODO: compute sparsity restricted to this DM's blocks, variables and sides.
1217 : Even fancier: compute the sparsity of the coupling of a contact secondary to the contact primary.
1218 : In any event, here we are in control of the matrix type and structure.
1219 : */
1220 0 : const auto & dof_map = *dmm->_dof_map;
1221 : PetscInt M, N, m, n;
1222 : MPI_Comm comm;
1223 0 : M = dof_map.n_dofs();
1224 0 : N = M;
1225 0 : m = static_cast<PetscInt>(dof_map.n_local_dofs());
1226 0 : n = m;
1227 0 : LibmeshPetscCallQ(PetscObjectGetComm((PetscObject)dm, &comm));
1228 0 : LibmeshPetscCallQ(MatCreate(comm, A));
1229 0 : LibmeshPetscCallQ(MatSetSizes(*A, m, n, M, N));
1230 0 : LibmeshPetscCallQ(MatSetType(*A, type));
1231 : /* TODO: set the prefix for *A and MatSetFromOptions(*A)? Might override the type and other
1232 : * settings made here. */
1233 0 : LibmeshPetscCallQ(MatSetUp(*A));
1234 0 : PetscFunctionReturn(PETSC_SUCCESS);
1235 : }
1236 :
1237 : static PetscErrorCode
1238 0 : DMView_Moose(DM dm, PetscViewer viewer)
1239 : {
1240 : PetscBool isascii;
1241 : const char *name, *prefix;
1242 0 : DM_Moose * dmm = (DM_Moose *)dm->data;
1243 :
1244 : PetscFunctionBegin;
1245 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1246 0 : if (isascii)
1247 : {
1248 0 : LibmeshPetscCallQ(PetscObjectGetName((PetscObject)dm, &name));
1249 0 : LibmeshPetscCallQ(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
1250 0 : LibmeshPetscCallQ(
1251 : PetscViewerASCIIPrintf(viewer, "DM Moose with name %s and prefix %s\n", name, prefix));
1252 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "variables:"));
1253 0 : for (const auto & vit : *(dmm->_var_ids))
1254 : {
1255 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "(%s,%u) ", vit.first.c_str(), vit.second));
1256 : }
1257 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1258 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "blocks:"));
1259 0 : for (const auto & bit : *(dmm->_block_ids))
1260 : {
1261 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "(%s,%d) ", bit.first.c_str(), bit.second));
1262 : }
1263 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1264 :
1265 0 : if (dmm->_side_ids->size())
1266 : {
1267 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "sides:"));
1268 0 : for (const auto & sit : *(dmm->_side_ids))
1269 : {
1270 0 : LibmeshPetscCallQ(
1271 : PetscViewerASCIIPrintf(viewer, "(%s,%d) ", sit.first.c_str(), sit.second));
1272 : }
1273 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1274 : }
1275 :
1276 0 : if (dmm->_unside_ids->size())
1277 : {
1278 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "unsides:"));
1279 0 : for (const auto & sit : *(dmm->_unside_ids))
1280 : {
1281 0 : LibmeshPetscCallQ(
1282 : PetscViewerASCIIPrintf(viewer, "(%s,%d) ", sit.first.c_str(), sit.second));
1283 : }
1284 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1285 : }
1286 :
1287 0 : if (dmm->_contact_names->size())
1288 : {
1289 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "contacts:"));
1290 0 : for (const auto & cit : *(dmm->_contact_names))
1291 : {
1292 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(
1293 : viewer, "(%s,%s,", cit.second.first.c_str(), cit.second.second.c_str()));
1294 0 : if ((*dmm->_contact_displaced)[cit.second])
1295 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "displaced) "));
1296 : else
1297 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "undisplaced) "));
1298 : }
1299 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1300 : }
1301 :
1302 0 : if (dmm->_uncontact_names->size())
1303 : {
1304 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "_uncontacts:"));
1305 0 : for (const auto & cit : *(dmm->_uncontact_names))
1306 : {
1307 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(
1308 : viewer, "(%s,%s,", cit.second.first.c_str(), cit.second.second.c_str()));
1309 0 : if ((*dmm->_uncontact_displaced)[cit.second])
1310 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "displaced) "));
1311 : else
1312 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "undisplaced) "));
1313 : }
1314 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1315 : }
1316 :
1317 0 : if (dmm->_splitlocs && dmm->_splitlocs->size())
1318 : {
1319 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "Field decomposition:"));
1320 : // FIX: decompositions might have different sizes and components on different ranks.
1321 0 : for (const auto & dit : *(dmm->_splitlocs))
1322 : {
1323 0 : std::string dname = dit.first;
1324 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, " %s", dname.c_str()));
1325 0 : }
1326 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1327 : }
1328 : }
1329 : else
1330 0 : SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Non-ASCII viewers are not supported");
1331 :
1332 0 : PetscFunctionReturn(PETSC_SUCCESS);
1333 : }
1334 :
1335 : static PetscErrorCode
1336 886 : DMMooseGetMeshBlocks_Private(DM dm, std::set<subdomain_id_type> & blocks)
1337 : {
1338 886 : DM_Moose * dmm = (DM_Moose *)(dm->data);
1339 :
1340 : PetscFunctionBegin;
1341 886 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
1342 886 : if (!dmm->_nl)
1343 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1344 :
1345 886 : const MeshBase & mesh = dmm->_system->get_mesh();
1346 : /* The following effectively is a verbatim copy of MeshBase::n_subdomains(). */
1347 : // This requires an inspection on every processor
1348 : libmesh_parallel_only(mesh.comm());
1349 723330 : for (const auto & elem : mesh.active_element_ptr_range())
1350 362108 : blocks.insert(elem->subdomain_id());
1351 : // Some subdomains may only live on other processors
1352 886 : mesh.comm().set_union(blocks);
1353 886 : PetscFunctionReturn(PETSC_SUCCESS);
1354 : }
1355 :
1356 : static PetscErrorCode
1357 443 : DMSetUp_Moose_Pre(DM dm)
1358 : {
1359 443 : DM_Moose * dmm = (DM_Moose *)(dm->data);
1360 :
1361 : PetscFunctionBegin;
1362 443 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
1363 443 : if (!dmm->_nl)
1364 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1365 :
1366 : /* Set up variables, blocks and sides. */
1367 443 : const auto & dofmap = *dmm->_dof_map;
1368 : /* libMesh mesh */
1369 443 : const MeshBase & mesh = dmm->_system->get_mesh();
1370 :
1371 : // Do sides
1372 443 : dmm->_nosides = PETSC_TRUE;
1373 443 : dmm->_side_ids->clear();
1374 443 : dmm->_side_names->clear();
1375 443 : if (dmm->_sides)
1376 : {
1377 40 : dmm->_nosides = PETSC_FALSE;
1378 116 : for (const auto & name : *(dmm->_sides))
1379 : {
1380 76 : boundary_id_type id = dmm->_nl->mesh().getBoundaryID(name);
1381 76 : dmm->_side_names->insert(std::make_pair(id, name));
1382 76 : dmm->_side_ids->insert(std::make_pair(name, id));
1383 : }
1384 40 : delete dmm->_sides;
1385 40 : dmm->_sides = LIBMESH_PETSC_NULLPTR;
1386 : }
1387 :
1388 : // Do unsides
1389 443 : dmm->_nounsides = PETSC_TRUE;
1390 443 : dmm->_unside_ids->clear();
1391 443 : dmm->_unside_names->clear();
1392 443 : if (dmm->_unsides)
1393 : {
1394 28 : dmm->_nounsides = PETSC_FALSE;
1395 84 : for (const auto & name : *(dmm->_unsides))
1396 : {
1397 56 : boundary_id_type id = dmm->_nl->mesh().getBoundaryID(name);
1398 56 : dmm->_unside_names->insert(std::make_pair(id, name));
1399 56 : dmm->_unside_ids->insert(std::make_pair(name, id));
1400 : }
1401 28 : delete dmm->_unsides;
1402 28 : dmm->_unsides = LIBMESH_PETSC_NULLPTR;
1403 : }
1404 :
1405 : // Do unside by var
1406 443 : dmm->_nounside_by_var = PETSC_TRUE;
1407 443 : dmm->_unside_by_var_set->clear();
1408 443 : if (dmm->_unside_by_var)
1409 : {
1410 12 : dmm->_nounside_by_var = PETSC_FALSE;
1411 36 : for (const auto & name : *(dmm->_unside_by_var))
1412 : {
1413 24 : const auto colon_pos = name.find(":");
1414 24 : auto unside_name = name.substr(0, colon_pos);
1415 24 : auto var_name = name.substr(colon_pos + 1);
1416 24 : boundary_id_type id = dmm->_nl->mesh().getBoundaryID(unside_name);
1417 24 : bool var_found = false;
1418 24 : for (unsigned int v = 0; v < dofmap.n_variables(); ++v)
1419 : {
1420 24 : const auto & vname = dofmap.variable(v).name();
1421 24 : if (vname == var_name)
1422 : {
1423 24 : dmm->_unside_by_var_set->insert(std::make_pair(id, v));
1424 24 : var_found = true;
1425 24 : break;
1426 : }
1427 : }
1428 24 : if (!var_found)
1429 0 : mooseError("No variable named '", var_name, "' found");
1430 24 : }
1431 12 : delete dmm->_unside_by_var;
1432 12 : dmm->_unside_by_var = LIBMESH_PETSC_NULLPTR;
1433 : }
1434 :
1435 443 : dmm->_nocontacts = PETSC_TRUE;
1436 :
1437 443 : if (dmm->_contacts)
1438 : {
1439 0 : dmm->_nocontacts = PETSC_FALSE;
1440 0 : for (const auto & cpair : *(dmm->_contacts))
1441 : {
1442 : try
1443 : {
1444 0 : if ((*dmm->_contact_displaced)[cpair])
1445 0 : dmm->_nl->feProblem().getDisplacedProblem()->geomSearchData().getPenetrationLocator(
1446 0 : cpair.first, cpair.second);
1447 : else
1448 0 : dmm->_nl->feProblem().geomSearchData().getPenetrationLocator(cpair.first, cpair.second);
1449 : }
1450 0 : catch (...)
1451 : {
1452 0 : std::ostringstream err;
1453 0 : err << "Problem retrieving contact for PenetrationLocator with primary " << cpair.first
1454 0 : << " and secondary " << cpair.second;
1455 0 : mooseError(err.str());
1456 0 : }
1457 0 : BoundaryID primary_id = dmm->_nl->mesh().getBoundaryID(cpair.first);
1458 0 : BoundaryID secondary_id = dmm->_nl->mesh().getBoundaryID(cpair.second);
1459 0 : DM_Moose::ContactID cid(primary_id, secondary_id);
1460 0 : dmm->_contact_names->insert(std::make_pair(cid, cpair));
1461 : }
1462 : }
1463 :
1464 443 : dmm->_nouncontacts = PETSC_TRUE;
1465 443 : if (dmm->_uncontacts)
1466 : {
1467 0 : dmm->_nouncontacts = PETSC_FALSE;
1468 0 : for (const auto & cpair : *(dmm->_uncontacts))
1469 : {
1470 : try
1471 : {
1472 0 : if ((*dmm->_uncontact_displaced)[cpair])
1473 0 : dmm->_nl->feProblem().getDisplacedProblem()->geomSearchData().getPenetrationLocator(
1474 0 : cpair.first, cpair.second);
1475 : else
1476 0 : dmm->_nl->feProblem().geomSearchData().getPenetrationLocator(cpair.first, cpair.second);
1477 : }
1478 0 : catch (...)
1479 : {
1480 0 : std::ostringstream err;
1481 0 : err << "Problem retrieving uncontact for PenetrationLocator with primary " << cpair.first
1482 0 : << " and secondary " << cpair.second;
1483 0 : mooseError(err.str());
1484 0 : }
1485 0 : BoundaryID primary_id = dmm->_nl->mesh().getBoundaryID(cpair.first);
1486 0 : BoundaryID secondary_id = dmm->_nl->mesh().getBoundaryID(cpair.second);
1487 0 : DM_Moose::ContactID cid(primary_id, secondary_id);
1488 0 : dmm->_uncontact_names->insert(std::make_pair(cid, cpair));
1489 : }
1490 : }
1491 :
1492 443 : dmm->_var_ids->clear();
1493 443 : dmm->_var_names->clear();
1494 : // FIXME: would be nice to invert this nested loop structure so we could iterate over the
1495 : // potentially smaller dmm->vars,
1496 : // but checking against dofmap.variable would still require a linear search, hence, no win. Would
1497 : // be nice to endow dofmap.variable
1498 : // with a fast search capability.
1499 1377 : for (unsigned int v = 0; v < dofmap.n_variables(); ++v)
1500 : {
1501 934 : std::string vname = dofmap.variable(v).name();
1502 934 : if (dmm->_vars && dmm->_vars->size() && dmm->_vars->find(vname) == dmm->_vars->end())
1503 322 : continue;
1504 612 : dmm->_var_ids->insert(std::pair<std::string, unsigned int>(vname, v));
1505 612 : dmm->_var_names->insert(std::pair<unsigned int, std::string>(v, vname));
1506 934 : }
1507 443 : if (dmm->_var_ids->size() == dofmap.n_variables())
1508 141 : dmm->_all_vars = PETSC_TRUE;
1509 : else
1510 302 : dmm->_all_vars = PETSC_FALSE;
1511 443 : if (dmm->_vars)
1512 : {
1513 302 : delete dmm->_vars;
1514 302 : dmm->_vars = LIBMESH_PETSC_NULLPTR;
1515 : }
1516 :
1517 443 : dmm->_block_ids->clear();
1518 443 : dmm->_block_names->clear();
1519 443 : std::set<subdomain_id_type> blocks;
1520 443 : LibmeshPetscCallQ(DMMooseGetMeshBlocks_Private(dm, blocks));
1521 443 : if (blocks.empty())
1522 0 : SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "No mesh blocks found.");
1523 :
1524 1048 : for (const auto & bid : blocks)
1525 : {
1526 605 : std::string bname = mesh.subdomain_name(bid);
1527 605 : if (!bname.length())
1528 : {
1529 : // Block names are currently implemented for Exodus meshes
1530 : // only, so we might have to make up our own block names and
1531 : // maintain our own mapping of block ids to names.
1532 551 : std::ostringstream ss;
1533 551 : ss << bid;
1534 551 : bname = ss.str();
1535 551 : }
1536 605 : if (dmm->_nosides && dmm->_nocontacts)
1537 : {
1538 : // If no sides and no contacts have been specified, by default (null or empty dmm->blocks) all
1539 : // blocks are included in the split Thus, skip this block only if it is explicitly excluded
1540 : // from a nonempty dmm->blocks.
1541 565 : if (dmm->_blocks && dmm->_blocks->size() &&
1542 0 : (dmm->_blocks->find(bname) ==
1543 0 : dmm->_blocks->end() && // We should allow users to use subdomain IDs
1544 565 : dmm->_blocks->find(std::to_string(bid)) == dmm->_blocks->end()))
1545 0 : continue;
1546 : }
1547 : else
1548 : {
1549 : // If sides or contacts have been specified, only the explicitly-specified blocks (those in
1550 : // dmm->blocks, if it's non-null) are in the split. Thus, include this block only if it is
1551 : // explicitly specified in a nonempty dmm->blocks. Equivalently, skip this block if
1552 : // dmm->blocks is dmm->blocks is null or empty or excludes this block.
1553 40 : if (!dmm->_blocks || !dmm->_blocks->size() ||
1554 0 : (dmm->_blocks->find(bname) ==
1555 0 : dmm->_blocks->end() // We should allow users to use subdomain IDs
1556 40 : && dmm->_blocks->find(std::to_string(bid)) == dmm->_blocks->end()))
1557 40 : continue;
1558 : }
1559 565 : dmm->_block_ids->insert(std::make_pair(bname, bid));
1560 565 : dmm->_block_names->insert(std::make_pair(bid, bname));
1561 605 : }
1562 :
1563 443 : if (dmm->_block_ids->size() == blocks.size())
1564 403 : dmm->_all_blocks = PETSC_TRUE;
1565 : else
1566 40 : dmm->_all_blocks = PETSC_FALSE;
1567 443 : if (dmm->_blocks)
1568 : {
1569 0 : delete dmm->_blocks;
1570 0 : dmm->_blocks = LIBMESH_PETSC_NULLPTR;
1571 : }
1572 :
1573 443 : std::string name = dmm->_system->name();
1574 443 : name += "_vars";
1575 1055 : for (const auto & vit : *(dmm->_var_names))
1576 612 : name += "_" + vit.second;
1577 :
1578 443 : name += "_blocks";
1579 :
1580 1008 : for (const auto & bit : *(dmm->_block_names))
1581 565 : name += "_" + bit.second;
1582 :
1583 443 : if (dmm->_side_names && dmm->_side_names->size())
1584 : {
1585 40 : name += "_sides";
1586 116 : for (const auto & sit : *(dmm->_side_names))
1587 76 : name += "_" + sit.second;
1588 : }
1589 443 : if (dmm->_unside_names && dmm->_unside_names->size())
1590 : {
1591 28 : name += "_unsides";
1592 84 : for (const auto & sit : *(dmm->_unside_names))
1593 56 : name += "_" + sit.second;
1594 : }
1595 443 : if (dmm->_contact_names && dmm->_contact_names->size())
1596 : {
1597 0 : name += "_contacts";
1598 0 : for (const auto & cit : *(dmm->_contact_names))
1599 0 : name += "_primary_" + cit.second.first + "_secondary_" + cit.second.second;
1600 : }
1601 443 : if (dmm->_uncontact_names && dmm->_uncontact_names->size())
1602 : {
1603 0 : name += "_uncontacts";
1604 0 : for (const auto & cit : *(dmm->_uncontact_names))
1605 0 : name += "_primary_" + cit.second.first + "_secondary_" + cit.second.second;
1606 : }
1607 443 : LibmeshPetscCallQ(PetscObjectSetName((PetscObject)dm, name.c_str()));
1608 443 : PetscFunctionReturn(PETSC_SUCCESS);
1609 443 : }
1610 :
1611 : PetscErrorCode
1612 0 : DMMooseReset(DM dm)
1613 : {
1614 : PetscBool ismoose;
1615 0 : DM_Moose * dmm = (DM_Moose *)(dm->data);
1616 :
1617 : PetscFunctionBegin;
1618 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose));
1619 0 : if (!ismoose)
1620 0 : PetscFunctionReturn(PETSC_SUCCESS);
1621 0 : if (!dmm->_nl)
1622 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1623 0 : LibmeshPetscCallQ(ISDestroy(&dmm->_embedding));
1624 0 : for (auto & it : *(dmm->_splits))
1625 : {
1626 0 : DM_Moose::SplitInfo & split = it.second;
1627 0 : LibmeshPetscCallQ(ISDestroy(&split._rembedding));
1628 0 : if (split._dm)
1629 0 : LibmeshPetscCallQ(DMMooseReset(split._dm));
1630 : }
1631 0 : dm->setupcalled = PETSC_FALSE;
1632 0 : PetscFunctionReturn(PETSC_SUCCESS);
1633 : }
1634 :
1635 : static PetscErrorCode
1636 443 : DMSetUp_Moose(DM dm)
1637 : {
1638 443 : DM_Moose * dmm = (DM_Moose *)(dm->data);
1639 :
1640 : PetscFunctionBegin;
1641 443 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
1642 443 : if (!dmm->_nl)
1643 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1644 443 : if (dmm->_print_embedding)
1645 : {
1646 : const char *name, *prefix;
1647 : IS embedding;
1648 :
1649 0 : LibmeshPetscCallQ(PetscObjectGetName((PetscObject)dm, &name));
1650 0 : LibmeshPetscCallQ(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
1651 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)dm)->comm),
1652 : "DM Moose with name %s and prefix %s\n",
1653 : name,
1654 : prefix));
1655 0 : if (dmm->_all_vars && dmm->_all_blocks && dmm->_nosides && dmm->_nounsides &&
1656 0 : dmm->_nocontacts && dmm->_nouncontacts)
1657 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)dm)->comm),
1658 0 : "\thas a trivial embedding\n"));
1659 : else
1660 : {
1661 0 : LibmeshPetscCallQ(DMMooseGetEmbedding_Private(dm, &embedding));
1662 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)dm)->comm),
1663 : "\thas embedding defined by IS:\n"));
1664 0 : LibmeshPetscCallQ(ISView(embedding, PETSC_VIEWER_STDOUT_(((PetscObject)dm)->comm)));
1665 0 : LibmeshPetscCallQ(ISDestroy(&embedding));
1666 : }
1667 : }
1668 : /*
1669 : Do not evaluate function, Jacobian or bounds for an embedded DM -- the subproblem might not have
1670 : enough information for that.
1671 : */
1672 443 : if (dmm->_all_vars && dmm->_all_blocks && dmm->_nosides && dmm->_nounsides && dmm->_nocontacts &&
1673 141 : dmm->_nouncontacts)
1674 : {
1675 141 : LibmeshPetscCallQ(DMSNESSetFunction(dm, SNESFunction_DMMoose, (void *)dm));
1676 141 : LibmeshPetscCallQ(DMSNESSetJacobian(dm, SNESJacobian_DMMoose, (void *)dm));
1677 141 : if (dmm->_nl->nonlinearSolver()->bounds || dmm->_nl->nonlinearSolver()->bounds_object)
1678 141 : LibmeshPetscCallQ(DMSetVariableBounds(dm, DMVariableBounds_Moose));
1679 : }
1680 443 : PetscFunctionReturn(PETSC_SUCCESS);
1681 : }
1682 : #if !PETSC_VERSION_LESS_THAN(3, 23, 0)
1683 : PetscErrorCode
1684 443 : DMSetFromOptions_Moose(DM dm, PetscOptionItems /*options*/)
1685 : #elif !PETSC_VERSION_LESS_THAN(3, 18, 0)
1686 : PetscErrorCode
1687 : DMSetFromOptions_Moose(DM dm, PetscOptionItems * /*options*/) // >= 3.18.0
1688 : #elif !PETSC_VERSION_LESS_THAN(3, 7, 0)
1689 : PetscErrorCode
1690 : DMSetFromOptions_Moose(PetscOptionItems * /*options*/, DM dm) // >= 3.7.0
1691 : #else
1692 : PetscErrorCode
1693 : DMSetFromOptions_Moose(PetscOptions * /*options*/, DM dm) // >= 3.6.0
1694 : #endif
1695 : {
1696 443 : DM_Moose * dmm = (DM_Moose *)dm->data;
1697 :
1698 : PetscFunctionBegin;
1699 443 : LibmeshPetscCallQ(DMMooseValidityCheck(dm));
1700 443 : if (!dmm->_nl)
1701 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1702 : // PETSc changed macro definitions in 3.18; the former correct usage
1703 : // is now a compiler error and the new usage is now a compiler
1704 : // warning.
1705 : #if !PETSC_VERSION_LESS_THAN(3, 18, 0)
1706 886 : PetscOptionsBegin(((PetscObject)dm)->comm, ((PetscObject)dm)->prefix, "DMMoose options", "DM");
1707 : #else
1708 : LibmeshPetscCallQ(PetscOptionsBegin(
1709 : ((PetscObject)dm)->comm, ((PetscObject)dm)->prefix, "DMMoose options", "DM"));
1710 : #endif
1711 443 : std::string opt, help;
1712 443 : PetscInt maxvars = dmm->_dof_map->n_variables();
1713 : char ** vars;
1714 443 : std::set<std::string> varset;
1715 443 : PetscInt nvars = maxvars;
1716 443 : LibmeshPetscCallQ(PetscMalloc(maxvars * sizeof(char *), &vars));
1717 443 : opt = "-dm_moose_vars";
1718 443 : help = "Variables in DMMoose";
1719 443 : LibmeshPetscCallQ(PetscOptionsStringArray(
1720 : opt.c_str(), help.c_str(), "DMMooseSetVars", vars, &nvars, LIBMESH_PETSC_NULLPTR));
1721 745 : for (PetscInt i = 0; i < nvars; ++i)
1722 : {
1723 302 : varset.insert(std::string(vars[i]));
1724 302 : LibmeshPetscCallQ(PetscFree(vars[i]));
1725 : }
1726 443 : LibmeshPetscCallQ(PetscFree(vars));
1727 443 : if (varset.size())
1728 302 : LibmeshPetscCallQ(DMMooseSetVariables(dm, varset));
1729 : //
1730 443 : std::set<subdomain_id_type> meshblocks;
1731 443 : LibmeshPetscCallQ(DMMooseGetMeshBlocks_Private(dm, meshblocks));
1732 443 : PetscInt maxblocks = meshblocks.size();
1733 : char ** blocks;
1734 443 : LibmeshPetscCallQ(PetscMalloc(maxblocks * sizeof(char *), &blocks));
1735 443 : std::set<std::string> blockset;
1736 443 : PetscInt nblocks = maxblocks;
1737 443 : opt = "-dm_moose_blocks";
1738 443 : help = "Blocks in DMMoose";
1739 443 : LibmeshPetscCallQ(PetscOptionsStringArray(
1740 : opt.c_str(), help.c_str(), "DMMooseSetBlocks", blocks, &nblocks, LIBMESH_PETSC_NULLPTR));
1741 443 : for (PetscInt i = 0; i < nblocks; ++i)
1742 : {
1743 0 : blockset.insert(std::string(blocks[i]));
1744 0 : LibmeshPetscCallQ(PetscFree(blocks[i]));
1745 : }
1746 443 : LibmeshPetscCallQ(PetscFree(blocks));
1747 443 : if (blockset.size())
1748 0 : LibmeshPetscCallQ(DMMooseSetBlocks(dm, blockset));
1749 443 : PetscInt maxsides = dmm->_system->get_mesh().get_boundary_info().get_global_boundary_ids().size();
1750 : char ** sides;
1751 443 : LibmeshPetscCallQ(PetscMalloc(maxsides * maxvars * sizeof(char *), &sides));
1752 443 : PetscInt nsides = maxsides;
1753 443 : std::set<std::string> sideset;
1754 :
1755 : // Do sides
1756 443 : opt = "-dm_moose_sides";
1757 443 : help = "Sides to include in DMMoose";
1758 443 : LibmeshPetscCallQ(PetscOptionsStringArray(
1759 : opt.c_str(), help.c_str(), "DMMooseSetSides", sides, &nsides, LIBMESH_PETSC_NULLPTR));
1760 519 : for (PetscInt i = 0; i < nsides; ++i)
1761 : {
1762 76 : sideset.insert(std::string(sides[i]));
1763 76 : LibmeshPetscCallQ(PetscFree(sides[i]));
1764 : }
1765 443 : if (sideset.size())
1766 40 : LibmeshPetscCallQ(DMMooseSetSides(dm, sideset));
1767 :
1768 : // Do unsides
1769 443 : opt = "-dm_moose_unsides";
1770 443 : help = "Sides to exclude from DMMoose";
1771 443 : nsides = maxsides;
1772 443 : LibmeshPetscCallQ(PetscOptionsStringArray(
1773 : opt.c_str(), help.c_str(), "DMMooseSetUnSides", sides, &nsides, LIBMESH_PETSC_NULLPTR));
1774 443 : sideset.clear();
1775 499 : for (PetscInt i = 0; i < nsides; ++i)
1776 : {
1777 56 : sideset.insert(std::string(sides[i]));
1778 56 : LibmeshPetscCallQ(PetscFree(sides[i]));
1779 : }
1780 443 : if (sideset.size())
1781 28 : LibmeshPetscCallQ(DMMooseSetUnSides(dm, sideset));
1782 :
1783 : // Do unsides by var
1784 443 : opt = "-dm_moose_unside_by_var";
1785 443 : help = "Sides to exclude from DMMoose on a by-var basis";
1786 443 : nsides = maxsides * maxvars;
1787 443 : LibmeshPetscCallQ(PetscOptionsStringArray(
1788 : opt.c_str(), help.c_str(), "DMMooseSetUnSideByVar", sides, &nsides, LIBMESH_PETSC_NULLPTR));
1789 443 : sideset.clear();
1790 467 : for (PetscInt i = 0; i < nsides; ++i)
1791 : {
1792 24 : sideset.insert(std::string(sides[i]));
1793 24 : LibmeshPetscCallQ(PetscFree(sides[i]));
1794 : }
1795 443 : if (sideset.size())
1796 12 : LibmeshPetscCallQ(DMMooseSetUnSideByVar(dm, sideset));
1797 :
1798 443 : LibmeshPetscCallQ(PetscFree(sides));
1799 443 : PetscInt maxcontacts = dmm->_nl->feProblem().geomSearchData()._penetration_locators.size();
1800 443 : std::shared_ptr<DisplacedProblem> displaced_problem = dmm->_nl->feProblem().getDisplacedProblem();
1801 443 : if (displaced_problem)
1802 0 : maxcontacts = PetscMax(
1803 : maxcontacts, (PetscInt)displaced_problem->geomSearchData()._penetration_locators.size());
1804 :
1805 443 : std::vector<DM_Moose::ContactName> contacts;
1806 443 : std::vector<PetscBool> contact_displaced;
1807 443 : PetscInt ncontacts = 0;
1808 443 : opt = "-dm_moose_ncontacts";
1809 : help =
1810 : "Number of contacts to include in DMMoose. For each <n> < "
1811 : "dm_moose_contacts\n\t-dm_moose_contact_<n> is a comma-separated <primary>,<secondary> pair "
1812 : "defining the contact surfaces"
1813 : "\t-dm_moose_contact_<n>_displaced <bool> determines whether the contact is defined on "
1814 443 : "the displaced mesh or not";
1815 443 : LibmeshPetscCallQ(PetscOptionsInt(opt.c_str(),
1816 : help.c_str(),
1817 : "DMMooseSetContacts",
1818 : ncontacts,
1819 : &ncontacts,
1820 : LIBMESH_PETSC_NULLPTR));
1821 443 : if (ncontacts > maxcontacts)
1822 0 : LIBMESH_SETERRQ2(((PetscObject)dm)->comm,
1823 : PETSC_ERR_ARG_SIZ,
1824 : "Number of requested contacts %" LIBMESH_PETSCINT_FMT
1825 : " exceeds the maximum number of contacts %" LIBMESH_PETSCINT_FMT,
1826 : ncontacts,
1827 : maxcontacts);
1828 443 : for (PetscInt i = 0; i < ncontacts; ++i)
1829 : {
1830 : {
1831 : char * primary_secondary[2];
1832 0 : PetscInt sz = 2;
1833 0 : std::ostringstream oopt, ohelp;
1834 0 : oopt << "-dm_moose_contact_" << i;
1835 0 : ohelp << "Primary and secondary for contact " << i;
1836 0 : LibmeshPetscCallQ(PetscOptionsStringArray(oopt.str().c_str(),
1837 : ohelp.str().c_str(),
1838 : "DMMooseSetContacts",
1839 : primary_secondary,
1840 : &sz,
1841 : LIBMESH_PETSC_NULLPTR));
1842 0 : if (sz != 2)
1843 0 : LIBMESH_SETERRQ2(
1844 : ((PetscObject)dm)->comm,
1845 : PETSC_ERR_ARG_SIZ,
1846 : "Expected 2 sideset IDs (primary & secondary) for contact %" LIBMESH_PETSCINT_FMT
1847 : ", got %" LIBMESH_PETSCINT_FMT " instead",
1848 : i,
1849 : sz);
1850 0 : contacts.push_back(DM_Moose::ContactName(std::string(primary_secondary[0]),
1851 0 : std::string(primary_secondary[1])));
1852 0 : LibmeshPetscCallQ(PetscFree(primary_secondary[0]));
1853 0 : LibmeshPetscCallQ(PetscFree(primary_secondary[1]));
1854 0 : }
1855 : {
1856 0 : PetscBool displaced = PETSC_FALSE;
1857 0 : std::ostringstream oopt, ohelp;
1858 0 : oopt << "-dm_moose_contact_" << i << "_displaced";
1859 0 : ohelp << "Whether contact " << i << " is determined using displaced mesh or not";
1860 0 : LibmeshPetscCallQ(PetscOptionsBool(oopt.str().c_str(),
1861 : ohelp.str().c_str(),
1862 : "DMMooseSetContacts",
1863 : PETSC_FALSE,
1864 : &displaced,
1865 : LIBMESH_PETSC_NULLPTR));
1866 0 : contact_displaced.push_back(displaced);
1867 0 : }
1868 : }
1869 443 : if (contacts.size())
1870 0 : LibmeshPetscCallQ(DMMooseSetContacts(dm, contacts, contact_displaced));
1871 : {
1872 443 : std::ostringstream oopt, ohelp;
1873 : PetscBool is_include_all_nodes;
1874 443 : oopt << "-dm_moose_includeAllContactNodes";
1875 443 : ohelp << "Whether to include all nodes on the contact surfaces into the subsolver";
1876 443 : LibmeshPetscCallQ(PetscOptionsBool(oopt.str().c_str(),
1877 : ohelp.str().c_str(),
1878 : "",
1879 : PETSC_FALSE,
1880 : &is_include_all_nodes,
1881 : LIBMESH_PETSC_NULLPTR));
1882 443 : dmm->_include_all_contact_nodes = is_include_all_nodes;
1883 443 : }
1884 443 : std::vector<DM_Moose::ContactName> uncontacts;
1885 443 : std::vector<PetscBool> uncontact_displaced;
1886 443 : PetscInt nuncontacts = 0;
1887 443 : opt = "-dm_moose_nuncontacts";
1888 : help =
1889 : "Number of contacts to exclude from DMMoose. For each <n> < "
1890 : "dm_moose_contacts\n\t-dm_moose_contact_<n> is a comma-separated <primary>,<secondary> pair "
1891 : "defining the contact surfaces"
1892 : "\t-dm_moose_contact_<n>_displaced <bool> determines whether the contact is defined on "
1893 443 : "the displaced mesh or not";
1894 443 : LibmeshPetscCallQ(PetscOptionsInt(opt.c_str(),
1895 : help.c_str(),
1896 : "DMMooseSetUnContacts",
1897 : nuncontacts,
1898 : &nuncontacts,
1899 : LIBMESH_PETSC_NULLPTR));
1900 443 : if (nuncontacts > maxcontacts)
1901 0 : LIBMESH_SETERRQ2(((PetscObject)dm)->comm,
1902 : PETSC_ERR_ARG_SIZ,
1903 : "Number of requested uncontacts %" LIBMESH_PETSCINT_FMT
1904 : " exceeds the maximum number of contacts %" LIBMESH_PETSCINT_FMT,
1905 : nuncontacts,
1906 : maxcontacts);
1907 443 : for (PetscInt i = 0; i < nuncontacts; ++i)
1908 : {
1909 : {
1910 : char * primary_secondary[2];
1911 0 : PetscInt sz = 2;
1912 0 : std::ostringstream oopt, ohelp;
1913 0 : oopt << "-dm_moose_uncontact_" << i;
1914 0 : ohelp << "Primary and secondary for uncontact " << i;
1915 0 : LibmeshPetscCallQ(PetscOptionsStringArray(oopt.str().c_str(),
1916 : ohelp.str().c_str(),
1917 : "DMMooseSetUnContacts",
1918 : primary_secondary,
1919 : &sz,
1920 : LIBMESH_PETSC_NULLPTR));
1921 0 : if (sz != 2)
1922 0 : LIBMESH_SETERRQ2(
1923 : ((PetscObject)dm)->comm,
1924 : PETSC_ERR_ARG_SIZ,
1925 : "Expected 2 sideset IDs (primary & secondary) for uncontact %" LIBMESH_PETSCINT_FMT
1926 : ", got %" LIBMESH_PETSCINT_FMT " instead",
1927 : i,
1928 : sz);
1929 0 : uncontacts.push_back(DM_Moose::ContactName(std::string(primary_secondary[0]),
1930 0 : std::string(primary_secondary[1])));
1931 0 : LibmeshPetscCallQ(PetscFree(primary_secondary[0]));
1932 0 : LibmeshPetscCallQ(PetscFree(primary_secondary[1]));
1933 0 : }
1934 : {
1935 0 : PetscBool displaced = PETSC_FALSE;
1936 0 : std::ostringstream oopt, ohelp;
1937 0 : oopt << "-dm_moose_uncontact_" << i << "_displaced";
1938 0 : ohelp << "Whether uncontact " << i << " is determined using displaced mesh or not";
1939 0 : LibmeshPetscCallQ(PetscOptionsBool(oopt.str().c_str(),
1940 : ohelp.str().c_str(),
1941 : "DMMooseSetUnContact",
1942 : PETSC_FALSE,
1943 : &displaced,
1944 : LIBMESH_PETSC_NULLPTR));
1945 0 : uncontact_displaced.push_back(displaced);
1946 0 : }
1947 : }
1948 443 : if (uncontacts.size())
1949 0 : LibmeshPetscCallQ(DMMooseSetUnContacts(dm, uncontacts, uncontact_displaced));
1950 :
1951 443 : PetscInt nsplits = 0;
1952 : /* Insert the usage of -dm_moose_fieldsplit_names into this help message, since the following
1953 : * if-clause might never fire, if -help is requested. */
1954 443 : const char * fdhelp = "Number of named fieldsplits defined by the DM.\n\
1955 : \tNames of fieldsplits are defined by -dm_moose_fieldsplit_names <splitname1> <splitname2> ...\n\
1956 : \tEach split can be configured with its own variables, blocks and sides, as any DMMoose";
1957 443 : LibmeshPetscCallQ(PetscOptionsInt(
1958 : "-dm_moose_nfieldsplits", fdhelp, "DMMooseSetSplitNames", nsplits, &nsplits, NULL));
1959 443 : if (nsplits)
1960 : {
1961 161 : PetscInt nnsplits = nsplits;
1962 161 : std::vector<std::string> split_names;
1963 : char ** splitnames;
1964 161 : LibmeshPetscCallQ(PetscMalloc(nsplits * sizeof(char *), &splitnames));
1965 161 : LibmeshPetscCallQ(PetscOptionsStringArray("-dm_moose_fieldsplit_names",
1966 : "Names of fieldsplits defined by the DM",
1967 : "DMMooseSetSplitNames",
1968 : splitnames,
1969 : &nnsplits,
1970 : LIBMESH_PETSC_NULLPTR));
1971 161 : if (!nnsplits)
1972 : {
1973 0 : for (PetscInt i = 0; i < nsplits; ++i)
1974 : {
1975 0 : std::ostringstream s;
1976 0 : s << i;
1977 0 : split_names.push_back(s.str());
1978 0 : }
1979 : }
1980 161 : else if (nsplits != nnsplits)
1981 0 : LIBMESH_SETERRQ2(((PetscObject)dm)->comm,
1982 : PETSC_ERR_ARG_SIZ,
1983 : "Expected %" LIBMESH_PETSCINT_FMT
1984 : " fieldsplit names, got %" LIBMESH_PETSCINT_FMT " instead",
1985 : nsplits,
1986 : nnsplits);
1987 : else
1988 : {
1989 483 : for (PetscInt i = 0; i < nsplits; ++i)
1990 : {
1991 322 : split_names.push_back(std::string(splitnames[i]));
1992 322 : LibmeshPetscCallQ(PetscFree(splitnames[i]));
1993 : }
1994 : }
1995 161 : LibmeshPetscCallQ(PetscFree(splitnames));
1996 161 : LibmeshPetscCallQ(DMMooseSetSplitNames(dm, split_names));
1997 161 : }
1998 443 : LibmeshPetscCallQ(PetscOptionsBool("-dm_moose_print_embedding",
1999 : "Print IS embedding DM's dofs",
2000 : "DMMoose",
2001 : dmm->_print_embedding,
2002 : &dmm->_print_embedding,
2003 : LIBMESH_PETSC_NULLPTR));
2004 443 : PetscOptionsEnd();
2005 443 : LibmeshPetscCallQ(DMSetUp_Moose_Pre(dm)); /* Need some preliminary set up because, strangely
2006 : enough, DMView() is called in DMSetFromOptions(). */
2007 443 : PetscFunctionReturn(PETSC_SUCCESS);
2008 : }
2009 :
2010 : static PetscErrorCode
2011 411 : DMDestroy_Moose(DM dm)
2012 : {
2013 411 : DM_Moose * dmm = (DM_Moose *)(dm->data);
2014 :
2015 : PetscFunctionBegin;
2016 411 : delete dmm->_name;
2017 411 : if (dmm->_vars)
2018 0 : delete dmm->_vars;
2019 411 : delete dmm->_var_ids;
2020 411 : delete dmm->_var_names;
2021 411 : if (dmm->_blocks)
2022 0 : delete dmm->_blocks;
2023 411 : delete dmm->_block_ids;
2024 411 : delete dmm->_block_names;
2025 411 : if (dmm->_sides)
2026 0 : delete dmm->_sides;
2027 411 : delete dmm->_side_ids;
2028 411 : delete dmm->_side_names;
2029 411 : if (dmm->_unsides)
2030 0 : delete dmm->_unsides;
2031 411 : delete dmm->_unside_ids;
2032 411 : delete dmm->_unside_names;
2033 411 : if (dmm->_unside_by_var)
2034 0 : delete dmm->_unside_by_var;
2035 411 : delete dmm->_unside_by_var_set;
2036 411 : if (dmm->_contacts)
2037 0 : delete dmm->_contacts;
2038 411 : delete dmm->_contact_names;
2039 411 : delete dmm->_contact_displaced;
2040 411 : if (dmm->_uncontacts)
2041 0 : delete dmm->_uncontacts;
2042 411 : delete dmm->_uncontact_names;
2043 411 : delete dmm->_uncontact_displaced;
2044 411 : if (dmm->_splits)
2045 : {
2046 701 : for (auto & sit : *(dmm->_splits))
2047 : {
2048 290 : LibmeshPetscCallQ(DMDestroy(&(sit.second._dm)));
2049 290 : LibmeshPetscCallQ(ISDestroy(&(sit.second._rembedding)));
2050 : }
2051 411 : delete dmm->_splits;
2052 : }
2053 411 : if (dmm->_splitlocs)
2054 145 : delete dmm->_splitlocs;
2055 411 : LibmeshPetscCallQ(ISDestroy(&dmm->_embedding));
2056 411 : LibmeshPetscCallQ(PetscFree(dm->data));
2057 411 : PetscFunctionReturn(PETSC_SUCCESS);
2058 : }
2059 :
2060 : PetscErrorCode
2061 443 : DMCreateMoose(MPI_Comm comm,
2062 : NonlinearSystemBase & nl,
2063 : const DofMapBase & dof_map,
2064 : const System & system,
2065 : const std::string & dm_name,
2066 : DM * dm)
2067 : {
2068 : PetscFunctionBegin;
2069 443 : LibmeshPetscCallQ(DMCreate(comm, dm));
2070 443 : LibmeshPetscCallQ(DMSetType(*dm, DMMOOSE));
2071 443 : LibmeshPetscCallQ(DMMooseSetNonlinearSystem(*dm, nl));
2072 443 : LibmeshPetscCallQ(DMMooseSetDofMap(*dm, dof_map));
2073 443 : LibmeshPetscCallQ(DMMooseSetSystem(*dm, system));
2074 443 : LibmeshPetscCallQ(DMMooseSetName(*dm, dm_name));
2075 443 : PetscFunctionReturn(PETSC_SUCCESS);
2076 : }
2077 :
2078 : EXTERN_C_BEGIN
2079 : PetscErrorCode
2080 443 : DMCreate_Moose(DM dm)
2081 : {
2082 : DM_Moose * dmm;
2083 :
2084 : PetscFunctionBegin;
2085 : PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
2086 : #if PETSC_RELEASE_LESS_THAN(3, 18, 0)
2087 : LibmeshPetscCallQ(PetscNewLog(dm, &dmm));
2088 : #else // PetscNewLog was deprecated
2089 443 : LibmeshPetscCallQ(PetscNew(&dmm));
2090 : #endif
2091 443 : dm->data = dmm;
2092 :
2093 443 : dmm->_name = new (std::string);
2094 443 : dmm->_var_ids = new (std::map<std::string, unsigned int>);
2095 443 : dmm->_block_ids = new (std::map<std::string, subdomain_id_type>);
2096 443 : dmm->_var_names = new (std::map<unsigned int, std::string>);
2097 443 : dmm->_block_names = new (std::map<unsigned int, std::string>);
2098 443 : dmm->_side_ids = new (std::map<std::string, BoundaryID>);
2099 443 : dmm->_side_names = new (std::map<BoundaryID, std::string>);
2100 443 : dmm->_unside_ids = new (std::map<std::string, BoundaryID>);
2101 443 : dmm->_unside_names = new (std::map<BoundaryID, std::string>);
2102 443 : dmm->_unside_by_var_set = new (std::set<std::pair<BoundaryID, unsigned int>>);
2103 443 : dmm->_contact_names = new (std::map<DM_Moose::ContactID, DM_Moose::ContactName>);
2104 443 : dmm->_uncontact_names = new (std::map<DM_Moose::ContactID, DM_Moose::ContactName>);
2105 443 : dmm->_contact_displaced = new (std::map<DM_Moose::ContactName, PetscBool>);
2106 443 : dmm->_uncontact_displaced = new (std::map<DM_Moose::ContactName, PetscBool>);
2107 :
2108 443 : dmm->_splits = new (std::map<std::string, DM_Moose::SplitInfo>);
2109 :
2110 443 : dmm->_print_embedding = PETSC_FALSE;
2111 :
2112 443 : dm->ops->createglobalvector = DMCreateGlobalVector_Moose;
2113 443 : dm->ops->createlocalvector = 0; // DMCreateLocalVector_Moose;
2114 443 : dm->ops->getcoloring = 0; // DMGetColoring_Moose;
2115 443 : dm->ops->creatematrix = DMCreateMatrix_Moose;
2116 443 : dm->ops->createinterpolation = 0; // DMCreateInterpolation_Moose;
2117 :
2118 443 : dm->ops->refine = 0; // DMRefine_Moose;
2119 443 : dm->ops->coarsen = 0; // DMCoarsen_Moose;
2120 : #if PETSC_RELEASE_LESS_THAN(3, 12, 0)
2121 : dm->ops->getinjection = 0; // DMGetInjection_Moose;
2122 : dm->ops->getaggregates = 0; // DMGetAggregates_Moose;
2123 : #else
2124 443 : dm->ops->createinjection = 0;
2125 : #endif
2126 :
2127 443 : dm->ops->createfielddecomposition = DMCreateFieldDecomposition_Moose;
2128 443 : dm->ops->createdomaindecomposition = DMCreateDomainDecomposition_Moose;
2129 :
2130 443 : dm->ops->destroy = DMDestroy_Moose;
2131 443 : dm->ops->view = DMView_Moose;
2132 443 : dm->ops->setfromoptions = DMSetFromOptions_Moose;
2133 443 : dm->ops->setup = DMSetUp_Moose;
2134 443 : PetscFunctionReturn(PETSC_SUCCESS);
2135 : }
2136 : EXTERN_C_END
2137 :
2138 : #undef __FUNCT__
2139 : #define __FUNCT__ "SNESUpdateDMMoose"
2140 : PetscErrorCode
2141 0 : SNESUpdateDMMoose(SNES snes, PetscInt iteration)
2142 : {
2143 : /* This is called any time the structure of the problem changes in a way that affects the Jacobian
2144 : sparsity pattern.
2145 : For example, this may happen when NodeFaceConstraints change Jacobian's sparsity pattern based
2146 : on newly-detected Penetration.
2147 : In that case certain preconditioners (e.g., PCASM) will not work, unless we tell them that the
2148 : sparsity pattern has changed.
2149 : For now we are rebuilding the whole KSP, when necessary.
2150 : */
2151 : DM dm;
2152 : KSP ksp;
2153 : const char * prefix;
2154 : MPI_Comm comm;
2155 : PC pc;
2156 :
2157 : PetscFunctionBegin;
2158 0 : if (iteration)
2159 : {
2160 : /* TODO: limit this only to situations when displaced (un)contact splits are present, as is
2161 : * DisplacedProblem(). */
2162 0 : LibmeshPetscCallQ(SNESGetDM(snes, &dm));
2163 0 : LibmeshPetscCallQ(DMMooseReset(dm));
2164 0 : LibmeshPetscCallQ(DMSetUp(dm));
2165 0 : LibmeshPetscCallQ(SNESGetKSP(snes, &ksp));
2166 : /* Should we rebuild the whole KSP? */
2167 0 : LibmeshPetscCallQ(PetscObjectGetOptionsPrefix((PetscObject)ksp, &prefix));
2168 0 : LibmeshPetscCallQ(PetscObjectGetComm((PetscObject)ksp, &comm));
2169 0 : LibmeshPetscCallQ(PCCreate(comm, &pc));
2170 0 : LibmeshPetscCallQ(PCSetDM(pc, dm));
2171 0 : LibmeshPetscCallQ(PCSetOptionsPrefix(pc, prefix));
2172 0 : LibmeshPetscCallQ(PCSetFromOptions(pc));
2173 0 : LibmeshPetscCallQ(KSPSetPC(ksp, pc));
2174 0 : LibmeshPetscCallQ(PCDestroy(&pc));
2175 : }
2176 0 : PetscFunctionReturn(PETSC_SUCCESS);
2177 : }
2178 :
2179 : PetscErrorCode
2180 129 : DMMooseRegisterAll()
2181 : {
2182 : static PetscBool DMMooseRegisterAllCalled = PETSC_FALSE;
2183 :
2184 : PetscFunctionBegin;
2185 129 : if (!DMMooseRegisterAllCalled)
2186 : {
2187 128 : LibmeshPetscCallQ(DMRegister(DMMOOSE, DMCreate_Moose));
2188 128 : DMMooseRegisterAllCalled = PETSC_TRUE;
2189 : }
2190 129 : PetscFunctionReturn(PETSC_SUCCESS);
2191 : }
|
