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