libMesh
petscdmlibmeshimpl.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2019 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 #include "libmesh/petsc_macro.h"
21 // This only works with petsc-3.3 and above.
22 
23 #if !PETSC_VERSION_LESS_THAN(3,3,0)
24 
25 // PETSc includes
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 // Local Includes
33 #include "libmesh/libmesh_common.h"
34 #include "libmesh/nonlinear_implicit_system.h"
35 #include "libmesh/petsc_nonlinear_solver.h"
36 #include "libmesh/petsc_linear_solver.h"
37 #include "libmesh/petsc_vector.h"
38 #include "libmesh/petsc_matrix.h"
39 #include "libmesh/petscdmlibmesh.h"
40 #include "libmesh/dof_map.h"
41 #include "libmesh/preconditioner.h"
42 #include "libmesh/elem.h"
43 #include "libmesh/parallel.h"
44 
45 
46 using namespace libMesh;
47 
48 
49 #define DMLIBMESH_NO_DECOMPOSITION 0
50 #define DMLIBMESH_FIELD_DECOMPOSITION 1
51 #define DMLIBMESH_DOMAIN_DECOMPOSITION 2
52 
53 #define DMLIBMESH_NO_EMBEDDING 0
54 #define DMLIBMESH_FIELD_EMBEDDING 1
55 #define DMLIBMESH_DOMAIN_EMBEDDING 2
56 
57 struct DM_libMesh
58 {
60  std::map<std::string, unsigned int> * varids;
61  std::map<unsigned int, std::string> * varnames;
62  std::map<std::string, unsigned int> * blockids;
63  std::map<unsigned int, std::string> * blocknames;
64  unsigned int decomposition_type;
65  std::vector<std::set<unsigned int>> * decomposition;
66  unsigned int embedding_type;
68  unsigned int vec_count;
69 };
70 
71 struct DMVec_libMesh {
72  std::string label;
73 };
74 
75 #undef __FUNCT__
76 #define __FUNCT__ "DMlibMeshGetVec_Private"
77 PetscErrorCode DMlibMeshGetVec_Private(DM /*dm*/, const char * /*name*/, Vec *)
78 {
79  PetscFunctionBegin;
80 
82 }
83 
84 
85 
86 PetscErrorCode DMlibMeshSetUpName_Private(DM dm);
87 
88 #undef __FUNCT__
89 #define __FUNCT__ "DMlibMeshSetSystem_libMesh"
91 {
92  const Parallel::Communicator & comm = sys.comm();
93 
94  PetscErrorCode ierr;
95  PetscFunctionBegin;
96  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
97  PetscBool islibmesh;
98  ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh);
99  if (!islibmesh) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
100 
101  if (dm->setupcalled) SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot reset the libMesh system after DM has been set up.");
102  DM_libMesh * dlm = (DM_libMesh *)(dm->data);
103  dlm->sys =&sys;
104  /* Initially populate the sets of active blockids and varids using all of the
105  existing blocks/variables (only variables are supported at the moment). */
106  DofMap & dofmap = dlm->sys->get_dof_map();
107  dlm->varids->clear();
108  dlm->varnames->clear();
109  for (auto v : IntRange<unsigned int>(0, dofmap.n_variables())) {
110  std::string vname = dofmap.variable(v).name();
111  dlm->varids->insert(std::pair<std::string,unsigned int>(vname,v));
112  dlm->varnames->insert(std::pair<unsigned int,std::string>(v,vname));
113  }
114  const MeshBase & mesh = dlm->sys->get_mesh();
115  dlm->blockids->clear();
116  dlm->blocknames->clear();
117  std::set<subdomain_id_type> blocks;
118  /* The following effectively is a verbatim copy of MeshBase::n_subdomains(). */
119  // This requires an inspection on every processor
120  libmesh_parallel_only(mesh.comm());
121  for (const auto & elem : mesh.active_element_ptr_range())
122  blocks.insert(elem->subdomain_id());
123  // Some subdomains may only live on other processors
124  comm.set_union(blocks);
125 
126  std::set<subdomain_id_type>::iterator bit = blocks.begin();
127  std::set<subdomain_id_type>::iterator bend = blocks.end();
128  if (bit == bend) SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "No mesh blocks found.");
129 
130  for (; bit != bend; ++bit) {
131  subdomain_id_type bid = *bit;
132  std::string bname = mesh.subdomain_name(bid);
133  if (!bname.length()) {
134  /* Block names are currently implemented for Exodus II meshes
135  only, so we might have to make up our own block names and
136  maintain our own mapping of block ids to names.
137  */
138  std::ostringstream ss;
139  ss << "dm" << bid;
140  bname = ss.str();
141  }
142  dlm->blockids->insert(std::pair<std::string,unsigned int>(bname,bid));
143  dlm->blocknames->insert(std::pair<unsigned int,std::string>(bid,bname));
144  }
147 }
148 
149 #undef __FUNCT__
150 #define __FUNCT__ "DMlibMeshGetSystem_libMesh"
152 {
153  PetscErrorCode ierr;
154 
155  PetscFunctionBegin;
156  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
157  PetscBool islibmesh;
158  ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh);CHKERRQ(ierr);
159  if (!islibmesh) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
160  DM_libMesh * dlm = (DM_libMesh *)(dm->data);
161  sys = dlm->sys;
163 }
164 
165 
166 #undef __FUNCT__
167 #define __FUNCT__ "DMlibMeshGetBlocks"
168 PetscErrorCode DMlibMeshGetBlocks(DM dm, PetscInt * n, char *** blocknames)
169 {
170  PetscErrorCode ierr;
171  PetscInt i;
172  PetscFunctionBegin;
173  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
174  PetscBool islibmesh;
175  ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh);
176  if (!islibmesh) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
177  DM_libMesh * dlm = (DM_libMesh *)(dm->data);
178  PetscValidPointer(n,2);
179  *n = cast_int<unsigned int>(dlm->blockids->size());
180  if (!blocknames) PetscFunctionReturn(0);
181  ierr = PetscMalloc(*n*sizeof(char *), blocknames); CHKERRQ(ierr);
182  i = 0;
183  for (const auto & pr : *(dlm->blockids))
184  {
185  ierr = PetscStrallocpy(pr.first.c_str(), *blocknames+i); CHKERRQ(ierr);
186  ++i;
187  }
189 }
190 
191 #undef __FUNCT__
192 #define __FUNCT__ "DMlibMeshGetVariables"
193 PetscErrorCode DMlibMeshGetVariables(DM dm, PetscInt * n, char *** varnames)
194 {
195  PetscErrorCode ierr;
196  PetscFunctionBegin;
197  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
198  PetscBool islibmesh;
199  PetscInt i;
200  ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh);
201  if (!islibmesh) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
202  DM_libMesh * dlm = (DM_libMesh *)(dm->data);
203  PetscValidPointer(n,2);
204  *n = cast_int<unsigned int>(dlm->varids->size());
205  if (!varnames) PetscFunctionReturn(0);
206  ierr = PetscMalloc(*n*sizeof(char *), varnames); CHKERRQ(ierr);
207  i = 0;
208  for (const auto & pr : *(dlm->varids))
209  {
210  ierr = PetscStrallocpy(pr.first.c_str(), *varnames+i); CHKERRQ(ierr);
211  ++i;
212  }
214 }
215 
216 #undef __FUNCT__
217 #define __FUNCT__ "DMlibMeshSetUpName_Private"
218 PetscErrorCode DMlibMeshSetUpName_Private(DM dm)
219 {
220  DM_libMesh * dlm = (DM_libMesh *)dm->data;
221  PetscErrorCode ierr;
222  PetscFunctionBegin;
223  std::string name = dlm->sys->name();
224  std::map<unsigned int, std::string> * dnames = PETSC_NULL, * enames = PETSC_NULL;
225  if (dlm->decomposition_type == DMLIBMESH_FIELD_DECOMPOSITION) {
226  name += ":dec:var:";
227  dnames = dlm->varnames;
228  }
229  if (dlm->decomposition_type == DMLIBMESH_DOMAIN_DECOMPOSITION) {
230  name += ":dec:block:";
231  dnames = dlm->blocknames;
232  }
233  if (dnames) {
234  for (auto decomp : *dlm->decomposition) {
235  for (std::set<unsigned int>::iterator dit_begin = decomp.begin(),
236  dit = dit_begin,
237  dit_end = decomp.end();
238  dit != dit_end; ++dit) {
239  unsigned int id = *dit;
240  if (dit != dit_begin)
241  name += ",";
242  name += (*dnames)[id];
243  }
244  name += ";";
245  }
246  }
247  if (dlm->embedding_type == DMLIBMESH_FIELD_EMBEDDING) {
248  name += ":emb:var:";
249  enames = dlm->varnames;
250  }
251  if (dlm->embedding_type == DMLIBMESH_DOMAIN_EMBEDDING) {
252  name += ":emb:block:";
253  enames = dlm->blocknames;
254  }
255  if (enames) {
256  for (auto eit = enames->begin(),
257  eit_end = enames->end(); eit != eit_end; ++eit) {
258  std::string & ename = eit->second;
259  if (eit != enames->begin())
260  name += ",";
261  name += ename;
262  }
263  name += ";";
264  }
265  ierr = PetscObjectSetName((PetscObject)dm, name.c_str()); CHKERRQ(ierr);
267 }
268 
269 
270 #undef __FUNCT__
271 #define __FUNCT__ "DMCreateFieldDecomposition_libMesh"
272 static PetscErrorCode DMCreateFieldDecomposition_libMesh(DM dm, PetscInt * len, char *** namelist, IS ** islist, DM ** dmlist)
273 {
274  PetscFunctionBegin;
275  PetscErrorCode ierr;
276  DM_libMesh * dlm = (DM_libMesh *)(dm->data);
278  IS emb;
279  if (dlm->decomposition_type != DMLIBMESH_FIELD_DECOMPOSITION) PetscFunctionReturn(0);
280 
281  *len = cast_int<unsigned int>(dlm->decomposition->size());
282  if (namelist) {ierr = PetscMalloc(*len*sizeof(char *), namelist); CHKERRQ(ierr);}
283  if (islist) {ierr = PetscMalloc(*len*sizeof(IS), islist); CHKERRQ(ierr);}
284  if (dmlist) {ierr = PetscMalloc(*len*sizeof(DM), dmlist); CHKERRQ(ierr);}
285  DofMap & dofmap = dlm->sys->get_dof_map();
286  for (auto d : index_range(*dlm->decomposition)) {
287  std::set<numeric_index_type> dindices;
288  std::string dname;
289  std::map<std::string, unsigned int> dvarids;
290  std::map<unsigned int, std::string> dvarnames;
291  unsigned int dvcount = 0;
292  for (const auto & v : (*dlm->decomposition)[d]) {
293  std::string vname = (*dlm->varnames)[v];
294  dvarids.insert(std::pair<std::string, unsigned int>(vname,v));
295  dvarnames.insert(std::pair<unsigned int,std::string>(v,vname));
296  if (!dvcount) dname = vname;
297  else dname += "_" + vname;
298  ++dvcount;
299  if (!islist) continue;
300  // Iterate only over this DM's blocks.
301  for (const auto & pr : *(dlm->blockids)) {
302  const subdomain_id_type sbd_id = cast_int<subdomain_id_type>(pr.second);
303  for (const auto & elem :
306  //unsigned int e_subdomain = elem->subdomain_id();
307  std::vector<numeric_index_type> evindices;
308  // Get the degree of freedom indices for the given variable off the current element.
309  dofmap.dof_indices(elem, evindices, v);
310  for (numeric_index_type dof : evindices) {
311  if (dof >= dofmap.first_dof() && dof < dofmap.end_dof()) // might want to use variable_first/last_local_dof instead
312  dindices.insert(dof);
313  }
314  }
315  }
316  }
317  if (namelist) {
318  ierr = PetscStrallocpy(dname.c_str(),(*namelist)+d); CHKERRQ(ierr);
319  }
320  if (islist) {
321  IS dis;
322  PetscInt * darray;
323  ierr = PetscMalloc(sizeof(PetscInt)*dindices.size(), &darray); CHKERRQ(ierr);
324  numeric_index_type i = 0;
325  for (const auto & id : dindices) {
326  darray[i] = id;
327  ++i;
328  }
329  ierr = ISCreateGeneral(((PetscObject)dm)->comm,
330  cast_int<PetscInt>(dindices.size()),
331  darray, PETSC_OWN_POINTER, &dis);
332  CHKERRQ(ierr);
333  if (dlm->embedding) {
334  /* Create a relative embedding into the parent's index space. */
335 #if PETSC_RELEASE_LESS_THAN(3,3,1)
336  ierr = ISMapFactorRight(dis,dlm->embedding, PETSC_TRUE, &emb); CHKERRQ(ierr);
337 #else
338  ierr = ISEmbed(dis,dlm->embedding, PETSC_TRUE, &emb); CHKERRQ(ierr);
339 #endif
340  PetscInt elen, dlen;
341  ierr = ISGetLocalSize(emb, &elen); CHKERRQ(ierr);
342  ierr = ISGetLocalSize(dis, &dlen); CHKERRQ(ierr);
343  if (elen != dlen) SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Failed to embed subdomain %D", d);
344  ierr = ISDestroy(&dis); CHKERRQ(ierr);
345  dis = emb;
346  }
347  else {
348  emb = dis;
349  }
350  (*islist)[d] = dis;
351  }
352  if (dmlist) {
353  DM ddm;
354  ierr = DMCreate(((PetscObject)dm)->comm, &ddm); CHKERRQ(ierr);
355  ierr = DMSetType(ddm, DMLIBMESH); CHKERRQ(ierr);
356  DM_libMesh * ddlm = (DM_libMesh *)(ddm->data);
357  ddlm->sys = dlm->sys;
358  /* copy over the block ids and names */
359  *ddlm->blockids = *dlm->blockids;
360  *ddlm->blocknames = *dlm->blocknames;
361  /* set the vars from the d-th part of the decomposition. */
362  *ddlm->varids = dvarids;
363  *ddlm->varnames = dvarnames;
364  ierr = PetscObjectReference((PetscObject)emb); CHKERRQ(ierr);
365  ddlm->embedding = emb;
366  ddlm->embedding_type = DMLIBMESH_FIELD_EMBEDDING;
367 
369  ierr = DMSetFromOptions(ddm); CHKERRQ(ierr);
370  (*dmlist)[d] = ddm;
371  }
372  }
374 }
375 
376 #undef __FUNCT__
377 #define __FUNCT__ "DMCreateDomainDecomposition_libMesh"
378 static PetscErrorCode DMCreateDomainDecomposition_libMesh(DM dm, PetscInt * len, char *** namelist, IS ** innerislist, IS ** outerislist, DM ** dmlist)
379 {
380  PetscFunctionBegin;
381  PetscErrorCode ierr;
382  DM_libMesh * dlm = (DM_libMesh *)(dm->data);
384  IS emb;
385  if (dlm->decomposition_type != DMLIBMESH_DOMAIN_DECOMPOSITION) PetscFunctionReturn(0);
386  *len = cast_int<unsigned int>(dlm->decomposition->size());
387  if (namelist) {ierr = PetscMalloc(*len*sizeof(char *), namelist); CHKERRQ(ierr);}
388  if (innerislist) {ierr = PetscMalloc(*len*sizeof(IS), innerislist); CHKERRQ(ierr);}
389  if (outerislist) *outerislist = PETSC_NULL; /* FIX: allow mesh-based overlap. */
390  if (dmlist) {ierr = PetscMalloc(*len*sizeof(DM), dmlist); CHKERRQ(ierr);}
391  for (auto d : index_range(*dlm->decomposition)) {
392  std::set<numeric_index_type> dindices;
393  std::string dname;
394  std::map<std::string, unsigned int> dblockids;
395  std::map<unsigned int,std::string> dblocknames;
396  unsigned int dbcount = 0;
397  for (const auto & b : (*dlm->decomposition)[d]) {
398  std::string bname = (*dlm->blocknames)[b];
399  dblockids.insert(std::pair<std::string, unsigned int>(bname,b));
400  dblocknames.insert(std::pair<unsigned int,std::string>(b,bname));
401  if (!dbcount) dname = bname;
402  else dname += "_" + bname;
403  ++dbcount;
404  if (!innerislist) continue;
405  const subdomain_id_type b_sbd_id = cast_int<subdomain_id_type>(b);
408  for ( ; el != end_el; ++el) {
409  const Elem * elem = *el;
410  std::vector<numeric_index_type> evindices;
411  // Iterate only over this DM's variables.
412  for (const auto & pr : *(dlm->varids)) {
413  // Get the degree of freedom indices for the given variable off the current element.
414  sys->get_dof_map().dof_indices(elem, evindices, pr.second);
415  for (const auto & dof : evindices) {
416  if (dof >= sys->get_dof_map().first_dof() && dof < sys->get_dof_map().end_dof()) // might want to use variable_first/last_local_dof instead
417  dindices.insert(dof);
418  }
419  }
420  }
421  }
422  if (namelist) {
423  ierr = PetscStrallocpy(dname.c_str(),(*namelist)+d); CHKERRQ(ierr);
424  }
425  if (innerislist) {
426  PetscInt * darray;
427  IS dis;
428  ierr = PetscMalloc(sizeof(PetscInt)*dindices.size(), &darray); CHKERRQ(ierr);
429  numeric_index_type i = 0;
430  for (const auto & id : dindices) {
431  darray[i] = id;
432  ++i;
433  }
434  ierr = ISCreateGeneral(((PetscObject)dm)->comm,
435  cast_int<PetscInt>(dindices.size()),
436  darray, PETSC_OWN_POINTER, &dis);
437  CHKERRQ(ierr);
438  if (dlm->embedding) {
439  /* Create a relative embedding into the parent's index space. */
440 #if PETSC_RELEASE_LESS_THAN(3,3,1)
441  ierr = ISMapFactorRight(dis,dlm->embedding, PETSC_TRUE, &emb); CHKERRQ(ierr);
442 #else
443  ierr = ISEmbed(dis,dlm->embedding, PETSC_TRUE, &emb); CHKERRQ(ierr);
444 #endif
445  PetscInt elen, dlen;
446  ierr = ISGetLocalSize(emb, &elen); CHKERRQ(ierr);
447  ierr = ISGetLocalSize(dis, &dlen); CHKERRQ(ierr);
448  if (elen != dlen) SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Failed to embed field %D", d);
449  ierr = ISDestroy(&dis); CHKERRQ(ierr);
450  dis = emb;
451  }
452  else {
453  emb = dis;
454  }
455  if (innerislist) {
456  ierr = PetscObjectReference((PetscObject)dis); CHKERRQ(ierr);
457  (*innerislist)[d] = dis;
458  }
459  ierr = ISDestroy(&dis); CHKERRQ(ierr);
460  }
461  if (dmlist) {
462  DM ddm;
463  ierr = DMCreate(((PetscObject)dm)->comm, &ddm); CHKERRQ(ierr);
464  ierr = DMSetType(ddm, DMLIBMESH); CHKERRQ(ierr);
465  DM_libMesh * ddlm = (DM_libMesh *)(ddm->data);
466  ddlm->sys = dlm->sys;
467  /* copy over the varids and varnames */
468  *ddlm->varids = *dlm->varids;
469  *ddlm->varnames = *dlm->varnames;
470  /* set the blocks from the d-th part of the decomposition. */
471  *ddlm->blockids = dblockids;
472  *ddlm->blocknames = dblocknames;
473  ierr = PetscObjectReference((PetscObject)emb); CHKERRQ(ierr);
474  ddlm->embedding = emb;
475  ddlm->embedding_type = DMLIBMESH_DOMAIN_EMBEDDING;
476 
478  ierr = DMSetFromOptions(ddm); CHKERRQ(ierr);
479  (*dmlist)[d] = ddm;
480  }
481  }
483 }
484 
485 
486 #undef __FUNCT__
487 #define __FUNCT__ "DMlibMeshCreateFieldDecompositionDM"
488 PetscErrorCode DMlibMeshCreateFieldDecompositionDM(DM dm, PetscInt dnumber, PetscInt * dsizes, char *** dvarlists, DM * ddm)
489 {
490  PetscErrorCode ierr;
491  PetscBool islibmesh;
492  PetscFunctionBegin;
493  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
494  ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh);
495  if (!islibmesh) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
496  if (dnumber < 0) SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative number %D of decomposition parts", dnumber);
497  PetscValidPointer(ddm,5);
498  DM_libMesh * dlm = (DM_libMesh *)(dm->data);
499  ierr = DMCreate(((PetscObject)dm)->comm, ddm); CHKERRQ(ierr);
500  ierr = DMSetType(*ddm, DMLIBMESH); CHKERRQ(ierr);
501  DM_libMesh * ddlm = (DM_libMesh *)((*ddm)->data);
502  ddlm->sys = dlm->sys;
503  ddlm->varids = dlm->varids;
504  ddlm->varnames = dlm->varnames;
505  ddlm->blockids = dlm->blockids;
506  ddlm->blocknames = dlm->blocknames;
507  ddlm->decomposition = new(std::vector<std::set<unsigned int>>);
508  ddlm->decomposition_type = DMLIBMESH_FIELD_DECOMPOSITION;
509  if (dnumber) {
510  for (PetscInt d = 0; d < dnumber; ++d) {
511  if (dsizes[d] < 0) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative size %D of decomposition part %D", dsizes[d],d);
512  ddlm->decomposition->push_back(std::set<unsigned int>());
513  for (PetscInt v = 0; v < dsizes[d]; ++v) {
514  std::string vname(dvarlists[d][v]);
515  std::map<std::string, unsigned int>::const_iterator vit = dlm->varids->find(vname);
516  if (vit == dlm->varids->end())
517  SETERRQ3(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Variable %D on the %D-th list with name %s is not owned by this DM", v, d, dvarlists[d][v]);
518  unsigned int vid = vit->second;
519  (*ddlm->decomposition)[d].insert(vid);
520  }
521  }
522  }
523  else { // Empty splits indicate default: split all variables with one per split.
524  PetscInt d = 0;
525  for (const auto & pr : (*ddlm->varids)) {
526  ddlm->decomposition->push_back(std::set<unsigned int>());
527  unsigned int vid = pr.second;
528  (*ddlm->decomposition)[d].insert(vid);
529  ++d;
530  }
531  }
533  ierr = DMSetFromOptions(*ddm); CHKERRQ(ierr);
534  ierr = DMSetUp(*ddm); CHKERRQ(ierr);
536 }
537 
538 #undef __FUNCT__
539 #define __FUNCT__ "DMlibMeshCreateDomainDecompositionDM"
540 PetscErrorCode DMlibMeshCreateDomainDecompositionDM(DM dm, PetscInt dnumber, PetscInt * dsizes, char *** dblocklists, DM * ddm)
541 {
542  PetscErrorCode ierr;
543  PetscBool islibmesh;
544  PetscFunctionBegin;
545  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
546  ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh);
547  if (!islibmesh) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
548  if (dnumber < 0) SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative number %D of decomposition parts", dnumber);
549  PetscValidPointer(ddm,5);
550  DM_libMesh * dlm = (DM_libMesh *)(dm->data);
551  ierr = DMCreate(((PetscObject)dm)->comm, ddm); CHKERRQ(ierr);
552  ierr = DMSetType(*ddm, DMLIBMESH); CHKERRQ(ierr);
553  DM_libMesh * ddlm = (DM_libMesh *)((*ddm)->data);
554  ddlm->sys = dlm->sys;
555  ddlm->varids = dlm->varids;
556  ddlm->varnames = dlm->varnames;
557  ddlm->blockids = dlm->blockids;
558  ddlm->blocknames = dlm->blocknames;
559  ddlm->decomposition = new(std::vector<std::set<unsigned int>>);
560  ddlm->decomposition_type = DMLIBMESH_DOMAIN_DECOMPOSITION;
561  if (dnumber) {
562  for (PetscInt d = 0; d < dnumber; ++d) {
563  if (dsizes[d] < 0) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative size %D of decomposition part %D", dsizes[d],d);
564  ddlm->decomposition->push_back(std::set<unsigned int>());
565  for (PetscInt b = 0; b < dsizes[d]; ++b) {
566  std::string bname(dblocklists[d][b]);
567  std::map<std::string, unsigned int>::const_iterator bit = dlm->blockids->find(bname);
568  if (bit == dlm->blockids->end())
569  SETERRQ3(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Block %D on the %D-th list with name %s is not owned by this DM", b, d, dblocklists[d][b]);
570  unsigned int bid = bit->second;
571  (*ddlm->decomposition)[d].insert(bid);
572  }
573  }
574  }
575  else { // Empty splits indicate default: split all blocks with one per split.
576  PetscInt d = 0;
577  for (const auto & pr : (*ddlm->blockids)) {
578  ddlm->decomposition->push_back(std::set<unsigned int>());
579  unsigned int bid = pr.second;
580  (*ddlm->decomposition)[d].insert(bid);
581  ++d;
582  }
583  }
585  ierr = DMSetFromOptions(*ddm); CHKERRQ(ierr);
586  ierr = DMSetUp(*ddm); CHKERRQ(ierr);
588 }
589 
590 struct token {
591  const char * s;
592  struct token * next;
593 };
594 
595 
596 // The following functions are only needed for older PETScs.
597 #if PETSC_RELEASE_LESS_THAN(3,3,1)
598 
599 #undef __FUNCT__
600 #define __FUNCT__ "DMlibMeshParseDecompositionDescriptor_Private"
601 static PetscErrorCode DMlibMeshParseDecompositionDescriptor_Private(DM dm, const char * ddesc, PetscInt * dtype, PetscInt * dcount, PetscInt ** dsizes, char **** dlists)
602 {
603  PetscFunctionBegin;
604  PetscErrorCode ierr;
605  PetscBool eq;
606  char * s0;
607  char * s;
608  char * ss;
609  struct token * llfirst = PETSC_NULL;
610  struct token * lllast = PETSC_NULL;
611  struct token * tok;
612  PetscInt stcount = 0, brcount = 0, d, i;
613  size_t len0, count;
614 
615  /*
616  Parse the decomposition descriptor.
617  Decomposition names could be of one of two forms:
618  var:v1,v2;v3,v4;v4,v5;
619  block:b1,b2;b3,b4;b4,b5;
620  resulting in an overlapping decomposition that groups
621  variables (v1,v2), (v3,v4), (v4,v5) or
622  blocks (b1,b2), (b3,b4), (b4,b5).
623  */
624  /* Copy the descriptor so that we can manipulate it in place. */
625  ierr = PetscStrallocpy(ddesc,&s0); CHKERRQ(ierr);
626  ierr = PetscStrlen(s0, &len0) ; CHKERRQ(ierr);
627  ierr = PetscStrstr(s0,":",&ss); CHKERRQ(ierr);
628  if (!ss) {
629  ss = s0+len0;
630  }
631  else {
632  *ss = 0;
633  }
634  ierr = PetscStrcmp(s0,"var",&eq); CHKERRQ(ierr);
635  if (eq) {
636  *dtype=DMLIBMESH_FIELD_DECOMPOSITION;
637  }
638  else {
639  ierr = PetscStrcmp(s0,"block",&eq);CHKERRQ(ierr);
640  if (!eq)
641  SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Could not determine decomposition type from descriptor: %s\n", ddesc); CHKERRQ(ierr);
642  *dtype=DMLIBMESH_DOMAIN_DECOMPOSITION;
643  }
644  ierr = PetscStrlen(s0,&count); CHKERRQ(ierr);
645  while (count < len0) {
646  struct token * st;
647  struct token * br;
648  ++ss; ++count;
649  s=ss;
650  while (*ss && *ss != ',' && *ss != ';') {
651  ++ss; ++count;
652  }
653  st = PETSC_NULL; br = PETSC_NULL;
654  if (*ss) {
655  /*
656  Found a separator, or a break.
657  Add an appropriate token to the list.
658  A token separator ',' produces no token.
659  */
660  if (*ss == ';') {
661  /* Create a break token: a token with a null string. */
662 #if PETSC_RELEASE_LESS_THAN(3,5,0)
663  ierr = PetscNew(struct token,&br);CHKERRQ(ierr);
664 #else
665  ierr = PetscNew(&br);CHKERRQ(ierr);
666 #endif
667  }
668  *ss = 0;
669  if (s != ss) {
670  /* A nonempty string. */
671 #if PETSC_RELEASE_LESS_THAN(3,5,0)
672  ierr = PetscNew(struct token, &st);CHKERRQ(ierr);
673 #else
674  ierr = PetscNew(&st);CHKERRQ(ierr);
675 #endif
676  st->s = s; /* The string will be properly copied below. */
677  }
678  /* Add the new tokens to the list. */
679  if (st) {
680  if (!lllast) {
681  llfirst = lllast = st;
682  }
683  else {
684  lllast->next = st; lllast = st;
685  }
686  }
687  if (br) {
688  if (!lllast) {
689  llfirst = lllast = br;
690  }
691  else {
692  lllast->next = br; lllast = br;
693  }
694  }
695  }
696  }
697  /* The result of parsing is in the linked list ll. */
698  /* Count up the strings and the breaks. */
699  tok = llfirst;
700  while (tok) {
701  if (tok->s)
702  ++stcount;
703  else
704  ++brcount;
705  tok = tok->next;
706  }
707  /* Allocate the space for the output. */
708  *dcount = brcount;
709  ierr = PetscMalloc(*dcount*sizeof(PetscInt), dsizes); CHKERRQ(ierr);
710  ierr = PetscMalloc(*dcount*sizeof(char **), dlists); CHKERRQ(ierr);
711  for (d = 0; d < *dcount; ++d) (*dsizes)[d] = 0;
712  tok = llfirst; d = 0;
713  while (tok) {
714  if (tok->s)
715  ++(*dsizes)[d];
716  else
717  ++d;
718  tok = tok->next;
719  }
720  for (d = 0; d < *dcount; ++d) {
721  ierr = PetscMalloc(sizeof(char **)*(*dsizes)[d], (* dlists)+d); CHKERRQ(ierr);
722  }
723  /* Now copy strings and destroy tokens. */
724  tok = llfirst; d = 0; i = 0;
725  while (tok) {
726  if (tok->s) {
727  ierr = PetscStrallocpy(tok->s, (*dlists)[d]+i); CHKERRQ(ierr);
728  ++i;
729  }
730  else {
731  ++d;
732  i = 0;
733  }
734  llfirst = tok;
735  tok = tok->next;
736  ierr = PetscFree(llfirst); CHKERRQ(ierr);
737  }
738  /* Deallocate workspace. */
739  ierr = PetscFree(s0); CHKERRQ(ierr);
741 }
742 
743 #undef __FUNCT__
744 #define __FUNCT__ "DMCreateFieldDecompositionDM_libMesh"
745 static PetscErrorCode DMCreateFieldDecompositionDM_libMesh(DM dm, const char * ddesc, DM * ddm)
746 {
747  PetscFunctionBegin;
748  PetscErrorCode ierr;
749  PetscInt dtype, dcount;
750  PetscInt * dsizes;
751  char *** dlists;
752  PetscFunctionBegin;
753  *ddm = PETSC_NULL;
754  ierr = DMlibMeshParseDecompositionDescriptor_Private(dm,ddesc,&dtype,&dcount,&dsizes,&dlists); CHKERRQ(ierr);
755  if (dtype == DMLIBMESH_FIELD_DECOMPOSITION){
756  ierr = DMlibMeshCreateFieldDecompositionDM(dm,dcount,dsizes,dlists,ddm); CHKERRQ(ierr);
757  }
758  else SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Unexpected unknown decomposition type for field decomposition descriptor %s", ddesc);
760 }
761 
762 #undef __FUNCT__
763 #define __FUNCT__ "DMCreateDomainDecompositionDM_libMesh"
764 static PetscErrorCode DMCreateDomainDecompositionDM_libMesh(DM dm, const char * ddesc, DM * ddm)
765 {
766  PetscFunctionBegin;
767  PetscErrorCode ierr;
768  PetscInt dtype, dcount;
769  PetscInt * dsizes;
770  char *** dlists;
771  PetscFunctionBegin;
772  *ddm = PETSC_NULL;
773  ierr = DMlibMeshParseDecompositionDescriptor_Private(dm,ddesc,&dtype,&dcount,&dsizes,&dlists); CHKERRQ(ierr);
774  if (dtype == DMLIBMESH_DOMAIN_DECOMPOSITION) {
775  ierr = DMlibMeshCreateDomainDecompositionDM(dm,dcount,dsizes,dlists,ddm); CHKERRQ(ierr);
776  }
777  else SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Unexpected unknown decomposition type for domain decomposition descriptor %s", ddesc);
779 }
780 
781 #endif
782 
783 
784 #undef __FUNCT__
785 #define __FUNCT__ "DMlibMeshFunction"
786 static PetscErrorCode DMlibMeshFunction(DM dm, Vec x, Vec r)
787 {
788  PetscErrorCode ierr;
789  PetscFunctionBegin;
790  libmesh_assert(x);
791  libmesh_assert(r);
792 
794  ierr = DMlibMeshGetSystem(dm, _sys);CHKERRQ(ierr);
795  NonlinearImplicitSystem & sys = *_sys;
796  PetscVector<Number> & X_sys = *libmesh_cast_ptr<PetscVector<Number> *>(sys.solution.get());
797  PetscVector<Number> & R_sys = *libmesh_cast_ptr<PetscVector<Number> *>(sys.rhs);
798  PetscVector<Number> X_global(x, _sys->comm()), R(r, _sys->comm());
799 
800  // Use the systems update() to get a good local version of the parallel solution
801  X_global.swap(X_sys);
802  R.swap(R_sys);
803 
805  _sys->update();
806 
807  // Swap back
808  X_global.swap(X_sys);
809  R.swap(R_sys);
810  R.zero();
811 
812  // if the user has provided both function pointers and objects only the pointer
813  // will be used, so catch that as an error
814  if (_sys->nonlinear_solver->residual && _sys->nonlinear_solver->residual_object)
815  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the Residual!");
816 
817  if (_sys->nonlinear_solver->matvec && _sys->nonlinear_solver->residual_and_jacobian_object)
818  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
819 
820  if (_sys->nonlinear_solver->residual != nullptr)
821  _sys->nonlinear_solver->residual(*(_sys->current_local_solution.get()), R, *_sys);
822 
823  else if (_sys->nonlinear_solver->residual_object != nullptr)
824  _sys->nonlinear_solver->residual_object->residual(*(_sys->current_local_solution.get()), R, *_sys);
825 
826  else if (_sys->nonlinear_solver->matvec != nullptr)
827  _sys->nonlinear_solver->matvec(*(_sys->current_local_solution.get()), &R, nullptr, *_sys);
828 
829  else if (_sys->nonlinear_solver->residual_and_jacobian_object != nullptr)
830  _sys->nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian(*(_sys->current_local_solution.get()), &R, nullptr, *_sys);
831 
832  else
833  libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
834 
835  R.close();
836  X_global.close();
838 }
839 
840 #if !PETSC_RELEASE_LESS_THAN(3,3,1)
841 #undef __FUNCT__
842 #define __FUNCT__ "SNESFunction_DMlibMesh"
843 static PetscErrorCode SNESFunction_DMlibMesh(SNES, Vec x, Vec r, void * ctx)
844 {
845  DM dm = (DM)ctx;
846  PetscErrorCode ierr;
847  PetscFunctionBegin;
848  ierr = DMlibMeshFunction(dm,x,r);CHKERRQ(ierr);
850 }
851 #endif
852 
853 
854 #undef __FUNCT__
855 #define __FUNCT__ "DMlibMeshJacobian"
856 static PetscErrorCode DMlibMeshJacobian(
857 #if PETSC_RELEASE_LESS_THAN(3,5,0)
858  DM dm, Vec x, Mat jac, Mat pc, MatStructure * msflag
859 #else
860  DM dm, Vec x, Mat jac, Mat pc
861 #endif
862  )
863 {
864  PetscErrorCode ierr;
865  PetscFunctionBegin;
867  ierr = DMlibMeshGetSystem(dm, _sys); CHKERRQ(ierr);
868  NonlinearImplicitSystem & sys = *_sys;
869 
870  PetscMatrix<Number> the_pc(pc,sys.comm());
871  PetscMatrix<Number> Jac(jac,sys.comm());
872  PetscVector<Number> & X_sys = *libmesh_cast_ptr<PetscVector<Number> *>(sys.solution.get());
873  PetscMatrix<Number> & Jac_sys = *libmesh_cast_ptr<PetscMatrix<Number> *>(sys.matrix);
874  PetscVector<Number> X_global(x, sys.comm());
875 
876  // Set the dof maps
877  the_pc.attach_dof_map(sys.get_dof_map());
878  Jac.attach_dof_map(sys.get_dof_map());
879 
880  // Use the systems update() to get a good local version of the parallel solution
881  X_global.swap(X_sys);
882  Jac.swap(Jac_sys);
883 
885  sys.update();
886 
887  X_global.swap(X_sys);
888  Jac.swap(Jac_sys);
889 
890  the_pc.zero();
891 
892  // if the user has provided both function pointers and objects only the pointer
893  // will be used, so catch that as an error
894  if (sys.nonlinear_solver->jacobian && sys.nonlinear_solver->jacobian_object)
895  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the Jacobian!");
896 
897  if (sys.nonlinear_solver->matvec && sys.nonlinear_solver->residual_and_jacobian_object)
898  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
899 
900  if (sys.nonlinear_solver->jacobian != nullptr)
901  sys.nonlinear_solver->jacobian(*(sys.current_local_solution.get()), the_pc, sys);
902 
903  else if (sys.nonlinear_solver->jacobian_object != nullptr)
904  sys.nonlinear_solver->jacobian_object->jacobian(*(sys.current_local_solution.get()), the_pc, sys);
905 
906  else if (sys.nonlinear_solver->matvec != nullptr)
907  sys.nonlinear_solver->matvec(*(sys.current_local_solution.get()), nullptr, &the_pc, sys);
908 
909  else if (sys.nonlinear_solver->residual_and_jacobian_object != nullptr)
910  sys.nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian(*(sys.current_local_solution.get()), nullptr, &the_pc, sys);
911 
912  else
913  libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
914 
915  the_pc.close();
916  Jac.close();
917  X_global.close();
918 #if PETSC_RELEASE_LESS_THAN(3,5,0)
919  *msflag = SAME_NONZERO_PATTERN;
920 #endif
922 }
923 
924 #if !PETSC_RELEASE_LESS_THAN(3,3,1)
925 #undef __FUNCT__
926 #define __FUNCT__ "SNESJacobian_DMlibMesh"
927 static PetscErrorCode SNESJacobian_DMlibMesh(
928 #if PETSC_RELEASE_LESS_THAN(3,5,0)
929  SNES, Vec x, Mat * jac, Mat * pc, MatStructure * flag, void * ctx
930 #else
931  SNES, Vec x, Mat jac, Mat pc, void * ctx
932 #endif
933  )
934 {
935  DM dm = (DM)ctx;
936  PetscErrorCode ierr;
937  PetscFunctionBegin;
938 #if PETSC_RELEASE_LESS_THAN(3,5,0)
939  ierr = DMlibMeshJacobian(dm,x,*jac,*pc,flag); CHKERRQ(ierr);
940 #else
941  ierr = DMlibMeshJacobian(dm,x,jac,pc); CHKERRQ(ierr);
942 #endif
944 }
945 #endif
946 
947 #undef __FUNCT__
948 #define __FUNCT__ "DMVariableBounds_libMesh"
949 static PetscErrorCode DMVariableBounds_libMesh(DM dm, Vec xl, Vec xu)
950 {
951  PetscErrorCode ierr;
953  ierr = DMlibMeshGetSystem(dm, _sys); CHKERRQ(ierr);
954  NonlinearImplicitSystem & sys = *_sys;
955  PetscVector<Number> XL(xl, sys.comm());
956  PetscVector<Number> XU(xu, sys.comm());
957  PetscFunctionBegin;
958 #if PETSC_VERSION_LESS_THAN(3,5,0) && PETSC_VERSION_RELEASE
959  ierr = VecSet(xl, SNES_VI_NINF);CHKERRQ(ierr);
960  ierr = VecSet(xu, SNES_VI_INF);CHKERRQ(ierr);
961 #else
962  // Workaround for nonstandard Q suffix warning with quad precision
963  const PetscReal petsc_inf = std::numeric_limits<PetscReal>::max() / 4;
964  ierr = VecSet(xl, -petsc_inf);CHKERRQ(ierr);
965  ierr = VecSet(xu, petsc_inf);CHKERRQ(ierr);
966 #endif
967  if (sys.nonlinear_solver->bounds != nullptr)
968  sys.nonlinear_solver->bounds(XL,XU,sys);
969  else if (sys.nonlinear_solver->bounds_object != nullptr)
970  sys.nonlinear_solver->bounds_object->bounds(XL,XU, sys);
971  else
972  SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "No bounds calculation in this libMesh object");
973 
975 }
976 
977 
978 #undef __FUNCT__
979 #define __FUNCT__ "DMCreateGlobalVector_libMesh"
980 static PetscErrorCode DMCreateGlobalVector_libMesh(DM dm, Vec *x)
981 {
982  PetscFunctionBegin;
983  PetscErrorCode ierr;
984  DM_libMesh * dlm = (DM_libMesh *)(dm->data);
985  PetscBool eq;
986 
987  ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH, &eq); CHKERRQ(ierr);
988 
989  if (!eq)
990  SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "DM of type %s, not of type %s", ((PetscObject)dm)->type, DMLIBMESH);
991 
992  if (!dlm->sys)
993  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No libMesh system set for DM_libMesh");
994 
996  PetscVector<Number> * pv = dynamic_cast<PetscVector<Number> *>(nv);
997  Vec v = pv->vec();
998  /* Unfortunately, currently this does not produce a ghosted vector, so nonlinear subproblem solves aren't going to be easily available.
999  Should work fine for getting vectors out for linear subproblem solvers. */
1000  if (dlm->embedding) {
1001  PetscInt n;
1002  ierr = VecCreate(((PetscObject)v)->comm, x); CHKERRQ(ierr);
1003  ierr = ISGetLocalSize(dlm->embedding, &n); CHKERRQ(ierr);
1004  ierr = VecSetSizes(*x,n,PETSC_DETERMINE); CHKERRQ(ierr);
1005  ierr = VecSetType(*x,((PetscObject)v)->type_name); CHKERRQ(ierr);
1006  ierr = VecSetFromOptions(*x); CHKERRQ(ierr);
1007  ierr = VecSetUp(*x); CHKERRQ(ierr);
1008  }
1009  else {
1010  ierr = VecDuplicate(v,x); CHKERRQ(ierr);
1011  }
1012  ierr = PetscObjectCompose((PetscObject)*x,"DM",(PetscObject)dm); CHKERRQ(ierr);
1014 }
1015 
1016 
1017 
1018 
1019 #undef __FUNCT__
1020 #define __FUNCT__ "DMCreateMatrix_libMesh"
1021 #if PETSC_VERSION_LT(3,5,0)
1022 static PetscErrorCode DMCreateMatrix_libMesh(DM dm, const MatType, Mat * A)
1023 #else
1024  static PetscErrorCode DMCreateMatrix_libMesh(DM dm, Mat * A)
1025 #endif
1027  PetscFunctionBegin;
1028  PetscErrorCode ierr;
1029  DM_libMesh * dlm = (DM_libMesh *)(dm->data);
1031 
1032  ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH, &eq); CHKERRQ(ierr);
1033 
1034  if (!eq)
1035  SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "DM of type %s, not of type %s", ((PetscObject)dm)->type, DMLIBMESH);
1036 
1037  if (!dlm->sys)
1038  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No libMesh system set for DM_libMesh");
1039 
1040  *A = (dynamic_cast<PetscMatrix<Number> *>(dlm->sys->matrix))->mat();
1041  ierr = PetscObjectReference((PetscObject)(*A)); CHKERRQ(ierr);
1043 }
1044 
1045 
1046 #undef __FUNCT__
1047 #define __FUNCT__ "DMView_libMesh"
1048 static PetscErrorCode DMView_libMesh(DM dm, PetscViewer viewer)
1049 {
1050  PetscErrorCode ierr;
1051  PetscBool isascii;
1052  const char * name, * prefix;
1053  DM_libMesh * dlm = (DM_libMesh *)dm->data;
1054  PetscFunctionBegin;
1055  ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii); CHKERRQ(ierr);
1056  if (isascii) {
1057  ierr = PetscObjectGetName((PetscObject)dm, &name); CHKERRQ(ierr);
1058  ierr = PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix); CHKERRQ(ierr);
1059  ierr = PetscViewerASCIIPrintf(viewer, "DM libMesh with name %s and prefix %s\n", name, prefix); CHKERRQ(ierr);
1060  ierr = PetscViewerASCIIPrintf(viewer, "blocks:", name, prefix); CHKERRQ(ierr);
1061  std::map<std::string,unsigned int>::iterator bit = dlm->blockids->begin();
1062  std::map<std::string,unsigned int>::const_iterator bend = dlm->blockids->end();
1063  for (; bit != bend; ++bit) {
1064  ierr = PetscViewerASCIIPrintf(viewer, "(%s,%D) ", bit->first.c_str(), bit->second); CHKERRQ(ierr);
1065  }
1066  ierr = PetscViewerASCIIPrintf(viewer, "\n"); CHKERRQ(ierr);
1067  ierr = PetscViewerASCIIPrintf(viewer, "variables:", name, prefix); CHKERRQ(ierr);
1068  std::map<std::string,unsigned int>::iterator vit = dlm->varids->begin();
1069  std::map<std::string,unsigned int>::const_iterator vend = dlm->varids->end();
1070  for (; vit != vend; ++vit) {
1071  ierr = PetscViewerASCIIPrintf(viewer, "(%s,%D) ", vit->first.c_str(), vit->second); CHKERRQ(ierr);
1072  }
1073  ierr = PetscViewerASCIIPrintf(viewer, "\n"); CHKERRQ(ierr);
1074  if (dlm->decomposition_type == DMLIBMESH_NO_DECOMPOSITION) {
1075  ierr = PetscViewerASCIIPrintf(viewer, "No decomposition\n"); CHKERRQ(ierr);
1076  }
1077  else {
1078  if (dlm->decomposition_type == DMLIBMESH_FIELD_DECOMPOSITION) {
1079  ierr = PetscViewerASCIIPrintf(viewer, "Field decomposition by variable: "); CHKERRQ(ierr);
1080  }
1081  else if (dlm->decomposition_type == DMLIBMESH_DOMAIN_DECOMPOSITION) {
1082  ierr = PetscViewerASCIIPrintf(viewer, "Domain decomposition by block: "); CHKERRQ(ierr);
1083  }
1084  else SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Unexpected decomposition type: %D", dlm->decomposition_type);
1085  /* FIX: decompositions might have different sizes and components on different ranks. */
1086  for (auto d : index_range(*dlm->decomposition)) {
1087  std::set<unsigned int>::iterator dbegin = (*dlm->decomposition)[d].begin();
1088  std::set<unsigned int>::iterator dit = (*dlm->decomposition)[d].begin();
1089  std::set<unsigned int>::iterator dend = (*dlm->decomposition)[d].end();
1090  for (; dit != dend; ++dit) {
1091  if (dit != dbegin) {
1092  ierr = PetscViewerASCIIPrintf(viewer, ","); CHKERRQ(ierr);
1093  }
1094  ierr = PetscViewerASCIIPrintf(viewer, "%D", *dit); CHKERRQ(ierr);
1095  }
1096  ierr = PetscViewerASCIIPrintf(viewer, ";"); CHKERRQ(ierr);
1097  }
1098  ierr = PetscViewerASCIIPrintf(viewer, "\n"); CHKERRQ(ierr);
1099  }
1100  }
1101 
1103 }
1104 
1105 #undef __FUNCT__
1106 #define __FUNCT__ "DMSetUp_libMesh"
1107 static PetscErrorCode DMSetUp_libMesh(DM dm)
1108 {
1109  PetscFunctionBegin;
1110  PetscErrorCode ierr;
1111  DM_libMesh * dlm = (DM_libMesh *)(dm->data);
1112  PetscBool eq;
1113 
1114  ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH, &eq); CHKERRQ(ierr);
1115 
1116  if (!eq)
1117  SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "DM of type %s, not of type %s", ((PetscObject)dm)->type, DMLIBMESH);
1118 
1119  if (!dlm->sys)
1120  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No libMesh system set for DM_libMesh");
1121  /*
1122  Do not evaluate function, Jacobian or bounds for an embedded DM -- the subproblem might not have enough information for that.
1123  */
1124  if (!dlm->embedding) {
1125 #if PETSC_RELEASE_LESS_THAN(3,3,1)
1126  ierr = DMSetFunction(dm, DMlibMeshFunction); CHKERRQ(ierr);
1127  ierr = DMSetJacobian(dm, DMlibMeshJacobian); CHKERRQ(ierr);
1128 #else
1129  ierr = DMSNESSetFunction(dm, SNESFunction_DMlibMesh, (void *)dm); CHKERRQ(ierr);
1130  ierr = DMSNESSetJacobian(dm, SNESJacobian_DMlibMesh, (void *)dm); CHKERRQ(ierr);
1131 #endif
1132  if (dlm->sys->nonlinear_solver->bounds || dlm->sys->nonlinear_solver->bounds_object)
1133  ierr = DMSetVariableBounds(dm, DMVariableBounds_libMesh); CHKERRQ(ierr);
1134  }
1135  else {
1136  /*
1137  Fow now we don't implement even these, although a linear "Dirichlet" subproblem is well-defined.
1138  Creating the submatrix, however, might require extracting the submatrix preallocation from an unassembled matrix.
1139  */
1140  dm->ops->createglobalvector = 0;
1141  dm->ops->creatematrix = 0;
1142  }
1144 }
1145 
1146 
1147 
1148 
1149 #undef __FUNCT__
1150 #define __FUNCT__ "DMDestroy_libMesh"
1151 static PetscErrorCode DMDestroy_libMesh(DM dm)
1152 {
1153  DM_libMesh * dlm = (DM_libMesh *)(dm->data);
1154  PetscErrorCode ierr;
1155  PetscFunctionBegin;
1156  delete dlm->varids;
1157  delete dlm->varnames;
1158  delete dlm->blockids;
1159  delete dlm->blocknames;
1160  delete dlm->decomposition;
1161  ierr = ISDestroy(&dlm->embedding); CHKERRQ(ierr);
1162  ierr = PetscFree(dm->data); CHKERRQ(ierr);
1163 
1165 }
1166 
1167 EXTERN_C_BEGIN
1168 #undef __FUNCT__
1169 #define __FUNCT__ "DMCreate_libMesh"
1170 PetscErrorCode DMCreate_libMesh(DM dm)
1171 {
1172  PetscErrorCode ierr;
1173  DM_libMesh * dlm;
1174 
1175  PetscFunctionBegin;
1176  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1177 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1178  ierr = PetscNewLog(dm,DM_libMesh,&dlm);CHKERRQ(ierr);
1179 #else
1180  ierr = PetscNewLog(dm,&dlm);CHKERRQ(ierr);
1181 #endif
1182  dm->data = dlm;
1183 
1184  /* DMlibMesh impl */
1185  dlm->varids = new(std::map<std::string, unsigned int>);
1186  dlm->blockids = new(std::map<std::string, unsigned int>);
1187  dlm->varnames = new(std::map<unsigned int, std::string>);
1188  dlm->blocknames = new(std::map<unsigned int, std::string>);
1189  dlm->decomposition = PETSC_NULL;
1190  dlm->decomposition_type = DMLIBMESH_NO_DECOMPOSITION;
1191 
1192  /* DM API */
1193  dm->ops->createglobalvector = DMCreateGlobalVector_libMesh;
1194  dm->ops->createlocalvector = 0; // DMCreateLocalVector_libMesh;
1195  dm->ops->getcoloring = 0; // DMGetColoring_libMesh;
1196  dm->ops->creatematrix = DMCreateMatrix_libMesh;
1197  dm->ops->createinterpolation= 0; // DMCreateInterpolation_libMesh;
1198 
1199  dm->ops->refine = 0; // DMRefine_libMesh;
1200  dm->ops->coarsen = 0; // DMCoarsen_libMesh;
1201 
1202  // * dm->ops->getinjection was renamed to dm->ops->createinjection in PETSc 5a84ad338 (5 Jul 2019)
1203  // * dm->ops-getaggregates was removed in PETSc 97779f9a (5 Jul 2019)
1204  // * Both changes were merged into PETSc master in 94aad3ce (7 Jul 2019).
1205 #if PETSC_RELEASE_LESS_THAN(3,12,0)
1206  dm->ops->getinjection = 0; // DMGetInjection_libMesh;
1207  dm->ops->getaggregates = 0; // DMGetAggregates_libMesh;
1208 #else
1209  dm->ops->createinjection = 0;
1210 #endif
1211 
1212 
1213 #if PETSC_RELEASE_LESS_THAN(3,3,1)
1214  dm->ops->createfielddecompositiondm = DMCreateFieldDecompositionDM_libMesh;
1215  dm->ops->createdomaindecompositiondm = DMCreateDomainDecompositionDM_libMesh;
1216 #endif
1217  dm->ops->createfielddecomposition = DMCreateFieldDecomposition_libMesh;
1218  dm->ops->createdomaindecomposition = DMCreateDomainDecomposition_libMesh;
1219 
1220  dm->ops->destroy = DMDestroy_libMesh;
1221  dm->ops->view = DMView_libMesh;
1222  dm->ops->setfromoptions = 0; // DMSetFromOptions_libMesh;
1223  dm->ops->setup = DMSetUp_libMesh;
1224 
1225  /* DMlibMesh API */
1226 #if PETSC_RELEASE_LESS_THAN(3,4,0)
1227  ierr = PetscObjectComposeFunction((PetscObject)dm,"DMlibMeshSetSystem_C",PETSC_NULL,(PetscVoidFunction)DMlibMeshSetSystem_libMesh);CHKERRQ(ierr);
1228  ierr = PetscObjectComposeFunction((PetscObject)dm,"DMlibMeshGetSystem_C",PETSC_NULL,(PetscVoidFunction)DMlibMeshGetSystem_libMesh);CHKERRQ(ierr);
1229 #else
1230  ierr = PetscObjectComposeFunction((PetscObject)dm,"DMlibMeshSetSystem_C",DMlibMeshSetSystem_libMesh);CHKERRQ(ierr);
1231  ierr = PetscObjectComposeFunction((PetscObject)dm,"DMlibMeshGetSystem_C",DMlibMeshGetSystem_libMesh);CHKERRQ(ierr);
1232 #endif
1233 
1235 }
1236 EXTERN_C_END
1237 
1238 
1239 #endif // #if !PETSC_VERSION_LESS_THAN(3,3,0)
DMCreateGlobalVector_libMesh
static PetscErrorCode DMCreateGlobalVector_libMesh(DM dm, Vec *x)
Definition: petscdmlibmeshimpl.C:980
libMesh::MeshBase::active_local_subdomain_elements_end
virtual element_iterator active_local_subdomain_elements_end(subdomain_id_type subdomain_id)=0
DMVec_libMesh
Definition: petscdmlibmeshimpl.C:71
dlm
DM_libMesh * dlm
Definition: petscdmlibmeshimpl.C:1029
DMView_libMesh
static PetscErrorCode DMView_libMesh(DM dm, PetscViewer viewer)
Definition: petscdmlibmeshimpl.C:1048
DM_libMesh::blockids
std::map< std::string, unsigned int > * blockids
Definition: petscdmlibmeshimpl.C:62
eq
PetscBool eq
Definition: petscdmlibmeshimpl.C:1030
DMCreateFieldDecomposition_libMesh
static PetscErrorCode DMCreateFieldDecomposition_libMesh(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
Definition: petscdmlibmeshimpl.C:272
libMesh::ExplicitSystem::rhs
NumericVector< Number > * rhs
The system matrix.
Definition: explicit_system.h:114
libMesh::DofMap::dof_indices
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1967
DMCreateDomainDecomposition_libMesh
static PetscErrorCode DMCreateDomainDecomposition_libMesh(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
Definition: petscdmlibmeshimpl.C:378
SNESJacobian_DMlibMesh
static PetscErrorCode SNESJacobian_DMlibMesh(#if PETSC_RELEASE_LESS_THAN(3, 5, 0) SNES, Vec x, Mat *jac, Mat *pc, MatStructure *flag, void *ctx #else SNES, Vec x, Mat jac, Mat pc, void *ctx #endif)
Definition: petscdmlibmeshimpl.C:927
DMlibMeshSetSystem_libMesh
PetscErrorCode DMlibMeshSetSystem_libMesh(DM dm, NonlinearImplicitSystem &sys)
Definition: petscdmlibmeshimpl.C:90
DM_libMesh::decomposition_type
unsigned int decomposition_type
Definition: petscdmlibmeshimpl.C:64
SNESFunction_DMlibMesh
static PetscErrorCode SNESFunction_DMlibMesh(SNES, Vec x, Vec r, void *ctx)
Definition: petscdmlibmeshimpl.C:843
DMVariableBounds_libMesh
static PetscErrorCode DMVariableBounds_libMesh(DM dm, Vec xl, Vec xu)
Definition: petscdmlibmeshimpl.C:949
libMesh::MeshBase::active_element_ptr_range
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
DMlibMeshSetUpName_Private
PetscErrorCode DMlibMeshSetUpName_Private(DM dm)
Definition: petscdmlibmeshimpl.C:218
libMesh::index_range
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:106
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
DM_libMesh::embedding
IS embedding
Definition: petscdmlibmeshimpl.C:67
end
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end.
Definition: variant_filter_iterator.h:343
libMesh::ParallelObject::comm
const Parallel::Communicator & comm() const
Definition: parallel_object.h:94
libMesh::DofMap::n_variables
unsigned int n_variables() const
Definition: dof_map.h:592
DMlibMeshJacobian
static PetscErrorCode DMlibMeshJacobian(#if PETSC_RELEASE_LESS_THAN(3, 5, 0) DM dm, Vec x, Mat jac, Mat pc, MatStructure *msflag #else DM dm, Vec x, Mat jac, Mat pc #endif)
Definition: petscdmlibmeshimpl.C:856
DMlibMeshGetBlocks
PetscErrorCode DMlibMeshGetBlocks(DM dm, PetscInt *n, char ***blocknames)
Definition: petscdmlibmeshimpl.C:168
DM_libMesh::decomposition
std::vector< std::set< unsigned int > > * decomposition
Definition: petscdmlibmeshimpl.C:65
mesh
MeshBase & mesh
Definition: mesh_communication.C:1257
DMSetUp_libMesh
static PetscErrorCode DMSetUp_libMesh(DM dm)
Definition: petscdmlibmeshimpl.C:1107
libMesh::DofMap::first_dof
dof_id_type first_dof(const processor_id_type proc) const
Definition: dof_map.h:650
token
Definition: petscdmlibmeshimpl.C:590
libMesh::ierr
ierr
Definition: petsc_dm_wrapper.C:72
libMesh::NonlinearImplicitSystem
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and non-linear solv...
Definition: nonlinear_implicit_system.h:54
DMlibMeshGetSystem_libMesh
PetscErrorCode DMlibMeshGetSystem_libMesh(DM dm, NonlinearImplicitSystem *&sys)
Definition: petscdmlibmeshimpl.C:151
PetscBool
PetscTruth PetscBool
Definition: petsc_macro.h:73
libMesh::NumericVector< Number >
libMesh::Variable::name
const std::string & name() const
Definition: variable.h:100
libMesh::libmesh_assert
libmesh_assert(ctx)
DMlibMeshCreateFieldDecompositionDM
PetscErrorCode DMlibMeshCreateFieldDecompositionDM(DM dm, PetscInt dnumber, PetscInt *dsizes, char ***dvarlists, DM *ddm)
Definition: petscdmlibmeshimpl.C:488
libMesh::IntRange
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
Definition: int_range.h:53
DM_libMesh::sys
NonlinearImplicitSystem * sys
Definition: petscdmlibmeshimpl.C:59
libMesh::MeshBase
This is the MeshBase class.
Definition: mesh_base.h:78
token::s
const char * s
Definition: petscdmlibmeshimpl.C:591
DM_libMesh::varids
std::map< std::string, unsigned int > * varids
Definition: petscdmlibmeshimpl.C:60
libMesh::DofMap::variable
const Variable & variable(const unsigned int c) const
Definition: dof_map.h:1894
libMesh::System::get_mesh
const MeshBase & get_mesh() const
Definition: system.h:2083
libMesh::DofMap::enforce_constraints_exactly
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=nullptr, bool homogeneous=false) const
Constrains the numeric vector v, which represents a solution defined on the mesh.
Definition: dof_map.h:2054
DMlibMeshGetSystem
PETSC_EXTERN PetscErrorCode DMlibMeshGetSystem(DM, libMesh::NonlinearImplicitSystem *&)
Definition: petscdmlibmesh.C:53
libMesh::System::current_local_solution
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1551
libMesh::ImplicitSystem::matrix
SparseMatrix< Number > * matrix
The system matrix.
Definition: implicit_system.h:393
DM_libMesh::vec_count
unsigned int vec_count
Definition: petscdmlibmeshimpl.C:68
DMlibMeshFunction
static PetscErrorCode DMlibMeshFunction(DM dm, Vec x, Vec r)
Definition: petscdmlibmeshimpl.C:786
DM_libMesh::blocknames
std::map< unsigned int, std::string > * blocknames
Definition: petscdmlibmeshimpl.C:63
libMesh::as_range
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
libMesh::PetscVector
This class provides a nice interface to PETSc's Vec object.
Definition: petsc_vector.h:72
A
static PetscErrorCode Mat * A
Definition: petscdmlibmeshimpl.C:1026
DMCreate_libMesh
PetscErrorCode DMCreate_libMesh(DM dm)
Definition: petscdmlibmeshimpl.C:1170
DM_libMesh::varnames
std::map< unsigned int, std::string > * varnames
Definition: petscdmlibmeshimpl.C:61
libMesh::MeshBase::const_element_iterator
The definition of the const_element_iterator struct.
Definition: mesh_base.h:1891
libMesh::numeric_index_type
dof_id_type numeric_index_type
Definition: id_types.h:99
libMesh::NonlinearImplicitSystem::nonlinear_solver
std::unique_ptr< NonlinearSolver< Number > > nonlinear_solver
The NonlinearSolver defines the default interface used to solve the nonlinear_implicit system.
Definition: nonlinear_implicit_system.h:261
DMCreateMatrix_libMesh
static PetscErrorCode DMCreateMatrix_libMesh(DM dm, const MatType, Mat *A) static PetscErrorCode DMCreateMatrix_libMesh(DM dm
DMlibMeshCreateDomainDecompositionDM
PetscErrorCode DMlibMeshCreateDomainDecompositionDM(DM dm, PetscInt dnumber, PetscInt *dsizes, char ***dblocklists, DM *ddm)
Definition: petscdmlibmeshimpl.C:540
PETSC_OWN_POINTER
Definition: petsc_macro.h:103
libMesh::System::solution
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1539
libMesh::ReferenceElem::get
const Elem & get(const ElemType type_in)
Definition: reference_elem.C:237
DM_libMesh
Definition: petscdmlibmeshimpl.C:57
DMCreateFieldDecompositionDM_libMesh
static PetscErrorCode DMCreateFieldDecompositionDM_libMesh(DM dm, const char *ddesc, DM *ddm)
Definition: petscdmlibmeshimpl.C:745
libMesh::PetscMatrix
This class provides a nice interface to the PETSc C-based data structures for parallel,...
Definition: petsc_matrix.h:87
libMesh::DofMap
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:176
libMesh::System::name
const std::string & name() const
Definition: system.h:2067
DMlibMeshGetVec_Private
PetscErrorCode DMlibMeshGetVec_Private(DM, const char *, Vec *)
Definition: petscdmlibmeshimpl.C:77
libMesh::MeshBase::subdomain_name
std::string & subdomain_name(subdomain_id_type id)
Definition: mesh_base.C:717
libMesh::ctx
void * ctx
Definition: petsc_dm_wrapper.C:71
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:100
DMlibMeshParseDecompositionDescriptor_Private
static PetscErrorCode DMlibMeshParseDecompositionDescriptor_Private(DM dm, const char *ddesc, PetscInt *dtype, PetscInt *dcount, PetscInt **dsizes, char ****dlists)
Definition: petscdmlibmeshimpl.C:601
DMDestroy_libMesh
static PetscErrorCode DMDestroy_libMesh(DM dm)
Definition: petscdmlibmeshimpl.C:1151
libMesh::PetscVector::vec
Vec vec()
Definition: petsc_vector.h:335
libMesh::System::get_dof_map
const DofMap & get_dof_map() const
Definition: system.h:2099
DMlibMeshGetVariables
PetscErrorCode DMlibMeshGetVariables(DM dm, PetscInt *n, char ***varnames)
Definition: petscdmlibmeshimpl.C:193
libMesh::TestClass
Definition: id_types.h:33
libMesh::CHKERRQ
CHKERRQ(ierr)
PETSC_ERR_ARG_WRONGSTATE
PETSC_ERR_ARG_WRONGSTATE
Definition: petscdmlibmeshimpl.C:1038
libMesh::PetscVector::swap
virtual void swap(NumericVector< T > &v) override
Swaps the contents of this with v.
Definition: petsc_vector.h:1194
DMCreateDomainDecompositionDM_libMesh
static PetscErrorCode DMCreateDomainDecompositionDM_libMesh(DM dm, const char *ddesc, DM *ddm)
Definition: petscdmlibmeshimpl.C:764
token::next
struct token * next
Definition: petscdmlibmeshimpl.C:592
libMesh::System::update
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:408
PetscFunctionReturn
PetscFunctionReturn(0)
DMVec_libMesh::label
std::string label
Definition: petscdmlibmeshimpl.C:72
DM_libMesh::embedding_type
unsigned int embedding_type
Definition: petscdmlibmeshimpl.C:66
libMesh::Quality::name
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
libMesh::MeshBase::active_local_subdomain_elements_begin
virtual element_iterator active_local_subdomain_elements_begin(subdomain_id_type subdomain_id)=0
libMesh::DofMap::end_dof
dof_id_type end_dof(const processor_id_type proc) const
Definition: dof_map.h:692