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