www.mooseframework.org
PetscDMMoose.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 // This only works with petsc-3.3 and above.
11 #include "libmesh/petsc_macro.h"
12 
13 #if defined(LIBMESH_HAVE_PETSC) && !PETSC_VERSION_LESS_THAN(3, 3, 0)
14 // Inside these guards we can use PETSC_VERSION_LT, which need not be
15 // modified upon transition from dev to a release.
16 
17 #include "PetscDMMoose.h"
18 
19 // PETSc includes
20 #include <petscerror.h>
21 #if !PETSC_RELEASE_LESS_THAN(3, 6, 0)
22 #include <petsc/private/dmimpl.h>
23 #else
24 #include <petsc-private/dmimpl.h>
25 #endif
26 
27 // MOOSE includes
28 #include "FEProblem.h"
29 #include "DisplacedProblem.h"
30 #include "MooseMesh.h"
31 #include "NonlinearSystem.h"
32 #include "PenetrationLocator.h"
33 #include "NearestNodeLocator.h"
34 #include "GeometricSearchData.h"
35 
36 #include "libmesh/nonlinear_implicit_system.h"
37 #include "libmesh/nonlinear_solver.h"
38 #include "libmesh/petsc_vector.h"
39 #include "libmesh/petsc_matrix.h"
40 #include "libmesh/dof_map.h"
41 #include "libmesh/preconditioner.h"
42 
43 struct DM_Moose
44 {
45  NonlinearSystemBase * _nl; // nonlinear system context
46  std::set<std::string> * _vars; // variables
47  std::map<std::string, unsigned int> * _var_ids;
48  std::map<unsigned int, std::string> * _var_names;
49  bool _all_vars; // whether all system variables are included
50  std::set<std::string> * _blocks; // mesh blocks
51  std::map<std::string, subdomain_id_type> * _block_ids;
52  std::map<unsigned int, std::string> * _block_names;
53  bool _all_blocks; // all blocks are included
54  std::set<std::string> * _sides; // mesh surfaces (edges in 2D)
55  std::map<BoundaryID, std::string> * _side_names;
56  std::map<std::string, BoundaryID> * _side_ids;
57  std::set<std::string> * _unsides; // excluded sides
58  std::map<std::string, BoundaryID> * _unside_ids;
59  std::map<BoundaryID, std::string> * _unside_names;
60  bool _nosides; // whether to include any sides
61  bool _nounsides; // whether to exclude any sides
62  typedef std::pair<std::string, std::string> ContactName;
63  typedef std::pair<BoundaryID, BoundaryID> ContactID;
64  std::set<ContactName> * _contacts;
65  std::map<ContactID, ContactName> * _contact_names;
66  std::set<ContactName> * _uncontacts;
67  std::map<ContactID, ContactName> * _uncontact_names;
68  std::map<ContactName, PetscBool> * _contact_displaced;
69  std::map<ContactName, PetscBool> * _uncontact_displaced;
73  // to locate splits without having to search, however,
74  // maintain a multimap from names to split locations (to enable
75  // the same split to appear in multiple spots (this might
76  // break the current implementation of PCFieldSplit, though).
77  std::multimap<std::string, unsigned int> * _splitlocs;
78  struct SplitInfo
79  {
80  DM _dm;
81  IS _rembedding; // relative embedding
82  };
83  std::map<std::string, SplitInfo> * _splits;
85  PetscBool _print_embedding;
86 };
87 
88 #undef __FUNCT__
89 #define __FUNCT__ "DMMooseGetContacts"
90 PetscErrorCode
92  std::vector<std::pair<std::string, std::string>> & contact_names,
93  std::vector<PetscBool> & displaced)
94 {
95  PetscErrorCode ierr;
96  PetscBool ismoose;
97 
99  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
100  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
101  CHKERRQ(ierr);
102  if (!ismoose)
103  SETERRQ2(((PetscObject)dm)->comm,
104  PETSC_ERR_ARG_WRONG,
105  "Got DM oftype %s, not of type %s",
106  ((PetscObject)dm)->type_name,
107  DMMOOSE);
108  DM_Moose * dmm = (DM_Moose *)dm->data;
109  for (const auto & it : *(dmm->_contact_names))
110  {
111  contact_names.push_back(it.second);
112  displaced.push_back((*dmm->_contact_displaced)[it.second]);
113  }
115 }
116 
117 #undef __FUNCT__
118 #define __FUNCT__ "DMMooseGetUnContacts"
119 PetscErrorCode
121  std::vector<std::pair<std::string, std::string>> & uncontact_names,
122  std::vector<PetscBool> & displaced)
123 {
124  PetscErrorCode ierr;
125  PetscBool ismoose;
126 
128  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
129  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
130  CHKERRQ(ierr);
131  if (!ismoose)
132  SETERRQ2(((PetscObject)dm)->comm,
133  PETSC_ERR_ARG_WRONG,
134  "Got DM oftype %s, not of type %s",
135  ((PetscObject)dm)->type_name,
136  DMMOOSE);
137  DM_Moose * dmm = (DM_Moose *)dm->data;
138  for (const auto & it : *(dmm->_uncontact_names))
139  {
140  uncontact_names.push_back(it.second);
141  displaced.push_back((*dmm->_uncontact_displaced)[it.second]);
142  }
144 }
145 
146 #undef __FUNCT__
147 #define __FUNCT__ "DMMooseGetSides"
148 PetscErrorCode
149 DMMooseGetSides(DM dm, std::vector<std::string> & side_names)
150 {
151  PetscErrorCode ierr;
152  PetscBool ismoose;
153 
155  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
156  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
157  CHKERRQ(ierr);
158  if (!ismoose)
159  SETERRQ2(((PetscObject)dm)->comm,
160  PETSC_ERR_ARG_WRONG,
161  "Got DM oftype %s, not of type %s",
162  ((PetscObject)dm)->type_name,
163  DMMOOSE);
164  DM_Moose * dmm = (DM_Moose *)dm->data;
165  for (const auto & it : *(dmm->_side_ids))
166  side_names.push_back(it.first);
168 }
169 
170 #undef __FUNCT__
171 #define __FUNCT__ "DMMooseGetUnSides"
172 PetscErrorCode
173 DMMooseGetUnSides(DM dm, std::vector<std::string> & side_names)
174 {
175  PetscErrorCode ierr;
176  PetscBool ismoose;
177 
179  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
180  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
181  CHKERRQ(ierr);
182  if (!ismoose)
183  SETERRQ2(((PetscObject)dm)->comm,
184  PETSC_ERR_ARG_WRONG,
185  "Got DM oftype %s, not of type %s",
186  ((PetscObject)dm)->type_name,
187  DMMOOSE);
188  DM_Moose * dmm = (DM_Moose *)dm->data;
189  for (const auto & it : *(dmm->_unside_ids))
190  side_names.push_back(it.first);
192 }
193 
194 #undef __FUNCT__
195 #define __FUNCT__ "DMMooseGetBlocks"
196 PetscErrorCode
197 DMMooseGetBlocks(DM dm, std::vector<std::string> & block_names)
198 {
199  PetscErrorCode ierr;
200  PetscBool ismoose;
201 
203  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
204  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
205  CHKERRQ(ierr);
206  if (!ismoose)
207  SETERRQ2(((PetscObject)dm)->comm,
208  PETSC_ERR_ARG_WRONG,
209  "Got DM oftype %s, not of type %s",
210  ((PetscObject)dm)->type_name,
211  DMMOOSE);
212  DM_Moose * dmm = (DM_Moose *)dm->data;
213  for (const auto & it : *(dmm->_block_ids))
214  block_names.push_back(it.first);
216 }
217 
218 #undef __FUNCT__
219 #define __FUNCT__ "DMMooseGetVariables"
220 PetscErrorCode
221 DMMooseGetVariables(DM dm, std::vector<std::string> & var_names)
222 {
223  PetscErrorCode ierr;
224  PetscBool ismoose;
225 
227  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
228  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
229  CHKERRQ(ierr);
230  if (!ismoose)
231  SETERRQ2(((PetscObject)dm)->comm,
232  PETSC_ERR_ARG_WRONG,
233  "Got DM oftype %s, not of type %s",
234  ((PetscObject)dm)->type_name,
235  DMMOOSE);
236  DM_Moose * dmm = (DM_Moose *)(dm->data);
237  for (const auto & it : *(dmm->_var_ids))
238  var_names.push_back(it.first);
240 }
241 
242 #undef __FUNCT__
243 #define __FUNCT__ "DMMooseSetNonlinearSystem"
244 PetscErrorCode
246 {
247  PetscErrorCode ierr;
248  PetscBool ismoose;
249 
251  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
252  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
253  CHKERRQ(ierr);
254  if (!ismoose)
255  SETERRQ2(((PetscObject)dm)->comm,
256  PETSC_ERR_ARG_WRONG,
257  "Got DM oftype %s, not of type %s",
258  ((PetscObject)dm)->type_name,
259  DMMOOSE);
260  if (dm->setupcalled)
261  SETERRQ(((PetscObject)dm)->comm,
263  "Cannot reset the NonlinearSystem after DM has been set up.");
264  DM_Moose * dmm = (DM_Moose *)(dm->data);
265  dmm->_nl = &nl;
267 }
268 
269 #undef __FUNCT__
270 #define __FUNCT__ "DMMooseSetVariables"
271 PetscErrorCode
272 DMMooseSetVariables(DM dm, const std::set<std::string> & vars)
273 {
274  PetscErrorCode ierr;
275  DM_Moose * dmm = (DM_Moose *)dm->data;
276  PetscBool ismoose;
277 
279  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
280  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
281  CHKERRQ(ierr);
282  if (!ismoose)
283  SETERRQ2(PETSC_COMM_SELF,
284  PETSC_ERR_ARG_WRONG,
285  "Got DM oftype %s, not of type %s",
286  ((PetscObject)dm)->type_name,
287  DMMOOSE);
288  if (dm->setupcalled)
289  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
290  if (dmm->_vars)
291  delete dmm->_vars;
292  dmm->_vars = new std::set<std::string>(vars);
294 }
295 
296 #undef __FUNCT__
297 #define __FUNCT__ "DMMooseSetBlocks"
298 PetscErrorCode
299 DMMooseSetBlocks(DM dm, const std::set<std::string> & blocks)
300 {
301  PetscErrorCode ierr;
302  DM_Moose * dmm = (DM_Moose *)dm->data;
303  PetscBool ismoose;
304 
306  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
307  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
308  CHKERRQ(ierr);
309  if (!ismoose)
310  SETERRQ2(PETSC_COMM_SELF,
311  PETSC_ERR_ARG_WRONG,
312  "Got DM oftype %s, not of type %s",
313  ((PetscObject)dm)->type_name,
314  DMMOOSE);
315  if (dm->setupcalled)
316  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
317  if (dmm->_blocks)
318  delete dmm->_blocks;
319  dmm->_blocks = new std::set<std::string>(blocks);
321 }
322 
323 #undef __FUNCT__
324 #define __FUNCT__ "DMMooseSetSides"
325 PetscErrorCode
326 DMMooseSetSides(DM dm, const std::set<std::string> & sides)
327 {
328  PetscErrorCode ierr;
329  DM_Moose * dmm = (DM_Moose *)dm->data;
330  PetscBool ismoose;
331 
333  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
334  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
335  CHKERRQ(ierr);
336  if (!ismoose)
337  SETERRQ2(PETSC_COMM_SELF,
338  PETSC_ERR_ARG_WRONG,
339  "Got DM oftype %s, not of type %s",
340  ((PetscObject)dm)->type_name,
341  DMMOOSE);
342  if (dm->setupcalled)
343  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
344  if (dmm->_sides)
345  delete dmm->_sides;
346  dmm->_sides = new std::set<std::string>(sides);
348 }
349 
350 #undef __FUNCT__
351 #define __FUNCT__ "DMMooseSetUnSides"
352 PetscErrorCode
353 DMMooseSetUnSides(DM dm, const std::set<std::string> & unsides)
354 {
355  PetscErrorCode ierr;
356  DM_Moose * dmm = (DM_Moose *)dm->data;
357  PetscBool ismoose;
358 
360  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
361  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
362  CHKERRQ(ierr);
363  if (!ismoose)
364  SETERRQ2(PETSC_COMM_SELF,
365  PETSC_ERR_ARG_WRONG,
366  "Got DM oftype %s, not of type %s",
367  ((PetscObject)dm)->type_name,
368  DMMOOSE);
369  if (dm->setupcalled)
370  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
371  if (dmm->_sides)
372  delete dmm->_sides;
373  dmm->_unsides = new std::set<std::string>(unsides);
375 }
376 
377 #undef __FUNCT__
378 #define __FUNCT__ "DMMooseSetContacts"
379 PetscErrorCode
381  const std::vector<std::pair<std::string, std::string>> & contacts,
382  const std::vector<PetscBool> & displaced)
383 {
384  PetscErrorCode ierr;
385  DM_Moose * dmm = (DM_Moose *)dm->data;
386  PetscBool ismoose;
387 
389  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
390  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
391  CHKERRQ(ierr);
392  if (!ismoose)
393  SETERRQ2(PETSC_COMM_SELF,
394  PETSC_ERR_ARG_WRONG,
395  "Got DM oftype %s, not of type %s",
396  ((PetscObject)dm)->type_name,
397  DMMOOSE);
398  if (dm->setupcalled)
399  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
400  if (contacts.size() != displaced.size())
401  SETERRQ2(PETSC_COMM_SELF,
402  PETSC_ERR_ARG_SIZ,
403  "Nonmatching sizes of the contact and displaced arrays: %D != %D",
404  contacts.size(),
405  displaced.size());
406  if (dmm->_contacts)
407  delete dmm->_contacts;
408  dmm->_contact_displaced->clear();
409  dmm->_contacts = new std::set<DM_Moose::ContactName>();
410  for (unsigned int i = 0; i < contacts.size(); ++i)
411  {
412  dmm->_contacts->insert(contacts[i]);
413  dmm->_contact_displaced->insert(std::make_pair(contacts[i], displaced[i]));
414  }
416 }
417 
418 #undef __FUNCT__
419 #define __FUNCT__ "DMMooseSetUnContacts"
420 PetscErrorCode
422  const std::vector<std::pair<std::string, std::string>> & uncontacts,
423  const std::vector<PetscBool> & displaced)
424 {
425  PetscErrorCode ierr;
426  DM_Moose * dmm = (DM_Moose *)dm->data;
427  PetscBool ismoose;
428 
430  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
431  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
432  CHKERRQ(ierr);
433  if (!ismoose)
434  SETERRQ2(PETSC_COMM_SELF,
435  PETSC_ERR_ARG_WRONG,
436  "Got DM oftype %s, not of type %s",
437  ((PetscObject)dm)->type_name,
438  DMMOOSE);
439  if (dm->setupcalled)
440  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for an already setup DM");
441  if (uncontacts.size() != displaced.size())
442  SETERRQ2(PETSC_COMM_SELF,
443  PETSC_ERR_ARG_SIZ,
444  "Nonmatching sizes of the uncontact and displaced arrays: %D != %D",
445  uncontacts.size(),
446  displaced.size());
447  if (dmm->_uncontacts)
448  delete dmm->_uncontacts;
449  dmm->_uncontact_displaced->clear();
450  dmm->_uncontacts = new std::set<DM_Moose::ContactName>();
451  for (unsigned int i = 0; i < uncontacts.size(); ++i)
452  {
453  dmm->_uncontacts->insert(uncontacts[i]);
454  dmm->_uncontact_displaced->insert(std::make_pair(uncontacts[i], displaced[i]));
455  }
457 }
458 
459 #undef __FUNCT__
460 #define __FUNCT__ "DMMooseGetNonlinearSystem"
461 PetscErrorCode
463 {
464  PetscErrorCode ierr;
465  PetscBool ismoose;
466 
468  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
469  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
470  CHKERRQ(ierr);
471  if (!ismoose)
472  SETERRQ2(((PetscObject)dm)->comm,
473  PETSC_ERR_ARG_WRONG,
474  "Got DM oftype %s, not of type %s",
475  ((PetscObject)dm)->type_name,
476  DMMOOSE);
477  DM_Moose * dmm = (DM_Moose *)(dm->data);
478  nl = dmm->_nl;
480 }
481 
482 #undef __FUNCT__
483 #define __FUNCT__ "DMMooseSetSplitNames"
484 PetscErrorCode
485 DMMooseSetSplitNames(DM dm, const std::vector<std::string> & split_names)
486 {
487  PetscErrorCode ierr;
488  PetscBool ismoose;
489 
491  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
492  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
493  CHKERRQ(ierr);
494  if (!ismoose)
495  SETERRQ2(((PetscObject)dm)->comm,
496  PETSC_ERR_ARG_WRONG,
497  "Got DM oftype %s, not of type %s",
498  ((PetscObject)dm)->type_name,
499  DMMOOSE);
500  DM_Moose * dmm = (DM_Moose *)(dm->data);
501 
502  if (dmm->_splits)
503  {
504  for (auto & it : *(dmm->_splits))
505  {
506  ierr = DMDestroy(&(it.second._dm));
507  CHKERRQ(ierr);
508  ierr = ISDestroy(&(it.second._rembedding));
509  CHKERRQ(ierr);
510  }
511  delete dmm->_splits;
512  dmm->_splits = PETSC_NULL;
513  }
514  if (dmm->_splitlocs)
515  {
516  delete dmm->_splitlocs;
517  dmm->_splitlocs = PETSC_NULL;
518  }
519  dmm->_splits = new std::map<std::string, DM_Moose::SplitInfo>();
520  dmm->_splitlocs = new std::multimap<std::string, unsigned int>();
521  for (unsigned int i = 0; i < split_names.size(); ++i)
522  {
523  DM_Moose::SplitInfo info;
524  info._dm = PETSC_NULL;
525  info._rembedding = PETSC_NULL;
526  std::string name = split_names[i];
527  (*dmm->_splits)[name] = info;
528  dmm->_splitlocs->insert(std::make_pair(name, i));
529  }
531 }
532 
533 #undef __FUNCT__
534 #define __FUNCT__ "DMMooseGetSplitNames"
535 PetscErrorCode
536 DMMooseGetSplitNames(DM dm, std::vector<std::string> & split_names)
537 {
538  PetscErrorCode ierr;
539  PetscBool ismoose;
540 
542  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
543  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
544  CHKERRQ(ierr);
545  if (!ismoose)
546  SETERRQ2(((PetscObject)dm)->comm,
547  PETSC_ERR_ARG_WRONG,
548  "Got DM oftype %s, not of type %s",
549  ((PetscObject)dm)->type_name,
550  DMMOOSE);
551  DM_Moose * dmm = (DM_Moose *)(dm->data);
552  if (!dm->setupcalled)
553  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "DM not set up");
554  split_names.clear();
555  split_names.reserve(dmm->_splitlocs->size());
556  if (dmm->_splitlocs && dmm->_splitlocs->size())
557  for (const auto & lit : *(dmm->_splitlocs))
558  {
559  std::string sname = lit.first;
560  unsigned int sloc = lit.second;
561  split_names[sloc] = sname;
562  }
564 }
565 
566 #undef __FUNCT__
567 #define __FUNCT__ "DMMooseGetEmbedding_Private"
568 static PetscErrorCode
569 DMMooseGetEmbedding_Private(DM dm, IS * embedding)
570 {
571  DM_Moose * dmm = (DM_Moose *)dm->data;
572  PetscErrorCode ierr;
573 
575  if (!embedding)
577  if (!dmm->_embedding)
578  {
579  // The rules interpreting the coexistence of blocks (un)sides/(un)contacts are these
580  // [sides and contacts behave similarly, so 'sides' means 'sides/contacts']
581  // ['ANY' means 'not NONE' and covers 'ALL' as well, unless there is a specific 'ALL' clause,
582  // which overrides 'ANY'; 'NOT ALL' means not ALL and not NONE]
583  // [there are always some blocks, since by default 'ALL' is assumed, unless it is overridden by
584  // a specific list, which implies ANY]
585  // In general,
586  // (1) ALL blocks and ANY sides are interpreted as the INTERSECTION of blocks and sides,
587  // equivalent to just the sides (since ALL blocks are assumed to be a cover).
588  // (2) NOT ALL blocks and ANY or NO sides are interpreted as the UNION of blocks and sides.
589  // (3a) ANY unsides and ANY blocks are interpreted as the DIFFERENCE of blocks and unsides.
590  // (3b) ANY unsides and ANY sides are interpreted as the DIFFERENCE of sides and unsides.
591  // (4) NO unsides means NO DIFFERENCE is needed.
592  // The result is easily computed by first computing the result of (1 & 2) followed by difference
593  // with the result of (3 & 4).
594  // To simply (1 & 2) observe the following:
595  // - The intersection is computed only if ALL blocks and ANY sides, and the result is the sides,
596  // so block dofs do not need to be computed.
597  // - Otherwise the union is computed, and initially consists of the blocks' dofs, to which the
598  // sides' dofs are added, if ANY.
599  // - The result is called 'indices'
600  // To satisfy (3 & 4) simply cmpute subtrahend set 'unindices' as all of the unsides' dofs:
601  // Then take the set difference of 'indices' and 'unindices', putting the result in 'dindices'.
602  if (!dmm->_all_vars || !dmm->_all_blocks || !dmm->_nosides || !dmm->_nounsides ||
604  {
605  DofMap & dofmap = dmm->_nl->system().get_dof_map();
606  std::set<dof_id_type> indices;
607  std::set<dof_id_type> unindices;
608  std::set<dof_id_type> cached_indices;
609  std::set<dof_id_type> cached_unindices;
610  const auto & node_to_elem_map = dmm->_nl->_fe_problem.mesh().nodeToElemMap();
611  for (const auto & vit : *(dmm->_var_ids))
612  {
613  unsigned int v = vit.second;
614  // Iterate only over this DM's blocks.
615  if (!dmm->_all_blocks || (dmm->_nosides && dmm->_nocontacts))
616  {
617  for (const auto & bit : *(dmm->_block_ids))
618  {
619  subdomain_id_type b = bit.second;
620  for (const auto & elem :
621  as_range(dmm->_nl->system().get_mesh().active_local_subdomain_elements_begin(b),
622  dmm->_nl->system().get_mesh().active_local_subdomain_elements_end(b)))
623  {
624  // Get the degree of freedom indices for the given variable off the current element.
625  std::vector<dof_id_type> evindices;
626  dofmap.dof_indices(elem, evindices, v);
627 
628  // might want to use variable_first/last_local_dof instead
629  for (const auto & dof : evindices)
630  if (dof >= dofmap.first_dof() && dof < dofmap.end_dof())
631  indices.insert(dof);
632  }
633 
634  // Sometime, we own nodes but do not own the elements the nodes connected to
635  {
636  bool is_on_current_block = false;
637  for (auto & node : dmm->_nl->system().get_mesh().local_node_ptr_range())
638  {
639  const unsigned int n_comp = node->n_comp(dmm->_nl->system().number(), v);
640 
641  // skip it if no dof
642  if (!n_comp)
643  continue;
644 
645  auto node_to_elem_pair = node_to_elem_map.find(node->id());
646  is_on_current_block = false;
647  for (const auto & elem_num : node_to_elem_pair->second)
648  {
649  // if one of incident elements belongs to a block, we consider
650  // the node lives in the block
651  Elem & neighbor_elem = dmm->_nl->system().get_mesh().elem_ref(elem_num);
652  if (neighbor_elem.subdomain_id() == b)
653  {
654  is_on_current_block = true;
655  break;
656  }
657  }
658  // we add indices for the current block only
659  if (!is_on_current_block)
660  continue;
661 
662  const dof_id_type index = node->dof_number(dmm->_nl->system().number(), v, 0);
663  if (index >= dofmap.first_dof() && index < dofmap.end_dof())
664  indices.insert(index);
665  }
666  }
667  }
668  }
669 
670  // Iterate over the sides from this split.
671  if (dmm->_side_ids->size())
672  {
673  // For some reason the following may return an empty node list
674  // std::vector<dof_id_type> snodes;
675  // std::vector<boundary_id_type> sides;
676  // dmm->nl->system().get_mesh().get_boundary_info().build_node_list(snodes, sides);
677  // // FIXME: make an array of (snode,side) pairs, sort on side and use std::lower_bound
678  // from <algorithm>
679  // for (dof_id_type i = 0; i < sides.size(); ++i) {
680  // boundary_id_type s = sides[i];
681  // if (!dmm->sidenames->count(s)) continue;
682  // const Node& node = dmm->nl->system().get_mesh().node_ref(snodes[i]);
683  // // determine v's dof on node and insert into indices
684  // }
686  for (const auto & bnode : bnodes)
687  {
688  BoundaryID boundary_id = bnode->_bnd_id;
689  if (dmm->_side_names->find(boundary_id) == dmm->_side_names->end())
690  continue;
691 
692  const Node * node = bnode->_node;
693  dof_id_type dof = node->dof_number(dmm->_nl->system().number(), v, 0);
694 
695  // might want to use variable_first/last_local_dof instead
696  if (dof >= dofmap.first_dof() && dof < dofmap.end_dof())
697  indices.insert(dof);
698  }
699  }
700 
701  // Iterate over the sides excluded from this split.
702  if (dmm->_unside_ids->size())
703  {
705  for (const auto & bnode : bnodes)
706  {
707  BoundaryID boundary_id = bnode->_bnd_id;
708  if (dmm->_unside_names->find(boundary_id) == dmm->_unside_names->end())
709  continue;
710  const Node * node = bnode->_node;
711 
712  // might want to use variable_first/last_local_dof instead
713  dof_id_type dof = node->dof_number(dmm->_nl->system().number(), v, 0);
714  if (dof >= dofmap.first_dof() && dof < dofmap.end_dof())
715  unindices.insert(dof);
716  }
717  }
718 
719  // Include all nodes on the contact surfaces
721  {
722  std::set<boundary_id_type> bc_id_set;
723  // loop over contacts
724  for (const auto & it : *(dmm->_contact_names))
725  {
726  bc_id_set.insert(it.first.first); // master
727  bc_id_set.insert(it.first.second); // slave
728  }
729  // loop over boundary elements
730  std::vector<dof_id_type> evindices;
732  for (const auto & belem : range)
733  {
734  const Elem * elem_bdry = belem->_elem;
735  BoundaryID boundary_id = belem->_bnd_id;
736 
737  if (bc_id_set.find(boundary_id) == bc_id_set.end())
738  continue;
739 
740  evindices.clear();
741  dofmap.dof_indices(elem_bdry, evindices, v);
742  for (const auto & edof : evindices)
743  if (edof >= dofmap.first_dof() && edof < dofmap.end_dof())
744  indices.insert(edof);
745  }
746  }
747 
748  // Iterate over the contacts included in this split.
749  if (dmm->_contact_names->size() && !(dmm->_include_all_contact_nodes))
750  {
751  std::vector<dof_id_type> evindices;
752  for (const auto & it : *(dmm->_contact_names))
753  {
754  PetscBool displaced = (*dmm->_contact_displaced)[it.second];
755  PenetrationLocator * locator;
756  if (displaced)
757  {
758  std::shared_ptr<DisplacedProblem> displaced_problem =
760  if (!displaced_problem)
761  {
762  std::ostringstream err;
763  err << "Cannot use a displaced contact (" << it.second.first << ","
764  << it.second.second << ") with an undisplaced problem";
765  mooseError(err.str());
766  }
767  locator = displaced_problem->geomSearchData()._penetration_locators[it.first];
768  }
769  else
770  locator = dmm->_nl->_fe_problem.geomSearchData()._penetration_locators[it.first];
771 
772  evindices.clear();
773  // penetration locator
774  auto lend = locator->_penetration_info.end();
775  for (auto lit = locator->_penetration_info.begin(); lit != lend; ++lit)
776  {
777  const dof_id_type slave_node_num = lit->first;
778  PenetrationInfo * pinfo = lit->second;
779  if (pinfo && pinfo->isCaptured())
780  {
781  Node & slave_node = dmm->_nl->system().get_mesh().node_ref(slave_node_num);
782  dof_id_type dof = slave_node.dof_number(dmm->_nl->system().number(), v, 0);
783  // might want to use variable_first/last_local_dof instead
784  if (dof >= dofmap.first_dof() && dof < dofmap.end_dof())
785  indices.insert(dof);
786  else
787  cached_indices.insert(dof); // cache nonlocal indices
788  // indices of slave elements
789  evindices.clear();
790 
791  auto node_to_elem_pair = node_to_elem_map.find(slave_node_num);
792  mooseAssert(node_to_elem_pair != node_to_elem_map.end(),
793  "Missing entry in node to elem map");
794  for (const auto & elem_num : node_to_elem_pair->second)
795  {
796  Elem & slave_elem = dmm->_nl->system().get_mesh().elem_ref(elem_num);
797  // Get the degree of freedom indices for the given variable off the current
798  // element.
799  evindices.clear();
800  dofmap.dof_indices(&slave_elem, evindices, v);
801  // might want to use variable_first/last_local_dof instead
802  for (const auto & edof : evindices)
803  if (edof >= dofmap.first_dof() && edof < dofmap.end_dof())
804  indices.insert(edof);
805  else
806  cached_indices.insert(edof);
807  }
808  // indices for master element
809  evindices.clear();
810  const Elem * master_elem = pinfo->_elem;
811  dofmap.dof_indices(master_elem, evindices, v);
812  for (const auto & edof : evindices)
813  if (edof >= dofmap.first_dof() && edof < dofmap.end_dof())
814  indices.insert(edof);
815  else
816  cached_indices.insert(edof);
817  } // if pinfo
818  } // for penetration
819  } // for contact names
820  } // if size of contact names
821 
823  {
824  std::set<boundary_id_type> bc_id_set;
825  // loop over contacts
826  for (const auto & it : *(dmm->_uncontact_names))
827  {
828  bc_id_set.insert(it.first.first);
829  bc_id_set.insert(it.first.second);
830  }
831  // loop over boundary elements
832  std::vector<dof_id_type> evindices;
834  for (const auto & belem : range)
835  {
836  const Elem * elem_bdry = belem->_elem;
837  unsigned short int side = belem->_side;
838  BoundaryID boundary_id = belem->_bnd_id;
839 
840  if (bc_id_set.find(boundary_id) == bc_id_set.end())
841  continue;
842 
843  std::unique_ptr<const Elem> side_bdry = elem_bdry->build_side_ptr(side, false);
844  evindices.clear();
845  dofmap.dof_indices(side_bdry.get(), evindices, v);
846  for (const auto & edof : evindices)
847  if (edof >= dofmap.first_dof() && edof < dofmap.end_dof())
848  unindices.insert(edof);
849  }
850  }
851 
852  // Iterate over the contacts excluded from this split.
854  {
855  std::vector<dof_id_type> evindices;
856  for (const auto & it : *(dmm->_uncontact_names))
857  {
858  PetscBool displaced = (*dmm->_uncontact_displaced)[it.second];
859  PenetrationLocator * locator;
860  if (displaced)
861  {
862  std::shared_ptr<DisplacedProblem> displaced_problem =
864  if (!displaced_problem)
865  {
866  std::ostringstream err;
867  err << "Cannot use a displaced uncontact (" << it.second.first << ","
868  << it.second.second << ") with an undisplaced problem";
869  mooseError(err.str());
870  }
871  locator = displaced_problem->geomSearchData()._penetration_locators[it.first];
872  }
873  else
874  locator = dmm->_nl->_fe_problem.geomSearchData()._penetration_locators[it.first];
875 
876  evindices.clear();
877  // penetration locator
878  auto lend = locator->_penetration_info.end();
879  for (auto lit = locator->_penetration_info.begin(); lit != lend; ++lit)
880  {
881  const dof_id_type slave_node_num = lit->first;
882  PenetrationInfo * pinfo = lit->second;
883  if (pinfo && pinfo->isCaptured())
884  {
885  Node & slave_node = dmm->_nl->system().get_mesh().node_ref(slave_node_num);
886  dof_id_type dof = slave_node.dof_number(dmm->_nl->system().number(), v, 0);
887  // might want to use variable_first/last_local_dof instead
888  if (dof >= dofmap.first_dof() && dof < dofmap.end_dof())
889  unindices.insert(dof);
890  else
891  cached_unindices.insert(dof);
892 
893  // indices for master element
894  evindices.clear();
895  const Elem * master_side = pinfo->_side;
896  dofmap.dof_indices(master_side, evindices, v);
897  // indices of master sides
898  for (const auto & edof : evindices)
899  if (edof >= dofmap.first_dof() && edof < dofmap.end_dof())
900  unindices.insert(edof);
901  else
902  cached_unindices.insert(edof);
903  } // if pinfo
904  } // for penetration
905  } // for uncontact names
906  } // if there exist uncontacts
907  } // variables
908 
909  std::vector<dof_id_type> local_vec_indices(cached_indices.size());
910  std::copy(cached_indices.begin(), cached_indices.end(), local_vec_indices.begin());
911  if (dmm->_contact_names->size() && !(dmm->_include_all_contact_nodes))
912  dmm->_nl->_fe_problem.mesh().comm().allgather(local_vec_indices, false);
913  // insert indices
914  for (const auto & dof : local_vec_indices)
915  if (dof >= dofmap.first_dof() && dof < dofmap.end_dof())
916  indices.insert(dof);
917 
918  local_vec_indices.clear();
919  local_vec_indices.resize(cached_unindices.size());
920  std::copy(cached_unindices.begin(), cached_unindices.end(), local_vec_indices.begin());
922  dmm->_nl->_fe_problem.mesh().comm().allgather(local_vec_indices, false);
923  // insert unindices
924  for (const auto & dof : local_vec_indices)
925  if (dof >= dofmap.first_dof() && dof < dofmap.end_dof())
926  unindices.insert(dof);
927 
928  std::set<dof_id_type> dindices;
929  std::set_difference(indices.begin(),
930  indices.end(),
931  unindices.begin(),
932  unindices.end(),
933  std::inserter(dindices, dindices.end()));
934  PetscInt * darray;
935  ierr = PetscMalloc(sizeof(PetscInt) * dindices.size(), &darray);
936  CHKERRQ(ierr);
937  dof_id_type i = 0;
938  for (const auto & dof : dindices)
939  {
940  darray[i] = dof;
941  ++i;
942  }
943  ierr = ISCreateGeneral(
944  ((PetscObject)dm)->comm, dindices.size(), darray, PETSC_OWN_POINTER, &dmm->_embedding);
945  CHKERRQ(ierr);
946  }
947  else
948  {
949  // if (dmm->allblocks && dmm->allvars && dmm->nosides && dmm->nounsides && dmm->nocontacts &&
950  // dmm->nouncontacts)
951  // DMCreateGlobalVector is defined()
952  Vec v;
953  PetscInt low, high;
954 
955  ierr = DMCreateGlobalVector(dm, &v);
956  CHKERRQ(ierr);
957  ierr = VecGetOwnershipRange(v, &low, &high);
958  CHKERRQ(ierr);
959  ierr = ISCreateStride(((PetscObject)dm)->comm, (high - low), low, 1, &dmm->_embedding);
960  CHKERRQ(ierr);
961  }
962  }
963  ierr = PetscObjectReference((PetscObject)(dmm->_embedding));
964  CHKERRQ(ierr);
965  *embedding = dmm->_embedding;
966 
968 }
969 
970 #undef __FUNCT__
971 #define __FUNCT__ "DMCreateFieldDecomposition_Moose"
972 static PetscErrorCode
974  DM dm, PetscInt * len, char *** namelist, IS ** islist, DM ** dmlist)
975 {
976  PetscErrorCode ierr;
977  DM_Moose * dmm = (DM_Moose *)(dm->data);
978 
980  /* Only called after DMSetUp(). */
981  if (!dmm->_splitlocs)
983  *len = dmm->_splitlocs->size();
984  if (namelist)
985  {
986  ierr = PetscMalloc(*len * sizeof(char *), namelist);
987  CHKERRQ(ierr);
988  }
989  if (islist)
990  {
991  ierr = PetscMalloc(*len * sizeof(IS), islist);
992  CHKERRQ(ierr);
993  }
994  if (dmlist)
995  {
996  ierr = PetscMalloc(*len * sizeof(DM), dmlist);
997  CHKERRQ(ierr);
998  }
999  for (const auto & dit : *(dmm->_splitlocs))
1000  {
1001  unsigned int d = dit.second;
1002  std::string dname = dit.first;
1003  DM_Moose::SplitInfo & dinfo = (*dmm->_splits)[dname];
1004  if (!dinfo._dm)
1005  {
1006  ierr = DMCreateMoose(((PetscObject)dm)->comm, *dmm->_nl, &dinfo._dm);
1007  CHKERRQ(ierr);
1008  ierr = PetscObjectSetOptionsPrefix((PetscObject)dinfo._dm, ((PetscObject)dm)->prefix);
1009  CHKERRQ(ierr);
1010  std::string suffix = std::string("fieldsplit_") + dname + "_";
1011  ierr = PetscObjectAppendOptionsPrefix((PetscObject)dinfo._dm, suffix.c_str());
1012  CHKERRQ(ierr);
1013  }
1014  ierr = DMSetFromOptions(dinfo._dm);
1015  CHKERRQ(ierr);
1016  ierr = DMSetUp(dinfo._dm);
1017  CHKERRQ(ierr);
1018  if (namelist)
1019  {
1020  ierr = PetscStrallocpy(dname.c_str(), (*namelist) + d);
1021  CHKERRQ(ierr);
1022  }
1023  if (islist)
1024  {
1025  if (!dinfo._rembedding)
1026  {
1027  IS dembedding, lembedding;
1028  ierr = DMMooseGetEmbedding_Private(dinfo._dm, &dembedding);
1029  CHKERRQ(ierr);
1030  if (dmm->_embedding)
1031  {
1032 /* Create a relative embedding into the parent's index space. */
1033 #if PETSC_VERSION_LT(3, 4, 0)
1034  ierr = ISMapFactorRight(dembedding, dmm->_embedding, PETSC_TRUE, &lembedding);
1035  CHKERRQ(ierr);
1036 #else
1037  ierr = ISEmbed(dembedding, dmm->_embedding, PETSC_TRUE, &lembedding);
1038  CHKERRQ(ierr);
1039 #endif
1040  const PetscInt * lindices;
1041  PetscInt len, dlen, llen, *rindices, off, i;
1042  ierr = ISGetLocalSize(dembedding, &dlen);
1043  CHKERRQ(ierr);
1044  ierr = ISGetLocalSize(lembedding, &llen);
1045  CHKERRQ(ierr);
1046  if (llen != dlen)
1047  SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Failed to embed split %D", d);
1048  ierr = ISDestroy(&dembedding);
1049  CHKERRQ(ierr);
1050  // Convert local embedding to global (but still relative) embedding
1051  ierr = PetscMalloc(llen * sizeof(PetscInt), &rindices);
1052  CHKERRQ(ierr);
1053  ierr = ISGetIndices(lembedding, &lindices);
1054  CHKERRQ(ierr);
1055  ierr = PetscMemcpy(rindices, lindices, llen * sizeof(PetscInt));
1056  CHKERRQ(ierr);
1057  ierr = ISDestroy(&lembedding);
1058  CHKERRQ(ierr);
1059  // We could get the index offset from a corresponding global vector, but subDMs don't yet
1060  // have global vectors
1061  ierr = ISGetLocalSize(dmm->_embedding, &len);
1062  CHKERRQ(ierr);
1063 
1064  ierr = MPI_Scan(&len,
1065  &off,
1066  1,
1067 #ifdef PETSC_USE_64BIT_INDICES
1068  MPI_LONG_LONG_INT,
1069 #else
1070  MPI_INT,
1071 #endif
1072  MPI_SUM,
1073  ((PetscObject)dm)->comm);
1074  CHKERRQ(ierr);
1075 
1076  off -= len;
1077  for (i = 0; i < llen; ++i)
1078  rindices[i] += off;
1079  ierr = ISCreateGeneral(
1080  ((PetscObject)dm)->comm, llen, rindices, PETSC_OWN_POINTER, &(dinfo._rembedding));
1081  CHKERRQ(ierr);
1082  }
1083  else
1084  {
1085  dinfo._rembedding = dembedding;
1086  }
1087  }
1088  ierr = PetscObjectReference((PetscObject)(dinfo._rembedding));
1089  CHKERRQ(ierr);
1090  (*islist)[d] = dinfo._rembedding;
1091  }
1092  if (dmlist)
1093  {
1094  ierr = PetscObjectReference((PetscObject)dinfo._dm);
1095  CHKERRQ(ierr);
1096  (*dmlist)[d] = dinfo._dm;
1097  }
1098  }
1100 }
1101 
1102 #undef __FUNCT__
1103 #define __FUNCT__ "DMCreateDomainDecomposition_Moose"
1104 static PetscErrorCode
1106  DM dm, PetscInt * len, char *** namelist, IS ** innerislist, IS ** outerislist, DM ** dmlist)
1107 {
1108  PetscErrorCode ierr;
1109 
1111  /* Use DMCreateFieldDecomposition_Moose() to obtain everything but outerislist, which is currently
1112  * PETSC_NULL. */
1113  if (outerislist)
1114  *outerislist = PETSC_NULL; /* FIX: allow mesh-based overlap. */
1115  ierr = DMCreateFieldDecomposition_Moose(dm, len, namelist, innerislist, dmlist);
1116  CHKERRQ(ierr);
1118 }
1119 
1120 #if PETSC_VERSION_LT(3, 4, 0)
1121 #undef __FUNCT__
1122 #define __FUNCT__ "DMCreateFieldDecompositionDM_Moose"
1123 PetscErrorCode
1124 DMCreateFieldDecompositionDM_Moose(DM dm, const char * /*name*/, DM * ddm)
1125 {
1126  PetscErrorCode ierr;
1127  PetscBool ismoose;
1128 
1130  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1131  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
1132  CHKERRQ(ierr);
1133  /* Return self. */
1134  if (*ddm)
1135  {
1136  ierr = PetscObjectReference((PetscObject)dm);
1137  CHKERRQ(ierr);
1138  *ddm = dm;
1139  }
1141 }
1142 
1143 #undef __FUNCT__
1144 #define __FUNCT__ "DMCreateDomainDecompositionDM_Moose"
1145 PetscErrorCode
1146 DMCreateDomainDecompositionDM_Moose(DM dm, const char * /*name*/, DM * ddm)
1147 {
1148  PetscErrorCode ierr;
1149  PetscBool ismoose;
1150 
1152  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1153  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
1154  CHKERRQ(ierr);
1155  /* Return self. */
1156  if (*ddm)
1157  {
1158  ierr = PetscObjectReference((PetscObject)dm);
1159  CHKERRQ(ierr);
1160  *ddm = dm;
1161  }
1163 }
1164 #endif
1165 
1166 #undef __FUNCT__
1167 #define __FUNCT__ "DMMooseFunction"
1168 static PetscErrorCode
1169 DMMooseFunction(DM dm, Vec x, Vec r)
1170 {
1171  PetscErrorCode ierr;
1172 
1174  libmesh_assert(x);
1175  libmesh_assert(r);
1176 
1177  NonlinearSystemBase * nl = NULL;
1179  CHKERRQ(ierr);
1180  PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(nl->system().solution.get());
1181  PetscVector<Number> X_global(x, nl->comm()), R(r, nl->comm());
1182 
1183  // Use the system's update() to get a good local version of the
1184  // parallel solution. system.update() does change the residual vector,
1185  // so there's no reason to swap PETSc's residual into the system for
1186  // this step.
1187  X_global.swap(X_sys);
1188  nl->system().update();
1189  X_global.swap(X_sys);
1190 
1191  // Enforce constraints (if any) exactly on the
1192  // current_local_solution. This is the solution vector that is
1193  // actually used in the computation of the residual below, and is
1194  // not locked by debug-enabled PETSc the way that "x" is.
1195  nl->system().get_dof_map().enforce_constraints_exactly(nl->system(),
1196  nl->system().current_local_solution.get());
1197 
1198  // Zero the residual vector before assembling
1199  R.zero();
1200 
1201  // if the user has provided both function pointers and objects only the pointer
1202  // will be used, so catch that as an error
1203  if (nl->nonlinearSolver()->residual && nl->nonlinearSolver()->residual_object)
1204  {
1205  std::ostringstream err;
1206  err << "ERROR: cannot specifiy both a function and object to compute the Residual!"
1207  << std::endl;
1208  mooseError(err.str());
1209  }
1210  if (nl->nonlinearSolver()->matvec && nl->nonlinearSolver()->residual_and_jacobian_object)
1211  {
1212  std::ostringstream err;
1213  err << "ERROR: cannot specifiy both a function and object to compute the combined Residual & "
1214  "Jacobian!"
1215  << std::endl;
1216  mooseError(err.str());
1217  }
1218  if (nl->nonlinearSolver()->residual != NULL)
1219  nl->nonlinearSolver()->residual(
1220  *(nl->system().current_local_solution.get()), R, nl->nonlinearSolver()->system());
1221  else if (nl->nonlinearSolver()->residual_object != NULL)
1222  nl->nonlinearSolver()->residual_object->residual(
1223  *(nl->system().current_local_solution.get()), R, nl->nonlinearSolver()->system());
1224  else if (nl->nonlinearSolver()->matvec != NULL)
1225  nl->nonlinearSolver()->matvec(
1226  *(nl->system().current_local_solution.get()), &R, NULL, nl->nonlinearSolver()->system());
1227  else if (nl->nonlinearSolver()->residual_and_jacobian_object != NULL)
1228  nl->nonlinearSolver()->residual_and_jacobian_object->residual_and_jacobian(
1229  *(nl->system().current_local_solution.get()), &R, NULL, nl->nonlinearSolver()->system());
1230  else
1231  {
1232  std::ostringstream err;
1233  err << "No suitable residual computation routine found";
1234  mooseError(err.str());
1235  }
1236  R.close();
1238 }
1239 
1240 #if !PETSC_VERSION_LT(3, 4, 0)
1241 #undef __FUNCT__
1242 #define __FUNCT__ "SNESFunction_DMMoose"
1243 static PetscErrorCode
1244 SNESFunction_DMMoose(SNES, Vec x, Vec r, void * ctx)
1245 {
1246  DM dm = (DM)ctx;
1247  PetscErrorCode ierr;
1248 
1250  ierr = DMMooseFunction(dm, x, r);
1251  CHKERRQ(ierr);
1253 }
1254 #endif
1255 
1256 #undef __FUNCT__
1257 #define __FUNCT__ "DMMooseJacobian"
1258 #if PETSC_VERSION_LT(3, 5, 0)
1259 static PetscErrorCode
1260 DMMooseJacobian(DM dm, Vec x, Mat jac, Mat pc, MatStructure * msflag)
1261 #else
1262 static PetscErrorCode
1263 DMMooseJacobian(DM dm, Vec x, Mat jac, Mat pc)
1264 #endif
1266  PetscErrorCode ierr;
1268 
1271  CHKERRQ(ierr);
1272 
1273  PetscMatrix<Number> the_pc(pc, nl->comm());
1274  PetscMatrix<Number> Jac(jac, nl->comm());
1275  PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(nl->system().solution.get());
1276  PetscVector<Number> X_global(x, nl->comm());
1277 
1278  // Set the dof maps
1279  the_pc.attach_dof_map(nl->system().get_dof_map());
1280  Jac.attach_dof_map(nl->system().get_dof_map());
1281 
1282  // Use the system's update() to get a good local version of the
1283  // parallel solution. system.update() does change the Jacobian, so
1284  // there's no reason to swap PETSc's Jacobian into the system for
1285  // this step.
1286  X_global.swap(X_sys);
1287  nl->system().update();
1288  X_global.swap(X_sys);
1289 
1290  // Enforce constraints (if any) exactly on the
1291  // current_local_solution. This is the solution vector that is
1292  // actually used in the computation of the Jacobian below, and is
1293  // not locked by debug-enabled PETSc the way that "x" is.
1294  nl->system().get_dof_map().enforce_constraints_exactly(nl->system(),
1295  nl->system().current_local_solution.get());
1296 
1297  // Zero out the preconditioner before computing the Jacobian.
1298  the_pc.zero();
1299 
1300  // if the user has provided both function pointers and objects only the pointer
1301  // will be used, so catch that as an error
1302  if (nl->nonlinearSolver()->jacobian && nl->nonlinearSolver()->jacobian_object)
1303  {
1304  std::ostringstream err;
1305  err << "ERROR: cannot specifiy both a function and object to compute the Jacobian!"
1306  << std::endl;
1307  mooseError(err.str());
1308  }
1309  if (nl->nonlinearSolver()->matvec && nl->nonlinearSolver()->residual_and_jacobian_object)
1310  {
1311  std::ostringstream err;
1312  err << "ERROR: cannot specifiy both a function and object to compute the combined Residual & "
1313  "Jacobian!"
1314  << std::endl;
1315  mooseError(err.str());
1316  }
1317  if (nl->nonlinearSolver()->jacobian != NULL)
1318  nl->nonlinearSolver()->jacobian(
1319  *(nl->system().current_local_solution.get()), the_pc, nl->nonlinearSolver()->system());
1320  else if (nl->nonlinearSolver()->jacobian_object != NULL)
1321  nl->nonlinearSolver()->jacobian_object->jacobian(
1322  *(nl->system().current_local_solution.get()), the_pc, nl->nonlinearSolver()->system());
1323  else if (nl->nonlinearSolver()->matvec != NULL)
1324  nl->nonlinearSolver()->matvec(*(nl->system().current_local_solution.get()),
1325  NULL,
1326  &the_pc,
1327  nl->nonlinearSolver()->system());
1328  else if (nl->nonlinearSolver()->residual_and_jacobian_object != NULL)
1329  nl->nonlinearSolver()->residual_and_jacobian_object->residual_and_jacobian(
1330  *(nl->system().current_local_solution.get()),
1331  NULL,
1332  &the_pc,
1333  nl->nonlinearSolver()->system());
1334  else
1335  {
1336  std::ostringstream err;
1337  err << "No suitable Jacobian routine or object";
1338  mooseError(err.str());
1339  }
1340  the_pc.close();
1341  Jac.close();
1342 #if PETSC_VERSION_LT(3, 5, 0)
1343  *msflag = SAME_NONZERO_PATTERN;
1344 #endif
1346 }
1347 
1348 #if !PETSC_VERSION_LT(3, 4, 0)
1349 #undef __FUNCT__
1350 #define __FUNCT__ "SNESJacobian_DMMoose"
1351 #if PETSC_VERSION_LT(3, 5, 0)
1352 static PetscErrorCode
1353 SNESJacobian_DMMoose(SNES, Vec x, Mat * jac, Mat * pc, MatStructure * flag, void * ctx)
1354 #else
1355 static PetscErrorCode
1356 SNESJacobian_DMMoose(SNES, Vec x, Mat jac, Mat pc, void * ctx)
1357 #endif
1359  DM dm = (DM)ctx;
1360  PetscErrorCode ierr;
1361 
1363 #if PETSC_VERSION_LT(3, 5, 0)
1364  ierr = DMMooseJacobian(dm, x, *jac, *pc, flag);
1365  CHKERRQ(ierr);
1366 #else
1367  ierr = DMMooseJacobian(dm, x, jac, pc);
1368  CHKERRQ(ierr);
1369 #endif
1371 }
1372 #endif
1373 
1374 #undef __FUNCT__
1375 #define __FUNCT__ "DMVariableBounds_Moose"
1376 static PetscErrorCode
1377 DMVariableBounds_Moose(DM dm, Vec xl, Vec xu)
1378 {
1379  PetscErrorCode ierr;
1380  NonlinearSystemBase * nl = NULL;
1381 
1384  CHKERRQ(ierr);
1385 
1386  PetscVector<Number> XL(xl, nl->comm());
1387  PetscVector<Number> XU(xu, nl->comm());
1388 
1389 #if PETSC_VERSION_LESS_THAN(3, 5, 0) && PETSC_VERSION_RELEASE
1390  ierr = VecSet(xl, SNES_VI_NINF);
1391  CHKERRQ(ierr);
1392  ierr = VecSet(xu, SNES_VI_INF);
1393  CHKERRQ(ierr);
1394 #else
1395  ierr = VecSet(xl, PETSC_NINFINITY);
1396  CHKERRQ(ierr);
1397  ierr = VecSet(xu, PETSC_INFINITY);
1398  CHKERRQ(ierr);
1399 #endif
1400  if (nl->nonlinearSolver()->bounds != NULL)
1401  nl->nonlinearSolver()->bounds(XL, XU, nl->nonlinearSolver()->system());
1402  else if (nl->nonlinearSolver()->bounds_object != NULL)
1403  nl->nonlinearSolver()->bounds_object->bounds(XL, XU, nl->nonlinearSolver()->system());
1404  else
1405  SETERRQ(
1406  ((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "No bounds calculation in this Moose object");
1408 }
1409 
1410 #undef __FUNCT__
1411 #define __FUNCT__ "DMCreateGlobalVector_Moose"
1412 static PetscErrorCode
1414 {
1415  PetscErrorCode ierr;
1416  DM_Moose * dmm = (DM_Moose *)(dm->data);
1417  PetscBool ismoose;
1418 
1420  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
1421  CHKERRQ(ierr);
1422  if (!ismoose)
1423  SETERRQ2(((PetscObject)dm)->comm,
1424  PETSC_ERR_ARG_WRONG,
1425  "DM of type %s, not of type %s",
1426  ((PetscObject)dm)->type,
1427  DMMOOSE);
1428  if (!dmm->_nl)
1429  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1430 
1431  NumericVector<Number> * nv = (dmm->_nl->system().solution).get();
1432  PetscVector<Number> * pv = dynamic_cast<PetscVector<Number> *>(nv);
1433  Vec v = pv->vec();
1434  /* Unfortunately, currently this does not produce a ghosted vector, so nonlinear subproblem solves
1435  aren't going to be easily available.
1436  Should work fine for getting vectors out for linear subproblem solvers. */
1437  if (dmm->_embedding)
1438  {
1439  PetscInt n;
1440  ierr = VecCreate(((PetscObject)v)->comm, x);
1441  CHKERRQ(ierr);
1442  ierr = ISGetLocalSize(dmm->_embedding, &n);
1443  CHKERRQ(ierr);
1444  ierr = VecSetSizes(*x, n, PETSC_DETERMINE);
1445  CHKERRQ(ierr);
1446  ierr = VecSetType(*x, ((PetscObject)v)->type_name);
1447  CHKERRQ(ierr);
1448  ierr = VecSetFromOptions(*x);
1449  CHKERRQ(ierr);
1450  ierr = VecSetUp(*x);
1451  CHKERRQ(ierr);
1452  }
1453  else
1454  {
1455  ierr = VecDuplicate(v, x);
1456  CHKERRQ(ierr);
1457  }
1458  ierr = PetscObjectCompose((PetscObject)*x, "DM", (PetscObject)dm);
1459  CHKERRQ(ierr);
1461 }
1462 
1463 #undef __FUNCT__
1464 #define __FUNCT__ "DMCreateMatrix_Moose"
1465 #if PETSC_VERSION_LT(3, 5, 0)
1466 static PetscErrorCode
1467 DMCreateMatrix_Moose(DM dm, const MatType type, Mat * A)
1468 #else
1469 static PetscErrorCode
1470 DMCreateMatrix_Moose(DM dm, Mat * A)
1471 #endif
1473  PetscErrorCode ierr;
1474  DM_Moose * dmm = (DM_Moose *)(dm->data);
1475  PetscBool ismoose;
1476 #if !PETSC_RELEASE_LESS_THAN(3, 5, 0)
1477  MatType type;
1478 #endif
1479 
1481  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
1482  CHKERRQ(ierr);
1483  if (!ismoose)
1484  SETERRQ2(((PetscObject)dm)->comm,
1485  PETSC_ERR_ARG_WRONG,
1486  "DM of type %s, not of type %s",
1487  ((PetscObject)dm)->type,
1488  DMMOOSE);
1489  if (!dmm->_nl)
1490  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1491 // No PETSC_VERSION_GE macro prior to petsc-3.4
1492 #if !PETSC_VERSION_LT(3, 5, 0)
1493  ierr = DMGetMatType(dm, &type);
1494  CHKERRQ(ierr);
1495 #endif
1496  /*
1497  The simplest thing for now: compute the sparsity_pattern using dof_map and init the matrix using
1498  that info.
1499  TODO: compute sparsity restricted to this DM's blocks, variables and sides.
1500  Even fancier: compute the sparsity of the coupling of a contact slave to the contact master.
1501  In any event, here we are in control of the matrix type and structure.
1502  */
1503  DofMap & dof_map = dmm->_nl->system().get_dof_map();
1504  PetscInt M, N, m, n;
1505  MPI_Comm comm;
1506  M = dof_map.n_dofs();
1507  N = M;
1508  m = static_cast<PetscInt>(dof_map.n_dofs_on_processor(dmm->_nl->system().processor_id()));
1509  n = m;
1510  ierr = PetscObjectGetComm((PetscObject)dm, &comm);
1511  CHKERRQ(ierr);
1512  ierr = MatCreate(comm, A);
1513  CHKERRQ(ierr);
1514  ierr = MatSetSizes(*A, m, n, M, N);
1515  CHKERRQ(ierr);
1516  ierr = MatSetType(*A, type);
1517  CHKERRQ(ierr);
1518  /* Set preallocation for the basic sparse matrix types (applies only if *A has the right type. */
1519  /* For now we ignore blocksize issues, since BAIJ doesn't play well with field decomposition by
1520  * variable. */
1521  const std::vector<numeric_index_type> & n_nz = dof_map.get_n_nz();
1522  const std::vector<numeric_index_type> & n_oz = dof_map.get_n_oz();
1523  ierr = MatSeqAIJSetPreallocation(*A, 0, (PetscInt *)(n_nz.empty() ? NULL : &n_nz[0]));
1524  CHKERRQ(ierr);
1525  ierr = MatMPIAIJSetPreallocation(*A,
1526  0,
1527  (PetscInt *)(n_nz.empty() ? NULL : &n_nz[0]),
1528  0,
1529  (PetscInt *)(n_oz.empty() ? NULL : &n_oz[0]));
1530  CHKERRQ(ierr);
1531  /* TODO: set the prefix for *A and MatSetFromOptions(*A)? Might override the type and other
1532  * settings made here. */
1533  ierr = MatSetUp(*A);
1534  CHKERRQ(ierr);
1536 }
1537 
1538 #undef __FUNCT__
1539 #define __FUNCT__ "DMView_Moose"
1540 static PetscErrorCode
1541 DMView_Moose(DM dm, PetscViewer viewer)
1542 {
1543  PetscErrorCode ierr;
1544  PetscBool isascii;
1545  const char *name, *prefix;
1546  DM_Moose * dmm = (DM_Moose *)dm->data;
1547 
1549  ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
1550  CHKERRQ(ierr);
1551  if (isascii)
1552  {
1553  ierr = PetscObjectGetName((PetscObject)dm, &name);
1554  CHKERRQ(ierr);
1555  ierr = PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
1556  CHKERRQ(ierr);
1557  ierr = PetscViewerASCIIPrintf(viewer, "DM Moose with name %s and prefix %s\n", name, prefix);
1558  CHKERRQ(ierr);
1559  ierr = PetscViewerASCIIPrintf(viewer, "variables:", name, prefix);
1560  CHKERRQ(ierr);
1561  for (const auto & vit : *(dmm->_var_ids))
1562  {
1563  ierr = PetscViewerASCIIPrintf(viewer, "(%s,%D) ", vit.first.c_str(), vit.second);
1564  CHKERRQ(ierr);
1565  }
1566  ierr = PetscViewerASCIIPrintf(viewer, "\n");
1567  CHKERRQ(ierr);
1568  ierr = PetscViewerASCIIPrintf(viewer, "blocks:");
1569  CHKERRQ(ierr);
1570  for (const auto & bit : *(dmm->_block_ids))
1571  {
1572  ierr = PetscViewerASCIIPrintf(viewer, "(%s,%D) ", bit.first.c_str(), bit.second);
1573  CHKERRQ(ierr);
1574  }
1575  ierr = PetscViewerASCIIPrintf(viewer, "\n");
1576  CHKERRQ(ierr);
1577 
1578  if (dmm->_side_ids->size())
1579  {
1580  ierr = PetscViewerASCIIPrintf(viewer, "sides:");
1581  CHKERRQ(ierr);
1582  for (const auto & sit : *(dmm->_side_ids))
1583  {
1584  ierr = PetscViewerASCIIPrintf(viewer, "(%s,%D) ", sit.first.c_str(), sit.second);
1585  CHKERRQ(ierr);
1586  }
1587  ierr = PetscViewerASCIIPrintf(viewer, "\n");
1588  CHKERRQ(ierr);
1589  }
1590 
1591  if (dmm->_unside_ids->size())
1592  {
1593  ierr = PetscViewerASCIIPrintf(viewer, "unsides:");
1594  CHKERRQ(ierr);
1595  for (const auto & sit : *(dmm->_unside_ids))
1596  {
1597  ierr = PetscViewerASCIIPrintf(viewer, "(%s,%D) ", sit.first.c_str(), sit.second);
1598  CHKERRQ(ierr);
1599  }
1600  ierr = PetscViewerASCIIPrintf(viewer, "\n");
1601  CHKERRQ(ierr);
1602  }
1603 
1604  if (dmm->_contact_names->size())
1605  {
1606  ierr = PetscViewerASCIIPrintf(viewer, "contacts:");
1607  CHKERRQ(ierr);
1608  for (const auto & cit : *(dmm->_contact_names))
1609  {
1610  ierr = PetscViewerASCIIPrintf(
1611  viewer, "(%s,%s,", cit.second.first.c_str(), cit.second.second.c_str());
1612  CHKERRQ(ierr);
1613  if ((*dmm->_contact_displaced)[cit.second])
1614  {
1615  ierr = PetscViewerASCIIPrintf(
1616  viewer, "displaced) ", cit.second.first.c_str(), cit.second.second.c_str());
1617  CHKERRQ(ierr);
1618  }
1619  else
1620  {
1621  ierr = PetscViewerASCIIPrintf(
1622  viewer, "undisplaced) ", cit.second.first.c_str(), cit.second.second.c_str());
1623  CHKERRQ(ierr);
1624  }
1625  }
1626  ierr = PetscViewerASCIIPrintf(viewer, "\n");
1627  CHKERRQ(ierr);
1628  }
1629 
1630  if (dmm->_uncontact_names->size())
1631  {
1632  ierr = PetscViewerASCIIPrintf(viewer, "_uncontacts:");
1633  CHKERRQ(ierr);
1634  for (const auto & cit : *(dmm->_uncontact_names))
1635  {
1636  ierr = PetscViewerASCIIPrintf(
1637  viewer, "(%s,%s,", cit.second.first.c_str(), cit.second.second.c_str());
1638  CHKERRQ(ierr);
1639  if ((*dmm->_uncontact_displaced)[cit.second])
1640  {
1641  ierr = PetscViewerASCIIPrintf(
1642  viewer, "displaced) ", cit.second.first.c_str(), cit.second.second.c_str());
1643  CHKERRQ(ierr);
1644  }
1645  else
1646  {
1647  ierr = PetscViewerASCIIPrintf(
1648  viewer, "undisplaced) ", cit.second.first.c_str(), cit.second.second.c_str());
1649  CHKERRQ(ierr);
1650  }
1651  }
1652  ierr = PetscViewerASCIIPrintf(viewer, "\n");
1653  CHKERRQ(ierr);
1654  }
1655 
1656  if (dmm->_splitlocs && dmm->_splitlocs->size())
1657  {
1658  ierr = PetscViewerASCIIPrintf(viewer, "Field decomposition:");
1659  CHKERRQ(ierr);
1660  // FIX: decompositions might have different sizes and components on different ranks.
1661  for (const auto & dit : *(dmm->_splitlocs))
1662  {
1663  std::string dname = dit.first;
1664  ierr = PetscViewerASCIIPrintf(viewer, " %s", dname.c_str());
1665  CHKERRQ(ierr);
1666  }
1667  ierr = PetscViewerASCIIPrintf(viewer, "\n");
1668  CHKERRQ(ierr);
1669  }
1670  }
1671  else
1672  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Non-ASCII viewers are not supported");
1673 
1675 }
1676 
1677 #undef __FUNCT__
1678 #define __FUNCT__ "DMMooseGetMeshBlocks_Private"
1679 static PetscErrorCode
1680 DMMooseGetMeshBlocks_Private(DM dm, std::set<subdomain_id_type> & blocks)
1681 {
1682  PetscErrorCode ierr;
1683  DM_Moose * dmm = (DM_Moose *)(dm->data);
1684  PetscBool ismoose;
1685 
1687  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
1688  CHKERRQ(ierr);
1689  if (!ismoose)
1690  SETERRQ2(((PetscObject)dm)->comm,
1691  PETSC_ERR_ARG_WRONG,
1692  "DM of type %s, not of type %s",
1693  ((PetscObject)dm)->type,
1694  DMMOOSE);
1695  if (!dmm->_nl)
1696  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1697 
1698  const MeshBase & mesh = dmm->_nl->system().get_mesh();
1699  /* The following effectively is a verbatim copy of MeshBase::n_subdomains(). */
1700  // This requires an inspection on every processor
1701  libmesh_parallel_only(mesh.comm());
1702  for (const auto & elem : mesh.active_element_ptr_range())
1703  blocks.insert(elem->subdomain_id());
1704  // Some subdomains may only live on other processors
1705  mesh.comm().set_union(blocks);
1707 }
1708 
1709 #undef __FUNCT__
1710 #define __FUNCT__ "DMSetUp_Moose_Pre"
1711 static PetscErrorCode
1713 {
1714  PetscErrorCode ierr;
1715  DM_Moose * dmm = (DM_Moose *)(dm->data);
1716  PetscBool ismoose;
1717 
1719  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
1720  CHKERRQ(ierr);
1721  if (!ismoose)
1722  SETERRQ2(((PetscObject)dm)->comm,
1723  PETSC_ERR_ARG_WRONG,
1724  "DM of type %s, not of type %s",
1725  ((PetscObject)dm)->type,
1726  DMMOOSE);
1727  if (!dmm->_nl)
1728  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1729 
1730  /* Set up variables, blocks and sides. */
1731  DofMap & dofmap = dmm->_nl->system().get_dof_map();
1732  /* libMesh mesh */
1733  const MeshBase & mesh = dmm->_nl->system().get_mesh();
1734 
1735  dmm->_nosides = PETSC_TRUE;
1736  dmm->_side_ids->clear();
1737  dmm->_side_names->clear();
1738  if (dmm->_sides)
1739  {
1740  dmm->_nosides = PETSC_FALSE;
1741  std::set<BoundaryID> ids;
1742  for (const auto & name : *(dmm->_sides))
1743  {
1744  boundary_id_type id = dmm->_nl->mesh().getBoundaryID(name);
1745  dmm->_side_names->insert(std::make_pair(id, name));
1746  dmm->_side_ids->insert(std::make_pair(name, id));
1747  }
1748  delete dmm->_sides;
1749  dmm->_sides = PETSC_NULL;
1750  }
1751  dmm->_nounsides = PETSC_TRUE;
1752  dmm->_unside_ids->clear();
1753  dmm->_unside_names->clear();
1754  if (dmm->_unsides)
1755  {
1756  dmm->_nounsides = PETSC_FALSE;
1757  std::set<BoundaryID> ids;
1758  for (const auto & name : *(dmm->_unsides))
1759  {
1760  boundary_id_type id = dmm->_nl->mesh().getBoundaryID(name);
1761  dmm->_unside_names->insert(std::make_pair(id, name));
1762  dmm->_unside_ids->insert(std::make_pair(name, id));
1763  }
1764  delete dmm->_unsides;
1765  dmm->_unsides = PETSC_NULL;
1766  }
1767  dmm->_nocontacts = PETSC_TRUE;
1768 
1769  if (dmm->_contacts)
1770  {
1771  dmm->_nocontacts = PETSC_FALSE;
1772  for (const auto & cpair : *(dmm->_contacts))
1773  {
1774  try
1775  {
1776  if ((*dmm->_contact_displaced)[cpair])
1777  dmm->_nl->_fe_problem.getDisplacedProblem()->geomSearchData().getPenetrationLocator(
1778  cpair.first, cpair.second);
1779  else
1780  dmm->_nl->_fe_problem.geomSearchData().getPenetrationLocator(cpair.first, cpair.second);
1781  }
1782  catch (...)
1783  {
1784  std::ostringstream err;
1785  err << "Problem retrieving contact for PenetrationLocator with master " << cpair.first
1786  << " and slave " << cpair.second;
1787  mooseError(err.str());
1788  }
1789  BoundaryID master_id = dmm->_nl->mesh().getBoundaryID(cpair.first);
1790  BoundaryID slave_id = dmm->_nl->mesh().getBoundaryID(cpair.second);
1791  DM_Moose::ContactID cid(master_id, slave_id);
1792  dmm->_contact_names->insert(std::make_pair(cid, cpair));
1793  }
1794  }
1795 
1796  dmm->_nouncontacts = PETSC_TRUE;
1797  if (dmm->_uncontacts)
1798  {
1799  dmm->_nouncontacts = PETSC_FALSE;
1800  for (const auto & cpair : *(dmm->_uncontacts))
1801  {
1802  try
1803  {
1804  if ((*dmm->_uncontact_displaced)[cpair])
1805  dmm->_nl->_fe_problem.getDisplacedProblem()->geomSearchData().getPenetrationLocator(
1806  cpair.first, cpair.second);
1807  else
1808  dmm->_nl->_fe_problem.geomSearchData().getPenetrationLocator(cpair.first, cpair.second);
1809  }
1810  catch (...)
1811  {
1812  std::ostringstream err;
1813  err << "Problem retrieving uncontact for PenetrationLocator with master " << cpair.first
1814  << " and slave " << cpair.second;
1815  mooseError(err.str());
1816  }
1817  BoundaryID master_id = dmm->_nl->mesh().getBoundaryID(cpair.first);
1818  BoundaryID slave_id = dmm->_nl->mesh().getBoundaryID(cpair.second);
1819  DM_Moose::ContactID cid(master_id, slave_id);
1820  dmm->_uncontact_names->insert(std::make_pair(cid, cpair));
1821  }
1822  }
1823 
1824  dmm->_var_ids->clear();
1825  dmm->_var_names->clear();
1826  // FIXME: would be nice to invert this nested loop structure so we could iterate over the
1827  // potentially smaller dmm->vars,
1828  // but checking against dofmap.variable would still require a linear search, hence, no win. Would
1829  // be nice to endow dofmap.variable
1830  // with a fast search capability.
1831  for (unsigned int v = 0; v < dofmap.n_variables(); ++v)
1832  {
1833  std::string vname = dofmap.variable(v).name();
1834  if (dmm->_vars && dmm->_vars->size() && dmm->_vars->find(vname) == dmm->_vars->end())
1835  continue;
1836  dmm->_var_ids->insert(std::pair<std::string, unsigned int>(vname, v));
1837  dmm->_var_names->insert(std::pair<unsigned int, std::string>(v, vname));
1838  }
1839  if (dmm->_var_ids->size() == dofmap.n_variables())
1840  dmm->_all_vars = PETSC_TRUE;
1841  else
1842  dmm->_all_vars = PETSC_FALSE;
1843  if (dmm->_vars)
1844  {
1845  delete dmm->_vars;
1846  dmm->_vars = PETSC_NULL;
1847  }
1848 
1849  dmm->_block_ids->clear();
1850  dmm->_block_names->clear();
1851  std::set<subdomain_id_type> blocks;
1852  ierr = DMMooseGetMeshBlocks_Private(dm, blocks);
1853  CHKERRQ(ierr);
1854  if (blocks.empty())
1855  SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "No mesh blocks found.");
1856 
1857  for (const auto & bid : blocks)
1858  {
1859  std::string bname = mesh.subdomain_name(bid);
1860  if (!bname.length())
1861  {
1862  // Block names are currently implemented for Exodus II meshes
1863  // only, so we might have to make up our own block names and
1864  // maintain our own mapping of block ids to names.
1865  std::ostringstream ss;
1866  ss << bid;
1867  bname = ss.str();
1868  }
1869  if (dmm->_nosides)
1870  {
1871  // If no sides have been specified, by default (null or empty dmm->blocks) all blocks are
1872  // included in the split
1873  // Thus, skip this block only if it is explicitly excluded from a nonempty dmm->blocks.
1874  if (dmm->_blocks && dmm->_blocks->size() && dmm->_blocks->find(bname) == dmm->_blocks->end())
1875  continue;
1876  }
1877  else
1878  {
1879  // If sides have been specified, only the explicitly-specified blocks (those in dmm->blocks,
1880  // if it's non-null) are in the split.
1881  // Thus, include this block only if it is explicitly specified in a nonempty dmm->blocks.
1882  // Equivalently, skip this block if dmm->blocks is dmm->blocks is null or empty or excludes
1883  // this block.
1884  if (!dmm->_blocks || !dmm->_blocks->size() ||
1885  dmm->_blocks->find(bname) == dmm->_blocks->end())
1886  continue;
1887  }
1888  dmm->_block_ids->insert(std::make_pair(bname, bid));
1889  dmm->_block_names->insert(std::make_pair(bid, bname));
1890  }
1891 
1892  if (dmm->_block_ids->size() == blocks.size())
1893  dmm->_all_blocks = PETSC_TRUE;
1894  else
1895  dmm->_all_blocks = PETSC_FALSE;
1896  if (dmm->_blocks)
1897  {
1898  delete dmm->_blocks;
1899  dmm->_blocks = PETSC_NULL;
1900  }
1901 
1902  std::string name = dmm->_nl->system().name();
1903  name += "_vars";
1904  for (const auto & vit : *(dmm->_var_names))
1905  name += "_" + vit.second;
1906 
1907  name += "_blocks";
1908 
1909  for (const auto & bit : *(dmm->_block_names))
1910  name += "_" + bit.second;
1911 
1912  if (dmm->_side_names && dmm->_side_names->size())
1913  {
1914  name += "_sides";
1915  for (const auto & sit : *(dmm->_side_names))
1916  name += "_" + sit.second;
1917  }
1918  if (dmm->_unside_names && dmm->_unside_names->size())
1919  {
1920  name += "_unsides";
1921  for (const auto & sit : *(dmm->_unside_names))
1922  name += "_" + sit.second;
1923  }
1924  if (dmm->_contact_names && dmm->_contact_names->size())
1925  {
1926  name += "_contacts";
1927  for (const auto & cit : *(dmm->_contact_names))
1928  name += "_master_" + cit.second.first + "_slave_" + cit.second.second;
1929  }
1930  if (dmm->_uncontact_names && dmm->_uncontact_names->size())
1931  {
1932  name += "_uncontacts";
1933  for (const auto & cit : *(dmm->_uncontact_names))
1934  name += "_master_" + cit.second.first + "_slave_" + cit.second.second;
1935  }
1936  ierr = PetscObjectSetName((PetscObject)dm, name.c_str());
1937  CHKERRQ(ierr);
1939 }
1940 
1941 #undef __FUNCT__
1942 #define __FUNCT__ "DMMooseReset"
1943 PetscErrorCode
1945 {
1946  PetscErrorCode ierr;
1947  DM_Moose * dmm = (DM_Moose *)(dm->data);
1948  PetscBool ismoose;
1949 
1951  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
1952  CHKERRQ(ierr);
1953  if (!ismoose)
1955  if (!dmm->_nl)
1956  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1957  ierr = ISDestroy(&dmm->_embedding);
1958  CHKERRQ(ierr);
1959  for (auto & it : *(dmm->_splits))
1960  {
1961  DM_Moose::SplitInfo & split = it.second;
1962  ierr = ISDestroy(&split._rembedding);
1963  CHKERRQ(ierr);
1964  if (split._dm)
1965  {
1966  ierr = DMMooseReset(split._dm);
1967  CHKERRQ(ierr);
1968  }
1969  }
1970  dm->setupcalled = PETSC_FALSE;
1972 }
1973 
1974 #undef __FUNCT__
1975 #define __FUNCT__ "DMSetUp_Moose"
1976 static PetscErrorCode
1978 {
1979  PetscErrorCode ierr;
1980  DM_Moose * dmm = (DM_Moose *)(dm->data);
1981  PetscBool ismoose;
1982 
1984  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
1985  CHKERRQ(ierr);
1986  if (!ismoose)
1987  SETERRQ2(((PetscObject)dm)->comm,
1988  PETSC_ERR_ARG_WRONG,
1989  "DM of type %s, not of type %s",
1990  ((PetscObject)dm)->type,
1991  DMMOOSE);
1992  if (!dmm->_nl)
1993  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
1994  if (dmm->_print_embedding)
1995  {
1996  const char *name, *prefix;
1997  IS embedding;
1998 
1999  ierr = PetscObjectGetName((PetscObject)dm, &name);
2000  CHKERRQ(ierr);
2001  ierr = PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
2002  CHKERRQ(ierr);
2003  ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)dm)->comm),
2004  "DM Moose with name %s and prefix %s\n",
2005  name,
2006  prefix);
2007  CHKERRQ(ierr);
2008  if (dmm->_all_vars && dmm->_all_blocks && dmm->_nosides && dmm->_nounsides &&
2010  {
2011  ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)dm)->comm),
2012  "\thas a trivial embedding\n");
2013  CHKERRQ(ierr);
2014  }
2015  else
2016  {
2017  ierr = DMMooseGetEmbedding_Private(dm, &embedding);
2018  CHKERRQ(ierr);
2019  ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)dm)->comm),
2020  "\thas embedding defined by IS:\n");
2021  CHKERRQ(ierr);
2022  ierr = ISView(embedding, PETSC_VIEWER_STDOUT_(((PetscObject)dm)->comm));
2023  CHKERRQ(ierr);
2024  ierr = ISDestroy(&embedding);
2025  CHKERRQ(ierr);
2026  }
2027  }
2028  /*
2029  Do not evaluate function, Jacobian or bounds for an embedded DM -- the subproblem might not have
2030  enough information for that.
2031  */
2033  dmm->_nouncontacts)
2034  {
2035 #if PETSC_VERSION_LT(3, 4, 0)
2036  ierr = DMSetFunction(dm, DMMooseFunction);
2037  CHKERRQ(ierr);
2038  ierr = DMSetJacobian(dm, DMMooseJacobian);
2039  CHKERRQ(ierr);
2040 #else
2041  ierr = DMSNESSetFunction(dm, SNESFunction_DMMoose, (void *)dm);
2042  CHKERRQ(ierr);
2043  ierr = DMSNESSetJacobian(dm, SNESJacobian_DMMoose, (void *)dm);
2044  CHKERRQ(ierr);
2045 #endif
2046  if (dmm->_nl->nonlinearSolver()->bounds || dmm->_nl->nonlinearSolver()->bounds_object)
2047  ierr = DMSetVariableBounds(dm, DMVariableBounds_Moose);
2048  CHKERRQ(ierr);
2049  }
2050  else
2051  {
2052  /*
2053  Fow now we don't implement even these, although a linear "Dirichlet" subproblem is
2054  well-defined.
2055  Creating the submatrix, however, might require extracting the submatrix preallocation from an
2056  unassembled matrix.
2057  */
2058  dm->ops->createglobalvector = 0;
2059  dm->ops->creatematrix = 0;
2060  }
2062 }
2063 
2064 #undef __FUNCT__
2065 #define __FUNCT__ "DMSetFromOptions_Moose"
2066 #if !PETSC_VERSION_LESS_THAN(3, 7, 0)
2067 PetscErrorCode DMSetFromOptions_Moose(PetscOptionItems * /*options*/, DM dm) // >= 3.7.0
2068 #elif !PETSC_RELEASE_LESS_THAN(3, 6, 0)
2069 PetscErrorCode DMSetFromOptions_Moose(PetscOptions * /*options*/, DM dm) // >= 3.6.0
2070 #else
2071 PetscErrorCode DMSetFromOptions_Moose(DM dm) // < 3.6.0
2072 #endif
2073 {
2074  PetscErrorCode ierr;
2075  PetscBool ismoose;
2076  DM_Moose * dmm = (DM_Moose *)dm->data;
2077 
2079  ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
2080  CHKERRQ(ierr);
2081  if (!ismoose)
2082  SETERRQ2(((PetscObject)dm)->comm,
2083  PETSC_ERR_ARG_WRONG,
2084  "DM of type %s, not of type %s",
2085  ((PetscObject)dm)->type,
2086  DMMOOSE);
2087  if (!dmm->_nl)
2088  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
2089  ierr = PetscOptionsBegin(
2090  ((PetscObject)dm)->comm, ((PetscObject)dm)->prefix, "DMMoose options", "DM");
2091  CHKERRQ(ierr);
2092  std::string opt, help;
2093  PetscInt maxvars = dmm->_nl->system().get_dof_map().n_variables();
2094  char ** vars;
2095  std::set<std::string> varset;
2096  PetscInt nvars = maxvars;
2097  ierr = PetscMalloc(maxvars * sizeof(char *), &vars);
2098  CHKERRQ(ierr);
2099  opt = "-dm_moose_vars";
2100  help = "Variables in DMMoose";
2101  ierr = PetscOptionsStringArray(
2102  opt.c_str(), help.c_str(), "DMMooseSetVars", vars, &nvars, PETSC_NULL);
2103  CHKERRQ(ierr);
2104  for (PetscInt i = 0; i < nvars; ++i)
2105  {
2106  varset.insert(std::string(vars[i]));
2107  ierr = PetscFree(vars[i]);
2108  CHKERRQ(ierr);
2109  }
2110  ierr = PetscFree(vars);
2111  CHKERRQ(ierr);
2112  if (varset.size())
2113  {
2114  ierr = DMMooseSetVariables(dm, varset);
2115  CHKERRQ(ierr);
2116  }
2117  //
2118  std::set<subdomain_id_type> meshblocks;
2119  ierr = DMMooseGetMeshBlocks_Private(dm, meshblocks);
2120  CHKERRQ(ierr);
2121  PetscInt maxblocks = meshblocks.size();
2122  char ** blocks;
2123  ierr = PetscMalloc(maxblocks * sizeof(char *), &blocks);
2124  CHKERRQ(ierr);
2125  std::set<std::string> blockset;
2126  PetscInt nblocks = maxblocks;
2127  opt = "-dm_moose_blocks";
2128  help = "Blocks in DMMoose";
2129  ierr = PetscOptionsStringArray(
2130  opt.c_str(), help.c_str(), "DMMooseSetBlocks", blocks, &nblocks, PETSC_NULL);
2131  CHKERRQ(ierr);
2132  for (PetscInt i = 0; i < nblocks; ++i)
2133  {
2134  blockset.insert(std::string(blocks[i]));
2135  ierr = PetscFree(blocks[i]);
2136  CHKERRQ(ierr);
2137  }
2138  ierr = PetscFree(blocks);
2139  CHKERRQ(ierr);
2140  if (blockset.size())
2141  {
2142  ierr = DMMooseSetBlocks(dm, blockset);
2143  CHKERRQ(ierr);
2144  }
2145  PetscInt maxsides = dmm->_nl->system().get_mesh().get_boundary_info().get_boundary_ids().size();
2146  char ** sides;
2147  ierr = PetscMalloc(maxsides * sizeof(char *), &sides);
2148  CHKERRQ(ierr);
2149  PetscInt nsides = maxsides;
2150  std::set<std::string> sideset;
2151  opt = "-dm_moose_sides";
2152  help = "Sides to include in DMMoose";
2153  ierr = PetscOptionsStringArray(
2154  opt.c_str(), help.c_str(), "DMMooseSetSides", sides, &nsides, PETSC_NULL);
2155  CHKERRQ(ierr);
2156  for (PetscInt i = 0; i < nsides; ++i)
2157  {
2158  sideset.insert(std::string(sides[i]));
2159  ierr = PetscFree(sides[i]);
2160  CHKERRQ(ierr);
2161  }
2162  if (sideset.size())
2163  {
2164  ierr = DMMooseSetSides(dm, sideset);
2165  CHKERRQ(ierr);
2166  }
2167  opt = "-dm_moose_unsides";
2168  help = "Sides to exclude from DMMoose";
2169  nsides = maxsides;
2170  ierr = PetscOptionsStringArray(
2171  opt.c_str(), help.c_str(), "DMMooseSetUnSides", sides, &nsides, PETSC_NULL);
2172  CHKERRQ(ierr);
2173  sideset.clear();
2174  for (PetscInt i = 0; i < nsides; ++i)
2175  {
2176  sideset.insert(std::string(sides[i]));
2177  ierr = PetscFree(sides[i]);
2178  CHKERRQ(ierr);
2179  }
2180  if (sideset.size())
2181  {
2182  ierr = DMMooseSetUnSides(dm, sideset);
2183  CHKERRQ(ierr);
2184  }
2185  ierr = PetscFree(sides);
2186  CHKERRQ(ierr);
2187  PetscInt maxcontacts = dmm->_nl->_fe_problem.geomSearchData()._penetration_locators.size();
2188  std::shared_ptr<DisplacedProblem> displaced_problem = dmm->_nl->_fe_problem.getDisplacedProblem();
2189  if (displaced_problem)
2190  maxcontacts = PetscMax(
2191  maxcontacts, (PetscInt)displaced_problem->geomSearchData()._penetration_locators.size());
2192 
2193  std::vector<DM_Moose::ContactName> contacts;
2194  std::vector<PetscBool> contact_displaced;
2195  PetscInt ncontacts = 0;
2196  opt = "-dm_moose_ncontacts";
2197  help = "Number of contacts to include in DMMoose. For each <n> < "
2198  "dm_moose_contacts\n\t-dm_moose_contact_<n> is a comma-separated <master>,<slave> pair "
2199  "defining the contact surfaces"
2200  "\t-dm_moose_contact_<n>_displaced <bool> determines whether the contact is defined on "
2201  "the displaced mesh or not";
2202  ierr = PetscOptionsInt(
2203  opt.c_str(), help.c_str(), "DMMooseSetContacts", ncontacts, &ncontacts, PETSC_NULL);
2204  CHKERRQ(ierr);
2205  if (ncontacts > maxcontacts)
2206  SETERRQ2(((PetscObject)dm)->comm,
2207  PETSC_ERR_ARG_SIZ,
2208  "Number of requested contacts %D exceeds the maximum number of contacts %D",
2209  ncontacts,
2210  maxcontacts);
2211  for (PetscInt i = 0; i < ncontacts; ++i)
2212  {
2213  {
2214  char * master_slave[2];
2215  PetscInt sz = 2;
2216  std::ostringstream oopt, ohelp;
2217  oopt << "-dm_moose_contact_" << i;
2218  ohelp << "Master and slave for contact " << i;
2219  ierr = PetscOptionsStringArray(oopt.str().c_str(),
2220  ohelp.str().c_str(),
2221  "DMMooseSetContacts",
2222  master_slave,
2223  &sz,
2224  PETSC_NULL);
2225  CHKERRQ(ierr);
2226  if (sz != 2)
2227  SETERRQ2(((PetscObject)dm)->comm,
2228  PETSC_ERR_ARG_SIZ,
2229  "Expected 2 sideset IDs (master & slave) for contact %D, got %D instead",
2230  i,
2231  sz);
2232  contacts.push_back(
2233  DM_Moose::ContactName(std::string(master_slave[0]), std::string(master_slave[1])));
2234  ierr = PetscFree(master_slave[0]);
2235  CHKERRQ(ierr);
2236  ierr = PetscFree(master_slave[1]);
2237  CHKERRQ(ierr);
2238  }
2239  {
2240  PetscBool displaced = PETSC_FALSE;
2241  std::ostringstream oopt, ohelp;
2242  oopt << "-dm_moose_contact_" << i << "_displaced";
2243  ohelp << "Whether contact " << i << " is determined using displaced mesh or not";
2244  ierr = PetscOptionsBool(oopt.str().c_str(),
2245  ohelp.str().c_str(),
2246  "DMMooseSetContacts",
2247  PETSC_FALSE,
2248  &displaced,
2249  PETSC_NULL);
2250  CHKERRQ(ierr);
2251  contact_displaced.push_back(displaced);
2252  }
2253  }
2254  if (contacts.size())
2255  {
2256  ierr = DMMooseSetContacts(dm, contacts, contact_displaced);
2257  CHKERRQ(ierr);
2258  }
2259  {
2260  std::ostringstream oopt, ohelp;
2261  PetscBool is_include_all_nodes;
2262  oopt << "-dm_moose_includeAllContactNodes";
2263  ohelp << "Whether to include all nodes on the contact surfaces into the subsolver";
2264  ierr = PetscOptionsBool(oopt.str().c_str(),
2265  ohelp.str().c_str(),
2266  "",
2267  PETSC_FALSE,
2268  &is_include_all_nodes,
2269  PETSC_NULL);
2270  CHKERRQ(ierr);
2271  dmm->_include_all_contact_nodes = is_include_all_nodes;
2272  }
2273  std::vector<DM_Moose::ContactName> uncontacts;
2274  std::vector<PetscBool> uncontact_displaced;
2275  PetscInt nuncontacts = 0;
2276  opt = "-dm_moose_nuncontacts";
2277  help = "Number of contacts to exclude from DMMoose. For each <n> < "
2278  "dm_moose_contacts\n\t-dm_moose_contact_<n> is a comma-separated <master>,<slave> pair "
2279  "defining the contact surfaces"
2280  "\t-dm_moose_contact_<n>_displaced <bool> determines whether the contact is defined on "
2281  "the displaced mesh or not";
2282  ierr = PetscOptionsInt(
2283  opt.c_str(), help.c_str(), "DMMooseSetUnContacts", nuncontacts, &nuncontacts, PETSC_NULL);
2284  CHKERRQ(ierr);
2285  if (nuncontacts > maxcontacts)
2286  SETERRQ2(((PetscObject)dm)->comm,
2287  PETSC_ERR_ARG_SIZ,
2288  "Number of requested uncontacts %D exceeds the maximum number of contacts %D",
2289  nuncontacts,
2290  maxcontacts);
2291  for (PetscInt i = 0; i < nuncontacts; ++i)
2292  {
2293  {
2294  char * master_slave[2];
2295  PetscInt sz = 2;
2296  std::ostringstream oopt, ohelp;
2297  oopt << "-dm_moose_uncontact_" << i;
2298  ohelp << "Master and slave for uncontact " << i;
2299  ierr = PetscOptionsStringArray(oopt.str().c_str(),
2300  ohelp.str().c_str(),
2301  "DMMooseSetUnContacts",
2302  master_slave,
2303  &sz,
2304  PETSC_NULL);
2305  CHKERRQ(ierr);
2306  if (sz != 2)
2307  SETERRQ2(((PetscObject)dm)->comm,
2308  PETSC_ERR_ARG_SIZ,
2309  "Expected 2 sideset IDs (master & slave) for uncontact %D, got %D instead",
2310  i,
2311  sz);
2312  uncontacts.push_back(
2313  DM_Moose::ContactName(std::string(master_slave[0]), std::string(master_slave[1])));
2314  ierr = PetscFree(master_slave[0]);
2315  CHKERRQ(ierr);
2316  ierr = PetscFree(master_slave[1]);
2317  CHKERRQ(ierr);
2318  }
2319  {
2320  PetscBool displaced = PETSC_FALSE;
2321  std::ostringstream oopt, ohelp;
2322  oopt << "-dm_moose_uncontact_" << i << "_displaced";
2323  ohelp << "Whether uncontact " << i << " is determined using displaced mesh or not";
2324  ierr = PetscOptionsBool(oopt.str().c_str(),
2325  ohelp.str().c_str(),
2326  "DMMooseSetUnContact",
2327  PETSC_FALSE,
2328  &displaced,
2329  PETSC_NULL);
2330  CHKERRQ(ierr);
2331  uncontact_displaced.push_back(displaced);
2332  }
2333  }
2334  if (uncontacts.size())
2335  {
2336  ierr = DMMooseSetUnContacts(dm, uncontacts, uncontact_displaced);
2337  CHKERRQ(ierr);
2338  }
2339 
2340  PetscInt nsplits = 0;
2341  /* Insert the usage of -dm_moose_fieldsplit_names into this help message, since the following
2342  * if-clause might never fire, if -help is requested. */
2343  const char * fdhelp = "Number of named fieldsplits defined by the DM.\n\
2344  \tNames of fieldsplits are defined by -dm_moose_fieldsplit_names <splitname1> <splitname2> ...\n\
2345  \tEach split can be configured with its own variables, blocks and sides, as any DMMoose";
2346  ierr = PetscOptionsInt(
2347  "-dm_moose_nfieldsplits", fdhelp, "DMMooseSetSplitNames", nsplits, &nsplits, NULL);
2348  CHKERRQ(ierr);
2349  if (nsplits)
2350  {
2351  PetscInt nnsplits = nsplits;
2352  std::vector<std::string> split_names;
2353  char ** splitnames;
2354  ierr = PetscMalloc(nsplits * sizeof(char *), &splitnames);
2355  CHKERRQ(ierr);
2356  ierr = PetscOptionsStringArray("-dm_moose_fieldsplit_names",
2357  "Names of fieldsplits defined by the DM",
2358  "DMMooseSetSplitNames",
2359  splitnames,
2360  &nnsplits,
2361  PETSC_NULL);
2362  CHKERRQ(ierr);
2363  if (!nnsplits)
2364  {
2365  for (PetscInt i = 0; i < nsplits; ++i)
2366  {
2367  std::ostringstream s;
2368  s << i;
2369  split_names.push_back(s.str());
2370  }
2371  }
2372  else if (nsplits != nnsplits)
2373  SETERRQ2(((PetscObject)dm)->comm,
2374  PETSC_ERR_ARG_SIZ,
2375  "Expected %D fieldsplit names, got %D instead",
2376  nsplits,
2377  nnsplits);
2378  else
2379  {
2380  for (PetscInt i = 0; i < nsplits; ++i)
2381  {
2382  split_names.push_back(std::string(splitnames[i]));
2383  ierr = PetscFree(splitnames[i]);
2384  CHKERRQ(ierr);
2385  }
2386  }
2387  ierr = PetscFree(splitnames);
2388  CHKERRQ(ierr);
2389  ierr = DMMooseSetSplitNames(dm, split_names);
2390  CHKERRQ(ierr);
2391  }
2392  ierr = PetscOptionsBool("-dm_moose_print_embedding",
2393  "Print IS embedding DM's dofs",
2394  "DMMoose",
2397  PETSC_NULL);
2398  CHKERRQ(ierr);
2399  ierr = PetscOptionsEnd();
2400  CHKERRQ(ierr);
2401  ierr = DMSetUp_Moose_Pre(dm);
2402  CHKERRQ(ierr); /* Need some preliminary set up because, strangely enough, DMView() is called in
2403  DMSetFromOptions(). */
2405 }
2406 
2407 #undef __FUNCT__
2408 #define __FUNCT__ "DMDestroy_Moose"
2409 static PetscErrorCode
2411 {
2412  DM_Moose * dmm = (DM_Moose *)(dm->data);
2413  PetscErrorCode ierr;
2414 
2416  if (dmm->_vars)
2417  delete dmm->_vars;
2418  delete dmm->_var_ids;
2419  delete dmm->_var_names;
2420  if (dmm->_blocks)
2421  delete dmm->_blocks;
2422  delete dmm->_block_ids;
2423  delete dmm->_block_names;
2424  if (dmm->_sides)
2425  delete dmm->_sides;
2426  delete dmm->_side_ids;
2427  delete dmm->_side_names;
2428  if (dmm->_unsides)
2429  delete dmm->_unsides;
2430  delete dmm->_unside_ids;
2431  delete dmm->_unside_names;
2432  if (dmm->_contacts)
2433  delete dmm->_contacts;
2434  delete dmm->_contact_names;
2435  delete dmm->_contact_displaced;
2436  if (dmm->_uncontacts)
2437  delete dmm->_uncontacts;
2438  delete dmm->_uncontact_names;
2439  delete dmm->_uncontact_displaced;
2440  if (dmm->_splits)
2441  {
2442  for (auto & sit : *(dmm->_splits))
2443  {
2444  ierr = DMDestroy(&(sit.second._dm));
2445  CHKERRQ(ierr);
2446  ierr = ISDestroy(&(sit.second._rembedding));
2447  CHKERRQ(ierr);
2448  }
2449  delete dmm->_splits;
2450  }
2451  if (dmm->_splitlocs)
2452  delete dmm->_splitlocs;
2453  ierr = ISDestroy(&dmm->_embedding);
2454  CHKERRQ(ierr);
2455  ierr = PetscFree(dm->data);
2456  CHKERRQ(ierr);
2458 }
2459 
2460 #undef __FUNCT__
2461 #define __FUNCT__ "DMCreateMoose"
2462 PetscErrorCode
2464 {
2465  PetscErrorCode ierr;
2466 
2468  ierr = DMCreate(comm, dm);
2469  CHKERRQ(ierr);
2470  ierr = DMSetType(*dm, DMMOOSE);
2471  CHKERRQ(ierr);
2473  CHKERRQ(ierr);
2475 }
2476 
2477 EXTERN_C_BEGIN
2478 #undef __FUNCT__
2479 #define __FUNCT__ "DMCreate_Moose"
2480 PetscErrorCode
2482 {
2483  PetscErrorCode ierr;
2484  DM_Moose * dmm;
2485 
2487  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
2488 #if PETSC_RELEASE_LESS_THAN(3, 5, 0)
2489  ierr = PetscNewLog(dm, DM_Moose, &dmm);
2490  CHKERRQ(ierr);
2491 #else
2492  ierr = PetscNewLog(dm, &dmm);
2493  CHKERRQ(ierr);
2494 #endif
2495  dm->data = dmm;
2496 
2497  dmm->_var_ids = new (std::map<std::string, unsigned int>);
2498  dmm->_block_ids = new (std::map<std::string, subdomain_id_type>);
2499  dmm->_var_names = new (std::map<unsigned int, std::string>);
2500  dmm->_block_names = new (std::map<unsigned int, std::string>);
2501  dmm->_side_ids = new (std::map<std::string, BoundaryID>);
2502  dmm->_side_names = new (std::map<BoundaryID, std::string>);
2503  dmm->_unside_ids = new (std::map<std::string, BoundaryID>);
2504  dmm->_unside_names = new (std::map<BoundaryID, std::string>);
2505  dmm->_contact_names = new (std::map<DM_Moose::ContactID, DM_Moose::ContactName>);
2506  dmm->_uncontact_names = new (std::map<DM_Moose::ContactID, DM_Moose::ContactName>);
2507  dmm->_contact_displaced = new (std::map<DM_Moose::ContactName, PetscBool>);
2508  dmm->_uncontact_displaced = new (std::map<DM_Moose::ContactName, PetscBool>);
2509 
2510  dmm->_splits = new (std::map<std::string, DM_Moose::SplitInfo>);
2511 
2512  dmm->_print_embedding = PETSC_FALSE;
2513 
2514  dm->ops->createglobalvector = DMCreateGlobalVector_Moose;
2515  dm->ops->createlocalvector = 0; // DMCreateLocalVector_Moose;
2516  dm->ops->getcoloring = 0; // DMGetColoring_Moose;
2517  dm->ops->creatematrix = DMCreateMatrix_Moose;
2518  dm->ops->createinterpolation = 0; // DMCreateInterpolation_Moose;
2519 
2520  dm->ops->refine = 0; // DMRefine_Moose;
2521  dm->ops->coarsen = 0; // DMCoarsen_Moose;
2522  dm->ops->getinjection = 0; // DMGetInjection_Moose;
2523  dm->ops->getaggregates = 0; // DMGetAggregates_Moose;
2524 
2525 #if PETSC_VERSION_LT(3, 4, 0)
2526  dm->ops->createfielddecompositiondm = DMCreateFieldDecompositionDM_Moose;
2527  dm->ops->createdomaindecompositiondm = DMCreateDomainDecompositionDM_Moose;
2528 #endif
2529  dm->ops->createfielddecomposition = DMCreateFieldDecomposition_Moose;
2530  dm->ops->createdomaindecomposition = DMCreateDomainDecomposition_Moose;
2531 
2532  dm->ops->destroy = DMDestroy_Moose;
2533  dm->ops->view = DMView_Moose;
2534  dm->ops->setfromoptions = DMSetFromOptions_Moose;
2535  dm->ops->setup = DMSetUp_Moose;
2537 }
2538 EXTERN_C_END
2539 
2540 #undef __FUNCT__
2541 #define __FUNCT__ "SNESUpdateDMMoose"
2542 PetscErrorCode
2543 SNESUpdateDMMoose(SNES snes, PetscInt iteration)
2544 {
2545  /* This is called any time the structure of the problem changes in a way that affects the Jacobian
2546  sparsity pattern.
2547  For example, this may happen when NodeFaceConstraints change Jacobian's sparsity pattern based
2548  on newly-detected Penetration.
2549  In that case certain preconditioners (e.g., PCASM) will not work, unless we tell them that the
2550  sparsity pattern has changed.
2551  For now we are rebuilding the whole KSP, when necessary.
2552  */
2553  PetscErrorCode ierr;
2554  DM dm;
2555  KSP ksp;
2556  const char * prefix;
2557  MPI_Comm comm;
2558  PC pc;
2559 
2561  if (iteration)
2562  {
2563  /* TODO: limit this only to situations when displaced (un)contact splits are present, as is
2564  * DisplacedProblem(). */
2565  ierr = SNESGetDM(snes, &dm);
2566  CHKERRQ(ierr);
2567  ierr = DMMooseReset(dm);
2568  CHKERRQ(ierr);
2569  ierr = DMSetUp(dm);
2570  CHKERRQ(ierr);
2571  ierr = SNESGetKSP(snes, &ksp);
2572  CHKERRQ(ierr);
2573  /* Should we rebuild the whole KSP? */
2574  ierr = PetscObjectGetOptionsPrefix((PetscObject)ksp, &prefix);
2575  CHKERRQ(ierr);
2576  ierr = PetscObjectGetComm((PetscObject)ksp, &comm);
2577  CHKERRQ(ierr);
2578  ierr = PCCreate(comm, &pc);
2579  CHKERRQ(ierr);
2580  ierr = PCSetDM(pc, dm);
2581  CHKERRQ(ierr);
2582  ierr = PCSetOptionsPrefix(pc, prefix);
2583  CHKERRQ(ierr);
2584  ierr = PCSetFromOptions(pc);
2585  CHKERRQ(ierr);
2586  ierr = KSPSetPC(ksp, pc);
2587  CHKERRQ(ierr);
2588  ierr = PCDestroy(&pc);
2589  CHKERRQ(ierr);
2590  }
2592 }
2593 
2594 #undef __FUNCT__
2595 #define __FUNCT__ "DMMooseRegisterAll"
2596 PetscErrorCode
2598 {
2599  static PetscBool DMMooseRegisterAllCalled = PETSC_FALSE;
2600  PetscErrorCode ierr;
2601 
2603  if (!DMMooseRegisterAllCalled)
2604  {
2605 #if PETSC_VERSION_LESS_THAN(3, 4, 0)
2606  ierr = DMRegister(DMMOOSE, PETSC_NULL, "DMCreate_Moose", DMCreate_Moose);
2607  CHKERRQ(ierr);
2608 #else
2609  ierr = DMRegister(DMMOOSE, DMCreate_Moose);
2610  CHKERRQ(ierr);
2611 #endif
2612  DMMooseRegisterAllCalled = PETSC_TRUE;
2613  }
2615 }
2616 #endif // #if defined(LIBMESH_HAVE_PETSC) && !PETSC_VERSION_LESS_THAN(3,3,0)
PetscErrorCode DMMooseRegisterAll()
static PetscErrorCode DMSetUp_Moose(DM dm)
PetscErrorCode DMSetFromOptions_Moose(PetscOptionItems *, DM dm) PetscErrorCode DMSetFromOptions_Moose(PetscOptions *
std::set< ContactName > * _uncontacts
Definition: PetscDMMoose.C:66
PetscErrorCode DMMooseGetBlocks(DM dm, std::vector< std::string > &block_names)
Definition: PetscDMMoose.C:197
DM_Moose * dmm
PetscErrorCode DMMooseSetContacts(DM dm, const std::vector< std::pair< std::string, std::string >> &contacts, const std::vector< PetscBool > &displaced)
Definition: PetscDMMoose.C:380
bool _all_vars
Definition: PetscDMMoose.C:49
static PetscErrorCode Vec Mat jac
PetscErrorCode DMMooseSetVariables(DM dm, const std::set< std::string > &vars)
Definition: PetscDMMoose.C:272
std::map< BoundaryID, std::string > * _side_names
Definition: PetscDMMoose.C:55
virtual NonlinearSolver< Number > * nonlinearSolver()=0
PetscInt N
if(nl->nonlinearSolver() ->matvec &&nl->nonlinearSolver() ->residual_and_jacobian_object)
PetscErrorCode DMMooseGetUnSides(DM dm, std::vector< std::string > &side_names)
Definition: PetscDMMoose.C:173
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:207
Data structure used to hold penetration information.
std::map< ContactName, PetscBool > * _contact_displaced
Definition: PetscDMMoose.C:68
PetscMatrix< Number > Jac(jac, nl->comm())
PetscErrorCode DMMooseGetNonlinearSystem(DM dm, NonlinearSystemBase *&nl)
Definition: PetscDMMoose.C:462
PenetrationLocator & getPenetrationLocator(const BoundaryName &master, const BoundaryName &slave, Order order=FIRST)
PetscVector< Number > & X_sys
bool _nounsides
Definition: PetscDMMoose.C:61
StoredRange< MooseMesh::const_bnd_elem_iterator, const BndElement * > ConstBndElemRange
Definition: MooseMesh.h:1259
PetscErrorCode DMMooseSetUnContacts(DM dm, const std::vector< std::pair< std::string, std::string >> &uncontacts, const std::vector< PetscBool > &displaced)
Definition: PetscDMMoose.C:421
bool _nouncontacts
Definition: PetscDMMoose.C:71
std::map< std::string, BoundaryID > * _side_ids
Definition: PetscDMMoose.C:56
bool _include_all_contact_nodes
Definition: PetscDMMoose.C:72
PetscErrorCode SNESUpdateDMMoose(SNES snes, PetscInt iteration)
PetscBool ismoose
static PetscErrorCode DMMooseGetMeshBlocks_Private(DM dm, std::set< subdomain_id_type > &blocks)
PetscFunctionBegin
std::map< ContactName, PetscBool > * _uncontact_displaced
Definition: PetscDMMoose.C:69
PetscErrorCode DMMooseReset(DM dm)
static PetscErrorCode Vec Mat Mat pc
PetscErrorCode DMMooseGetUnContacts(DM dm, std::vector< std::pair< std::string, std::string >> &uncontact_names, std::vector< PetscBool > &displaced)
Definition: PetscDMMoose.C:120
static PetscErrorCode Vec x
PetscErrorCode DMMooseSetNonlinearSystem(DM dm, NonlinearSystemBase &nl)
Definition: PetscDMMoose.C:245
PetscVector< Number > X_global(x, nl->comm())
PETSC_ERR_ARG_WRONGSTATE
static PetscErrorCode DMSetUp_Moose_Pre(DM dm)
std::set< std::string > * _sides
Definition: PetscDMMoose.C:54
virtual GeometricSearchData & geomSearchData() override
PetscErrorCode DMCreateDomainDecompositionDM_Moose(DM dm, const char *, DM *ddm)
bool _nocontacts
Definition: PetscDMMoose.C:70
static PetscErrorCode DMVariableBounds_Moose(DM dm, Vec xl, Vec xu)
std::map< std::string, unsigned int > * _var_ids
Definition: PetscDMMoose.C:47
std::map< unsigned int, std::string > * _var_names
Definition: PetscDMMoose.C:48
NonlinearSystemBase * nl
Nonlinear system to be solved.
std::vector< std::string > split(const std::string &str, const std::string &delimiter)
Python like split function for strings.
Definition: MooseUtils.C:736
PetscErrorCode DMMooseGetContacts(DM dm, std::vector< std::pair< std::string, std::string >> &contact_names, std::vector< PetscBool > &displaced)
Definition: PetscDMMoose.C:91
FEProblemBase & _fe_problem
PetscMatrix< Number > the_pc(pc, nl->comm())
virtual std::shared_ptr< DisplacedProblem > getDisplacedProblem()
static PetscErrorCode DMView_Moose(DM dm, PetscViewer viewer)
std::map< ContactID, ContactName > * _uncontact_names
Definition: PetscDMMoose.C:67
static PetscErrorCode Mat * A
PetscInt m
const Elem * _elem
boundary_id_type BoundaryID
bool isCaptured() const
DofMap & dof_map
static PetscErrorCode DMMooseJacobian(DM dm, Vec x, Mat jac, Mat pc, MatStructure *msflag) static PetscErrorCode DMMooseJacobian(DM dm
const std::vector< numeric_index_type > & n_nz
std::set< ContactName > * _contacts
Definition: PetscDMMoose.C:64
PetscErrorCode DMMooseGetSplitNames(DM dm, std::vector< std::string > &split_names)
Definition: PetscDMMoose.C:536
PetscErrorCode DMCreateMoose(MPI_Comm comm, NonlinearSystemBase &nl, DM *dm)
IS _embedding
Definition: PetscDMMoose.C:84
bool _all_blocks
Definition: PetscDMMoose.C:53
const Elem * _side
std::map< std::string, SplitInfo > * _splits
Definition: PetscDMMoose.C:83
* msflag
static PetscErrorCode DMCreateFieldDecomposition_Moose(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
Definition: PetscDMMoose.C:973
PetscFunctionReturn(0)
std::map< std::pair< unsigned int, unsigned int >, PenetrationLocator * > _penetration_locators
std::pair< BoundaryID, BoundaryID > ContactID
Definition: PetscDMMoose.C:63
CHKERRQ(ierr)
static PetscErrorCode DMCreateMatrix_Moose(DM dm, const MatType type, Mat *A) static PetscErrorCode DMCreateMatrix_Moose(DM dm
std::set< std::string > * _unsides
Definition: PetscDMMoose.C:57
std::set< std::string > * _vars
Definition: PetscDMMoose.C:46
static PetscErrorCode SNESJacobian_DMMoose(SNES, Vec x, Mat *jac, Mat *pc, MatStructure *flag, void *ctx) static PetscErrorCode SNESJacobian_DMMoose(SNES
static PetscErrorCode DMMooseGetEmbedding_Private(DM dm, IS *embedding)
Definition: PetscDMMoose.C:569
MatType type
std::map< std::string, BoundaryID > * _unside_ids
Definition: PetscDMMoose.C:58
PetscErrorCode DMCreate_Moose(DM dm)
PetscErrorCode DMMooseGetSides(DM dm, std::vector< std::string > &side_names)
Definition: PetscDMMoose.C:149
PetscBool _print_embedding
Definition: PetscDMMoose.C:85
NonlinearSystemBase * _nl
Definition: PetscDMMoose.C:45
PetscInt M
static PetscErrorCode SNESFunction_DMMoose(SNES, Vec x, Vec r, void *ctx)
virtual System & system() override
Get the reference to the libMesh system.
PetscInt n
virtual MooseMesh & mesh()
Definition: SystemBase.h:103
PetscErrorCode DMMooseSetSides(DM dm, const std::set< std::string > &sides)
Definition: PetscDMMoose.C:326
MPI_Comm comm
virtual MooseMesh & mesh() override
std::map< std::string, subdomain_id_type > * _block_ids
Definition: PetscDMMoose.C:51
static PetscErrorCode DMMooseFunction(DM dm, Vec x, Vec r)
ierr
StoredRange< MooseMesh::const_bnd_elem_iterator, const BndElement * > * getBoundaryElementRange()
Definition: MooseMesh.C:801
StoredRange< MooseMesh::const_bnd_node_iterator, const BndNode * > ConstBndNodeRange
Some useful StoredRange typedefs.
Definition: MooseMesh.h:1258
std::set< std::string > * _blocks
Definition: PetscDMMoose.C:50
PetscErrorCode DMCreateFieldDecompositionDM_Moose(DM dm, const char *, DM *ddm)
static PetscErrorCode DMCreateGlobalVector_Moose(DM dm, Vec *x)
PetscErrorCode DMMooseSetBlocks(DM dm, const std::set< std::string > &blocks)
Definition: PetscDMMoose.C:299
std::map< ContactID, ContactName > * _contact_names
Definition: PetscDMMoose.C:65
bool _nosides
Definition: PetscDMMoose.C:60
std::multimap< std::string, unsigned int > * _splitlocs
Definition: PetscDMMoose.C:77
std::pair< std::string, std::string > ContactName
Definition: PetscDMMoose.C:62
static PetscErrorCode DMDestroy_Moose(DM dm)
static PetscErrorCode DMCreateDomainDecomposition_Moose(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
PetscErrorCode DMMooseSetSplitNames(DM dm, const std::vector< std::string > &split_names)
Definition: PetscDMMoose.C:485
static PetscErrorCode Vec Mat Mat void * ctx
StoredRange< MooseMesh::const_bnd_node_iterator, const BndNode * > * getBoundaryNodeRange()
Definition: MooseMesh.C:788
const std::vector< numeric_index_type > & n_oz
PetscErrorCode DMMooseSetUnSides(DM dm, const std::set< std::string > &unsides)
Definition: PetscDMMoose.C:353
BoundaryID getBoundaryID(const BoundaryName &boundary_name) const
Get the associated BoundaryID for the boundary name.
Definition: MooseMesh.C:1007
PetscErrorCode DMMooseGetVariables(DM dm, std::vector< std::string > &var_names)
Definition: PetscDMMoose.C:221
std::map< BoundaryID, std::string > * _unside_names
Definition: PetscDMMoose.C:59
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:690
std::map< unsigned int, std::string > * _block_names
Definition: PetscDMMoose.C:52