Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 :
22 : #ifdef LIBMESH_HAVE_PETSC
23 :
24 : #include "libmesh/ignore_warnings.h"
25 :
26 : // PETSc includes
27 : #if !PETSC_VERSION_LESS_THAN(3,6,0)
28 : # include <petsc/private/dmimpl.h>
29 : #else
30 : # include <petsc-private/dmimpl.h>
31 : #endif
32 :
33 : #include "libmesh/restore_warnings.h"
34 :
35 : // Local Includes
36 : #include "libmesh/libmesh_common.h"
37 : #include "libmesh/nonlinear_implicit_system.h"
38 : #include "libmesh/petsc_nonlinear_solver.h"
39 : #include "libmesh/petsc_linear_solver.h"
40 : #include "libmesh/petsc_vector.h"
41 : #include "libmesh/petsc_matrix_base.h"
42 : #include "libmesh/petscdmlibmesh.h"
43 : #include "libmesh/dof_map.h"
44 : #include "libmesh/preconditioner.h"
45 : #include "libmesh/elem.h"
46 : #include "libmesh/parallel.h"
47 :
48 :
49 : using namespace libMesh;
50 :
51 :
52 : #define DMLIBMESH_NO_DECOMPOSITION 0
53 : #define DMLIBMESH_FIELD_DECOMPOSITION 1
54 : #define DMLIBMESH_DOMAIN_DECOMPOSITION 2
55 :
56 : #define DMLIBMESH_NO_EMBEDDING 0
57 : #define DMLIBMESH_FIELD_EMBEDDING 1
58 : #define DMLIBMESH_DOMAIN_EMBEDDING 2
59 :
60 : struct DM_libMesh
61 : {
62 : NonlinearImplicitSystem * sys;
63 : std::map<std::string, unsigned int> * varids;
64 : std::map<unsigned int, std::string> * varnames;
65 : std::map<std::string, unsigned int> * blockids;
66 : std::map<unsigned int, std::string> * blocknames;
67 : unsigned int decomposition_type;
68 : std::vector<std::set<unsigned int>> * decomposition;
69 : unsigned int embedding_type;
70 : IS embedding;
71 : unsigned int vec_count;
72 : };
73 :
74 : struct DMVec_libMesh {
75 : std::string label;
76 : };
77 :
78 : #undef __FUNCT__
79 : #define __FUNCT__ "DMlibMeshGetVec_Private"
80 0 : PetscErrorCode DMlibMeshGetVec_Private(DM /*dm*/, const char * /*name*/, Vec *)
81 : {
82 : PetscFunctionBegin;
83 :
84 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
85 : }
86 :
87 :
88 :
89 : PetscErrorCode DMlibMeshSetUpName_Private(DM dm);
90 :
91 : #undef __FUNCT__
92 : #define __FUNCT__ "DMlibMeshSetSystem_libMesh"
93 29400 : PetscErrorCode DMlibMeshSetSystem_libMesh(DM dm, NonlinearImplicitSystem & sys)
94 : {
95 1680 : const Parallel::Communicator & comm = sys.comm();
96 :
97 : PetscFunctionBegin;
98 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
99 : PetscBool islibmesh;
100 29400 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh));
101 29400 : if (!islibmesh) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
102 :
103 29400 : if (dm->setupcalled) SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot reset the libMesh system after DM has been set up.");
104 29400 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
105 29400 : dlm->sys =&sys;
106 : /* Initially populate the sets of active blockids and varids using all of the
107 : existing blocks/variables (only variables are supported at the moment). */
108 840 : DofMap & dofmap = dlm->sys->get_dof_map();
109 29400 : dlm->varids->clear();
110 29400 : dlm->varnames->clear();
111 63140 : for (auto v : make_range(dofmap.n_variables())) {
112 33740 : std::string vname = dofmap.variable(v).name();
113 33740 : dlm->varids->insert(std::pair<std::string,unsigned int>(vname,v));
114 66516 : dlm->varnames->insert(std::pair<unsigned int,std::string>(v,vname));
115 : }
116 29400 : const MeshBase & mesh = dlm->sys->get_mesh();
117 29400 : dlm->blockids->clear();
118 29400 : dlm->blocknames->clear();
119 1680 : std::set<subdomain_id_type> blocks;
120 : /* The following effectively is a verbatim copy of MeshBase::n_subdomains(). */
121 : // This requires an inspection on every processor
122 840 : libmesh_parallel_only(mesh.comm());
123 60931982 : for (const auto & elem : mesh.active_element_ptr_range())
124 31361063 : blocks.insert(elem->subdomain_id());
125 : // Some subdomains may only live on other processors
126 29400 : comm.set_union(blocks);
127 :
128 840 : std::set<subdomain_id_type>::iterator bit = blocks.begin();
129 840 : std::set<subdomain_id_type>::iterator bend = blocks.end();
130 29400 : if (bit == bend) SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "No mesh blocks found.");
131 :
132 59360 : for (; bit != bend; ++bit) {
133 29960 : subdomain_id_type bid = *bit;
134 29960 : std::string bname = mesh.subdomain_name(bid);
135 29960 : if (!bname.length()) {
136 : /* Block names are currently implemented for Exodus II meshes
137 : only, so we might have to make up our own block names and
138 : maintain our own mapping of block ids to names.
139 : */
140 29960 : std::ostringstream ss;
141 29104 : ss << "dm" << bid;
142 30816 : bname = ss.str();
143 28248 : }
144 29960 : dlm->blockids->insert(std::pair<std::string,unsigned int>(bname,bid));
145 59064 : dlm->blocknames->insert(std::pair<unsigned int,std::string>(bid,bname));
146 : }
147 29400 : LibmeshPetscCallQ(DMlibMeshSetUpName_Private(dm));
148 840 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
149 : }
150 :
151 : #undef __FUNCT__
152 : #define __FUNCT__ "DMlibMeshGetSystem_libMesh"
153 14140 : PetscErrorCode DMlibMeshGetSystem_libMesh(DM dm, NonlinearImplicitSystem *& sys)
154 : {
155 : PetscFunctionBegin;
156 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
157 : PetscBool islibmesh;
158 14140 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh));
159 14140 : if (!islibmesh) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
160 14140 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
161 14140 : sys = dlm->sys;
162 14140 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
163 : }
164 :
165 :
166 : #undef __FUNCT__
167 : #define __FUNCT__ "DMlibMeshGetBlocks"
168 0 : PetscErrorCode DMlibMeshGetBlocks(DM dm, PetscInt * n, char *** blocknames)
169 : {
170 : PetscInt i;
171 : PetscFunctionBegin;
172 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
173 : PetscBool islibmesh;
174 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh));
175 0 : if (!islibmesh) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
176 0 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
177 : #if PETSC_RELEASE_GREATER_EQUALS(3,20,0)
178 : PetscAssertPointer(n,2);
179 : #else
180 : PetscValidPointer(n,2);
181 : #endif
182 0 : *n = cast_int<unsigned int>(dlm->blockids->size());
183 0 : if (!blocknames) PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
184 0 : LibmeshPetscCallQ(PetscMalloc(*n*sizeof(char *), blocknames));
185 0 : i = 0;
186 0 : for (const auto & pr : *(dlm->blockids))
187 : {
188 0 : LibmeshPetscCallQ(PetscStrallocpy(pr.first.c_str(), *blocknames+i));
189 0 : ++i;
190 : }
191 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
192 : }
193 :
194 : #undef __FUNCT__
195 : #define __FUNCT__ "DMlibMeshGetVariables"
196 0 : PetscErrorCode DMlibMeshGetVariables(DM dm, PetscInt * n, char *** varnames)
197 : {
198 : PetscFunctionBegin;
199 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
200 : PetscBool islibmesh;
201 : PetscInt i;
202 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh));
203 0 : if (!islibmesh) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
204 0 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
205 : #if PETSC_RELEASE_GREATER_EQUALS(3,20,0)
206 : PetscAssertPointer(n,2);
207 : #else
208 : PetscValidPointer(n,2);
209 : #endif
210 0 : *n = cast_int<unsigned int>(dlm->varids->size());
211 0 : if (!varnames) PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
212 0 : LibmeshPetscCallQ(PetscMalloc(*n*sizeof(char *), varnames));
213 0 : i = 0;
214 0 : for (const auto & pr : *(dlm->varids))
215 : {
216 0 : LibmeshPetscCallQ(PetscStrallocpy(pr.first.c_str(), *varnames+i));
217 0 : ++i;
218 : }
219 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
220 : }
221 :
222 : #undef __FUNCT__
223 : #define __FUNCT__ "DMlibMeshSetUpName_Private"
224 29400 : PetscErrorCode DMlibMeshSetUpName_Private(DM dm)
225 : {
226 29400 : DM_libMesh * dlm = (DM_libMesh *)dm->data;
227 : PetscFunctionBegin;
228 29400 : std::string name = dlm->sys->name();
229 840 : std::map<unsigned int, std::string> * dnames = LIBMESH_PETSC_NULLPTR,
230 840 : * enames = LIBMESH_PETSC_NULLPTR;
231 29400 : if (dlm->decomposition_type == DMLIBMESH_FIELD_DECOMPOSITION) {
232 0 : name += ":dec:var:";
233 0 : dnames = dlm->varnames;
234 : }
235 29400 : if (dlm->decomposition_type == DMLIBMESH_DOMAIN_DECOMPOSITION) {
236 0 : name += ":dec:block:";
237 0 : dnames = dlm->blocknames;
238 : }
239 29400 : if (dnames) {
240 0 : for (auto decomp : *dlm->decomposition) {
241 0 : for (std::set<unsigned int>::iterator dit_begin = decomp.begin(),
242 0 : dit = dit_begin,
243 0 : dit_end = decomp.end();
244 0 : dit != dit_end; ++dit) {
245 0 : unsigned int id = *dit;
246 0 : if (dit != dit_begin)
247 0 : name += ",";
248 0 : name += (*dnames)[id];
249 : }
250 0 : name += ";";
251 : }
252 : }
253 29400 : if (dlm->embedding_type == DMLIBMESH_FIELD_EMBEDDING) {
254 0 : name += ":emb:var:";
255 0 : enames = dlm->varnames;
256 : }
257 29400 : if (dlm->embedding_type == DMLIBMESH_DOMAIN_EMBEDDING) {
258 0 : name += ":emb:block:";
259 0 : enames = dlm->blocknames;
260 : }
261 29400 : if (enames) {
262 0 : for (auto eit = enames->begin(),
263 0 : eit_end = enames->end(); eit != eit_end; ++eit) {
264 0 : std::string & ename = eit->second;
265 0 : if (eit != enames->begin())
266 0 : name += ",";
267 0 : name += ename;
268 : }
269 0 : name += ";";
270 : }
271 29400 : LibmeshPetscCallQ(PetscObjectSetName((PetscObject)dm, name.c_str()));
272 30240 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
273 : }
274 :
275 :
276 : #undef __FUNCT__
277 : #define __FUNCT__ "DMCreateFieldDecomposition_libMesh"
278 0 : static PetscErrorCode DMCreateFieldDecomposition_libMesh(DM dm, PetscInt * len, char *** namelist, IS ** islist, DM ** dmlist)
279 : {
280 : PetscFunctionBegin;
281 0 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
282 0 : NonlinearImplicitSystem * sys = dlm->sys;
283 : IS emb;
284 0 : if (dlm->decomposition_type != DMLIBMESH_FIELD_DECOMPOSITION) PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
285 :
286 0 : *len = cast_int<unsigned int>(dlm->decomposition->size());
287 0 : if (namelist) {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(char *), namelist));}
288 0 : if (islist) {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(IS), islist));}
289 0 : if (dmlist) {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(DM), dmlist));}
290 0 : DofMap & dofmap = dlm->sys->get_dof_map();
291 0 : for (auto d : index_range(*dlm->decomposition)) {
292 0 : std::set<numeric_index_type> dindices;
293 0 : std::string dname;
294 0 : std::map<std::string, unsigned int> dvarids;
295 0 : std::map<unsigned int, std::string> dvarnames;
296 0 : unsigned int dvcount = 0;
297 0 : for (const auto & v : (*dlm->decomposition)[d]) {
298 0 : std::string vname = (*dlm->varnames)[v];
299 0 : dvarids.insert(std::pair<std::string, unsigned int>(vname,v));
300 0 : dvarnames.insert(std::pair<unsigned int,std::string>(v,vname));
301 0 : if (!dvcount) dname = vname;
302 0 : else dname += "_" + vname;
303 0 : ++dvcount;
304 0 : if (!islist) continue;
305 : // Iterate only over this DM's blocks.
306 0 : for (const auto & pr : *(dlm->blockids)) {
307 0 : const subdomain_id_type sbd_id = cast_int<subdomain_id_type>(pr.second);
308 0 : for (const auto & elem :
309 0 : as_range(sys->get_mesh().active_local_subdomain_elements_begin(sbd_id),
310 0 : sys->get_mesh().active_local_subdomain_elements_end(sbd_id))) {
311 : //unsigned int e_subdomain = elem->subdomain_id();
312 0 : std::vector<numeric_index_type> evindices;
313 : // Get the degree of freedom indices for the given variable off the current element.
314 0 : dofmap.dof_indices(elem, evindices, v);
315 0 : for (numeric_index_type dof : evindices) {
316 0 : if (dof >= dofmap.first_dof() && dof < dofmap.end_dof()) // might want to use variable_first/last_local_dof instead
317 0 : dindices.insert(dof);
318 : }
319 0 : }
320 : }
321 : }
322 0 : if (namelist) {
323 0 : LibmeshPetscCallQ(PetscStrallocpy(dname.c_str(),(*namelist)+d));
324 : }
325 0 : if (islist) {
326 : IS dis;
327 : PetscInt * darray;
328 0 : LibmeshPetscCallQ(PetscMalloc(sizeof(PetscInt)*dindices.size(), &darray));
329 0 : numeric_index_type i = 0;
330 0 : for (const auto & id : dindices) {
331 0 : darray[i] = id;
332 0 : ++i;
333 : }
334 0 : LibmeshPetscCallQ(ISCreateGeneral(((PetscObject)dm)->comm,
335 : cast_int<PetscInt>(dindices.size()),
336 : darray, PETSC_OWN_POINTER, &dis));
337 0 : if (dlm->embedding) {
338 : /* Create a relative embedding into the parent's index space. */
339 0 : LibmeshPetscCallQ(ISEmbed(dis,dlm->embedding, PETSC_TRUE, &emb));
340 : PetscInt elen, dlen;
341 0 : LibmeshPetscCallQ(ISGetLocalSize(emb, &elen));
342 0 : LibmeshPetscCallQ(ISGetLocalSize(dis, &dlen));
343 0 : if (elen != dlen) LIBMESH_SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Failed to embed subdomain %zu", d);
344 0 : LibmeshPetscCallQ(ISDestroy(&dis));
345 0 : dis = emb;
346 : }
347 : else {
348 0 : emb = dis;
349 : }
350 0 : (*islist)[d] = dis;
351 : }
352 0 : if (dmlist) {
353 : DM ddm;
354 0 : LibmeshPetscCallQ(DMCreate(((PetscObject)dm)->comm, &ddm));
355 0 : LibmeshPetscCallQ(DMSetType(ddm, DMLIBMESH));
356 0 : DM_libMesh * ddlm = (DM_libMesh *)(ddm->data);
357 0 : ddlm->sys = dlm->sys;
358 : /* copy over the block ids and names */
359 0 : *ddlm->blockids = *dlm->blockids;
360 0 : *ddlm->blocknames = *dlm->blocknames;
361 : /* set the vars from the d-th part of the decomposition. */
362 0 : *ddlm->varids = dvarids;
363 0 : *ddlm->varnames = dvarnames;
364 0 : LibmeshPetscCallQ(PetscObjectReference((PetscObject)emb));
365 0 : ddlm->embedding = emb;
366 0 : ddlm->embedding_type = DMLIBMESH_FIELD_EMBEDDING;
367 :
368 0 : LibmeshPetscCallQ(DMlibMeshSetUpName_Private(ddm));
369 0 : LibmeshPetscCallQ(DMSetFromOptions(ddm));
370 0 : (*dmlist)[d] = ddm;
371 : }
372 : }
373 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
374 : }
375 :
376 : #undef __FUNCT__
377 : #define __FUNCT__ "DMCreateDomainDecomposition_libMesh"
378 0 : static PetscErrorCode DMCreateDomainDecomposition_libMesh(DM dm, PetscInt * len, char *** namelist, IS ** innerislist, IS ** outerislist, DM ** dmlist)
379 : {
380 : PetscFunctionBegin;
381 0 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
382 0 : NonlinearImplicitSystem * sys = dlm->sys;
383 : IS emb;
384 0 : if (dlm->decomposition_type != DMLIBMESH_DOMAIN_DECOMPOSITION) PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
385 0 : *len = cast_int<unsigned int>(dlm->decomposition->size());
386 0 : if (namelist) {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(char *), namelist));}
387 0 : if (innerislist) {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(IS), innerislist));}
388 0 : if (outerislist) *outerislist = LIBMESH_PETSC_NULLPTR; /* FIX: allow mesh-based overlap. */
389 0 : if (dmlist) {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(DM), dmlist));}
390 0 : for (auto d : index_range(*dlm->decomposition)) {
391 0 : std::set<numeric_index_type> dindices;
392 0 : std::string dname;
393 0 : std::map<std::string, unsigned int> dblockids;
394 0 : std::map<unsigned int,std::string> dblocknames;
395 0 : unsigned int dbcount = 0;
396 0 : for (const auto & b : (*dlm->decomposition)[d]) {
397 0 : std::string bname = (*dlm->blocknames)[b];
398 0 : dblockids.insert(std::pair<std::string, unsigned int>(bname,b));
399 0 : dblocknames.insert(std::pair<unsigned int,std::string>(b,bname));
400 0 : if (!dbcount) dname = bname;
401 0 : else dname += "_" + bname;
402 0 : ++dbcount;
403 0 : if (!innerislist) continue;
404 0 : const subdomain_id_type b_sbd_id = cast_int<subdomain_id_type>(b);
405 0 : MeshBase::const_element_iterator el = sys->get_mesh().active_local_subdomain_elements_begin(b_sbd_id);
406 0 : const MeshBase::const_element_iterator end_el = sys->get_mesh().active_local_subdomain_elements_end(b_sbd_id);
407 0 : for ( ; el != end_el; ++el) {
408 0 : const Elem * elem = *el;
409 0 : std::vector<numeric_index_type> evindices;
410 : // Iterate only over this DM's variables.
411 0 : for (const auto & pr : *(dlm->varids)) {
412 : // Get the degree of freedom indices for the given variable off the current element.
413 0 : sys->get_dof_map().dof_indices(elem, evindices, pr.second);
414 0 : for (const auto & dof : evindices) {
415 0 : 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
416 0 : dindices.insert(dof);
417 : }
418 : }
419 : }
420 : }
421 0 : if (namelist) {
422 0 : LibmeshPetscCallQ(PetscStrallocpy(dname.c_str(),(*namelist)+d));
423 : }
424 0 : if (innerislist) {
425 : PetscInt * darray;
426 : IS dis;
427 0 : LibmeshPetscCallQ(PetscMalloc(sizeof(PetscInt)*dindices.size(), &darray));
428 0 : numeric_index_type i = 0;
429 0 : for (const auto & id : dindices) {
430 0 : darray[i] = id;
431 0 : ++i;
432 : }
433 0 : LibmeshPetscCallQ(ISCreateGeneral(((PetscObject)dm)->comm,
434 : cast_int<PetscInt>(dindices.size()),
435 : darray, PETSC_OWN_POINTER, &dis));
436 0 : if (dlm->embedding) {
437 : /* Create a relative embedding into the parent's index space. */
438 0 : LibmeshPetscCallQ(ISEmbed(dis,dlm->embedding, PETSC_TRUE, &emb));
439 : PetscInt elen, dlen;
440 0 : LibmeshPetscCallQ(ISGetLocalSize(emb, &elen));
441 0 : LibmeshPetscCallQ(ISGetLocalSize(dis, &dlen));
442 0 : if (elen != dlen) LIBMESH_SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Failed to embed field %zu" , d);
443 0 : LibmeshPetscCallQ(ISDestroy(&dis));
444 0 : dis = emb;
445 : }
446 : else {
447 0 : emb = dis;
448 : }
449 0 : if (innerislist) {
450 0 : LibmeshPetscCallQ(PetscObjectReference((PetscObject)dis));
451 0 : (*innerislist)[d] = dis;
452 : }
453 0 : LibmeshPetscCallQ(ISDestroy(&dis));
454 : }
455 0 : if (dmlist) {
456 : DM ddm;
457 0 : LibmeshPetscCallQ(DMCreate(((PetscObject)dm)->comm, &ddm));
458 0 : LibmeshPetscCallQ(DMSetType(ddm, DMLIBMESH));
459 0 : DM_libMesh * ddlm = (DM_libMesh *)(ddm->data);
460 0 : ddlm->sys = dlm->sys;
461 : /* copy over the varids and varnames */
462 0 : *ddlm->varids = *dlm->varids;
463 0 : *ddlm->varnames = *dlm->varnames;
464 : /* set the blocks from the d-th part of the decomposition. */
465 0 : *ddlm->blockids = dblockids;
466 0 : *ddlm->blocknames = dblocknames;
467 0 : LibmeshPetscCallQ(PetscObjectReference((PetscObject)emb));
468 0 : ddlm->embedding = emb;
469 0 : ddlm->embedding_type = DMLIBMESH_DOMAIN_EMBEDDING;
470 :
471 0 : LibmeshPetscCallQ(DMlibMeshSetUpName_Private(ddm));
472 0 : LibmeshPetscCallQ(DMSetFromOptions(ddm));
473 0 : (*dmlist)[d] = ddm;
474 : }
475 : }
476 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
477 : }
478 :
479 :
480 : #undef __FUNCT__
481 : #define __FUNCT__ "DMlibMeshCreateFieldDecompositionDM"
482 0 : PetscErrorCode DMlibMeshCreateFieldDecompositionDM(DM dm, PetscInt dnumber, PetscInt * dsizes, char *** dvarlists, DM * ddm)
483 : {
484 : PetscBool islibmesh;
485 : PetscFunctionBegin;
486 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
487 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh));
488 0 : if (!islibmesh) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
489 0 : if (dnumber < 0) LIBMESH_SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative number %" LIBMESH_PETSCINT_FMT " of decomposition parts", dnumber);
490 : #if PETSC_RELEASE_GREATER_EQUALS(3,20,0)
491 : PetscAssertPointer(ddm,5);
492 : #else
493 : PetscValidPointer(ddm,5);
494 : #endif
495 0 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
496 0 : LibmeshPetscCallQ(DMCreate(((PetscObject)dm)->comm, ddm));
497 0 : LibmeshPetscCallQ(DMSetType(*ddm, DMLIBMESH));
498 0 : DM_libMesh * ddlm = (DM_libMesh *)((*ddm)->data);
499 0 : ddlm->sys = dlm->sys;
500 0 : ddlm->varids = dlm->varids;
501 0 : ddlm->varnames = dlm->varnames;
502 0 : ddlm->blockids = dlm->blockids;
503 0 : ddlm->blocknames = dlm->blocknames;
504 0 : ddlm->decomposition = new(std::vector<std::set<unsigned int>>);
505 0 : ddlm->decomposition_type = DMLIBMESH_FIELD_DECOMPOSITION;
506 0 : if (dnumber) {
507 0 : for (PetscInt d = 0; d < dnumber; ++d) {
508 0 : if (dsizes[d] < 0) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative size %" LIBMESH_PETSCINT_FMT " of decomposition part %" LIBMESH_PETSCINT_FMT, dsizes[d],d);
509 0 : ddlm->decomposition->push_back(std::set<unsigned int>());
510 0 : for (PetscInt v = 0; v < dsizes[d]; ++v) {
511 0 : std::string vname(dvarlists[d][v]);
512 0 : std::map<std::string, unsigned int>::const_iterator vit = dlm->varids->find(vname);
513 0 : if (vit == dlm->varids->end())
514 0 : LIBMESH_SETERRQ3(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Variable %" LIBMESH_PETSCINT_FMT " on the %" LIBMESH_PETSCINT_FMT "-th list with name %s is not owned by this DM", v, d, dvarlists[d][v]);
515 0 : unsigned int vid = vit->second;
516 0 : (*ddlm->decomposition)[d].insert(vid);
517 : }
518 : }
519 : }
520 : else { // Empty splits indicate default: split all variables with one per split.
521 0 : PetscInt d = 0;
522 0 : for (const auto & pr : (*ddlm->varids)) {
523 0 : ddlm->decomposition->push_back(std::set<unsigned int>());
524 0 : unsigned int vid = pr.second;
525 0 : (*ddlm->decomposition)[d].insert(vid);
526 0 : ++d;
527 : }
528 : }
529 0 : LibmeshPetscCallQ(DMlibMeshSetUpName_Private(*ddm));
530 0 : LibmeshPetscCallQ(DMSetFromOptions(*ddm));
531 0 : LibmeshPetscCallQ(DMSetUp(*ddm));
532 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
533 : }
534 :
535 : #undef __FUNCT__
536 : #define __FUNCT__ "DMlibMeshCreateDomainDecompositionDM"
537 0 : PetscErrorCode DMlibMeshCreateDomainDecompositionDM(DM dm, PetscInt dnumber, PetscInt * dsizes, char *** dblocklists, DM * ddm)
538 : {
539 : PetscBool islibmesh;
540 : PetscFunctionBegin;
541 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
542 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh));
543 0 : if (!islibmesh) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
544 0 : if (dnumber < 0) LIBMESH_SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative number %" LIBMESH_PETSCINT_FMT " of decomposition parts", dnumber);
545 : #if PETSC_RELEASE_GREATER_EQUALS(3,20,0)
546 : PetscAssertPointer(ddm,5);
547 : #else
548 : PetscValidPointer(ddm,5);
549 : #endif
550 0 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
551 0 : LibmeshPetscCallQ(DMCreate(((PetscObject)dm)->comm, ddm));
552 0 : LibmeshPetscCallQ(DMSetType(*ddm, DMLIBMESH));
553 0 : DM_libMesh * ddlm = (DM_libMesh *)((*ddm)->data);
554 0 : ddlm->sys = dlm->sys;
555 0 : ddlm->varids = dlm->varids;
556 0 : ddlm->varnames = dlm->varnames;
557 0 : ddlm->blockids = dlm->blockids;
558 0 : ddlm->blocknames = dlm->blocknames;
559 0 : ddlm->decomposition = new(std::vector<std::set<unsigned int>>);
560 0 : ddlm->decomposition_type = DMLIBMESH_DOMAIN_DECOMPOSITION;
561 0 : if (dnumber) {
562 0 : for (PetscInt d = 0; d < dnumber; ++d) {
563 0 : if (dsizes[d] < 0) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative size %" LIBMESH_PETSCINT_FMT " of decomposition part %" LIBMESH_PETSCINT_FMT, dsizes[d],d);
564 0 : ddlm->decomposition->push_back(std::set<unsigned int>());
565 0 : for (PetscInt b = 0; b < dsizes[d]; ++b) {
566 0 : std::string bname(dblocklists[d][b]);
567 0 : std::map<std::string, unsigned int>::const_iterator bit = dlm->blockids->find(bname);
568 0 : if (bit == dlm->blockids->end())
569 0 : LIBMESH_SETERRQ3(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Block %" LIBMESH_PETSCINT_FMT " on the %" LIBMESH_PETSCINT_FMT "-th list with name %s is not owned by this DM", b, d, dblocklists[d][b]);
570 0 : unsigned int bid = bit->second;
571 0 : (*ddlm->decomposition)[d].insert(bid);
572 : }
573 : }
574 : }
575 : else { // Empty splits indicate default: split all blocks with one per split.
576 0 : PetscInt d = 0;
577 0 : for (const auto & pr : (*ddlm->blockids)) {
578 0 : ddlm->decomposition->push_back(std::set<unsigned int>());
579 0 : unsigned int bid = pr.second;
580 0 : (*ddlm->decomposition)[d].insert(bid);
581 0 : ++d;
582 : }
583 : }
584 0 : LibmeshPetscCallQ(DMlibMeshSetUpName_Private(*ddm));
585 0 : LibmeshPetscCallQ(DMSetFromOptions(*ddm));
586 0 : LibmeshPetscCallQ(DMSetUp(*ddm));
587 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
588 : }
589 :
590 : struct token {
591 : const char * s;
592 : struct token * next;
593 : };
594 :
595 :
596 :
597 : #undef __FUNCT__
598 : #define __FUNCT__ "DMlibMeshFunction"
599 0 : static PetscErrorCode DMlibMeshFunction(DM dm, Vec x, Vec r)
600 : {
601 : PetscFunctionBegin;
602 0 : libmesh_assert(x);
603 0 : libmesh_assert(r);
604 :
605 : NonlinearImplicitSystem * _sys;
606 0 : LibmeshPetscCallQ(DMlibMeshGetSystem(dm, _sys));
607 0 : NonlinearImplicitSystem & sys = *_sys;
608 0 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
609 0 : PetscVector<Number> & R_sys = *cast_ptr<PetscVector<Number> *>(sys.rhs);
610 0 : PetscVector<Number> X_global(x, _sys->comm()), R(r, _sys->comm());
611 :
612 : // Use the systems update() to get a good local version of the parallel solution
613 0 : X_global.swap(X_sys);
614 0 : R.swap(R_sys);
615 :
616 0 : _sys->get_dof_map().enforce_constraints_exactly(*_sys);
617 0 : _sys->update();
618 :
619 : // Swap back
620 0 : X_global.swap(X_sys);
621 0 : R.swap(R_sys);
622 0 : R.zero();
623 :
624 : // if the user has provided both function pointers and objects only the pointer
625 : // will be used, so catch that as an error
626 0 : libmesh_error_msg_if(_sys->nonlinear_solver->residual && _sys->nonlinear_solver->residual_object,
627 : "ERROR: cannot specify both a function and object to compute the Residual!");
628 :
629 0 : libmesh_error_msg_if(_sys->nonlinear_solver->matvec && _sys->nonlinear_solver->residual_and_jacobian_object,
630 : "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
631 :
632 0 : if (_sys->nonlinear_solver->residual != nullptr)
633 0 : _sys->nonlinear_solver->residual(*(_sys->current_local_solution.get()), R, *_sys);
634 :
635 0 : else if (_sys->nonlinear_solver->residual_object != nullptr)
636 0 : _sys->nonlinear_solver->residual_object->residual(*(_sys->current_local_solution.get()), R, *_sys);
637 :
638 0 : else if (_sys->nonlinear_solver->matvec != nullptr)
639 0 : _sys->nonlinear_solver->matvec(*(_sys->current_local_solution.get()), &R, nullptr, *_sys);
640 :
641 0 : else if (_sys->nonlinear_solver->residual_and_jacobian_object != nullptr)
642 0 : _sys->nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian(*(_sys->current_local_solution.get()), &R, nullptr, *_sys);
643 :
644 : else
645 0 : libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
646 :
647 0 : R.close();
648 0 : X_global.close();
649 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
650 0 : }
651 :
652 : #undef __FUNCT__
653 : #define __FUNCT__ "SNESFunction_DMlibMesh"
654 0 : static PetscErrorCode SNESFunction_DMlibMesh(SNES, Vec x, Vec r, void * ctx)
655 : {
656 0 : DM dm = (DM)ctx;
657 : PetscFunctionBegin;
658 0 : LibmeshPetscCallQ(DMlibMeshFunction(dm,x,r));
659 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
660 : }
661 :
662 :
663 : #undef __FUNCT__
664 : #define __FUNCT__ "DMlibMeshJacobian"
665 0 : static PetscErrorCode DMlibMeshJacobian(DM dm, Vec x, Mat jac, Mat pc)
666 : {
667 : PetscFunctionBegin;
668 : NonlinearImplicitSystem * _sys;
669 0 : LibmeshPetscCallQ(DMlibMeshGetSystem(dm, _sys));
670 0 : NonlinearImplicitSystem & sys = *_sys;
671 :
672 0 : libmesh_assert(pc);
673 0 : libmesh_assert(jac);
674 0 : PetscMatrixBase<Number> & the_pc = *PetscMatrixBase<Number>::get_context(pc, sys.comm());
675 0 : PetscMatrixBase<Number> & Jac = *PetscMatrixBase<Number>::get_context(jac, sys.comm());
676 0 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
677 0 : PetscMatrixBase<Number> & Jac_sys = *cast_ptr<PetscMatrixBase<Number> *>(sys.matrix);
678 0 : PetscVector<Number> X_global(x, sys.comm());
679 :
680 : // Set the dof maps
681 0 : the_pc.attach_dof_map(sys.get_dof_map());
682 0 : Jac.attach_dof_map(sys.get_dof_map());
683 :
684 : // Use the systems update() to get a good local version of the parallel solution
685 0 : X_global.swap(X_sys);
686 0 : Jac.swap(Jac_sys);
687 :
688 0 : sys.get_dof_map().enforce_constraints_exactly(sys);
689 0 : sys.update();
690 :
691 0 : X_global.swap(X_sys);
692 0 : Jac.swap(Jac_sys);
693 :
694 0 : the_pc.zero();
695 :
696 : // if the user has provided both function pointers and objects only the pointer
697 : // will be used, so catch that as an error
698 0 : libmesh_error_msg_if(sys.nonlinear_solver->jacobian && sys.nonlinear_solver->jacobian_object,
699 : "ERROR: cannot specify both a function and object to compute the Jacobian!");
700 :
701 0 : libmesh_error_msg_if(sys.nonlinear_solver->matvec && sys.nonlinear_solver->residual_and_jacobian_object,
702 : "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
703 :
704 0 : if (sys.nonlinear_solver->jacobian != nullptr)
705 0 : sys.nonlinear_solver->jacobian(*(sys.current_local_solution.get()), the_pc, sys);
706 :
707 0 : else if (sys.nonlinear_solver->jacobian_object != nullptr)
708 0 : sys.nonlinear_solver->jacobian_object->jacobian(*(sys.current_local_solution.get()), the_pc, sys);
709 :
710 0 : else if (sys.nonlinear_solver->matvec != nullptr)
711 0 : sys.nonlinear_solver->matvec(*(sys.current_local_solution.get()), nullptr, &the_pc, sys);
712 :
713 0 : else if (sys.nonlinear_solver->residual_and_jacobian_object != nullptr)
714 0 : sys.nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian(*(sys.current_local_solution.get()), nullptr, &the_pc, sys);
715 :
716 : else
717 0 : libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
718 :
719 0 : the_pc.close();
720 0 : Jac.close();
721 0 : X_global.close();
722 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
723 0 : }
724 :
725 : #undef __FUNCT__
726 : #define __FUNCT__ "SNESJacobian_DMlibMesh"
727 0 : static PetscErrorCode SNESJacobian_DMlibMesh(SNES, Vec x, Mat jac, Mat pc, void * ctx)
728 : {
729 0 : DM dm = (DM)ctx;
730 : PetscFunctionBegin;
731 0 : LibmeshPetscCallQ(DMlibMeshJacobian(dm,x,jac,pc));
732 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
733 : }
734 :
735 : #undef __FUNCT__
736 : #define __FUNCT__ "DMVariableBounds_libMesh"
737 14140 : static PetscErrorCode DMVariableBounds_libMesh(DM dm, Vec xl, Vec xu)
738 : {
739 : NonlinearImplicitSystem * _sys;
740 14140 : LibmeshPetscCallQ(DMlibMeshGetSystem(dm, _sys));
741 14140 : NonlinearImplicitSystem & sys = *_sys;
742 14948 : PetscVector<Number> XL(xl, sys.comm());
743 14948 : PetscVector<Number> XU(xu, sys.comm());
744 : PetscFunctionBegin;
745 : // Workaround for nonstandard Q suffix warning with quad precision
746 404 : const PetscReal petsc_inf = std::numeric_limits<PetscReal>::max() / 4;
747 14140 : LibmeshPetscCallQ(VecSet(xl, -petsc_inf));
748 14140 : LibmeshPetscCallQ(VecSet(xu, petsc_inf));
749 14140 : if (sys.nonlinear_solver->bounds != nullptr)
750 0 : sys.nonlinear_solver->bounds(XL,XU,sys);
751 14140 : else if (sys.nonlinear_solver->bounds_object != nullptr)
752 14140 : sys.nonlinear_solver->bounds_object->bounds(XL,XU, sys);
753 : else
754 0 : SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "No bounds calculation in this libMesh object");
755 :
756 404 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
757 13332 : }
758 :
759 :
760 : #undef __FUNCT__
761 : #define __FUNCT__ "DMCreateGlobalVector_libMesh"
762 14140 : static PetscErrorCode DMCreateGlobalVector_libMesh(DM dm, Vec *x)
763 : {
764 : PetscFunctionBegin;
765 14140 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
766 : PetscBool eq;
767 :
768 14140 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH, &eq));
769 :
770 14140 : if (!eq)
771 0 : LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
772 :
773 14140 : if (!dlm->sys)
774 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No libMesh system set for DM_libMesh");
775 :
776 404 : NumericVector<Number> * nv = (dlm->sys->solution).get();
777 14140 : PetscVector<Number> * pv = dynamic_cast<PetscVector<Number> *>(nv);
778 808 : Vec v = pv->vec();
779 : /* Unfortunately, currently this does not produce a ghosted vector, so nonlinear subproblem solves aren't going to be easily available.
780 : Should work fine for getting vectors out for linear subproblem solvers. */
781 14140 : if (dlm->embedding) {
782 : PetscInt n;
783 0 : LibmeshPetscCallQ(VecCreate(((PetscObject)v)->comm, x));
784 0 : LibmeshPetscCallQ(ISGetLocalSize(dlm->embedding, &n));
785 0 : LibmeshPetscCallQ(VecSetSizes(*x,n,PETSC_DETERMINE));
786 0 : LibmeshPetscCallQ(VecSetType(*x,((PetscObject)v)->type_name));
787 0 : LibmeshPetscCallQ(VecSetFromOptions(*x));
788 0 : LibmeshPetscCallQ(VecSetUp(*x));
789 : }
790 : else {
791 14140 : LibmeshPetscCallQ(VecDuplicate(v,x));
792 : }
793 :
794 : #if PETSC_VERSION_LESS_THAN(3,13,0)
795 808 : LibmeshPetscCallQ(PetscObjectCompose((PetscObject)*x,"DM",(PetscObject)dm));
796 : #else
797 13332 : LibmeshPetscCallQ(VecSetDM(*x, dm));
798 : #endif
799 :
800 404 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
801 : }
802 :
803 :
804 :
805 :
806 : #undef __FUNCT__
807 : #define __FUNCT__ "DMCreateMatrix_libMesh"
808 0 : static PetscErrorCode DMCreateMatrix_libMesh(DM dm, Mat * A)
809 : {
810 : PetscFunctionBegin;
811 0 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
812 : PetscBool eq;
813 :
814 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH, &eq));
815 :
816 0 : if (!eq)
817 0 : LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
818 :
819 0 : if (!dlm->sys)
820 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No libMesh system set for DM_libMesh");
821 :
822 0 : *A = (dynamic_cast<PetscMatrixBase<Number> *>(dlm->sys->matrix))->mat();
823 0 : LibmeshPetscCallQ(PetscObjectReference((PetscObject)(*A)));
824 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
825 : }
826 :
827 :
828 : #undef __FUNCT__
829 : #define __FUNCT__ "DMView_libMesh"
830 0 : static PetscErrorCode DMView_libMesh(DM dm, PetscViewer viewer)
831 : {
832 : PetscBool isascii;
833 : const char * name, * prefix;
834 0 : DM_libMesh * dlm = (DM_libMesh *)dm->data;
835 : PetscFunctionBegin;
836 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
837 0 : if (isascii) {
838 0 : LibmeshPetscCallQ(PetscObjectGetName((PetscObject)dm, &name));
839 0 : LibmeshPetscCallQ(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
840 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "DM libMesh with name %s and prefix %s\n", name, prefix));
841 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "blocks:"));
842 0 : std::map<std::string,unsigned int>::iterator bit = dlm->blockids->begin();
843 0 : std::map<std::string,unsigned int>::const_iterator bend = dlm->blockids->end();
844 0 : for (; bit != bend; ++bit) {
845 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "(%s,%d) ", bit->first.c_str(), bit->second));
846 : }
847 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
848 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "variables:"));
849 0 : std::map<std::string,unsigned int>::iterator vit = dlm->varids->begin();
850 0 : std::map<std::string,unsigned int>::const_iterator vend = dlm->varids->end();
851 0 : for (; vit != vend; ++vit) {
852 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "(%s,%d) ", vit->first.c_str(), vit->second));
853 : }
854 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
855 0 : if (dlm->decomposition_type == DMLIBMESH_NO_DECOMPOSITION) {
856 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "No decomposition\n"));
857 : }
858 : else {
859 0 : if (dlm->decomposition_type == DMLIBMESH_FIELD_DECOMPOSITION) {
860 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "Field decomposition by variable: "));
861 : }
862 0 : else if (dlm->decomposition_type == DMLIBMESH_DOMAIN_DECOMPOSITION) {
863 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "Domain decomposition by block: "));
864 : }
865 0 : else LIBMESH_SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Unexpected decomposition type: %d", dlm->decomposition_type);
866 : /* FIX: decompositions might have different sizes and components on different ranks. */
867 0 : for (auto d : index_range(*dlm->decomposition)) {
868 0 : std::set<unsigned int>::iterator dbegin = (*dlm->decomposition)[d].begin();
869 0 : std::set<unsigned int>::iterator dit = (*dlm->decomposition)[d].begin();
870 0 : std::set<unsigned int>::iterator dend = (*dlm->decomposition)[d].end();
871 0 : for (; dit != dend; ++dit) {
872 0 : if (dit != dbegin) {
873 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, ","));
874 : }
875 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "%u", *dit));
876 : }
877 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, ";"));
878 : }
879 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
880 : }
881 : }
882 :
883 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
884 : }
885 :
886 : #undef __FUNCT__
887 : #define __FUNCT__ "DMSetUp_libMesh"
888 29400 : static PetscErrorCode DMSetUp_libMesh(DM dm)
889 : {
890 : PetscFunctionBegin;
891 29400 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
892 : PetscBool eq;
893 :
894 29400 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH, &eq));
895 :
896 29400 : if (!eq)
897 0 : LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
898 :
899 29400 : if (!dlm->sys)
900 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No libMesh system set for DM_libMesh");
901 : /*
902 : Do not evaluate function, Jacobian or bounds for an embedded DM -- the subproblem might not have enough information for that.
903 : */
904 29400 : if (!dlm->embedding) {
905 29400 : LibmeshPetscCallQ(DMSNESSetFunction(dm, SNESFunction_DMlibMesh, (void *)dm));
906 29400 : LibmeshPetscCallQ(DMSNESSetJacobian(dm, SNESJacobian_DMlibMesh, (void *)dm));
907 30240 : if (dlm->sys->nonlinear_solver->bounds || dlm->sys->nonlinear_solver->bounds_object)
908 28280 : LibmeshPetscCallQ(DMSetVariableBounds(dm, DMVariableBounds_libMesh));
909 : }
910 : else {
911 : /*
912 : Fow now we don't implement even these, although a linear "Dirichlet" subproblem is well-defined.
913 : Creating the submatrix, however, might require extracting the submatrix preallocation from an unassembled matrix.
914 : */
915 0 : dm->ops->createglobalvector = 0;
916 0 : dm->ops->creatematrix = 0;
917 : }
918 840 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
919 : }
920 :
921 :
922 :
923 :
924 : #undef __FUNCT__
925 : #define __FUNCT__ "DMDestroy_libMesh"
926 29400 : static PetscErrorCode DMDestroy_libMesh(DM dm)
927 : {
928 29400 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
929 : PetscFunctionBegin;
930 57960 : delete dlm->varids;
931 57960 : delete dlm->varnames;
932 57960 : delete dlm->blockids;
933 57960 : delete dlm->blocknames;
934 29400 : delete dlm->decomposition;
935 29400 : LibmeshPetscCallQ(ISDestroy(&dlm->embedding));
936 29400 : LibmeshPetscCallQ(PetscFree(dm->data));
937 :
938 29400 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
939 : }
940 :
941 : EXTERN_C_BEGIN
942 : #undef __FUNCT__
943 : #define __FUNCT__ "DMCreate_libMesh"
944 29400 : PetscErrorCode DMCreate_libMesh(DM dm)
945 : {
946 : DM_libMesh * dlm;
947 :
948 : PetscFunctionBegin;
949 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
950 29400 : LibmeshPetscCallQ(PetscNew(&dlm));
951 29400 : dm->data = dlm;
952 :
953 : /* DMlibMesh impl */
954 29400 : dlm->varids = new(std::map<std::string, unsigned int>);
955 29400 : dlm->blockids = new(std::map<std::string, unsigned int>);
956 29400 : dlm->varnames = new(std::map<unsigned int, std::string>);
957 29400 : dlm->blocknames = new(std::map<unsigned int, std::string>);
958 29400 : dlm->decomposition = LIBMESH_PETSC_NULLPTR;
959 29400 : dlm->decomposition_type = DMLIBMESH_NO_DECOMPOSITION;
960 :
961 : /* DM API */
962 29400 : dm->ops->createglobalvector = DMCreateGlobalVector_libMesh;
963 29400 : dm->ops->createlocalvector = 0; // DMCreateLocalVector_libMesh;
964 29400 : dm->ops->getcoloring = 0; // DMGetColoring_libMesh;
965 29400 : dm->ops->creatematrix = DMCreateMatrix_libMesh;
966 29400 : dm->ops->createinterpolation= 0; // DMCreateInterpolation_libMesh;
967 :
968 29400 : dm->ops->refine = 0; // DMRefine_libMesh;
969 29400 : dm->ops->coarsen = 0; // DMCoarsen_libMesh;
970 :
971 : // * dm->ops->getinjection was renamed to dm->ops->createinjection in PETSc 5a84ad338 (5 Jul 2019)
972 : // * dm->ops-getaggregates was removed in PETSc 97779f9a (5 Jul 2019)
973 : // * Both changes were merged into PETSc master in 94aad3ce (7 Jul 2019).
974 : #if PETSC_VERSION_LESS_THAN(3,12,0)
975 1680 : dm->ops->getinjection = 0; // DMGetInjection_libMesh;
976 1680 : dm->ops->getaggregates = 0; // DMGetAggregates_libMesh;
977 : #else
978 27720 : dm->ops->createinjection = 0;
979 : #endif
980 :
981 :
982 29400 : dm->ops->createfielddecomposition = DMCreateFieldDecomposition_libMesh;
983 29400 : dm->ops->createdomaindecomposition = DMCreateDomainDecomposition_libMesh;
984 :
985 29400 : dm->ops->destroy = DMDestroy_libMesh;
986 29400 : dm->ops->view = DMView_libMesh;
987 29400 : dm->ops->setfromoptions = 0; // DMSetFromOptions_libMesh;
988 29400 : dm->ops->setup = DMSetUp_libMesh;
989 :
990 : /* DMlibMesh API */
991 29400 : LibmeshPetscCallQ(PetscObjectComposeFunction((PetscObject)dm,"DMlibMeshSetSystem_C",DMlibMeshSetSystem_libMesh));
992 29400 : LibmeshPetscCallQ(PetscObjectComposeFunction((PetscObject)dm,"DMlibMeshGetSystem_C",DMlibMeshGetSystem_libMesh));
993 :
994 29400 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
995 : }
996 : EXTERN_C_END
997 :
998 : #endif // LIBMESH_HAVE_PETSC
|
