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