20 #include "libmesh/petsc_macro.h"
23 #if !PETSC_VERSION_LESS_THAN(3,3,0)
26 #if !PETSC_RELEASE_LESS_THAN(3,6,0)
27 # include <petsc/private/dmimpl.h>
29 # include <petsc-private/dmimpl.h>
33 #include "libmesh/libmesh_common.h"
34 #include "libmesh/nonlinear_implicit_system.h"
35 #include "libmesh/petsc_nonlinear_solver.h"
36 #include "libmesh/petsc_linear_solver.h"
37 #include "libmesh/petsc_vector.h"
38 #include "libmesh/petsc_matrix.h"
39 #include "libmesh/petscdmlibmesh.h"
40 #include "libmesh/dof_map.h"
41 #include "libmesh/preconditioner.h"
42 #include "libmesh/elem.h"
43 #include "libmesh/parallel.h"
49 #define DMLIBMESH_NO_DECOMPOSITION 0
50 #define DMLIBMESH_FIELD_DECOMPOSITION 1
51 #define DMLIBMESH_DOMAIN_DECOMPOSITION 2
53 #define DMLIBMESH_NO_EMBEDDING 0
54 #define DMLIBMESH_FIELD_EMBEDDING 1
55 #define DMLIBMESH_DOMAIN_EMBEDDING 2
60 std::map<std::string, unsigned int> *
varids;
61 std::map<unsigned int, std::string> *
varnames;
62 std::map<std::string, unsigned int> *
blockids;
76 #define __FUNCT__ "DMlibMeshGetVec_Private"
89 #define __FUNCT__ "DMlibMeshSetSystem_libMesh"
92 const Parallel::Communicator & comm = sys.
comm();
96 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
98 ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh);
99 if (!islibmesh) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
101 if (dm->setupcalled) SETERRQ(((PetscObject)dm)->comm,
PETSC_ERR_ARG_WRONGSTATE,
"Cannot reset the libMesh system after DM has been set up.");
111 dlm->
varids->insert(std::pair<std::string,unsigned int>(vname,v));
112 dlm->
varnames->insert(std::pair<unsigned int,std::string>(v,vname));
117 std::set<subdomain_id_type> blocks;
122 blocks.insert(elem->subdomain_id());
124 comm.set_union(blocks);
126 std::set<subdomain_id_type>::iterator bit = blocks.begin();
127 std::set<subdomain_id_type>::iterator bend = blocks.end();
128 if (bit == bend) SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_PLIB,
"No mesh blocks found.");
130 for (; bit != bend; ++bit) {
133 if (!bname.length()) {
138 std::ostringstream ss;
142 dlm->
blockids->insert(std::pair<std::string,unsigned int>(bname,bid));
143 dlm->
blocknames->insert(std::pair<unsigned int,std::string>(bid,bname));
150 #define __FUNCT__ "DMlibMeshGetSystem_libMesh"
156 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
158 ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh);
CHKERRQ(
ierr);
159 if (!islibmesh) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
167 #define __FUNCT__ "DMlibMeshGetBlocks"
173 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
175 ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh);
176 if (!islibmesh) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
178 PetscValidPointer(n,2);
179 *n = cast_int<unsigned int>(
dlm->
blockids->size());
185 ierr = PetscStrallocpy(pr.first.c_str(), *blocknames+i);
CHKERRQ(
ierr);
192 #define __FUNCT__ "DMlibMeshGetVariables"
197 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
200 ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh);
201 if (!islibmesh) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
203 PetscValidPointer(n,2);
204 *n = cast_int<unsigned int>(
dlm->
varids->size());
217 #define __FUNCT__ "DMlibMeshSetUpName_Private"
224 std::map<unsigned int, std::string> * dnames = PETSC_NULL, * enames = PETSC_NULL;
230 name +=
":dec:block:";
235 for (std::set<unsigned int>::iterator dit_begin = decomp.begin(),
237 dit_end = decomp.end();
238 dit != dit_end; ++dit) {
239 unsigned int id = *dit;
240 if (dit != dit_begin)
242 name += (*dnames)[id];
252 name +=
":emb:block:";
256 for (
auto eit = enames->begin(),
257 eit_end = enames->end(); eit != eit_end; ++eit) {
258 std::string & ename = eit->second;
259 if (eit != enames->begin())
271 #define __FUNCT__ "DMCreateFieldDecomposition_libMesh"
282 if (namelist) {
ierr = PetscMalloc(*len*
sizeof(
char *), namelist);
CHKERRQ(
ierr);}
283 if (islist) {
ierr = PetscMalloc(*len*
sizeof(IS), islist);
CHKERRQ(
ierr);}
284 if (dmlist) {
ierr = PetscMalloc(*len*
sizeof(DM), dmlist);
CHKERRQ(
ierr);}
287 std::set<numeric_index_type> dindices;
289 std::map<std::string, unsigned int> dvarids;
290 std::map<unsigned int, std::string> dvarnames;
291 unsigned int dvcount = 0;
294 dvarids.insert(std::pair<std::string, unsigned int>(vname,v));
295 dvarnames.insert(std::pair<unsigned int,std::string>(v,vname));
296 if (!dvcount) dname = vname;
297 else dname +=
"_" + vname;
299 if (!islist)
continue;
303 for (
const auto & elem :
307 std::vector<numeric_index_type> evindices;
312 dindices.insert(dof);
323 ierr = PetscMalloc(
sizeof(PetscInt)*dindices.size(), &darray);
CHKERRQ(
ierr);
325 for (
const auto &
id : dindices) {
329 ierr = ISCreateGeneral(((PetscObject)dm)->comm,
330 cast_int<PetscInt>(dindices.size()),
335 #if PETSC_RELEASE_LESS_THAN(3,3,1)
343 if (elen != dlen) SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB,
"Failed to embed subdomain %D", d);
377 #define __FUNCT__ "DMCreateDomainDecomposition_libMesh"
387 if (namelist) {
ierr = PetscMalloc(*len*
sizeof(
char *), namelist);
CHKERRQ(
ierr);}
388 if (innerislist) {
ierr = PetscMalloc(*len*
sizeof(IS), innerislist);
CHKERRQ(
ierr);}
389 if (outerislist) *outerislist = PETSC_NULL;
390 if (dmlist) {
ierr = PetscMalloc(*len*
sizeof(DM), dmlist);
CHKERRQ(
ierr);}
392 std::set<numeric_index_type> dindices;
394 std::map<std::string, unsigned int> dblockids;
395 std::map<unsigned int,std::string> dblocknames;
396 unsigned int dbcount = 0;
399 dblockids.insert(std::pair<std::string, unsigned int>(bname,b));
400 dblocknames.insert(std::pair<unsigned int,std::string>(b,bname));
401 if (!dbcount) dname = bname;
402 else dname +=
"_" + bname;
404 if (!innerislist)
continue;
408 for ( ; el != end_el; ++el) {
409 const Elem * elem = *el;
410 std::vector<numeric_index_type> evindices;
415 for (
const auto & dof : evindices) {
417 dindices.insert(dof);
428 ierr = PetscMalloc(
sizeof(PetscInt)*dindices.size(), &darray);
CHKERRQ(
ierr);
430 for (
const auto &
id : dindices) {
434 ierr = ISCreateGeneral(((PetscObject)dm)->comm,
435 cast_int<PetscInt>(dindices.size()),
440 #if PETSC_RELEASE_LESS_THAN(3,3,1)
448 if (elen != dlen) SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB,
"Failed to embed field %D", d);
457 (*innerislist)[d] = dis;
487 #define __FUNCT__ "DMlibMeshCreateFieldDecompositionDM"
493 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
494 ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh);
495 if (!islibmesh) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
496 if (dnumber < 0) SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"Negative number %D of decomposition parts", dnumber);
497 PetscValidPointer(ddm,5);
507 ddlm->
decomposition =
new(std::vector<std::set<unsigned int>>);
510 for (PetscInt d = 0; d < dnumber; ++d) {
511 if (dsizes[d] < 0) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"Negative size %D of decomposition part %D", dsizes[d],d);
513 for (PetscInt v = 0; v < dsizes[d]; ++v) {
514 std::string vname(dvarlists[d][v]);
515 std::map<std::string, unsigned int>::const_iterator vit =
dlm->
varids->find(vname);
517 SETERRQ3(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"Variable %D on the %D-th list with name %s is not owned by this DM", v, d, dvarlists[d][v]);
518 unsigned int vid = vit->second;
525 for (
const auto & pr : (*ddlm->
varids)) {
527 unsigned int vid = pr.second;
539 #define __FUNCT__ "DMlibMeshCreateDomainDecompositionDM"
545 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
546 ierr = PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh);
547 if (!islibmesh) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
548 if (dnumber < 0) SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"Negative number %D of decomposition parts", dnumber);
549 PetscValidPointer(ddm,5);
559 ddlm->
decomposition =
new(std::vector<std::set<unsigned int>>);
562 for (PetscInt d = 0; d < dnumber; ++d) {
563 if (dsizes[d] < 0) SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"Negative size %D of decomposition part %D", dsizes[d],d);
565 for (PetscInt b = 0; b < dsizes[d]; ++b) {
566 std::string bname(dblocklists[d][b]);
567 std::map<std::string, unsigned int>::const_iterator bit =
dlm->
blockids->find(bname);
569 SETERRQ3(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"Block %D on the %D-th list with name %s is not owned by this DM", b, d, dblocklists[d][b]);
570 unsigned int bid = bit->second;
577 for (
const auto & pr : (*ddlm->
blockids)) {
579 unsigned int bid = pr.second;
597 #if PETSC_RELEASE_LESS_THAN(3,3,1)
600 #define __FUNCT__ "DMlibMeshParseDecompositionDescriptor_Private"
609 struct token * llfirst = PETSC_NULL;
610 struct token * lllast = PETSC_NULL;
612 PetscInt stcount = 0, brcount = 0, d, i;
636 *dtype=DMLIBMESH_FIELD_DECOMPOSITION;
641 SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"Could not determine decomposition type from descriptor: %s\n", ddesc);
CHKERRQ(
ierr);
642 *dtype=DMLIBMESH_DOMAIN_DECOMPOSITION;
645 while (count < len0) {
650 while (*ss && *ss !=
',' && *ss !=
';') {
653 st = PETSC_NULL; br = PETSC_NULL;
662 #if PETSC_RELEASE_LESS_THAN(3,5,0)
671 #if PETSC_RELEASE_LESS_THAN(3,5,0)
681 llfirst = lllast = st;
684 lllast->
next = st; lllast = st;
689 llfirst = lllast = br;
692 lllast->
next = br; lllast = br;
711 for (d = 0; d < *dcount; ++d) (*dsizes)[d] = 0;
712 tok = llfirst; d = 0;
720 for (d = 0; d < *dcount; ++d) {
721 ierr = PetscMalloc(
sizeof(
char **)*(*dsizes)[d], (* dlists)+d);
CHKERRQ(
ierr);
724 tok = llfirst; d = 0; i = 0;
744 #define __FUNCT__ "DMCreateFieldDecompositionDM_libMesh"
749 PetscInt dtype, dcount;
755 if (dtype == DMLIBMESH_FIELD_DECOMPOSITION){
758 else SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB,
"Unexpected unknown decomposition type for field decomposition descriptor %s", ddesc);
763 #define __FUNCT__ "DMCreateDomainDecompositionDM_libMesh"
768 PetscInt dtype, dcount;
774 if (dtype == DMLIBMESH_DOMAIN_DECOMPOSITION) {
777 else SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB,
"Unexpected unknown decomposition type for domain decomposition descriptor %s", ddesc);
785 #define __FUNCT__ "DMlibMeshFunction"
801 X_global.
swap(X_sys);
808 X_global.swap(X_sys);
815 libmesh_error_msg(
"ERROR: cannot specify both a function and object to compute the Residual!");
818 libmesh_error_msg(
"ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
833 libmesh_error_msg(
"Error! Unable to compute residual and/or Jacobian!");
840 #if !PETSC_RELEASE_LESS_THAN(3,3,1)
842 #define __FUNCT__ "SNESFunction_DMlibMesh"
855 #define __FUNCT__ "DMlibMeshJacobian"
857 #
if PETSC_RELEASE_LESS_THAN(3,5,0)
858 DM dm, Vec x, Mat jac, Mat pc, MatStructure * msflag
860 DM dm, Vec x, Mat jac, Mat pc
881 X_global.
swap(X_sys);
887 X_global.swap(X_sys);
895 libmesh_error_msg(
"ERROR: cannot specify both a function and object to compute the Jacobian!");
898 libmesh_error_msg(
"ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
913 libmesh_error_msg(
"Error! Unable to compute residual and/or Jacobian!");
918 #if PETSC_RELEASE_LESS_THAN(3,5,0)
919 *msflag = SAME_NONZERO_PATTERN;
924 #if !PETSC_RELEASE_LESS_THAN(3,3,1)
926 #define __FUNCT__ "SNESJacobian_DMlibMesh"
928 #
if PETSC_RELEASE_LESS_THAN(3,5,0)
929 SNES, Vec x, Mat * jac, Mat * pc, MatStructure * flag,
void *
ctx
931 SNES, Vec x, Mat jac, Mat pc,
void *
ctx
938 #if PETSC_RELEASE_LESS_THAN(3,5,0)
948 #define __FUNCT__ "DMVariableBounds_libMesh"
958 #if PETSC_VERSION_LESS_THAN(3,5,0) && PETSC_VERSION_RELEASE
963 const PetscReal petsc_inf = std::numeric_limits<PetscReal>::max() / 4;
972 SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"No bounds calculation in this libMesh object");
979 #define __FUNCT__ "DMCreateGlobalVector_libMesh"
990 SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"DM of type %s, not of type %s", ((PetscObject)dm)->type, DMLIBMESH);
1012 ierr = PetscObjectCompose((PetscObject)*x,
"DM",(PetscObject)dm);
CHKERRQ(
ierr);
1020 #define __FUNCT__ "DMCreateMatrix_libMesh"
1021 #if PETSC_VERSION_LT(3,5,0)
1035 SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"DM of type %s, not of type %s", ((PetscObject)dm)->type, DMLIBMESH);
1047 #define __FUNCT__ "DMView_libMesh"
1050 PetscErrorCode
ierr;
1052 const char *
name, * prefix;
1055 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
CHKERRQ(
ierr);
1058 ierr = PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
CHKERRQ(
ierr);
1059 ierr = PetscViewerASCIIPrintf(viewer,
"DM libMesh with name %s and prefix %s\n",
name, prefix);
CHKERRQ(
ierr);
1061 std::map<std::string,unsigned int>::iterator bit =
dlm->
blockids->begin();
1062 std::map<std::string,unsigned int>::const_iterator bend =
dlm->
blockids->end();
1063 for (; bit != bend; ++bit) {
1064 ierr = PetscViewerASCIIPrintf(viewer,
"(%s,%D) ", bit->first.c_str(), bit->second);
CHKERRQ(
ierr);
1068 std::map<std::string,unsigned int>::iterator vit =
dlm->
varids->begin();
1069 std::map<std::string,unsigned int>::const_iterator vend =
dlm->
varids->end();
1070 for (; vit != vend; ++vit) {
1071 ierr = PetscViewerASCIIPrintf(viewer,
"(%s,%D) ", vit->first.c_str(), vit->second);
CHKERRQ(
ierr);
1075 ierr = PetscViewerASCIIPrintf(viewer,
"No decomposition\n");
CHKERRQ(
ierr);
1079 ierr = PetscViewerASCIIPrintf(viewer,
"Field decomposition by variable: ");
CHKERRQ(
ierr);
1082 ierr = PetscViewerASCIIPrintf(viewer,
"Domain decomposition by block: ");
CHKERRQ(
ierr);
1084 else SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB,
"Unexpected decomposition type: %D",
dlm->
decomposition_type);
1087 std::set<unsigned int>::iterator dbegin = (*
dlm->
decomposition)[d].begin();
1090 for (; dit != dend; ++dit) {
1091 if (dit != dbegin) {
1106 #define __FUNCT__ "DMSetUp_libMesh"
1110 PetscErrorCode
ierr;
1117 SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG,
"DM of type %s, not of type %s", ((PetscObject)dm)->type, DMLIBMESH);
1125 #if PETSC_RELEASE_LESS_THAN(3,3,1)
1140 dm->ops->createglobalvector = 0;
1141 dm->ops->creatematrix = 0;
1150 #define __FUNCT__ "DMDestroy_libMesh"
1154 PetscErrorCode
ierr;
1169 #define __FUNCT__ "DMCreate_libMesh"
1172 PetscErrorCode
ierr;
1176 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1177 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1185 dlm->
varids =
new(std::map<std::string, unsigned int>);
1186 dlm->
blockids =
new(std::map<std::string, unsigned int>);
1187 dlm->
varnames =
new(std::map<unsigned int, std::string>);
1194 dm->ops->createlocalvector = 0;
1195 dm->ops->getcoloring = 0;
1197 dm->ops->createinterpolation= 0;
1199 dm->ops->refine = 0;
1200 dm->ops->coarsen = 0;
1205 #if PETSC_RELEASE_LESS_THAN(3,12,0)
1206 dm->ops->getinjection = 0;
1207 dm->ops->getaggregates = 0;
1209 dm->ops->createinjection = 0;
1213 #if PETSC_RELEASE_LESS_THAN(3,3,1)
1222 dm->ops->setfromoptions = 0;
1226 #if PETSC_RELEASE_LESS_THAN(3,4,0)
1239 #endif // #if !PETSC_VERSION_LESS_THAN(3,3,0)