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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo