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