Actual source code: dmmbmg.cxx
1: #include <petsc/private/dmmbimpl.h>
2: #include <petscdmmoab.h>
3: #include <MBTagConventions.hpp>
4: #include <moab/NestedRefine.hpp>
6: /*@C
7: DMMoabGenerateHierarchy - Generate a multi-level uniform refinement hierarchy
8: by succesively refining a coarse mesh, already defined in the DM object
9: provided by the user.
11: Collective
13: Input Parameter:
14: . dmb - The DMMoab object
16: Output Parameters:
17: + nlevels - The number of levels of refinement needed to generate the hierarchy
18: - ldegrees - The degree of refinement at each level in the hierarchy
20: Level: beginner
22: @*/
23: PetscErrorCode DMMoabGenerateHierarchy(DM dm, PetscInt nlevels, PetscInt *ldegrees)
24: {
25: DM_Moab *dmmoab;
26: moab::ErrorCode merr;
27: PetscInt *pdegrees, ilevel;
28: std::vector<moab::EntityHandle> hsets;
31: dmmoab = (DM_Moab *)(dm)->data;
33: if (!ldegrees) {
34: PetscMalloc1(nlevels, &pdegrees);
35: for (ilevel = 0; ilevel < nlevels; ilevel++) pdegrees[ilevel] = 2; /* default = Degree 2 refinement */
36: } else pdegrees = ldegrees;
38: /* initialize set level refinement data for hierarchy */
39: dmmoab->nhlevels = nlevels;
41: /* Instantiate the nested refinement class */
42: #ifdef MOAB_HAVE_MPI
43: dmmoab->hierarchy = new moab::NestedRefine(dynamic_cast<moab::Core *>(dmmoab->mbiface), dmmoab->pcomm, dmmoab->fileset);
44: #else
45: dmmoab->hierarchy = new moab::NestedRefine(dynamic_cast<moab::Core *>(dmmoab->mbiface), NULL, dmmoab->fileset);
46: #endif
48: PetscMalloc1(nlevels + 1, &dmmoab->hsets);
50: /* generate the mesh hierarchy */
51: merr = dmmoab->hierarchy->generate_mesh_hierarchy(nlevels, pdegrees, hsets, false);
52: MBERRNM(merr);
54: #ifdef MOAB_HAVE_MPI
55: if (dmmoab->pcomm->size() > 1) {
56: merr = dmmoab->hierarchy->exchange_ghosts(hsets, dmmoab->nghostrings);
57: MBERRNM(merr);
58: }
59: #endif
61: /* copy the mesh sets for nested refinement hierarchy */
62: dmmoab->hsets[0] = hsets[0];
63: for (ilevel = 1; ilevel <= nlevels; ilevel++) {
64: dmmoab->hsets[ilevel] = hsets[ilevel];
66: #ifdef MOAB_HAVE_MPI
67: merr = dmmoab->pcomm->assign_global_ids(hsets[ilevel], dmmoab->dim, 0, false, true, false);
68: MBERRNM(merr);
69: #endif
71: /* Update material and other geometric tags from parent to child sets */
72: merr = dmmoab->hierarchy->update_special_tags(ilevel, hsets[ilevel]);
73: MBERRNM(merr);
74: }
76: hsets.clear();
77: if (!ldegrees) PetscFree(pdegrees);
78: return 0;
79: }
81: /*@C
82: DMRefineHierarchy_Moab - Generate a multi-level DM hierarchy
83: by succesively refining a coarse mesh.
85: Collective
87: Input Parameter:
88: . dm - The DMMoab object
90: Output Parameters:
91: + nlevels - The number of levels of refinement needed to generate the hierarchy
92: - dmf - The DM objects after successive refinement of the hierarchy
94: Level: beginner
96: @*/
97: PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmf[])
98: {
99: PetscInt i;
102: DMRefine(dm, PetscObjectComm((PetscObject)dm), &dmf[0]);
103: for (i = 1; i < nlevels; i++) DMRefine(dmf[i - 1], PetscObjectComm((PetscObject)dm), &dmf[i]);
104: return 0;
105: }
107: /*@C
108: DMCoarsenHierarchy_Moab - Generate a multi-level DM hierarchy
109: by succesively coarsening a refined mesh.
111: Collective
113: Input Parameter:
114: . dm - The DMMoab object
116: Output Parameters:
117: + nlevels - The number of levels of refinement needed to generate the hierarchy
118: - dmc - The DM objects after successive coarsening of the hierarchy
120: Level: beginner
122: @*/
123: PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmc[])
124: {
125: PetscInt i;
128: DMCoarsen(dm, PetscObjectComm((PetscObject)dm), &dmc[0]);
129: for (i = 1; i < nlevels; i++) DMCoarsen(dmc[i - 1], PetscObjectComm((PetscObject)dm), &dmc[i]);
130: return 0;
131: }
133: PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM, PetscInt *, PetscInt *, PetscInt *, PetscInt *, PetscBool);
135: /*@C
136: DMCreateInterpolation_Moab - Generate the interpolation operators to transform
137: operators (matrices, vectors) from parent level to child level as defined by
138: the DM inputs provided by the user.
140: Collective
142: Input Parameters:
143: + dm1 - The DMMoab object
144: - dm2 - the second, finer DMMoab object
146: Output Parameters:
147: + interpl - The interpolation operator for transferring data between the levels
148: - vec - The scaling vector (optional)
150: Level: developer
152: @*/
153: PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmp, DM dmc, Mat *interpl, Vec *vec)
154: {
155: DM_Moab *dmbp, *dmbc;
156: moab::ErrorCode merr;
157: PetscInt dim;
158: PetscReal factor;
159: PetscInt innz, *nnz, ionz, *onz;
160: PetscInt nlsizp, nlsizc, nlghsizp, ngsizp, ngsizc;
161: const PetscBool use_consistent_bases = PETSC_TRUE;
165: dmbp = (DM_Moab *)(dmp)->data;
166: dmbc = (DM_Moab *)(dmc)->data;
167: nlsizp = dmbp->nloc; // *dmb1->numFields;
168: nlsizc = dmbc->nloc; // *dmb2->numFields;
169: ngsizp = dmbp->n; // *dmb1->numFields;
170: ngsizc = dmbc->n; // *dmb2->numFields;
171: nlghsizp = (dmbp->nloc + dmbp->nghost); // *dmb1->numFields;
173: // Columns = Parent DoFs ; Rows = Child DoFs
174: // Interpolation matrix: \sum_{i=1}^P Owned(Child) * (Owned(Parent) + Ghosted(Parent))
175: // Size: nlsizc * nlghsizp
176: PetscInfo(NULL, "Creating interpolation matrix %" PetscInt_FMT " X %" PetscInt_FMT " to apply transformation between levels %" PetscInt_FMT " -> %" PetscInt_FMT ".\n", ngsizc, nlghsizp, dmbp->hlevel, dmbc->hlevel);
178: DMGetDimension(dmp, &dim);
180: /* allocate the nnz, onz arrays based on block size and local nodes */
181: PetscCalloc2(nlsizc, &nnz, nlsizc, &onz);
183: /* Loop through the local elements and compute the relation between the current parent and the refined_level. */
184: for (moab::Range::iterator iter = dmbc->vowned->begin(); iter != dmbc->vowned->end(); iter++) {
185: const moab::EntityHandle vhandle = *iter;
186: /* define local variables */
187: moab::EntityHandle parent;
188: std::vector<moab::EntityHandle> adjs;
189: moab::Range found;
191: /* store the vertex DoF number */
192: const int ldof = dmbc->lidmap[vhandle - dmbc->seqstart];
194: /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects
195: to the current vertex. We can then decipher if a vertex is ghosted or not and compute the
196: non-zero pattern accordingly. */
197: merr = dmbc->hierarchy->get_adjacencies(vhandle, dmbc->dim, adjs);
198: MBERRNM(merr);
200: /* loop over vertices and update the number of connectivity */
201: for (unsigned jter = 0; jter < adjs.size(); jter++) {
202: const moab::EntityHandle jhandle = adjs[jter];
204: /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */
205: merr = dmbc->hierarchy->child_to_parent(jhandle, dmbc->hlevel, dmbp->hlevel, &parent);
206: MBERRNM(merr);
208: /* Get connectivity information in canonical ordering for the local element */
209: std::vector<moab::EntityHandle> connp;
210: merr = dmbp->hierarchy->get_connectivity(parent, dmbp->hlevel, connp);
211: MBERRNM(merr);
213: for (unsigned ic = 0; ic < connp.size(); ++ic) {
214: /* loop over each element connected to the adjacent vertex and update as needed */
215: /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */
216: if (found.find(connp[ic]) != found.end()) continue; /* make sure we don't double count shared vertices */
217: if (dmbp->vghost->find(connp[ic]) != dmbp->vghost->end()) onz[ldof]++; /* update out-of-proc onz */
218: else nnz[ldof]++; /* else local vertex */
219: found.insert(connp[ic]);
220: }
221: }
222: }
224: for (int i = 0; i < nlsizc; i++) nnz[i] += 1; /* self count the node */
226: ionz = onz[0];
227: innz = nnz[0];
228: for (int tc = 0; tc < nlsizc; tc++) {
229: // check for maximum allowed sparsity = fully dense
230: nnz[tc] = std::min(nlsizp, nnz[tc]);
231: onz[tc] = std::min(ngsizp - nlsizp, onz[tc]);
233: PetscInfo(NULL, " %d: NNZ = %d, ONZ = %d\n", tc, nnz[tc], onz[tc]);
235: innz = (innz < nnz[tc] ? nnz[tc] : innz);
236: ionz = (ionz < onz[tc] ? onz[tc] : ionz);
237: }
239: /* create interpolation matrix */
240: MatCreate(PetscObjectComm((PetscObject)dmc), interpl);
241: MatSetSizes(*interpl, nlsizc, nlsizp, ngsizc, ngsizp);
242: MatSetType(*interpl, MATAIJ);
243: MatSetFromOptions(*interpl);
245: MatSeqAIJSetPreallocation(*interpl, innz, nnz);
246: MatMPIAIJSetPreallocation(*interpl, innz, nnz, ionz, onz);
248: /* clean up temporary memory */
249: PetscFree2(nnz, onz);
251: /* set up internal matrix data-structures */
252: MatSetUp(*interpl);
254: /* Define variables for assembly */
255: std::vector<moab::EntityHandle> children;
256: std::vector<moab::EntityHandle> connp, connc;
257: std::vector<PetscReal> pcoords, ccoords, values_phi;
259: if (use_consistent_bases) {
260: const moab::EntityHandle ehandle = dmbp->elocal->front();
262: merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children);
263: MBERRNM(merr);
265: /* Get connectivity and coordinates of the parent vertices */
266: merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp);
267: MBERRNM(merr);
268: merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc);
269: MBERRNM(merr);
271: std::vector<PetscReal> natparam(3 * connc.size(), 0.0);
272: pcoords.resize(connp.size() * 3);
273: ccoords.resize(connc.size() * 3);
274: values_phi.resize(connp.size() * connc.size());
275: /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */
276: merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]);
277: MBERRNM(merr);
278: merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]);
279: MBERRNM(merr);
281: /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */
282: for (unsigned tc = 0; tc < connc.size(); tc++) {
283: const PetscInt offset = tc * 3;
285: /* Scale ccoords relative to pcoords */
286: DMMoabPToRMapping(dim, connp.size(), &pcoords[0], &ccoords[offset], &natparam[offset], &values_phi[connp.size() * tc]);
287: }
288: } else {
289: factor = std::pow(2.0 /*degree_P_for_refinement*/, (dmbc->hlevel - dmbp->hlevel) * dmbp->dim * 1.0);
290: }
292: /* TODO: Decipher the correct non-zero pattern. There is still some issue with onz allocation */
293: MatSetOption(*interpl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
295: /* Loop through the remaining vertices. These vertices appear only on the current refined_level. */
296: for (moab::Range::iterator iter = dmbp->elocal->begin(); iter != dmbp->elocal->end(); iter++) {
297: const moab::EntityHandle ehandle = *iter;
299: /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */
300: children.clear();
301: connc.clear();
302: merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children);
303: MBERRNM(merr);
305: /* Get connectivity and coordinates of the parent vertices */
306: merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp);
307: MBERRNM(merr);
308: merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc);
309: MBERRNM(merr);
311: pcoords.resize(connp.size() * 3);
312: ccoords.resize(connc.size() * 3);
313: /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */
314: merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]);
315: MBERRNM(merr);
316: merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]);
317: MBERRNM(merr);
319: std::vector<int> dofsp(connp.size()), dofsc(connc.size());
320: /* TODO: specific to scalar system - use GetDofs */
321: DMMoabGetDofsBlocked(dmp, connp.size(), &connp[0], &dofsp[0]);
322: DMMoabGetDofsBlocked(dmc, connc.size(), &connc[0], &dofsc[0]);
324: /* Compute the actual interpolation weights when projecting solution/residual between levels */
325: if (use_consistent_bases) {
326: /* Use the cached values of natural parameteric coordinates and basis pre-evaluated.
327: We are making an assumption here that UMR used in GMG to generate the hierarchy uses
328: the same template for all elements; This will fail for mixed element meshes (TRI/QUAD).
330: TODO: Fix the above assumption by caching data for families (especially for Tets and mixed meshes)
331: */
333: /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */
334: for (unsigned tc = 0; tc < connc.size(); tc++) {
335: /* TODO: Check if we should be using INSERT_VALUES instead */
336: MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi[connp.size() * tc], ADD_VALUES);
337: }
338: } else {
339: /* Compute the interpolation weights by determining distance of 1-ring
340: neighbor vertices from current vertex
342: This should be used only when FEM basis is not used for the discretization.
343: Else, the consistent interface to compute the basis function for interpolation
344: between the levels should be evaluated correctly to preserve convergence of GMG.
345: Shephard's basis will be terrible for any unsmooth problems.
346: */
347: values_phi.resize(connp.size());
348: for (unsigned tc = 0; tc < connc.size(); tc++) {
349: PetscReal normsum = 0.0;
350: for (unsigned tp = 0; tp < connp.size(); tp++) {
351: values_phi[tp] = 0.0;
352: for (unsigned k = 0; k < 3; k++) values_phi[tp] += std::pow(pcoords[tp * 3 + k] - ccoords[k + tc * 3], dim);
353: if (values_phi[tp] < 1e-12) {
354: values_phi[tp] = 1e12;
355: } else {
356: //values_phi[tp] = std::pow(values_phi[tp], -1.0/dim);
357: values_phi[tp] = std::pow(values_phi[tp], -1.0);
358: normsum += values_phi[tp];
359: }
360: }
361: for (unsigned tp = 0; tp < connp.size(); tp++) {
362: if (values_phi[tp] > 1e11) values_phi[tp] = factor * 0.5 / connp.size();
363: else values_phi[tp] = factor * values_phi[tp] * 0.5 / (connp.size() * normsum);
364: }
365: MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi[0], ADD_VALUES);
366: }
367: }
368: }
369: if (vec) *vec = NULL;
370: MatAssemblyBegin(*interpl, MAT_FINAL_ASSEMBLY);
371: MatAssemblyEnd(*interpl, MAT_FINAL_ASSEMBLY);
372: return 0;
373: }
375: /*@C
376: DMCreateInjection_Moab - Generate a multi-level uniform refinement hierarchy
377: by succesively refining a coarse mesh, already defined in the DM object
378: provided by the user.
380: Collective
382: Input Parameter:
383: . dmb - The DMMoab object
385: Output Parameters:
386: + nlevels - The number of levels of refinement needed to generate the hierarchy
387: - ldegrees - The degree of refinement at each level in the hierarchy
389: Level: beginner
391: @*/
392: PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dm1, DM dm2, VecScatter *ctx)
393: {
394: //DM_Moab *dmmoab;
398: //dmmoab = (DM_Moab*)(dm1)->data;
400: PetscPrintf(PETSC_COMM_WORLD, "[DMCreateInjection_Moab] :: Placeholder\n");
401: return 0;
402: }
404: static PetscErrorCode DMMoab_UMR_Private(DM dm, MPI_Comm comm, PetscBool refine, DM *dmref)
405: {
406: PetscInt i, dim;
407: DM dm2;
408: moab::ErrorCode merr;
409: DM_Moab *dmb = (DM_Moab *)dm->data, *dd2;
414: if ((dmb->hlevel == dmb->nhlevels && refine) || (dmb->hlevel == 0 && !refine)) {
415: if (dmb->hlevel + 1 > dmb->nhlevels && refine) PetscInfo(NULL, "Invalid multigrid refinement hierarchy level specified (%" PetscInt_FMT "). MOAB UMR max levels = %" PetscInt_FMT ". Creating a NULL object.\n", dmb->hlevel + 1, dmb->nhlevels);
416: if (dmb->hlevel - 1 < 0 && !refine) PetscInfo(NULL, "Invalid multigrid coarsen hierarchy level specified (%" PetscInt_FMT "). Creating a NULL object.\n", dmb->hlevel - 1);
417: *dmref = PETSC_NULL;
418: return 0;
419: }
421: DMMoabCreate(PetscObjectComm((PetscObject)dm), &dm2);
422: dd2 = (DM_Moab *)dm2->data;
424: dd2->mbiface = dmb->mbiface;
425: #ifdef MOAB_HAVE_MPI
426: dd2->pcomm = dmb->pcomm;
427: #endif
428: dd2->icreatedinstance = PETSC_FALSE;
429: dd2->nghostrings = dmb->nghostrings;
431: /* set the new level based on refinement/coarsening */
432: if (refine) {
433: dd2->hlevel = dmb->hlevel + 1;
434: } else {
435: dd2->hlevel = dmb->hlevel - 1;
436: }
438: /* Copy the multilevel hierarchy pointers in MOAB */
439: dd2->hierarchy = dmb->hierarchy;
440: dd2->nhlevels = dmb->nhlevels;
441: PetscMalloc1(dd2->nhlevels + 1, &dd2->hsets);
442: for (i = 0; i <= dd2->nhlevels; i++) dd2->hsets[i] = dmb->hsets[i];
443: dd2->fileset = dd2->hsets[dd2->hlevel];
445: /* do the remaining initializations for DMMoab */
446: dd2->bs = dmb->bs;
447: dd2->numFields = dmb->numFields;
448: dd2->rw_dbglevel = dmb->rw_dbglevel;
449: dd2->partition_by_rank = dmb->partition_by_rank;
450: PetscStrcpy(dd2->extra_read_options, dmb->extra_read_options);
451: PetscStrcpy(dd2->extra_write_options, dmb->extra_write_options);
452: dd2->read_mode = dmb->read_mode;
453: dd2->write_mode = dmb->write_mode;
455: /* set global ID tag handle */
456: DMMoabSetLocalToGlobalTag(dm2, dmb->ltog_tag);
458: merr = dd2->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dd2->material_tag);
459: MBERRNM(merr);
461: DMSetOptionsPrefix(dm2, ((PetscObject)dm)->prefix);
462: DMGetDimension(dm, &dim);
463: DMSetDimension(dm2, dim);
465: /* allow overloaded (user replaced) operations to be inherited by refinement clones */
466: dm2->ops->creatematrix = dm->ops->creatematrix;
468: /* copy fill information if given */
469: DMMoabSetBlockFills(dm2, dmb->dfill, dmb->ofill);
471: /* copy vector type information */
472: DMSetMatType(dm2, dm->mattype);
473: DMSetVecType(dm2, dm->vectype);
474: dd2->numFields = dmb->numFields;
475: if (dmb->numFields) DMMoabSetFieldNames(dm2, dmb->numFields, dmb->fieldNames);
477: DMSetFromOptions(dm2);
479: /* recreate Dof numbering for the refined DM and make sure the distribution is correctly populated */
480: DMSetUp(dm2);
482: *dmref = dm2;
483: return 0;
484: }
486: /*@C
487: DMRefine_Moab - Generate a multi-level uniform refinement hierarchy
488: by succesively refining a coarse mesh, already defined in the DM object
489: provided by the user.
491: Collective on dm
493: Input Parameters:
494: + dm - The DMMoab object
495: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
497: Output Parameter:
498: . dmf - the refined DM, or NULL
500: Note: If no refinement was done, the return value is NULL
502: Level: developer
504: @*/
505: PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM *dmf)
506: {
509: DMMoab_UMR_Private(dm, comm, PETSC_TRUE, dmf);
510: return 0;
511: }
513: /*@C
514: DMCoarsen_Moab - Generate a multi-level uniform refinement hierarchy
515: by succesively refining a coarse mesh, already defined in the DM object
516: provided by the user.
518: Collective on dm
520: Input Parameters:
521: + dm - The DMMoab object
522: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
524: Output Parameter:
525: . dmf - the coarsened DM, or NULL
527: Note: If no coarsening was done, the return value is NULL
529: Level: developer
531: @*/
532: PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM *dmc)
533: {
535: DMMoab_UMR_Private(dm, comm, PETSC_FALSE, dmc);
536: return 0;
537: }