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