Actual source code: dmmbmg.cxx
petsc-3.9.4 2018-09-11
1: #include <petsc/private/dmmbimpl.h>
2: #include <petscdmmoab.h>
3: #include <MBTagConventions.hpp>
4: #include <moab/NestedRefine.hpp>
6: /*@
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 on MPI_Comm
13: Input Parameter:
14: + dmb - The DMMoab object
16: Output Parameter:
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: .keywords: DMMoab, create, refinement
23: @*/
24: PetscErrorCode DMMoabGenerateHierarchy(DM dm, PetscInt nlevels, PetscInt *ldegrees)
25: {
26: DM_Moab *dmmoab;
27: PetscErrorCode ierr;
28: moab::ErrorCode merr;
29: PetscInt *pdegrees, ilevel;
30: std::vector<moab::EntityHandle> hsets;
34: dmmoab = (DM_Moab*)(dm)->data;
36: if (!ldegrees) {
37: PetscMalloc1(nlevels, &pdegrees);
38: for (ilevel = 0; ilevel < nlevels; ilevel++) pdegrees[ilevel] = 2; /* default = Degree 2 refinement */
39: }
40: else pdegrees = ldegrees;
42: /* initialize set level refinement data for hierarchy */
43: dmmoab->nhlevels = nlevels;
45: /* Instantiate the nested refinement class */
46: #ifdef MOAB_HAVE_MPI
47: dmmoab->hierarchy = new moab::NestedRefine(dynamic_cast<moab::Core*>(dmmoab->mbiface), dmmoab->pcomm, dmmoab->fileset);
48: #else
49: dmmoab->hierarchy = new moab::NestedRefine(dynamic_cast<moab::Core*>(dmmoab->mbiface), NULL, dmmoab->fileset);
50: #endif
52: PetscMalloc1(nlevels + 1, &dmmoab->hsets);
54: /* generate the mesh hierarchy */
55: merr = dmmoab->hierarchy->generate_mesh_hierarchy(nlevels, pdegrees, hsets, false);MBERRNM(merr);
57: #ifdef MOAB_HAVE_MPI
58: if (dmmoab->pcomm->size() > 1) {
59: merr = dmmoab->hierarchy->exchange_ghosts(hsets, dmmoab->nghostrings);MBERRNM(merr);
60: }
61: #endif
63: /* copy the mesh sets for nested refinement hierarchy */
64: dmmoab->hsets[0] = hsets[0];
65: for (ilevel = 1; ilevel <= nlevels; ilevel++)
66: {
67: dmmoab->hsets[ilevel] = hsets[ilevel];
69: #ifdef MOAB_HAVE_MPI
70: merr = dmmoab->pcomm->assign_global_ids(hsets[ilevel], dmmoab->dim, 0, false, true, false);MBERRNM(merr);
71: #endif
73: /* Update material and other geometric tags from parent to child sets */
74: merr = dmmoab->hierarchy->update_special_tags(ilevel, hsets[ilevel]);MBERRNM(merr);
75: }
77: hsets.clear();
78: if (!ldegrees) {
79: PetscFree(pdegrees);
80: }
81: return(0);
82: }
84: /*@
85: DMRefineHierarchy_Moab - Generate a multi-level DM hierarchy
86: by succesively refining a coarse mesh.
88: Collective on MPI_Comm
90: Input Parameter:
91: + dm - The DMMoab object
93: Output Parameter:
94: + nlevels - The number of levels of refinement needed to generate the hierarchy
95: . dmf - The DM objects after successive refinement of the hierarchy
97: Level: beginner
99: .keywords: DMMoab, generate, refinement
100: @*/
101: PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmf[])
102: {
103: PetscErrorCode ierr;
104: PetscInt i;
108: DMRefine(dm, PetscObjectComm((PetscObject)dm), &dmf[0]);
109: for (i = 1; i < nlevels; i++) {
110: DMRefine(dmf[i - 1], PetscObjectComm((PetscObject)dm), &dmf[i]);
111: }
112: return(0);
113: }
115: /*@
116: DMCoarsenHierarchy_Moab - Generate a multi-level DM hierarchy
117: by succesively coarsening a refined mesh.
119: Collective on MPI_Comm
121: Input Parameter:
122: + dm - The DMMoab object
124: Output Parameter:
125: + nlevels - The number of levels of refinement needed to generate the hierarchy
126: . dmc - The DM objects after successive coarsening of the hierarchy
128: Level: beginner
130: .keywords: DMMoab, generate, coarsening
131: @*/
132: PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmc[])
133: {
134: PetscErrorCode ierr;
135: PetscInt i;
139: DMCoarsen(dm, PetscObjectComm((PetscObject)dm), &dmc[0]);
140: for (i = 1; i < nlevels; i++) {
141: DMCoarsen(dmc[i - 1], PetscObjectComm((PetscObject)dm), &dmc[i]);
142: }
143: return(0);
144: }
146: PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM, PetscInt*, PetscInt*, PetscInt*, PetscInt*, PetscBool);
148: /*@
149: DMCreateInterpolation_Moab - Generate the interpolation operators to transform
150: operators (matrices, vectors) from parent level to child level as defined by
151: the DM inputs provided by the user.
153: Collective on MPI_Comm
155: Input Parameter:
156: + dm1 - The DMMoab object
157: - dm2 - the second, finer DMMoab object
159: Output Parameter:
160: + interpl - The interpolation operator for transferring data between the levels
161: - vec - The scaling vector (optional)
163: Level: developer
165: .keywords: DMMoab, create, refinement
166: @*/
167: PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmp, DM dmc, Mat* interpl, Vec* vec)
168: {
169: DM_Moab *dmbp, *dmbc;
170: PetscErrorCode ierr;
171: moab::ErrorCode merr;
172: PetscInt dim;
173: PetscReal factor;
174: PetscInt innz, *nnz, ionz, *onz;
175: PetscInt nlsizp, nlsizc, nlghsizp, ngsizp, ngsizc;
176: const PetscBool use_consistent_bases=PETSC_TRUE;
181: dmbp = (DM_Moab*)(dmp)->data;
182: dmbc = (DM_Moab*)(dmc)->data;
183: nlsizp = dmbp->nloc;// *dmb1->numFields;
184: nlsizc = dmbc->nloc;// *dmb2->numFields;
185: ngsizp = dmbp->n; // *dmb1->numFields;
186: ngsizc = dmbc->n; // *dmb2->numFields;
187: nlghsizp = (dmbp->nloc + dmbp->nghost); // *dmb1->numFields;
189: // Columns = Parent DoFs ; Rows = Child DoFs
190: // Interpolation matrix: \sum_{i=1}^P Owned(Child) * (Owned(Parent) + Ghosted(Parent))
191: // Size: nlsizc * nlghsizp
192: PetscInfo4(NULL, "Creating interpolation matrix %D X %D to apply transformation between levels %D -> %D.\n", ngsizc, nlghsizp, dmbp->hlevel, dmbc->hlevel);
194: DMGetDimension(dmp, &dim);
196: /* allocate the nnz, onz arrays based on block size and local nodes */
197: PetscCalloc2(nlsizc, &nnz, nlsizc, &onz);
199: /* Loop through the local elements and compute the relation between the current parent and the refined_level. */
200: for (moab::Range::iterator iter = dmbc->vowned->begin(); iter != dmbc->vowned->end(); iter++) {
202: const moab::EntityHandle vhandle = *iter;
203: /* define local variables */
204: moab::EntityHandle parent;
205: std::vector<moab::EntityHandle> adjs;
206: moab::Range found;
208: /* store the vertex DoF number */
209: const int ldof = dmbc->lidmap[vhandle - dmbc->seqstart];
211: /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects
212: to the current vertex. We can then decipher if a vertex is ghosted or not and compute the
213: non-zero pattern accordingly. */
214: merr = dmbc->hierarchy->get_adjacencies(vhandle, dmbc->dim, adjs);MBERRNM(merr);
216: /* loop over vertices and update the number of connectivity */
217: for (unsigned jter = 0; jter < adjs.size(); jter++) {
219: const moab::EntityHandle jhandle = adjs[jter];
221: /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */
222: merr = dmbc->hierarchy->child_to_parent(jhandle, dmbc->hlevel, dmbp->hlevel, &parent);MBERRNM(merr);
224: /* Get connectivity information in canonical ordering for the local element */
225: std::vector<moab::EntityHandle> connp;
226: merr = dmbp->hierarchy->get_connectivity(parent, dmbp->hlevel, connp);MBERRNM(merr);
228: for (unsigned ic=0; ic < connp.size(); ++ic) {
230: /* loop over each element connected to the adjacent vertex and update as needed */
231: /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */
232: if (found.find(connp[ic]) != found.end()) continue; /* make sure we don't double count shared vertices */
233: if (dmbp->vghost->find(connp[ic]) != dmbp->vghost->end()) onz[ldof]++; /* update out-of-proc onz */
234: else nnz[ldof]++; /* else local vertex */
235: found.insert(connp[ic]);
236: }
237: }
238: }
240: for (int i = 0; i < nlsizc; i++)
241: nnz[i] += 1; /* self count the node */
243: ionz = onz[0];
244: innz = nnz[0];
245: for (int tc = 0; tc < nlsizc; tc++) {
246: // check for maximum allowed sparsity = fully dense
247: nnz[tc] = std::min(nlsizp, nnz[tc]);
248: onz[tc] = std::min(ngsizp - nlsizp, onz[tc]);
250: PetscInfo3(NULL, " %d: NNZ = %d, ONZ = %d\n", tc, nnz[tc], onz[tc]);
252: innz = (innz < nnz[tc] ? nnz[tc] : innz);
253: ionz = (ionz < onz[tc] ? onz[tc] : ionz);
254: }
256: /* create interpolation matrix */
257: MatCreate(PetscObjectComm((PetscObject)dmc), interpl);
258: MatSetSizes(*interpl, nlsizc, nlsizp, ngsizc, ngsizp);
259: MatSetType(*interpl, MATAIJ);
260: MatSetFromOptions(*interpl);
262: MatSeqAIJSetPreallocation(*interpl, innz, nnz);
263: MatMPIAIJSetPreallocation(*interpl, innz, nnz, ionz, onz);
265: /* clean up temporary memory */
266: PetscFree2(nnz, onz);
268: /* set up internal matrix data-structures */
269: MatSetUp(*interpl);
271: /* Define variables for assembly */
272: std::vector<moab::EntityHandle> children;
273: std::vector<moab::EntityHandle> connp, connc;
274: std::vector<PetscReal> pcoords, ccoords, values_phi;
276: if (use_consistent_bases) {
277: const moab::EntityHandle ehandle = dmbp->elocal->front();
279: merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children);MBERRNM(merr);
281: /* Get connectivity and coordinates of the parent vertices */
282: merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp);MBERRNM(merr);
283: merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc);MBERRNM(merr);
285: std::vector<PetscReal> natparam(3*connc.size(), 0.0);
286: pcoords.resize(connp.size() * 3);
287: ccoords.resize(connc.size() * 3);
288: values_phi.resize(connp.size()*connc.size());
289: /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */
290: merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]);MBERRNM(merr);
291: merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]);MBERRNM(merr);
293: /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */
294: for (unsigned tc = 0; tc < connc.size(); tc++) {
295: const PetscInt offset = tc * 3;
297: /* Scale ccoords relative to pcoords */
298: DMMoabPToRMapping(dim, connp.size(), &pcoords[0], &ccoords[offset], &natparam[offset], &values_phi[connp.size()*tc]);
299: }
300: }
301: else {
302: factor = std::pow(2.0 /*degree_P_for_refinement*/, (dmbc->hlevel - dmbp->hlevel) * dmbp->dim * 1.0);
303: }
305: /* TODO: Decipher the correct non-zero pattern. There is still some issue with onz allocation */
306: MatSetOption(*interpl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
308: /* Loop through the remaining vertices. These vertices appear only on the current refined_level. */
309: for (moab::Range::iterator iter = dmbp->elocal->begin(); iter != dmbp->elocal->end(); iter++) {
311: const moab::EntityHandle ehandle = *iter;
313: /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */
314: children.clear();
315: connc.clear();
316: merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children);MBERRNM(merr);
318: /* Get connectivity and coordinates of the parent vertices */
319: merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp);MBERRNM(merr);
320: merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc);MBERRNM(merr);
322: pcoords.resize(connp.size() * 3);
323: ccoords.resize(connc.size() * 3);
324: /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */
325: merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]);MBERRNM(merr);
326: merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]);MBERRNM(merr);
328: std::vector<int> dofsp(connp.size()), dofsc(connc.size());
329: /* TODO: specific to scalar system - use GetDofs */
330: DMMoabGetDofsBlocked(dmp, connp.size(), &connp[0], &dofsp[0]);
331: DMMoabGetDofsBlocked(dmc, connc.size(), &connc[0], &dofsc[0]);
333: /* Compute the actual interpolation weights when projecting solution/residual between levels */
334: if (use_consistent_bases) {
336: /* Use the cached values of natural parameteric coordinates and basis pre-evaluated.
337: We are making an assumption here that UMR used in GMG to generate the hierarchy uses
338: the same template for all elements; This will fail for mixed element meshes (TRI/QUAD).
340: TODO: Fix the above assumption by caching data for families (especially for Tets and mixed meshes)
341: */
343: /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */
344: for (unsigned tc = 0; tc < connc.size(); tc++) {
345: /* TODO: Check if we should be using INSERT_VALUES instead */
346: MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi[connp.size()*tc], ADD_VALUES);
347: }
348: }
349: else {
350: /* Compute the interpolation weights by determining distance of 1-ring
351: neighbor vertices from current vertex
353: This should be used only when FEM basis is not used for the discretization.
354: Else, the consistent interface to compute the basis function for interpolation
355: between the levels should be evaluated correctly to preserve convergence of GMG.
356: Shephard's basis will be terrible for any unsmooth problems.
357: */
358: values_phi.resize(connp.size());
359: for (unsigned tc = 0; tc < connc.size(); tc++) {
361: PetscReal normsum = 0.0;
362: for (unsigned tp = 0; tp < connp.size(); tp++) {
363: values_phi[tp] = 0.0;
364: for (unsigned k = 0; k < 3; k++)
365: values_phi[tp] += std::pow(pcoords[tp * 3 + k] - ccoords[k + tc * 3], dim);
366: if (values_phi[tp] < 1e-12) {
367: values_phi[tp] = 1e12;
368: }
369: else {
370: //values_phi[tp] = std::pow(values_phi[tp], -1.0/dim);
371: values_phi[tp] = std::pow(values_phi[tp], -1.0);
372: normsum += values_phi[tp];
373: }
374: }
375: for (unsigned tp = 0; tp < connp.size(); tp++) {
376: if (values_phi[tp] > 1e11)
377: values_phi[tp] = factor * 0.5 / connp.size();
378: else
379: values_phi[tp] = factor * values_phi[tp] * 0.5 / (connp.size() * normsum);
380: }
381: MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi[0], ADD_VALUES);
382: }
383: }
384: }
385: if (vec) *vec = NULL;
386: MatAssemblyBegin(*interpl, MAT_FINAL_ASSEMBLY);
387: MatAssemblyEnd(*interpl, MAT_FINAL_ASSEMBLY);
388: return(0);
389: }
391: /*@
392: DMCreateInjection_Moab - Generate a multi-level uniform refinement hierarchy
393: by succesively refining a coarse mesh, already defined in the DM object
394: provided by the user.
396: Collective on MPI_Comm
398: Input Parameter:
399: . dmb - The DMMoab object
401: Output Parameter:
402: . nlevels - The number of levels of refinement needed to generate the hierarchy
403: + ldegrees - The degree of refinement at each level in the hierarchy
405: Level: beginner
407: .keywords: DMMoab, create, refinement
408: @*/
409: PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dm1, DM dm2, VecScatter* ctx)
410: {
411: //DM_Moab *dmmoab;
416: //dmmoab = (DM_Moab*)(dm1)->data;
418: PetscPrintf(PETSC_COMM_WORLD, "[DMCreateInjection_Moab] :: Placeholder\n");
419: return(0);
420: }
422: static PetscErrorCode DMMoab_UMR_Private(DM dm, MPI_Comm comm, PetscBool refine, DM *dmref)
423: {
424: PetscErrorCode ierr;
425: PetscInt i, dim;
426: DM dm2;
427: moab::ErrorCode merr;
428: DM_Moab *dmb = (DM_Moab*)dm->data, *dd2;
434: if ( (dmb->hlevel == dmb->nhlevels && refine) || (dmb->hlevel == 0 && !refine) ) {
435: if (dmb->hlevel + 1 > dmb->nhlevels && refine) PetscInfo2(NULL, "Invalid multigrid refinement hierarchy level specified (%D). MOAB UMR max levels = %D. Creating a NULL object.\n", dmb->hlevel + 1, dmb->nhlevels);
436: if (dmb->hlevel - 1 < 0 && !refine) PetscInfo1(NULL, "Invalid multigrid coarsen hierarchy level specified (%D). Creating a NULL object.\n", dmb->hlevel - 1);
437: *dmref = PETSC_NULL;
438: return(0);
439: }
441: DMMoabCreate(PetscObjectComm((PetscObject)dm), &dm2);
442: dd2 = (DM_Moab*)dm2->data;
444: dd2->mbiface = dmb->mbiface;
445: #ifdef MOAB_HAVE_MPI
446: dd2->pcomm = dmb->pcomm;
447: #endif
448: dd2->icreatedinstance = PETSC_FALSE;
449: dd2->nghostrings = dmb->nghostrings;
451: /* set the new level based on refinement/coarsening */
452: if (refine) {
453: dd2->hlevel = dmb->hlevel + 1;
454: }
455: else {
456: dd2->hlevel = dmb->hlevel - 1;
457: }
459: /* Copy the multilevel hierarchy pointers in MOAB */
460: dd2->hierarchy = dmb->hierarchy;
461: dd2->nhlevels = dmb->nhlevels;
462: PetscMalloc1(dd2->nhlevels + 1, &dd2->hsets);
463: for (i = 0; i <= dd2->nhlevels; i++) {
464: dd2->hsets[i] = dmb->hsets[i];
465: }
466: dd2->fileset = dd2->hsets[dd2->hlevel];
468: /* do the remaining initializations for DMMoab */
469: dd2->bs = dmb->bs;
470: dd2->numFields = dmb->numFields;
471: dd2->rw_dbglevel = dmb->rw_dbglevel;
472: dd2->partition_by_rank = dmb->partition_by_rank;
473: PetscStrcpy(dd2->extra_read_options, dmb->extra_read_options);
474: PetscStrcpy(dd2->extra_write_options, dmb->extra_write_options);
475: dd2->read_mode = dmb->read_mode;
476: dd2->write_mode = dmb->write_mode;
478: /* set global ID tag handle */
479: DMMoabSetLocalToGlobalTag(dm2, dmb->ltog_tag);
481: merr = dd2->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dd2->material_tag);MBERRNM(merr);
483: DMSetOptionsPrefix(dm2, ((PetscObject)dm)->prefix);
484: DMGetDimension(dm, &dim);
485: DMSetDimension(dm2, dim);
487: /* allow overloaded (user replaced) operations to be inherited by refinement clones */
488: dm2->ops->creatematrix = dm->ops->creatematrix;
490: /* copy fill information if given */
491: DMMoabSetBlockFills(dm2, dmb->dfill, dmb->ofill);
493: /* copy vector type information */
494: DMSetMatType(dm2, dm->mattype);
495: DMSetVecType(dm2, dm->vectype);
496: dd2->numFields = dmb->numFields;
497: if (dmb->numFields) {
498: DMMoabSetFieldNames(dm2, dmb->numFields, dmb->fieldNames);
499: }
501: DMSetFromOptions(dm2);
503: /* recreate Dof numbering for the refined DM and make sure the distribution is correctly populated */
504: DMSetUp(dm2);
506: *dmref = dm2;
507: return(0);
508: }
511: /*@
512: DMRefine_Moab - Generate a multi-level uniform refinement hierarchy
513: by succesively refining a coarse mesh, already defined in the DM object
514: provided by the user.
516: Collective on DM
518: Input Parameter:
519: + dm - The DMMoab object
520: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
522: Output Parameter:
523: . dmf - the refined DM, or NULL
525: Note: If no refinement was done, the return value is NULL
527: Level: developer
529: .keywords: DMMoab, create, refinement
530: @*/
531: PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM* dmf)
532: {
533: PetscErrorCode ierr;
538: DMMoab_UMR_Private(dm, comm, PETSC_TRUE, dmf);
539: return(0);
540: }
542: /*@
543: DMCoarsen_Moab - Generate a multi-level uniform refinement hierarchy
544: by succesively refining a coarse mesh, already defined in the DM object
545: provided by the user.
547: Collective on DM
549: Input Parameter:
550: . dm - The DMMoab object
551: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
553: Output Parameter:
554: . dmf - the coarsened DM, or NULL
556: Note: If no coarsening was done, the return value is NULL
558: Level: developer
560: .keywords: DMMoab, create, coarsening
561: @*/
562: PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM* dmc)
563: {
564: PetscErrorCode ierr;
569: DMMoab_UMR_Private(dm, comm, PETSC_FALSE, dmc);
570: return(0);
571: }