https://mooseframework.inl.gov
PetscDMMoose.C
Go to the documentation of this file.
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 checkSize(const std::string & split_name, const I1 split_size, const I2 size_expected_by_parent)
40 {
41  if (libMesh::cast_int<libMesh::numeric_index_type>(split_size) !=
42  libMesh::cast_int<libMesh::numeric_index_type>(size_expected_by_parent))
43  mooseError("Split '",
44  split_name,
45  "' has size ",
46  libMesh::cast_int<libMesh::numeric_index_type>(split_size),
47  " but the parent split expected size ",
48  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 }
52 
53 struct DM_Moose
54 {
55  NonlinearSystemBase * _nl; // nonlinear system context
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
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;
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 
102  PetscBool _print_embedding;
103 
105  std::string * _name;
106 
110  void checkChildSize(DM child, PetscInt child_size, const std::string & child_name);
111 };
112 
113 void
114 DM_Moose::checkChildSize(DM child, const PetscInt child_size, const std::string & child_name)
115 {
116  for (const auto & split : *_splits)
117  if (split.second._dm == child)
118  {
119  mooseAssert(split.first == child_name, "These should match");
120  PetscInt parent_expected_size;
121  auto ierr = ISGetLocalSize(split.second._rembedding, &parent_expected_size);
122  if (ierr)
123  mooseError("Unable to get size");
124  checkSize(child_name, child_size, parent_expected_size);
125  return;
126  }
127 
128  mooseError("No child DM match");
129 }
130 
131 PetscErrorCode
133 {
134  PetscBool ismoose;
135 
137  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
138  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose));
139  if (!ismoose)
140  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  PetscFunctionReturn(PETSC_SUCCESS);
146 }
147 
148 PetscErrorCode
150  std::vector<std::pair<std::string, std::string>> & contact_names,
151  std::vector<PetscBool> & displaced)
152 {
155  DM_Moose * dmm = (DM_Moose *)dm->data;
156  for (const auto & it : *(dmm->_contact_names))
157  {
158  contact_names.push_back(it.second);
159  displaced.push_back((*dmm->_contact_displaced)[it.second]);
160  }
161  PetscFunctionReturn(PETSC_SUCCESS);
162 }
163 
164 PetscErrorCode
166  std::vector<std::pair<std::string, std::string>> & uncontact_names,
167  std::vector<PetscBool> & displaced)
168 {
171  DM_Moose * dmm = (DM_Moose *)dm->data;
172  for (const auto & it : *(dmm->_uncontact_names))
173  {
174  uncontact_names.push_back(it.second);
175  displaced.push_back((*dmm->_uncontact_displaced)[it.second]);
176  }
177  PetscFunctionReturn(PETSC_SUCCESS);
178 }
179 
180 PetscErrorCode
181 DMMooseGetSides(DM dm, std::vector<std::string> & side_names)
182 {
185  DM_Moose * dmm = (DM_Moose *)dm->data;
186  for (const auto & it : *(dmm->_side_ids))
187  side_names.push_back(it.first);
188  PetscFunctionReturn(PETSC_SUCCESS);
189 }
190 
191 PetscErrorCode
192 DMMooseGetUnSides(DM dm, std::vector<std::string> & side_names)
193 {
196  DM_Moose * dmm = (DM_Moose *)dm->data;
197  for (const auto & it : *(dmm->_unside_ids))
198  side_names.push_back(it.first);
199  PetscFunctionReturn(PETSC_SUCCESS);
200 }
201 
202 PetscErrorCode
203 DMMooseGetBlocks(DM dm, std::vector<std::string> & block_names)
204 {
207  DM_Moose * dmm = (DM_Moose *)dm->data;
208  for (const auto & it : *(dmm->_block_ids))
209  block_names.push_back(it.first);
210  PetscFunctionReturn(PETSC_SUCCESS);
211 }
212 
213 PetscErrorCode
214 DMMooseGetVariables(DM dm, std::vector<std::string> & var_names)
215 {
218  DM_Moose * dmm = (DM_Moose *)(dm->data);
219  for (const auto & it : *(dmm->_var_ids))
220  var_names.push_back(it.first);
221  PetscFunctionReturn(PETSC_SUCCESS);
222 }
223 
224 PetscErrorCode
226 {
229  if (dm->setupcalled)
230  SETERRQ(((PetscObject)dm)->comm,
232  "Cannot reset the NonlinearSystem after DM has been set up.");
233  DM_Moose * dmm = (DM_Moose *)(dm->data);
234  dmm->_nl = &nl;
235  PetscFunctionReturn(PETSC_SUCCESS);
236 }
237 
238 PetscErrorCode
239 DMMooseSetDofMap(DM dm, const DofMapBase & dof_map)
240 {
243  if (dm->setupcalled)
244  SETERRQ(((PetscObject)dm)->comm,
246  "Cannot reset the degree of freedom map after DM has been set up.");
247  DM_Moose * dmm = (DM_Moose *)(dm->data);
248  dmm->_dof_map = &dof_map;
249  PetscFunctionReturn(PETSC_SUCCESS);
250 }
251 
252 PetscErrorCode
253 DMMooseSetSystem(DM dm, const System & system)
254 {
257  if (dm->setupcalled)
258  SETERRQ(((PetscObject)dm)->comm,
260  "Cannot reset the degree of freedom map after DM has been set up.");
261  DM_Moose * dmm = (DM_Moose *)(dm->data);
262  dmm->_system = &system;
263  PetscFunctionReturn(PETSC_SUCCESS);
264 }
265 
266 PetscErrorCode
267 DMMooseSetName(DM dm, const std::string & dm_name)
268 {
271  if (dm->setupcalled)
272  SETERRQ(((PetscObject)dm)->comm,
274  "Cannot reset the MOOSE DM name after DM has been set up.");
275  DM_Moose * dmm = (DM_Moose *)(dm->data);
276  *dmm->_name = dm_name;
277  PetscFunctionReturn(PETSC_SUCCESS);
278 }
279 
280 PetscErrorCode
282 {
285  if (dm->setupcalled)
286  SETERRQ(((PetscObject)dm)->comm,
288  "Cannot reset the parent DM after the child DM has been set up.");
289 
290  DM_Moose * dmm = (DM_Moose *)(dm->data);
291  dmm->_parent = parent;
292  PetscFunctionReturn(PETSC_SUCCESS);
293 }
294 
295 PetscErrorCode
296 DMMooseSetVariables(DM dm, const std::set<std::string> & vars)
297 {
298  DM_Moose * dmm = (DM_Moose *)dm->data;
299 
302  if (dm->setupcalled)
303  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
304  if (dmm->_vars)
305  delete dmm->_vars;
306  std::set<std::string> processed_vars;
307  for (const auto & var_name : vars)
308  {
309  const auto * const var =
310  dmm->_nl->hasVariable(var_name)
311  ? static_cast<MooseVariableBase *>(&dmm->_nl->getVariable(0, var_name))
312  : static_cast<MooseVariableBase *>(&dmm->_nl->getScalarVariable(0, var_name));
313  if (var->isArray())
314  for (const auto i : make_range(var->count()))
315  processed_vars.insert(var->arrayVariableComponent(i));
316  else
317  processed_vars.insert(var_name);
318  }
319 
320  dmm->_vars = new std::set<std::string>(std::move(processed_vars));
321  PetscFunctionReturn(PETSC_SUCCESS);
322 }
323 
324 PetscErrorCode
325 DMMooseSetBlocks(DM dm, const std::set<std::string> & blocks)
326 {
327  DM_Moose * dmm = (DM_Moose *)dm->data;
328 
331  if (dm->setupcalled)
332  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
333  if (dmm->_blocks)
334  delete dmm->_blocks;
335  dmm->_blocks = new std::set<std::string>(blocks);
336  PetscFunctionReturn(PETSC_SUCCESS);
337 }
338 
339 PetscErrorCode
340 DMMooseSetSides(DM dm, const std::set<std::string> & sides)
341 {
342  DM_Moose * dmm = (DM_Moose *)dm->data;
343 
346  if (dm->setupcalled)
347  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
348  if (dmm->_sides)
349  delete dmm->_sides;
350  dmm->_sides = new std::set<std::string>(sides);
351  PetscFunctionReturn(PETSC_SUCCESS);
352 }
353 
354 PetscErrorCode
355 DMMooseSetUnSides(DM dm, const std::set<std::string> & unsides)
356 {
357  DM_Moose * dmm = (DM_Moose *)dm->data;
358 
361  if (dm->setupcalled)
362  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
363  if (dmm->_unsides)
364  delete dmm->_unsides;
365  dmm->_unsides = new std::set<std::string>(unsides);
366  PetscFunctionReturn(PETSC_SUCCESS);
367 }
368 
369 PetscErrorCode
370 DMMooseSetUnSideByVar(DM dm, const std::set<std::string> & unside_by_var)
371 {
372  DM_Moose * dmm = (DM_Moose *)dm->data;
373 
376  if (dm->setupcalled)
377  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
378  if (dmm->_unside_by_var)
379  delete dmm->_unside_by_var;
380  dmm->_unside_by_var = new std::set<std::string>(unside_by_var);
381  PetscFunctionReturn(PETSC_SUCCESS);
382 }
383 
384 PetscErrorCode
386  const std::vector<std::pair<std::string, std::string>> & contacts,
387  const std::vector<PetscBool> & displaced)
388 {
389  DM_Moose * dmm = (DM_Moose *)dm->data;
390 
393  if (dm->setupcalled)
394  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
395  if (contacts.size() != displaced.size())
396  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  if (dmm->_contacts)
403  delete dmm->_contacts;
404  dmm->_contact_displaced->clear();
405  dmm->_contacts = new std::set<DM_Moose::ContactName>();
406  for (unsigned int i = 0; i < contacts.size(); ++i)
407  {
408  dmm->_contacts->insert(contacts[i]);
409  dmm->_contact_displaced->insert(std::make_pair(contacts[i], displaced[i]));
410  }
411  PetscFunctionReturn(PETSC_SUCCESS);
412 }
413 
414 PetscErrorCode
416  const std::vector<std::pair<std::string, std::string>> & uncontacts,
417  const std::vector<PetscBool> & displaced)
418 {
419  DM_Moose * dmm = (DM_Moose *)dm->data;
420 
423  if (dm->setupcalled)
424  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
425  if (uncontacts.size() != displaced.size())
426  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  if (dmm->_uncontacts)
434  delete dmm->_uncontacts;
435  dmm->_uncontact_displaced->clear();
436  dmm->_uncontacts = new std::set<DM_Moose::ContactName>();
437  for (unsigned int i = 0; i < uncontacts.size(); ++i)
438  {
439  dmm->_uncontacts->insert(uncontacts[i]);
440  dmm->_uncontact_displaced->insert(std::make_pair(uncontacts[i], displaced[i]));
441  }
442  PetscFunctionReturn(PETSC_SUCCESS);
443 }
444 
445 PetscErrorCode
447 {
450  DM_Moose * dmm = (DM_Moose *)(dm->data);
451  nl = dmm->_nl;
452  PetscFunctionReturn(PETSC_SUCCESS);
453 }
454 
455 PetscErrorCode
456 DMMooseSetSplitNames(DM dm, const std::vector<std::string> & split_names)
457 {
460  DM_Moose * dmm = (DM_Moose *)(dm->data);
461 
462  if (dmm->_splits)
463  {
464  for (auto & it : *(dmm->_splits))
465  {
466  LibmeshPetscCallQ(DMDestroy(&(it.second._dm)));
467  LibmeshPetscCallQ(ISDestroy(&(it.second._rembedding)));
468  }
469  delete dmm->_splits;
470  dmm->_splits = LIBMESH_PETSC_NULLPTR;
471  }
472  if (dmm->_splitlocs)
473  {
474  delete dmm->_splitlocs;
475  dmm->_splitlocs = LIBMESH_PETSC_NULLPTR;
476  }
477  dmm->_splits = new std::map<std::string, DM_Moose::SplitInfo>();
478  dmm->_splitlocs = new std::multimap<std::string, unsigned int>();
479  for (unsigned int i = 0; i < split_names.size(); ++i)
480  {
482  info._dm = LIBMESH_PETSC_NULLPTR;
483  info._rembedding = LIBMESH_PETSC_NULLPTR;
484  std::string name = split_names[i];
485  (*dmm->_splits)[name] = info;
486  dmm->_splitlocs->insert(std::make_pair(name, i));
487  }
488  PetscFunctionReturn(PETSC_SUCCESS);
489 }
490 
491 PetscErrorCode
492 DMMooseGetSplitNames(DM dm, std::vector<std::string> & split_names)
493 {
496  DM_Moose * dmm = (DM_Moose *)(dm->data);
497  if (!dm->setupcalled)
498  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "DM not set up");
499  split_names.clear();
500  split_names.reserve(dmm->_splitlocs->size());
501  if (dmm->_splitlocs && dmm->_splitlocs->size())
502  for (const auto & lit : *(dmm->_splitlocs))
503  {
504  std::string sname = lit.first;
505  unsigned int sloc = lit.second;
506  split_names[sloc] = sname;
507  }
508  PetscFunctionReturn(PETSC_SUCCESS);
509 }
510 
511 static PetscErrorCode
512 DMMooseGetEmbedding_Private(DM dm, IS * embedding)
513 {
514  DM_Moose * dmm = (DM_Moose *)dm->data;
515 
517  if (!embedding)
518  PetscFunctionReturn(PETSC_SUCCESS);
519  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  if (!dmm->_all_vars || !dmm->_all_blocks || !dmm->_nosides || !dmm->_nounsides ||
545  !dmm->_nounside_by_var || !dmm->_nocontacts || !dmm->_nouncontacts)
546  {
547  auto & dofmap = *dmm->_dof_map;
548  // Put this outside the lambda scope to avoid constant memory reallocation
549  std::vector<dof_id_type> node_indices;
550  auto process_nodal_dof_indices =
551  [&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  dofmap.dof_indices(&node, node_indices, var_num);
557  for (const auto index : node_indices)
558  {
559  if (index >= dofmap.first_dof() && index < dofmap.end_dof())
560  local_indices.insert(index);
561  else if (nonlocal_indices)
562  nonlocal_indices->insert(index);
563  }
564  };
565 
566  auto process_elem_dof_indices =
567  [&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  for (const auto index : elem_indices)
572  {
573  if (index >= dofmap.first_dof() && index < dofmap.end_dof())
574  local_indices.insert(index);
575  else if (nonlocal_indices)
576  nonlocal_indices->insert(index);
577  }
578  };
579 
580  std::set<dof_id_type> indices;
581  std::set<dof_id_type> unindices;
582  std::set<dof_id_type> cached_indices;
583  std::set<dof_id_type> cached_unindices;
584  auto & lm_mesh = dmm->_system->get_mesh();
585  const auto & node_to_elem_map = dmm->_nl->feProblem().mesh().nodeToElemMap();
586  for (const auto & vit : *(dmm->_var_ids))
587  {
588  unsigned int v = vit.second;
589  // Iterate only over this DM's blocks.
590  if (!dmm->_all_blocks || (dmm->_nosides && dmm->_nocontacts))
591  for (const auto & bit : *(dmm->_block_ids))
592  {
593  subdomain_id_type b = bit.second;
594  for (const auto & elem : as_range(lm_mesh.active_local_subdomain_elements_begin(b),
595  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  std::vector<dof_id_type> evindices;
599  dofmap.dof_indices(elem, evindices, v);
600  process_elem_dof_indices(evindices, indices);
601  }
602 
603  // Sometime, we own nodes but do not own the elements the nodes are connected to
604  {
605  bool is_on_current_block = false;
606  for (auto & node : lm_mesh.local_node_ptr_range())
607  {
608  const unsigned int n_comp = node->n_comp(dmm->_system->number(), v);
609 
610  // skip it if no dof
611  if (!n_comp)
612  continue;
613 
614  auto node_to_elem_pair = node_to_elem_map.find(node->id());
615  is_on_current_block = false;
616  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  const Elem & neighbor_elem = lm_mesh.elem_ref(elem_num);
621  if (neighbor_elem.subdomain_id() == b)
622  {
623  is_on_current_block = true;
624  break;
625  }
626  }
627  // we add indices for the current block only
628  if (!is_on_current_block)
629  continue;
630 
631  process_nodal_dof_indices(*node, v, indices);
632  }
633  }
634  }
635 
636  // Iterate over the sides from this split.
637  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  ConstBndNodeRange & bnodes = *dmm->_nl->mesh().getBoundaryNodeRange();
652  for (const auto & bnode : bnodes)
653  {
654  BoundaryID boundary_id = bnode->_bnd_id;
655  if (dmm->_side_names->find(boundary_id) == dmm->_side_names->end())
656  continue;
657 
658  const Node * node = bnode->_node;
659  process_nodal_dof_indices(*node, v, indices);
660  }
661  }
662 
663  // Iterate over the sides excluded from this split.
664  if (dmm->_unside_ids->size())
665  {
666  ConstBndNodeRange & bnodes = *dmm->_nl->mesh().getBoundaryNodeRange();
667  for (const auto & bnode : bnodes)
668  {
669  BoundaryID boundary_id = bnode->_bnd_id;
670  if (dmm->_unside_names->find(boundary_id) == dmm->_unside_names->end())
671  continue;
672  const Node * node = bnode->_node;
673  process_nodal_dof_indices(*node, v, unindices);
674  }
675  }
676  if (dmm->_unside_by_var_set->size())
677  {
678  std::set<BoundaryID> eligible_bids;
679  for (const auto & [bid, var] : *(dmm->_unside_by_var_set))
680  if (var == v)
681  eligible_bids.insert(bid);
682 
683  ConstBndNodeRange & bnodes = *dmm->_nl->mesh().getBoundaryNodeRange();
684  for (const auto & bnode : bnodes)
685  {
686  BoundaryID boundary_id = bnode->_bnd_id;
687  if (eligible_bids.count(boundary_id))
688  {
689  const Node * node = bnode->_node;
690  process_nodal_dof_indices(*node, v, unindices);
691  }
692  }
693  }
694 
695  auto process_contact_all_nodes =
696  [dmm, process_nodal_dof_indices, v](const auto & contact_names,
697  auto & indices_to_insert_to)
698  {
699  std::set<boundary_id_type> bc_id_set;
700  // loop over contacts
701  for (const auto & [contact_bid_pair, contact_bname_pair] : contact_names)
702  {
703  libmesh_ignore(contact_bname_pair);
704  bc_id_set.insert(contact_bid_pair.first); // primary
705  bc_id_set.insert(contact_bid_pair.second); // secondary
706  }
707  // loop over boundary elements
709  for (const auto & belem : range)
710  {
711  const Elem * elem_bdry = belem->_elem;
712  const auto side = belem->_side;
713  BoundaryID boundary_id = belem->_bnd_id;
714 
715  if (bc_id_set.find(boundary_id) == bc_id_set.end())
716  continue;
717 
718  for (const auto node_idx : elem_bdry->node_index_range())
719  if (elem_bdry->is_node_on_side(node_idx, side))
720  process_nodal_dof_indices(elem_bdry->node_ref(node_idx), v, indices_to_insert_to);
721  }
722  };
723 
724  auto process_contact_some_nodes =
725  [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  std::vector<dof_id_type> evindices;
731  for (const auto & it : contact_names)
732  {
733  PetscBool displaced = (*dmm->_uncontact_displaced)[it.second];
734  PenetrationLocator * locator;
735  if (displaced)
736  {
737  std::shared_ptr<DisplacedProblem> displaced_problem =
739  if (!displaced_problem)
740  {
741  std::ostringstream err;
742  err << "Cannot use a displaced uncontact (" << it.second.first << ","
743  << it.second.second << ") with an undisplaced problem";
744  mooseError(err.str());
745  }
746  locator = displaced_problem->geomSearchData()._penetration_locators[it.first];
747  }
748  else
749  locator = dmm->_nl->feProblem().geomSearchData()._penetration_locators[it.first];
750 
751  evindices.clear();
752  // penetration locator
753  auto lend = locator->_penetration_info.end();
754  for (auto lit = locator->_penetration_info.begin(); lit != lend; ++lit)
755  {
756  const dof_id_type secondary_node_num = lit->first;
757  PenetrationInfo * pinfo = lit->second;
758  if (pinfo && pinfo->isCaptured())
759  {
760  const Node & secondary_node = lm_mesh.node_ref(secondary_node_num);
761  process_nodal_dof_indices(
762  secondary_node, v, indices_to_insert_to, &nonlocal_indices_to_insert_to);
763 
764  // indices for primary element
765  evindices.clear();
766  const Elem * primary_side = pinfo->_side;
767  dofmap.dof_indices(primary_side, evindices, v);
768  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  };
774 
775  // Include all nodes on the contact surfaces
776  if (dmm->_contact_names->size() && dmm->_include_all_contact_nodes)
777  process_contact_all_nodes(*dmm->_contact_names, indices);
778 
779  // Iterate over the contacts included in this split.
780  if (dmm->_contact_names->size() && !(dmm->_include_all_contact_nodes))
781  process_contact_some_nodes(*dmm->_contact_names, indices, cached_indices);
782 
783  // Exclude all nodes on the contact surfaces
784  if (dmm->_uncontact_names->size() && dmm->_include_all_contact_nodes)
785  process_contact_all_nodes(*dmm->_uncontact_names, unindices);
786 
787  // Iterate over the contacts excluded from this split.
788  if (dmm->_uncontact_names->size() && !(dmm->_include_all_contact_nodes))
789  process_contact_some_nodes(*dmm->_uncontact_names, unindices, cached_unindices);
790  } // variables
791 
792  std::vector<dof_id_type> local_vec_indices(cached_indices.size());
793  std::copy(cached_indices.begin(), cached_indices.end(), local_vec_indices.begin());
794  if (dmm->_contact_names->size() && !(dmm->_include_all_contact_nodes))
795  dmm->_nl->feProblem().mesh().comm().allgather(local_vec_indices, false);
796  // insert indices
797  for (const auto & dof : local_vec_indices)
798  if (dof >= dofmap.first_dof() && dof < dofmap.end_dof())
799  indices.insert(dof);
800 
801  local_vec_indices.clear();
802  local_vec_indices.resize(cached_unindices.size());
803  std::copy(cached_unindices.begin(), cached_unindices.end(), local_vec_indices.begin());
804  if (dmm->_uncontact_names->size() && !(dmm->_include_all_contact_nodes))
805  dmm->_nl->feProblem().mesh().comm().allgather(local_vec_indices, false);
806  // insert unindices
807  for (const auto & dof : local_vec_indices)
808  if (dof >= dofmap.first_dof() && dof < dofmap.end_dof())
809  unindices.insert(dof);
810 
811  std::set<dof_id_type> dindices;
812  std::set_difference(indices.begin(),
813  indices.end(),
814  unindices.begin(),
815  unindices.end(),
816  std::inserter(dindices, dindices.end()));
817  PetscInt * darray;
818  LibmeshPetscCallQ(PetscMalloc(sizeof(PetscInt) * dindices.size(), &darray));
819  dof_id_type i = 0;
820  for (const auto & dof : dindices)
821  {
822  darray[i] = dof;
823  ++i;
824  }
825  LibmeshPetscCallQ(ISCreateGeneral(
826  ((PetscObject)dm)->comm, dindices.size(), darray, PETSC_OWN_POINTER, &dmm->_embedding));
827  }
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  LibmeshPetscCallQ(DMCreateGlobalVector(dm, &v));
837  LibmeshPetscCallQ(VecGetOwnershipRange(v, &low, &high));
839  ISCreateStride(((PetscObject)dm)->comm, (high - low), low, 1, &dmm->_embedding));
840  }
841  }
842  LibmeshPetscCallQ(PetscObjectReference((PetscObject)(dmm->_embedding)));
843  *embedding = dmm->_embedding;
844 
845  PetscFunctionReturn(PETSC_SUCCESS);
846 }
847 
848 static PetscErrorCode
850  DM dm, PetscInt * len, char *** namelist, IS ** islist, DM ** dmlist)
851 {
852  DM_Moose * dmm = (DM_Moose *)(dm->data);
853 
855 
856  PetscInt split_size_sum = 0;
857 
858  /* Only called after DMSetUp(). */
859  if (!dmm->_splitlocs)
860  PetscFunctionReturn(PETSC_SUCCESS);
861  *len = dmm->_splitlocs->size();
862  if (namelist)
863  LibmeshPetscCallQ(PetscMalloc(*len * sizeof(char *), namelist));
864  if (islist)
865  LibmeshPetscCallQ(PetscMalloc(*len * sizeof(IS), islist));
866  if (dmlist)
867  LibmeshPetscCallQ(PetscMalloc(*len * sizeof(DM), dmlist));
868  for (const auto & dit : *(dmm->_splitlocs))
869  {
870  unsigned int d = dit.second;
871  std::string dname = dit.first;
872  DM_Moose::SplitInfo & dinfo = (*dmm->_splits)[dname];
873  if (!dinfo._dm)
874  {
876  ((PetscObject)dm)->comm, *dmm->_nl, *dmm->_dof_map, *dmm->_system, dname, &dinfo._dm));
878  PetscObjectSetOptionsPrefix((PetscObject)dinfo._dm, ((PetscObject)dm)->prefix));
879  std::string suffix = std::string("fieldsplit_") + dname + "_";
880  LibmeshPetscCallQ(PetscObjectAppendOptionsPrefix((PetscObject)dinfo._dm, suffix.c_str()));
882  }
883  LibmeshPetscCallQ(DMSetFromOptions(dinfo._dm));
884  LibmeshPetscCallQ(DMSetUp(dinfo._dm));
885  if (namelist)
886  LibmeshPetscCallQ(PetscStrallocpy(dname.c_str(), (*namelist) + d));
887  if (islist)
888  {
889  if (!dinfo._rembedding)
890  {
891  IS dembedding, lembedding;
893  if (dmm->_embedding)
894  {
895  // Create a relative embedding into the parent's index space.
896  LibmeshPetscCallQ(ISEmbed(dembedding, dmm->_embedding, PETSC_TRUE, &lembedding));
897  const PetscInt * lindices;
898  PetscInt len, dlen, llen, *rindices, off, i;
899  LibmeshPetscCallQ(ISGetLocalSize(dembedding, &dlen));
900  LibmeshPetscCallQ(ISGetLocalSize(lembedding, &llen));
901  if (llen != dlen)
902  LIBMESH_SETERRQ1(
903  ((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Failed to embed split %u", d);
904  LibmeshPetscCallQ(ISDestroy(&dembedding));
905  // Convert local embedding to global (but still relative) embedding
906  LibmeshPetscCallQ(PetscMalloc(llen * sizeof(PetscInt), &rindices));
907  LibmeshPetscCallQ(ISGetIndices(lembedding, &lindices));
908  LibmeshPetscCallQ(PetscMemcpy(rindices, lindices, llen * sizeof(PetscInt)));
909  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  LibmeshPetscCallQ(ISGetLocalSize(dmm->_embedding, &len));
913 
914  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  off -= len;
926  for (i = 0; i < llen; ++i)
927  rindices[i] += off;
928  LibmeshPetscCallQ(ISCreateGeneral(
929  ((PetscObject)dm)->comm, llen, rindices, PETSC_OWN_POINTER, &(dinfo._rembedding)));
930  }
931  else
932  {
933  dinfo._rembedding = dembedding;
934  }
935  }
936  LibmeshPetscCallQ(PetscObjectReference((PetscObject)(dinfo._rembedding)));
937  (*islist)[d] = dinfo._rembedding;
938  PetscInt is_size;
939  LibmeshPetscCallQ(ISGetLocalSize(dinfo._rembedding, &is_size));
940  split_size_sum += is_size;
941  }
942  if (dmlist)
943  {
944  LibmeshPetscCallQ(PetscObjectReference((PetscObject)dinfo._dm));
945  (*dmlist)[d] = dinfo._dm;
946  }
947  }
948 
949  mooseAssert(islist, "What does it even mean if this is NULL?");
950 
951  if (dmm->_parent)
952  dmm->_parent->checkChildSize(dm, split_size_sum, *dmm->_name);
953  else
954  checkSize(*dmm->_name, split_size_sum, dmm->_dof_map->n_local_dofs());
955 
956  PetscFunctionReturn(PETSC_SUCCESS);
957 }
958 
959 static PetscErrorCode
961  DM dm, PetscInt * len, char *** namelist, IS ** innerislist, IS ** outerislist, DM ** dmlist)
962 {
964  /* Use DMCreateFieldDecomposition_Moose() to obtain everything but outerislist, which is currently
965  * LIBMESH_PETSC_NULLPTR. */
966  if (outerislist)
967  *outerislist = LIBMESH_PETSC_NULLPTR; /* FIX: allow mesh-based overlap. */
968  LibmeshPetscCallQ(DMCreateFieldDecomposition_Moose(dm, len, namelist, innerislist, dmlist));
969  PetscFunctionReturn(PETSC_SUCCESS);
970 }
971 
972 static PetscErrorCode
973 DMMooseFunction(DM dm, Vec x, Vec r)
974 {
976  libmesh_assert(x);
977  libmesh_assert(r);
978 
979  NonlinearSystemBase * nl = NULL;
981  PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(nl->system().solution.get());
982  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  X_global.swap(X_sys);
989  nl->system().update();
990  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.
997  nl->system().current_local_solution.get());
998 
999  // Zero the residual vector before assembling
1000  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  if (nl->nonlinearSolver()->residual && nl->nonlinearSolver()->residual_object)
1005  {
1006  std::ostringstream err;
1007  err << "ERROR: cannot specifiy both a function and object to compute the Residual!"
1008  << std::endl;
1009  mooseError(err.str());
1010  }
1011  if (nl->nonlinearSolver()->matvec && nl->nonlinearSolver()->residual_and_jacobian_object)
1012  {
1013  std::ostringstream err;
1014  err << "ERROR: cannot specifiy both a function and object to compute the combined Residual & "
1015  "Jacobian!"
1016  << std::endl;
1017  mooseError(err.str());
1018  }
1019  if (nl->nonlinearSolver()->residual != NULL)
1020  nl->nonlinearSolver()->residual(
1021  *(nl->system().current_local_solution.get()), R, nl->nonlinearSolver()->system());
1022  else if (nl->nonlinearSolver()->residual_object != NULL)
1023  nl->nonlinearSolver()->residual_object->residual(
1024  *(nl->system().current_local_solution.get()), R, nl->nonlinearSolver()->system());
1025  else if (nl->nonlinearSolver()->matvec != NULL)
1026  nl->nonlinearSolver()->matvec(
1027  *(nl->system().current_local_solution.get()), &R, NULL, nl->nonlinearSolver()->system());
1028  else if (nl->nonlinearSolver()->residual_and_jacobian_object != NULL)
1029  nl->nonlinearSolver()->residual_and_jacobian_object->residual_and_jacobian(
1030  *(nl->system().current_local_solution.get()), &R, NULL, nl->nonlinearSolver()->system());
1031  else
1032  {
1033  std::ostringstream err;
1034  err << "No suitable residual computation routine found";
1035  mooseError(err.str());
1036  }
1037  R.close();
1038  PetscFunctionReturn(PETSC_SUCCESS);
1039 }
1040 
1041 static PetscErrorCode
1042 SNESFunction_DMMoose(SNES, Vec x, Vec r, void * ctx)
1043 {
1044  DM dm = (DM)ctx;
1045 
1048  PetscFunctionReturn(PETSC_SUCCESS);
1049 }
1050 
1051 static PetscErrorCode
1052 DMMooseJacobian(DM dm, Vec x, Mat jac, Mat pc)
1053 {
1054  NonlinearSystemBase * nl = NULL;
1055 
1058 
1059  PetscMatrix<Number> the_pc(pc, nl->comm());
1060  PetscMatrix<Number> Jac(jac, nl->comm());
1061  PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(nl->system().solution.get());
1062  PetscVector<Number> X_global(x, nl->comm());
1063 
1064  // Set the dof maps
1065  the_pc.attach_dof_map(nl->system().get_dof_map());
1066  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  X_global.swap(X_sys);
1073  nl->system().update();
1074  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.
1081  nl->system().current_local_solution.get());
1082 
1083  // Zero out the preconditioner before computing the Jacobian.
1084  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  if (nl->nonlinearSolver()->jacobian && nl->nonlinearSolver()->jacobian_object)
1089  {
1090  std::ostringstream err;
1091  err << "ERROR: cannot specifiy both a function and object to compute the Jacobian!"
1092  << std::endl;
1093  mooseError(err.str());
1094  }
1095  if (nl->nonlinearSolver()->matvec && nl->nonlinearSolver()->residual_and_jacobian_object)
1096  {
1097  std::ostringstream err;
1098  err << "ERROR: cannot specifiy both a function and object to compute the combined Residual & "
1099  "Jacobian!"
1100  << std::endl;
1101  mooseError(err.str());
1102  }
1103  if (nl->nonlinearSolver()->jacobian != NULL)
1104  nl->nonlinearSolver()->jacobian(
1105  *(nl->system().current_local_solution.get()), the_pc, nl->nonlinearSolver()->system());
1106  else if (nl->nonlinearSolver()->jacobian_object != NULL)
1107  nl->nonlinearSolver()->jacobian_object->jacobian(
1108  *(nl->system().current_local_solution.get()), the_pc, nl->nonlinearSolver()->system());
1109  else if (nl->nonlinearSolver()->matvec != NULL)
1110  nl->nonlinearSolver()->matvec(*(nl->system().current_local_solution.get()),
1111  NULL,
1112  &the_pc,
1113  nl->nonlinearSolver()->system());
1114  else if (nl->nonlinearSolver()->residual_and_jacobian_object != NULL)
1115  nl->nonlinearSolver()->residual_and_jacobian_object->residual_and_jacobian(
1116  *(nl->system().current_local_solution.get()),
1117  NULL,
1118  &the_pc,
1119  nl->nonlinearSolver()->system());
1120  else
1121  {
1122  std::ostringstream err;
1123  err << "No suitable Jacobian routine or object";
1124  mooseError(err.str());
1125  }
1126  the_pc.close();
1127  Jac.close();
1128  PetscFunctionReturn(PETSC_SUCCESS);
1129 }
1130 
1131 static PetscErrorCode
1132 SNESJacobian_DMMoose(SNES, Vec x, Mat jac, Mat pc, void * ctx)
1133 {
1134  DM dm = (DM)ctx;
1135 
1137  LibmeshPetscCallQ(DMMooseJacobian(dm, x, jac, pc));
1138  PetscFunctionReturn(PETSC_SUCCESS);
1139 }
1140 
1141 static PetscErrorCode
1142 DMVariableBounds_Moose(DM dm, Vec xl, Vec xu)
1143 {
1144  NonlinearSystemBase * nl = NULL;
1145 
1148 
1149  PetscVector<Number> XL(xl, nl->comm());
1150  PetscVector<Number> XU(xu, nl->comm());
1151 
1152  LibmeshPetscCallQ(VecSet(xl, PETSC_NINFINITY));
1153  LibmeshPetscCallQ(VecSet(xu, PETSC_INFINITY));
1154  if (nl->nonlinearSolver()->bounds != NULL)
1155  nl->nonlinearSolver()->bounds(XL, XU, nl->nonlinearSolver()->system());
1156  else if (nl->nonlinearSolver()->bounds_object != NULL)
1157  nl->nonlinearSolver()->bounds_object->bounds(XL, XU, nl->nonlinearSolver()->system());
1158  else
1159  SETERRQ(
1160  ((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "No bounds calculation in this Moose object");
1161  PetscFunctionReturn(PETSC_SUCCESS);
1162 }
1163 
1164 static PetscErrorCode
1166 {
1167  DM_Moose * dmm = (DM_Moose *)(dm->data);
1168 
1171  if (!dmm->_nl)
1172  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1173 
1174  NumericVector<Number> * nv = (dmm->_system->solution).get();
1175  PetscVector<Number> * pv = dynamic_cast<PetscVector<Number> *>(nv);
1176  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  if (dmm->_embedding)
1181  {
1182  PetscInt n;
1183  LibmeshPetscCallQ(VecCreate(((PetscObject)v)->comm, x));
1184  LibmeshPetscCallQ(ISGetLocalSize(dmm->_embedding, &n));
1185  LibmeshPetscCallQ(VecSetSizes(*x, n, PETSC_DETERMINE));
1186  LibmeshPetscCallQ(VecSetType(*x, ((PetscObject)v)->type_name));
1187  LibmeshPetscCallQ(VecSetFromOptions(*x));
1188  LibmeshPetscCallQ(VecSetUp(*x));
1189  }
1190  else
1191  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  LibmeshPetscCallQ(VecSetDM(*x, dm));
1197 #endif
1198  PetscFunctionReturn(PETSC_SUCCESS);
1199 }
1200 
1201 static PetscErrorCode
1203 {
1204  DM_Moose * dmm = (DM_Moose *)(dm->data);
1205  MatType type;
1206 
1209  if (!dmm->_nl)
1210  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1211  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  const auto & dof_map = *dmm->_dof_map;
1221  PetscInt M, N, m, n;
1222  MPI_Comm comm;
1223  M = dof_map.n_dofs();
1224  N = M;
1225  m = static_cast<PetscInt>(dof_map.n_local_dofs());
1226  n = m;
1227  LibmeshPetscCallQ(PetscObjectGetComm((PetscObject)dm, &comm));
1228  LibmeshPetscCallQ(MatCreate(comm, A));
1229  LibmeshPetscCallQ(MatSetSizes(*A, m, n, M, N));
1230  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  LibmeshPetscCallQ(MatSetUp(*A));
1234  PetscFunctionReturn(PETSC_SUCCESS);
1235 }
1236 
1237 static PetscErrorCode
1238 DMView_Moose(DM dm, PetscViewer viewer)
1239 {
1240  PetscBool isascii;
1241  const char *name, *prefix;
1242  DM_Moose * dmm = (DM_Moose *)dm->data;
1243 
1245  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1246  if (isascii)
1247  {
1248  LibmeshPetscCallQ(PetscObjectGetName((PetscObject)dm, &name));
1249  LibmeshPetscCallQ(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
1251  PetscViewerASCIIPrintf(viewer, "DM Moose with name %s and prefix %s\n", name, prefix));
1252  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "variables:"));
1253  for (const auto & vit : *(dmm->_var_ids))
1254  {
1255  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "(%s,%u) ", vit.first.c_str(), vit.second));
1256  }
1257  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1258  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "blocks:"));
1259  for (const auto & bit : *(dmm->_block_ids))
1260  {
1261  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "(%s,%d) ", bit.first.c_str(), bit.second));
1262  }
1263  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1264 
1265  if (dmm->_side_ids->size())
1266  {
1267  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "sides:"));
1268  for (const auto & sit : *(dmm->_side_ids))
1269  {
1271  PetscViewerASCIIPrintf(viewer, "(%s,%d) ", sit.first.c_str(), sit.second));
1272  }
1273  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1274  }
1275 
1276  if (dmm->_unside_ids->size())
1277  {
1278  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "unsides:"));
1279  for (const auto & sit : *(dmm->_unside_ids))
1280  {
1282  PetscViewerASCIIPrintf(viewer, "(%s,%d) ", sit.first.c_str(), sit.second));
1283  }
1284  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1285  }
1286 
1287  if (dmm->_contact_names->size())
1288  {
1289  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "contacts:"));
1290  for (const auto & cit : *(dmm->_contact_names))
1291  {
1292  LibmeshPetscCallQ(PetscViewerASCIIPrintf(
1293  viewer, "(%s,%s,", cit.second.first.c_str(), cit.second.second.c_str()));
1294  if ((*dmm->_contact_displaced)[cit.second])
1295  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "displaced) "));
1296  else
1297  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "undisplaced) "));
1298  }
1299  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1300  }
1301 
1302  if (dmm->_uncontact_names->size())
1303  {
1304  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "_uncontacts:"));
1305  for (const auto & cit : *(dmm->_uncontact_names))
1306  {
1307  LibmeshPetscCallQ(PetscViewerASCIIPrintf(
1308  viewer, "(%s,%s,", cit.second.first.c_str(), cit.second.second.c_str()));
1309  if ((*dmm->_uncontact_displaced)[cit.second])
1310  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "displaced) "));
1311  else
1312  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "undisplaced) "));
1313  }
1314  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1315  }
1316 
1317  if (dmm->_splitlocs && dmm->_splitlocs->size())
1318  {
1319  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "Field decomposition:"));
1320  // FIX: decompositions might have different sizes and components on different ranks.
1321  for (const auto & dit : *(dmm->_splitlocs))
1322  {
1323  std::string dname = dit.first;
1324  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, " %s", dname.c_str()));
1325  }
1326  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
1327  }
1328  }
1329  else
1330  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Non-ASCII viewers are not supported");
1331 
1332  PetscFunctionReturn(PETSC_SUCCESS);
1333 }
1334 
1335 static PetscErrorCode
1336 DMMooseGetMeshBlocks_Private(DM dm, std::set<subdomain_id_type> & blocks)
1337 {
1338  DM_Moose * dmm = (DM_Moose *)(dm->data);
1339 
1342  if (!dmm->_nl)
1343  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1344 
1345  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  for (const auto & elem : mesh.active_element_ptr_range())
1350  blocks.insert(elem->subdomain_id());
1351  // Some subdomains may only live on other processors
1353  PetscFunctionReturn(PETSC_SUCCESS);
1354 }
1355 
1356 static PetscErrorCode
1358 {
1359  DM_Moose * dmm = (DM_Moose *)(dm->data);
1360 
1363  if (!dmm->_nl)
1364  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1365 
1366  /* Set up variables, blocks and sides. */
1367  const auto & dofmap = *dmm->_dof_map;
1368  /* libMesh mesh */
1369  const MeshBase & mesh = dmm->_system->get_mesh();
1370 
1371  // Do sides
1372  dmm->_nosides = PETSC_TRUE;
1373  dmm->_side_ids->clear();
1374  dmm->_side_names->clear();
1375  if (dmm->_sides)
1376  {
1377  dmm->_nosides = PETSC_FALSE;
1378  for (const auto & name : *(dmm->_sides))
1379  {
1380  boundary_id_type id = dmm->_nl->mesh().getBoundaryID(name);
1381  dmm->_side_names->insert(std::make_pair(id, name));
1382  dmm->_side_ids->insert(std::make_pair(name, id));
1383  }
1384  delete dmm->_sides;
1385  dmm->_sides = LIBMESH_PETSC_NULLPTR;
1386  }
1387 
1388  // Do unsides
1389  dmm->_nounsides = PETSC_TRUE;
1390  dmm->_unside_ids->clear();
1391  dmm->_unside_names->clear();
1392  if (dmm->_unsides)
1393  {
1394  dmm->_nounsides = PETSC_FALSE;
1395  for (const auto & name : *(dmm->_unsides))
1396  {
1397  boundary_id_type id = dmm->_nl->mesh().getBoundaryID(name);
1398  dmm->_unside_names->insert(std::make_pair(id, name));
1399  dmm->_unside_ids->insert(std::make_pair(name, id));
1400  }
1401  delete dmm->_unsides;
1402  dmm->_unsides = LIBMESH_PETSC_NULLPTR;
1403  }
1404 
1405  // Do unside by var
1406  dmm->_nounside_by_var = PETSC_TRUE;
1407  dmm->_unside_by_var_set->clear();
1408  if (dmm->_unside_by_var)
1409  {
1410  dmm->_nounside_by_var = PETSC_FALSE;
1411  for (const auto & name : *(dmm->_unside_by_var))
1412  {
1413  const auto colon_pos = name.find(":");
1414  auto unside_name = name.substr(0, colon_pos);
1415  auto var_name = name.substr(colon_pos + 1);
1416  boundary_id_type id = dmm->_nl->mesh().getBoundaryID(unside_name);
1417  bool var_found = false;
1418  for (unsigned int v = 0; v < dofmap.n_variables(); ++v)
1419  {
1420  const auto & vname = dofmap.variable(v).name();
1421  if (vname == var_name)
1422  {
1423  dmm->_unside_by_var_set->insert(std::make_pair(id, v));
1424  var_found = true;
1425  break;
1426  }
1427  }
1428  if (!var_found)
1429  mooseError("No variable named '", var_name, "' found");
1430  }
1431  delete dmm->_unside_by_var;
1432  dmm->_unside_by_var = LIBMESH_PETSC_NULLPTR;
1433  }
1434 
1435  dmm->_nocontacts = PETSC_TRUE;
1436 
1437  if (dmm->_contacts)
1438  {
1439  dmm->_nocontacts = PETSC_FALSE;
1440  for (const auto & cpair : *(dmm->_contacts))
1441  {
1442  try
1443  {
1444  if ((*dmm->_contact_displaced)[cpair])
1445  dmm->_nl->feProblem().getDisplacedProblem()->geomSearchData().getPenetrationLocator(
1446  cpair.first, cpair.second);
1447  else
1448  dmm->_nl->feProblem().geomSearchData().getPenetrationLocator(cpair.first, cpair.second);
1449  }
1450  catch (...)
1451  {
1452  std::ostringstream err;
1453  err << "Problem retrieving contact for PenetrationLocator with primary " << cpair.first
1454  << " and secondary " << cpair.second;
1455  mooseError(err.str());
1456  }
1457  BoundaryID primary_id = dmm->_nl->mesh().getBoundaryID(cpair.first);
1458  BoundaryID secondary_id = dmm->_nl->mesh().getBoundaryID(cpair.second);
1459  DM_Moose::ContactID cid(primary_id, secondary_id);
1460  dmm->_contact_names->insert(std::make_pair(cid, cpair));
1461  }
1462  }
1463 
1464  dmm->_nouncontacts = PETSC_TRUE;
1465  if (dmm->_uncontacts)
1466  {
1467  dmm->_nouncontacts = PETSC_FALSE;
1468  for (const auto & cpair : *(dmm->_uncontacts))
1469  {
1470  try
1471  {
1472  if ((*dmm->_uncontact_displaced)[cpair])
1473  dmm->_nl->feProblem().getDisplacedProblem()->geomSearchData().getPenetrationLocator(
1474  cpair.first, cpair.second);
1475  else
1476  dmm->_nl->feProblem().geomSearchData().getPenetrationLocator(cpair.first, cpair.second);
1477  }
1478  catch (...)
1479  {
1480  std::ostringstream err;
1481  err << "Problem retrieving uncontact for PenetrationLocator with primary " << cpair.first
1482  << " and secondary " << cpair.second;
1483  mooseError(err.str());
1484  }
1485  BoundaryID primary_id = dmm->_nl->mesh().getBoundaryID(cpair.first);
1486  BoundaryID secondary_id = dmm->_nl->mesh().getBoundaryID(cpair.second);
1487  DM_Moose::ContactID cid(primary_id, secondary_id);
1488  dmm->_uncontact_names->insert(std::make_pair(cid, cpair));
1489  }
1490  }
1491 
1492  dmm->_var_ids->clear();
1493  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  for (unsigned int v = 0; v < dofmap.n_variables(); ++v)
1500  {
1501  std::string vname = dofmap.variable(v).name();
1502  if (dmm->_vars && dmm->_vars->size() && dmm->_vars->find(vname) == dmm->_vars->end())
1503  continue;
1504  dmm->_var_ids->insert(std::pair<std::string, unsigned int>(vname, v));
1505  dmm->_var_names->insert(std::pair<unsigned int, std::string>(v, vname));
1506  }
1507  if (dmm->_var_ids->size() == dofmap.n_variables())
1508  dmm->_all_vars = PETSC_TRUE;
1509  else
1510  dmm->_all_vars = PETSC_FALSE;
1511  if (dmm->_vars)
1512  {
1513  delete dmm->_vars;
1514  dmm->_vars = LIBMESH_PETSC_NULLPTR;
1515  }
1516 
1517  dmm->_block_ids->clear();
1518  dmm->_block_names->clear();
1519  std::set<subdomain_id_type> blocks;
1521  if (blocks.empty())
1522  SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "No mesh blocks found.");
1523 
1524  for (const auto & bid : blocks)
1525  {
1526  std::string bname = mesh.subdomain_name(bid);
1527  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  std::ostringstream ss;
1533  ss << bid;
1534  bname = ss.str();
1535  }
1536  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  if (dmm->_blocks && dmm->_blocks->size() &&
1542  (dmm->_blocks->find(bname) ==
1543  dmm->_blocks->end() && // We should allow users to use subdomain IDs
1544  dmm->_blocks->find(std::to_string(bid)) == dmm->_blocks->end()))
1545  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  if (!dmm->_blocks || !dmm->_blocks->size() ||
1554  (dmm->_blocks->find(bname) ==
1555  dmm->_blocks->end() // We should allow users to use subdomain IDs
1556  && dmm->_blocks->find(std::to_string(bid)) == dmm->_blocks->end()))
1557  continue;
1558  }
1559  dmm->_block_ids->insert(std::make_pair(bname, bid));
1560  dmm->_block_names->insert(std::make_pair(bid, bname));
1561  }
1562 
1563  if (dmm->_block_ids->size() == blocks.size())
1564  dmm->_all_blocks = PETSC_TRUE;
1565  else
1566  dmm->_all_blocks = PETSC_FALSE;
1567  if (dmm->_blocks)
1568  {
1569  delete dmm->_blocks;
1570  dmm->_blocks = LIBMESH_PETSC_NULLPTR;
1571  }
1572 
1573  std::string name = dmm->_system->name();
1574  name += "_vars";
1575  for (const auto & vit : *(dmm->_var_names))
1576  name += "_" + vit.second;
1577 
1578  name += "_blocks";
1579 
1580  for (const auto & bit : *(dmm->_block_names))
1581  name += "_" + bit.second;
1582 
1583  if (dmm->_side_names && dmm->_side_names->size())
1584  {
1585  name += "_sides";
1586  for (const auto & sit : *(dmm->_side_names))
1587  name += "_" + sit.second;
1588  }
1589  if (dmm->_unside_names && dmm->_unside_names->size())
1590  {
1591  name += "_unsides";
1592  for (const auto & sit : *(dmm->_unside_names))
1593  name += "_" + sit.second;
1594  }
1595  if (dmm->_contact_names && dmm->_contact_names->size())
1596  {
1597  name += "_contacts";
1598  for (const auto & cit : *(dmm->_contact_names))
1599  name += "_primary_" + cit.second.first + "_secondary_" + cit.second.second;
1600  }
1601  if (dmm->_uncontact_names && dmm->_uncontact_names->size())
1602  {
1603  name += "_uncontacts";
1604  for (const auto & cit : *(dmm->_uncontact_names))
1605  name += "_primary_" + cit.second.first + "_secondary_" + cit.second.second;
1606  }
1607  LibmeshPetscCallQ(PetscObjectSetName((PetscObject)dm, name.c_str()));
1608  PetscFunctionReturn(PETSC_SUCCESS);
1609 }
1610 
1611 PetscErrorCode
1613 {
1614  PetscBool ismoose;
1615  DM_Moose * dmm = (DM_Moose *)(dm->data);
1616 
1618  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose));
1619  if (!ismoose)
1620  PetscFunctionReturn(PETSC_SUCCESS);
1621  if (!dmm->_nl)
1622  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1623  LibmeshPetscCallQ(ISDestroy(&dmm->_embedding));
1624  for (auto & it : *(dmm->_splits))
1625  {
1626  DM_Moose::SplitInfo & split = it.second;
1627  LibmeshPetscCallQ(ISDestroy(&split._rembedding));
1628  if (split._dm)
1630  }
1631  dm->setupcalled = PETSC_FALSE;
1632  PetscFunctionReturn(PETSC_SUCCESS);
1633 }
1634 
1635 static PetscErrorCode
1637 {
1638  DM_Moose * dmm = (DM_Moose *)(dm->data);
1639 
1642  if (!dmm->_nl)
1643  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1644  if (dmm->_print_embedding)
1645  {
1646  const char *name, *prefix;
1647  IS embedding;
1648 
1649  LibmeshPetscCallQ(PetscObjectGetName((PetscObject)dm, &name));
1650  LibmeshPetscCallQ(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
1651  LibmeshPetscCallQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)dm)->comm),
1652  "DM Moose with name %s and prefix %s\n",
1653  name,
1654  prefix));
1655  if (dmm->_all_vars && dmm->_all_blocks && dmm->_nosides && dmm->_nounsides &&
1656  dmm->_nocontacts && dmm->_nouncontacts)
1657  LibmeshPetscCallQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)dm)->comm),
1658  "\thas a trivial embedding\n"));
1659  else
1660  {
1662  LibmeshPetscCallQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)dm)->comm),
1663  "\thas embedding defined by IS:\n"));
1664  LibmeshPetscCallQ(ISView(embedding, PETSC_VIEWER_STDOUT_(((PetscObject)dm)->comm)));
1665  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  if (dmm->_all_vars && dmm->_all_blocks && dmm->_nosides && dmm->_nounsides && dmm->_nocontacts &&
1673  dmm->_nouncontacts)
1674  {
1675  LibmeshPetscCallQ(DMSNESSetFunction(dm, SNESFunction_DMMoose, (void *)dm));
1676  LibmeshPetscCallQ(DMSNESSetJacobian(dm, SNESJacobian_DMMoose, (void *)dm));
1677  if (dmm->_nl->nonlinearSolver()->bounds || dmm->_nl->nonlinearSolver()->bounds_object)
1678  LibmeshPetscCallQ(DMSetVariableBounds(dm, DMVariableBounds_Moose));
1679  }
1680  PetscFunctionReturn(PETSC_SUCCESS);
1681 }
1682 #if !PETSC_VERSION_LESS_THAN(3, 23, 0)
1683 PetscErrorCode
1684 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
1696  DM_Moose * dmm = (DM_Moose *)dm->data;
1697 
1700  if (!dmm->_nl)
1701  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  PetscOptionsBegin(((PetscObject)dm)->comm, ((PetscObject)dm)->prefix, "DMMoose options", "DM");
1707 #else
1709  ((PetscObject)dm)->comm, ((PetscObject)dm)->prefix, "DMMoose options", "DM"));
1710 #endif
1711  std::string opt, help;
1712  PetscInt maxvars = dmm->_dof_map->n_variables();
1713  char ** vars;
1714  std::set<std::string> varset;
1715  PetscInt nvars = maxvars;
1716  LibmeshPetscCallQ(PetscMalloc(maxvars * sizeof(char *), &vars));
1717  opt = "-dm_moose_vars";
1718  help = "Variables in DMMoose";
1719  LibmeshPetscCallQ(PetscOptionsStringArray(
1720  opt.c_str(), help.c_str(), "DMMooseSetVars", vars, &nvars, LIBMESH_PETSC_NULLPTR));
1721  for (PetscInt i = 0; i < nvars; ++i)
1722  {
1723  varset.insert(std::string(vars[i]));
1724  LibmeshPetscCallQ(PetscFree(vars[i]));
1725  }
1726  LibmeshPetscCallQ(PetscFree(vars));
1727  if (varset.size())
1729  //
1730  std::set<subdomain_id_type> meshblocks;
1732  PetscInt maxblocks = meshblocks.size();
1733  char ** blocks;
1734  LibmeshPetscCallQ(PetscMalloc(maxblocks * sizeof(char *), &blocks));
1735  std::set<std::string> blockset;
1736  PetscInt nblocks = maxblocks;
1737  opt = "-dm_moose_blocks";
1738  help = "Blocks in DMMoose";
1739  LibmeshPetscCallQ(PetscOptionsStringArray(
1740  opt.c_str(), help.c_str(), "DMMooseSetBlocks", blocks, &nblocks, LIBMESH_PETSC_NULLPTR));
1741  for (PetscInt i = 0; i < nblocks; ++i)
1742  {
1743  blockset.insert(std::string(blocks[i]));
1744  LibmeshPetscCallQ(PetscFree(blocks[i]));
1745  }
1746  LibmeshPetscCallQ(PetscFree(blocks));
1747  if (blockset.size())
1750  char ** sides;
1751  LibmeshPetscCallQ(PetscMalloc(maxsides * maxvars * sizeof(char *), &sides));
1752  PetscInt nsides = maxsides;
1753  std::set<std::string> sideset;
1754 
1755  // Do sides
1756  opt = "-dm_moose_sides";
1757  help = "Sides to include in DMMoose";
1758  LibmeshPetscCallQ(PetscOptionsStringArray(
1759  opt.c_str(), help.c_str(), "DMMooseSetSides", sides, &nsides, LIBMESH_PETSC_NULLPTR));
1760  for (PetscInt i = 0; i < nsides; ++i)
1761  {
1762  sideset.insert(std::string(sides[i]));
1763  LibmeshPetscCallQ(PetscFree(sides[i]));
1764  }
1765  if (sideset.size())
1767 
1768  // Do unsides
1769  opt = "-dm_moose_unsides";
1770  help = "Sides to exclude from DMMoose";
1771  nsides = maxsides;
1772  LibmeshPetscCallQ(PetscOptionsStringArray(
1773  opt.c_str(), help.c_str(), "DMMooseSetUnSides", sides, &nsides, LIBMESH_PETSC_NULLPTR));
1774  sideset.clear();
1775  for (PetscInt i = 0; i < nsides; ++i)
1776  {
1777  sideset.insert(std::string(sides[i]));
1778  LibmeshPetscCallQ(PetscFree(sides[i]));
1779  }
1780  if (sideset.size())
1782 
1783  // Do unsides by var
1784  opt = "-dm_moose_unside_by_var";
1785  help = "Sides to exclude from DMMoose on a by-var basis";
1786  nsides = maxsides * maxvars;
1787  LibmeshPetscCallQ(PetscOptionsStringArray(
1788  opt.c_str(), help.c_str(), "DMMooseSetUnSideByVar", sides, &nsides, LIBMESH_PETSC_NULLPTR));
1789  sideset.clear();
1790  for (PetscInt i = 0; i < nsides; ++i)
1791  {
1792  sideset.insert(std::string(sides[i]));
1793  LibmeshPetscCallQ(PetscFree(sides[i]));
1794  }
1795  if (sideset.size())
1797 
1798  LibmeshPetscCallQ(PetscFree(sides));
1800  std::shared_ptr<DisplacedProblem> displaced_problem = dmm->_nl->feProblem().getDisplacedProblem();
1802  maxcontacts = PetscMax(
1803  maxcontacts, (PetscInt)displaced_problem->geomSearchData()._penetration_locators.size());
1804 
1805  std::vector<DM_Moose::ContactName> contacts;
1806  std::vector<PetscBool> contact_displaced;
1807  PetscInt ncontacts = 0;
1808  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  "the displaced mesh or not";
1815  LibmeshPetscCallQ(PetscOptionsInt(opt.c_str(),
1816  help.c_str(),
1817  "DMMooseSetContacts",
1818  ncontacts,
1819  &ncontacts,
1820  LIBMESH_PETSC_NULLPTR));
1821  if (ncontacts > maxcontacts)
1822  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  for (PetscInt i = 0; i < ncontacts; ++i)
1829  {
1830  {
1831  char * primary_secondary[2];
1832  PetscInt sz = 2;
1833  std::ostringstream oopt, ohelp;
1834  oopt << "-dm_moose_contact_" << i;
1835  ohelp << "Primary and secondary for contact " << i;
1836  LibmeshPetscCallQ(PetscOptionsStringArray(oopt.str().c_str(),
1837  ohelp.str().c_str(),
1838  "DMMooseSetContacts",
1839  primary_secondary,
1840  &sz,
1841  LIBMESH_PETSC_NULLPTR));
1842  if (sz != 2)
1843  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  contacts.push_back(DM_Moose::ContactName(std::string(primary_secondary[0]),
1851  std::string(primary_secondary[1])));
1852  LibmeshPetscCallQ(PetscFree(primary_secondary[0]));
1853  LibmeshPetscCallQ(PetscFree(primary_secondary[1]));
1854  }
1855  {
1856  PetscBool displaced = PETSC_FALSE;
1857  std::ostringstream oopt, ohelp;
1858  oopt << "-dm_moose_contact_" << i << "_displaced";
1859  ohelp << "Whether contact " << i << " is determined using displaced mesh or not";
1860  LibmeshPetscCallQ(PetscOptionsBool(oopt.str().c_str(),
1861  ohelp.str().c_str(),
1862  "DMMooseSetContacts",
1863  PETSC_FALSE,
1864  &displaced,
1865  LIBMESH_PETSC_NULLPTR));
1866  contact_displaced.push_back(displaced);
1867  }
1868  }
1869  if (contacts.size())
1871  {
1872  std::ostringstream oopt, ohelp;
1873  PetscBool is_include_all_nodes;
1874  oopt << "-dm_moose_includeAllContactNodes";
1875  ohelp << "Whether to include all nodes on the contact surfaces into the subsolver";
1876  LibmeshPetscCallQ(PetscOptionsBool(oopt.str().c_str(),
1877  ohelp.str().c_str(),
1878  "",
1879  PETSC_FALSE,
1881  LIBMESH_PETSC_NULLPTR));
1883  }
1884  std::vector<DM_Moose::ContactName> uncontacts;
1885  std::vector<PetscBool> uncontact_displaced;
1886  PetscInt nuncontacts = 0;
1887  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  "the displaced mesh or not";
1894  LibmeshPetscCallQ(PetscOptionsInt(opt.c_str(),
1895  help.c_str(),
1896  "DMMooseSetUnContacts",
1897  nuncontacts,
1898  &nuncontacts,
1899  LIBMESH_PETSC_NULLPTR));
1901  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  for (PetscInt i = 0; i < nuncontacts; ++i)
1908  {
1909  {
1910  char * primary_secondary[2];
1911  PetscInt sz = 2;
1912  std::ostringstream oopt, ohelp;
1913  oopt << "-dm_moose_uncontact_" << i;
1914  ohelp << "Primary and secondary for uncontact " << i;
1915  LibmeshPetscCallQ(PetscOptionsStringArray(oopt.str().c_str(),
1916  ohelp.str().c_str(),
1917  "DMMooseSetUnContacts",
1918  primary_secondary,
1919  &sz,
1920  LIBMESH_PETSC_NULLPTR));
1921  if (sz != 2)
1922  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  uncontacts.push_back(DM_Moose::ContactName(std::string(primary_secondary[0]),
1930  std::string(primary_secondary[1])));
1931  LibmeshPetscCallQ(PetscFree(primary_secondary[0]));
1932  LibmeshPetscCallQ(PetscFree(primary_secondary[1]));
1933  }
1934  {
1935  PetscBool displaced = PETSC_FALSE;
1936  std::ostringstream oopt, ohelp;
1937  oopt << "-dm_moose_uncontact_" << i << "_displaced";
1938  ohelp << "Whether uncontact " << i << " is determined using displaced mesh or not";
1939  LibmeshPetscCallQ(PetscOptionsBool(oopt.str().c_str(),
1940  ohelp.str().c_str(),
1941  "DMMooseSetUnContact",
1942  PETSC_FALSE,
1943  &displaced,
1944  LIBMESH_PETSC_NULLPTR));
1945  uncontact_displaced.push_back(displaced);
1946  }
1947  }
1948  if (uncontacts.size())
1950 
1951  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  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  LibmeshPetscCallQ(PetscOptionsInt(
1958  "-dm_moose_nfieldsplits", fdhelp, "DMMooseSetSplitNames", nsplits, &nsplits, NULL));
1959  if (nsplits)
1960  {
1961  PetscInt nnsplits = nsplits;
1962  std::vector<std::string> split_names;
1963  char ** splitnames;
1964  LibmeshPetscCallQ(PetscMalloc(nsplits * sizeof(char *), &splitnames));
1965  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  if (!nnsplits)
1972  {
1973  for (PetscInt i = 0; i < nsplits; ++i)
1974  {
1975  std::ostringstream s;
1976  s << i;
1977  split_names.push_back(s.str());
1978  }
1979  }
1980  else if (nsplits != nnsplits)
1981  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  for (PetscInt i = 0; i < nsplits; ++i)
1990  {
1991  split_names.push_back(std::string(splitnames[i]));
1992  LibmeshPetscCallQ(PetscFree(splitnames[i]));
1993  }
1994  }
1995  LibmeshPetscCallQ(PetscFree(splitnames));
1996  LibmeshPetscCallQ(DMMooseSetSplitNames(dm, split_names));
1997  }
1998  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  PetscOptionsEnd();
2005  LibmeshPetscCallQ(DMSetUp_Moose_Pre(dm)); /* Need some preliminary set up because, strangely
2006  enough, DMView() is called in DMSetFromOptions(). */
2007  PetscFunctionReturn(PETSC_SUCCESS);
2008 }
2009 
2010 static PetscErrorCode
2012 {
2013  DM_Moose * dmm = (DM_Moose *)(dm->data);
2014 
2016  delete dmm->_name;
2017  if (dmm->_vars)
2018  delete dmm->_vars;
2019  delete dmm->_var_ids;
2020  delete dmm->_var_names;
2021  if (dmm->_blocks)
2022  delete dmm->_blocks;
2023  delete dmm->_block_ids;
2024  delete dmm->_block_names;
2025  if (dmm->_sides)
2026  delete dmm->_sides;
2027  delete dmm->_side_ids;
2028  delete dmm->_side_names;
2029  if (dmm->_unsides)
2030  delete dmm->_unsides;
2031  delete dmm->_unside_ids;
2032  delete dmm->_unside_names;
2033  if (dmm->_unside_by_var)
2034  delete dmm->_unside_by_var;
2035  delete dmm->_unside_by_var_set;
2036  if (dmm->_contacts)
2037  delete dmm->_contacts;
2038  delete dmm->_contact_names;
2039  delete dmm->_contact_displaced;
2040  if (dmm->_uncontacts)
2041  delete dmm->_uncontacts;
2042  delete dmm->_uncontact_names;
2043  delete dmm->_uncontact_displaced;
2044  if (dmm->_splits)
2045  {
2046  for (auto & sit : *(dmm->_splits))
2047  {
2048  LibmeshPetscCallQ(DMDestroy(&(sit.second._dm)));
2049  LibmeshPetscCallQ(ISDestroy(&(sit.second._rembedding)));
2050  }
2051  delete dmm->_splits;
2052  }
2053  if (dmm->_splitlocs)
2054  delete dmm->_splitlocs;
2055  LibmeshPetscCallQ(ISDestroy(&dmm->_embedding));
2056  LibmeshPetscCallQ(PetscFree(dm->data));
2057  PetscFunctionReturn(PETSC_SUCCESS);
2058 }
2059 
2060 PetscErrorCode
2061 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 {
2069  LibmeshPetscCallQ(DMCreate(comm, dm));
2070  LibmeshPetscCallQ(DMSetType(*dm, DMMOOSE));
2074  LibmeshPetscCallQ(DMMooseSetName(*dm, dm_name));
2075  PetscFunctionReturn(PETSC_SUCCESS);
2076 }
2077 
2078 EXTERN_C_BEGIN
2079 PetscErrorCode
2081 {
2082  DM_Moose * dmm;
2083 
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  LibmeshPetscCallQ(PetscNew(&dmm));
2090 #endif
2091  dm->data = dmm;
2092 
2093  dmm->_name = new (std::string);
2094  dmm->_var_ids = new (std::map<std::string, unsigned int>);
2095  dmm->_block_ids = new (std::map<std::string, subdomain_id_type>);
2096  dmm->_var_names = new (std::map<unsigned int, std::string>);
2097  dmm->_block_names = new (std::map<unsigned int, std::string>);
2098  dmm->_side_ids = new (std::map<std::string, BoundaryID>);
2099  dmm->_side_names = new (std::map<BoundaryID, std::string>);
2100  dmm->_unside_ids = new (std::map<std::string, BoundaryID>);
2101  dmm->_unside_names = new (std::map<BoundaryID, std::string>);
2102  dmm->_unside_by_var_set = new (std::set<std::pair<BoundaryID, unsigned int>>);
2103  dmm->_contact_names = new (std::map<DM_Moose::ContactID, DM_Moose::ContactName>);
2104  dmm->_uncontact_names = new (std::map<DM_Moose::ContactID, DM_Moose::ContactName>);
2105  dmm->_contact_displaced = new (std::map<DM_Moose::ContactName, PetscBool>);
2106  dmm->_uncontact_displaced = new (std::map<DM_Moose::ContactName, PetscBool>);
2107 
2108  dmm->_splits = new (std::map<std::string, DM_Moose::SplitInfo>);
2109 
2110  dmm->_print_embedding = PETSC_FALSE;
2111 
2112  dm->ops->createglobalvector = DMCreateGlobalVector_Moose;
2113  dm->ops->createlocalvector = 0; // DMCreateLocalVector_Moose;
2114  dm->ops->getcoloring = 0; // DMGetColoring_Moose;
2115  dm->ops->creatematrix = DMCreateMatrix_Moose;
2116  dm->ops->createinterpolation = 0; // DMCreateInterpolation_Moose;
2117 
2118  dm->ops->refine = 0; // DMRefine_Moose;
2119  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  dm->ops->createinjection = 0;
2125 #endif
2126 
2127  dm->ops->createfielddecomposition = DMCreateFieldDecomposition_Moose;
2128  dm->ops->createdomaindecomposition = DMCreateDomainDecomposition_Moose;
2129 
2130  dm->ops->destroy = DMDestroy_Moose;
2131  dm->ops->view = DMView_Moose;
2132  dm->ops->setfromoptions = DMSetFromOptions_Moose;
2133  dm->ops->setup = DMSetUp_Moose;
2134  PetscFunctionReturn(PETSC_SUCCESS);
2135 }
2136 EXTERN_C_END
2137 
2138 #undef __FUNCT__
2139 #define __FUNCT__ "SNESUpdateDMMoose"
2140 PetscErrorCode
2141 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 
2158  if (iteration)
2159  {
2160  /* TODO: limit this only to situations when displaced (un)contact splits are present, as is
2161  * DisplacedProblem(). */
2162  LibmeshPetscCallQ(SNESGetDM(snes, &dm));
2164  LibmeshPetscCallQ(DMSetUp(dm));
2165  LibmeshPetscCallQ(SNESGetKSP(snes, &ksp));
2166  /* Should we rebuild the whole KSP? */
2167  LibmeshPetscCallQ(PetscObjectGetOptionsPrefix((PetscObject)ksp, &prefix));
2168  LibmeshPetscCallQ(PetscObjectGetComm((PetscObject)ksp, &comm));
2169  LibmeshPetscCallQ(PCCreate(comm, &pc));
2170  LibmeshPetscCallQ(PCSetDM(pc, dm));
2171  LibmeshPetscCallQ(PCSetOptionsPrefix(pc, prefix));
2172  LibmeshPetscCallQ(PCSetFromOptions(pc));
2173  LibmeshPetscCallQ(KSPSetPC(ksp, pc));
2174  LibmeshPetscCallQ(PCDestroy(&pc));
2175  }
2176  PetscFunctionReturn(PETSC_SUCCESS);
2177 }
2178 
2179 PetscErrorCode
2181 {
2182  static PetscBool DMMooseRegisterAllCalled = PETSC_FALSE;
2183 
2185  if (!DMMooseRegisterAllCalled)
2186  {
2187  LibmeshPetscCallQ(DMRegister(DMMOOSE, DMCreate_Moose));
2188  DMMooseRegisterAllCalled = PETSC_TRUE;
2189  }
2190  PetscFunctionReturn(PETSC_SUCCESS);
2191 }
PetscErrorCode DMMooseRegisterAll()
std::string name(const ElemQuality q)
PetscInt nvars
static PetscErrorCode DMSetUp_Moose(DM dm)
PetscOptionsEnd()
MooseMesh & mesh()
Definition: SystemBase.h:99
std::set< ContactName > * _uncontacts
Definition: PetscDMMoose.C:82
void allgather(const T &send_data, std::vector< T, A > &recv_data) const
PetscErrorCode DMMooseGetBlocks(DM dm, std::vector< std::string > &block_names)
Definition: PetscDMMoose.C:203
std::map< std::pair< BoundaryID, BoundaryID >, PenetrationLocator * > _penetration_locators
PetscErrorCode DMMooseSetContacts(DM dm, const std::vector< std::pair< std::string, std::string >> &contacts, const std::vector< PetscBool > &displaced)
Definition: PetscDMMoose.C:385
bool _all_vars
Definition: PetscDMMoose.C:62
const DofMapBase * _dof_map
Definition: PetscDMMoose.C:56
std::set< std::pair< BoundaryID, unsigned int > > * _unside_by_var_set
Definition: PetscDMMoose.C:74
PetscErrorCode DMMooseSetVariables(DM dm, const std::set< std::string > &vars)
Definition: PetscDMMoose.C:296
std::map< BoundaryID, std::string > * _side_names
Definition: PetscDMMoose.C:68
std::shared_ptr< DisplacedProblem > displaced_problem
PetscErrorCode DMCreateMoose(MPI_Comm comm, NonlinearSystemBase &nl, const DofMapBase &dof_map, const System &system, const std::string &dm_name, DM *dm)
Create a MOOSE DM.
PetscInt maxsides
PetscInt maxvars
std::set< std::string > varset
MPI_Info info
char ** blocks
std::set< std::string > sideset
PetscErrorCode DMMooseGetUnSides(DM dm, std::vector< std::string > &side_names)
Definition: PetscDMMoose.C:192
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:333
Data structure used to hold penetration information.
virtual unsigned int n_variables() const=0
std::map< ContactName, PetscBool > * _contact_displaced
Definition: PetscDMMoose.C:84
std::vector< DM_Moose::ContactName > uncontacts
PetscErrorCode DMMooseGetNonlinearSystem(DM dm, NonlinearSystemBase *&nl)
Definition: PetscDMMoose.C:446
PetscBool is_include_all_nodes
bool _nounsides
Definition: PetscDMMoose.C:76
PetscInt nsides
DM_Moose * _parent
Definition: PetscDMMoose.C:58
char ** vars
PetscErrorCode DMMooseSetUnContacts(DM dm, const std::vector< std::pair< std::string, std::string >> &uncontacts, const std::vector< PetscBool > &displaced)
Definition: PetscDMMoose.C:415
bool _nouncontacts
Definition: PetscDMMoose.C:87
PetscErrorCode DMSetFromOptions_Moose(DM dm, PetscOptionItems) PetscErrorCode DMSetFromOptions_Moose(DM dm
std::string opt
std::map< std::string, BoundaryID > * _side_ids
Definition: PetscDMMoose.C:69
bool _include_all_contact_nodes
Definition: PetscDMMoose.C:88
std::string * _name
The name of this DM.
Definition: PetscDMMoose.C:105
std::set< std::string > blockset
if(subdm)
PetscErrorCode SNESUpdateDMMoose(SNES snes, PetscInt iteration)
MeshBase & mesh
static PetscErrorCode DMMooseGetMeshBlocks_Private(DM dm, std::set< subdomain_id_type > &blocks)
PetscErrorCode PetscOptionItems *PetscErrorCode DM dm
PetscFunctionBegin
void checkSize(const std::string &split_name, const I1 split_size, const I2 size_expected_by_parent)
Definition: PetscDMMoose.C:39
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const=0
std::map< ContactName, PetscBool > * _uncontact_displaced
Definition: PetscDMMoose.C:85
const Parallel::Communicator & comm() const
virtual libMesh::NonlinearSolver< Number > * nonlinearSolver()=0
std::vector< PetscBool > uncontact_displaced
PetscErrorCode DMMooseReset(DM dm)
PetscErrorCode DMMooseGetUnContacts(DM dm, std::vector< std::pair< std::string, std::string >> &uncontact_names, std::vector< PetscBool > &displaced)
Definition: PetscDMMoose.C:165
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
PetscErrorCode DMMooseSetNonlinearSystem(DM dm, NonlinearSystemBase &nl)
Definition: PetscDMMoose.C:225
const BoundaryInfo & get_boundary_info() const
const MeshBase & get_mesh() const
PETSC_ERR_ARG_WRONGSTATE
PetscInt nblocks
static PetscErrorCode DMSetUp_Moose_Pre(DM dm)
std::set< std::string > * _sides
Definition: PetscDMMoose.C:67
virtual GeometricSearchData & geomSearchData() override
bool _nocontacts
Definition: PetscDMMoose.C:86
static PetscErrorCode DMVariableBounds_Moose(DM dm, Vec xl, Vec xu)
std::map< std::string, unsigned int > * _var_ids
Definition: PetscDMMoose.C:60
virtual void swap(NumericVector< T > &v) override
std::map< unsigned int, std::string > * _var_names
Definition: PetscDMMoose.C:61
EXTERN_C_BEGIN PetscErrorCode DMCreate_Moose(DM dm)
Nonlinear system to be solved.
static PetscErrorCode DMMooseJacobian(DM dm, Vec x, Mat jac, Mat pc)
PetscErrorCode DMMooseGetContacts(DM dm, std::vector< std::pair< std::string, std::string >> &contact_names, std::vector< PetscBool > &displaced)
Definition: PetscDMMoose.C:149
void libmesh_ignore(const Args &...)
FEProblemBase & feProblem()
Definition: SystemBase.h:103
static PetscErrorCode SNESJacobian_DMMoose(SNES, Vec x, Mat jac, Mat pc, void *ctx)
unsigned int number() const
int8_t boundary_id_type
const Node & node_ref(const unsigned int i) const
static PetscErrorCode DMView_Moose(DM dm, PetscViewer viewer)
std::map< ContactID, ContactName > * _uncontact_names
Definition: PetscDMMoose.C:83
char ** sides
boundary_id_type BoundaryID
PetscInt nsplits
bool isCaptured() const
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
std::set< ContactName > * _contacts
Definition: PetscDMMoose.C:80
std::unique_ptr< NumericVector< Number > > solution
libmesh_assert(ctx)
PetscErrorCode DMMooseGetSplitNames(DM dm, std::vector< std::string > &split_names)
Definition: PetscDMMoose.C:492
PetscErrorCode DMMooseSetParentDM(DM dm, DM_Moose *parent)
Definition: PetscDMMoose.C:281
std::set< std::string > * _unside_by_var
Definition: PetscDMMoose.C:73
bool _all_blocks
Definition: PetscDMMoose.C:66
std::string & subdomain_name(subdomain_id_type id)
OStreamProxy err(std::cerr)
PetscInt maxblocks
const Elem * _side
std::map< std::string, SplitInfo > * _splits
Definition: PetscDMMoose.C:99
static PetscErrorCode DMCreateFieldDecomposition_Moose(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
Definition: PetscDMMoose.C:849
virtual bool hasVariable(const std::string &var_name) const
Query a system for a variable.
Definition: SystemBase.C:843
std::pair< BoundaryID, BoundaryID > ContactID
Definition: PetscDMMoose.C:79
std::set< std::string > * _unsides
Definition: PetscDMMoose.C:70
std::set< std::string > * _vars
Definition: PetscDMMoose.C:59
std::string help
tbb::split split
virtual MooseVariableScalar & getScalarVariable(THREAD_ID tid, const std::string &var_name) const
Gets a reference to a scalar variable with specified number.
Definition: SystemBase.C:145
static PetscErrorCode DMMooseGetEmbedding_Private(DM dm, IS *embedding)
Definition: PetscDMMoose.C:512
LibmeshPetscCallQ(DMShellGetContext(dm, &ctx))
virtual void update()
std::map< std::string, BoundaryID > * _unside_ids
Definition: PetscDMMoose.C:71
virtual std::shared_ptr< const DisplacedProblem > getDisplacedProblem() const
PetscErrorCode DMMooseGetSides(DM dm, std::vector< std::string > &side_names)
Definition: PetscDMMoose.C:181
PetscBool _print_embedding
Definition: PetscDMMoose.C:102
subdomain_id_type subdomain_id() const
NonlinearSystemBase * _nl
Definition: PetscDMMoose.C:55
PetscErrorCode DMMooseSetName(DM dm, const std::string &dm_name)
Definition: PetscDMMoose.C:267
static PetscErrorCode SNESFunction_DMMoose(SNES, Vec x, Vec r, void *ctx)
contact_displaced
PetscErrorCode DMMooseValidityCheck(DM dm)
Definition: PetscDMMoose.C:132
const System * _system
Definition: PetscDMMoose.C:57
PetscInt nuncontacts
const std::set< boundary_id_type > & get_global_boundary_ids() const
PetscErrorCode DMMooseSetSides(DM dm, const std::set< std::string > &sides)
Definition: PetscDMMoose.C:340
IntRange< T > make_range(T beg, T end)
virtual MooseMesh & mesh() override
dof_id_type n_local_dofs() const
std::unique_ptr< NumericVector< Number > > current_local_solution
std::map< std::string, subdomain_id_type > * _block_ids
Definition: PetscDMMoose.C:64
IntRange< unsigned short > node_index_range() const
PetscInt maxcontacts
void checkChildSize(DM child, PetscInt child_size, const std::string &child_name)
Check whether the size of the child matches the size we expect.
Definition: PetscDMMoose.C:114
PetscErrorCode DMMooseSetDofMap(DM dm, const DofMapBase &dof_map)
Definition: PetscDMMoose.C:239
PetscErrorCode DMMooseSetUnSideByVar(DM dm, const std::set< std::string > &unside_by_var)
Definition: PetscDMMoose.C:370
void * ctx
static PetscErrorCode DMMooseFunction(DM dm, Vec x, Vec r)
Definition: PetscDMMoose.C:973
libMesh::StoredRange< MooseMesh::const_bnd_elem_iterator, const BndElement * > * getBoundaryElementRange()
Definition: MooseMesh.C:1303
MooseVariableFieldBase & getVariable(THREAD_ID tid, const std::string &var_name) const
Gets a reference to a variable of with specified name.
Definition: SystemBase.C:90
std::set< subdomain_id_type > meshblocks
const std::string & name() const
std::set< std::string > * _blocks
Definition: PetscDMMoose.C:63
const char * fdhelp
PetscFunctionReturn(LIBMESH_PETSC_SUCCESS)
static PetscErrorCode DMCreateGlobalVector_Moose(DM dm, Vec *x)
PetscErrorCode DMMooseSetBlocks(DM dm, const std::set< std::string > &blocks)
Definition: PetscDMMoose.C:325
std::map< ContactID, ContactName > * _contact_names
Definition: PetscDMMoose.C:81
PetscErrorCode DMMooseSetSystem(DM dm, const System &system)
Definition: PetscDMMoose.C:253
bool _nosides
Definition: PetscDMMoose.C:75
std::multimap< std::string, unsigned int > * _splitlocs
Definition: PetscDMMoose.C:93
std::pair< std::string, std::string > ContactName
Definition: PetscDMMoose.C:78
const DofMap & get_dof_map() const
static PetscErrorCode DMDestroy_Moose(DM dm)
static PetscErrorCode DMCreateDomainDecomposition_Moose(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
Definition: PetscDMMoose.C:960
PetscErrorCode DMMooseSetSplitNames(DM dm, const std::vector< std::string > &split_names)
Definition: PetscDMMoose.C:456
static PetscErrorCode DMCreateMatrix_Moose(DM dm, Mat *A)
contacts
for(PetscInt i=0;i< nvars;++i)
libMesh::StoredRange< MooseMesh::const_bnd_node_iterator, const BndNode * > * getBoundaryNodeRange()
Definition: MooseMesh.C:1289
PetscOptionsBegin(((PetscObject) dm) ->comm,((PetscObject) dm) ->prefix, "DMMoose options", "DM")
PetscErrorCode DMMooseSetUnSides(DM dm, const std::set< std::string > &unsides)
Definition: PetscDMMoose.C:355
bool _nounside_by_var
Definition: PetscDMMoose.C:77
Base variable class.
BoundaryID getBoundaryID(const BoundaryName &boundary_name) const
Get the associated BoundaryID for the boundary name.
Definition: MooseMesh.C:1692
PetscErrorCode DMMooseGetVariables(DM dm, std::vector< std::string > &var_names)
Definition: PetscDMMoose.C:214
uint8_t dof_id_type
std::map< BoundaryID, std::string > * _unside_names
Definition: PetscDMMoose.C:72
PenetrationLocator & getPenetrationLocator(const BoundaryName &primary, const BoundaryName &secondary, libMesh::Order order=libMesh::FIRST)
const std::map< dof_id_type, std::vector< dof_id_type > > & nodeToElemMap()
If not already created, creates a map from every node to all elements to which they are connected...
Definition: MooseMesh.C:1178
virtual void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di, const unsigned int vn, int p_level=-12345) const=0
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=nullptr, bool homogeneous=false) const
void set_union(T &data, const unsigned int root_id) const
std::map< unsigned int, std::string > * _block_names
Definition: PetscDMMoose.C:65
virtual libMesh::System & system() override
Get the reference to the libMesh system.