xref: /petsc/src/dm/impls/moab/dmmbmg.cxx (revision a4e35b1925eceef64945ea472b84f2bf06a67b5e)
1755f3dfbSVijay Mahadevan #include <petsc/private/dmmbimpl.h>
2b117cd09SVijay Mahadevan #include <petscdmmoab.h>
3b117cd09SVijay Mahadevan #include <MBTagConventions.hpp>
4b117cd09SVijay Mahadevan #include <moab/NestedRefine.hpp>
5b117cd09SVijay Mahadevan 
6cab5ea25SPierre Jolivet /*@C
7b117cd09SVijay Mahadevan   DMMoabGenerateHierarchy - Generate a multi-level uniform refinement hierarchy
820f4b53cSBarry Smith   by succesively refining a coarse mesh, already defined in the `DM` object
9b117cd09SVijay Mahadevan   provided by the user.
10b117cd09SVijay Mahadevan 
11d083f849SBarry Smith   Collective
12b117cd09SVijay Mahadevan 
13b117cd09SVijay Mahadevan   Input Parameter:
1460225df5SJacob Faibussowitsch . dm - The `DMMOAB` object
15b117cd09SVijay Mahadevan 
16d8d19677SJose E. Roman   Output Parameters:
17b117cd09SVijay Mahadevan + nlevels  - The number of levels of refinement needed to generate the hierarchy
18a2b725a8SWilliam Gropp - ldegrees - The degree of refinement at each level in the hierarchy
19b117cd09SVijay Mahadevan 
20b117cd09SVijay Mahadevan   Level: beginner
21b117cd09SVijay Mahadevan 
22*a4e35b19SJacob Faibussowitsch .seealso: `DMMoabCreate()`
23b117cd09SVijay Mahadevan @*/
24d71ae5a4SJacob Faibussowitsch PetscErrorCode DMMoabGenerateHierarchy(DM dm, PetscInt nlevels, PetscInt *ldegrees)
25d71ae5a4SJacob Faibussowitsch {
26b117cd09SVijay Mahadevan   DM_Moab                        *dmmoab;
27b117cd09SVijay Mahadevan   moab::ErrorCode                 merr;
282417220eSVijay Mahadevan   PetscInt                       *pdegrees, ilevel;
29e882eb38SVijay Mahadevan   std::vector<moab::EntityHandle> hsets;
30b117cd09SVijay Mahadevan 
31b117cd09SVijay Mahadevan   PetscFunctionBegin;
32b117cd09SVijay Mahadevan   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
33b117cd09SVijay Mahadevan   dmmoab = (DM_Moab *)(dm)->data;
34b117cd09SVijay Mahadevan 
35b117cd09SVijay Mahadevan   if (!ldegrees) {
369566063dSJacob Faibussowitsch     PetscCall(PetscMalloc1(nlevels, &pdegrees));
372417220eSVijay Mahadevan     for (ilevel = 0; ilevel < nlevels; ilevel++) pdegrees[ilevel] = 2; /* default = Degree 2 refinement */
389371c9d4SSatish Balay   } else pdegrees = ldegrees;
39b117cd09SVijay Mahadevan 
40b117cd09SVijay Mahadevan   /* initialize set level refinement data for hierarchy */
41b117cd09SVijay Mahadevan   dmmoab->nhlevels = nlevels;
42b117cd09SVijay Mahadevan 
43b117cd09SVijay Mahadevan   /* Instantiate the nested refinement class */
449daf19fdSVijay Mahadevan #ifdef MOAB_HAVE_MPI
453f1c6e43SVijay Mahadevan   dmmoab->hierarchy = new moab::NestedRefine(dynamic_cast<moab::Core *>(dmmoab->mbiface), dmmoab->pcomm, dmmoab->fileset);
469daf19fdSVijay Mahadevan #else
479daf19fdSVijay Mahadevan   dmmoab->hierarchy = new moab::NestedRefine(dynamic_cast<moab::Core *>(dmmoab->mbiface), NULL, dmmoab->fileset);
489daf19fdSVijay Mahadevan #endif
49b117cd09SVijay Mahadevan 
509566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(nlevels + 1, &dmmoab->hsets));
51e882eb38SVijay Mahadevan 
52b117cd09SVijay Mahadevan   /* generate the mesh hierarchy */
539371c9d4SSatish Balay   merr = dmmoab->hierarchy->generate_mesh_hierarchy(nlevels, pdegrees, hsets, false);
549371c9d4SSatish Balay   MBERRNM(merr);
55e882eb38SVijay Mahadevan 
569daf19fdSVijay Mahadevan #ifdef MOAB_HAVE_MPI
57755f3dfbSVijay Mahadevan   if (dmmoab->pcomm->size() > 1) {
589371c9d4SSatish Balay     merr = dmmoab->hierarchy->exchange_ghosts(hsets, dmmoab->nghostrings);
599371c9d4SSatish Balay     MBERRNM(merr);
60755f3dfbSVijay Mahadevan   }
619daf19fdSVijay Mahadevan #endif
6249d66b22SVijay Mahadevan 
6364e1c140SVijay Mahadevan   /* copy the mesh sets for nested refinement hierarchy */
64e92d1c7cSVijay Mahadevan   dmmoab->hsets[0] = hsets[0];
659371c9d4SSatish Balay   for (ilevel = 1; ilevel <= nlevels; ilevel++) {
662417220eSVijay Mahadevan     dmmoab->hsets[ilevel] = hsets[ilevel];
67e92d1c7cSVijay Mahadevan 
689c368985SVijay Mahadevan #ifdef MOAB_HAVE_MPI
699371c9d4SSatish Balay     merr = dmmoab->pcomm->assign_global_ids(hsets[ilevel], dmmoab->dim, 0, false, true, false);
709371c9d4SSatish Balay     MBERRNM(merr);
719c368985SVijay Mahadevan #endif
72e92d1c7cSVijay Mahadevan 
73e92d1c7cSVijay Mahadevan     /* Update material and other geometric tags from parent to child sets */
749371c9d4SSatish Balay     merr = dmmoab->hierarchy->update_special_tags(ilevel, hsets[ilevel]);
759371c9d4SSatish Balay     MBERRNM(merr);
762417220eSVijay Mahadevan   }
7764e1c140SVijay Mahadevan 
7864e1c140SVijay Mahadevan   hsets.clear();
7948a46eb9SPierre Jolivet   if (!ldegrees) PetscCall(PetscFree(pdegrees));
803ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
81b117cd09SVijay Mahadevan }
82b117cd09SVijay Mahadevan 
83*a4e35b19SJacob Faibussowitsch // PetscClangLinter pragma ignore: -fdoc-*
84*a4e35b19SJacob Faibussowitsch /*
8520f4b53cSBarry Smith   DMRefineHierarchy_Moab - Generate a multi-level `DM` hierarchy
86e882eb38SVijay Mahadevan   by succesively refining a coarse mesh.
87b117cd09SVijay Mahadevan 
88d083f849SBarry Smith   Collective
89b117cd09SVijay Mahadevan 
90b117cd09SVijay Mahadevan   Input Parameter:
9120f4b53cSBarry Smith . dm - The `DMMOAB` object
92b117cd09SVijay Mahadevan 
93d8d19677SJose E. Roman   Output Parameters:
94e882eb38SVijay Mahadevan + nlevels - The number of levels of refinement needed to generate the hierarchy
95a2b725a8SWilliam Gropp - dmf     - The DM objects after successive refinement of the hierarchy
96b117cd09SVijay Mahadevan 
97b117cd09SVijay Mahadevan   Level: beginner
98*a4e35b19SJacob Faibussowitsch */
99d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmf[])
100d71ae5a4SJacob Faibussowitsch {
101b117cd09SVijay Mahadevan   PetscInt i;
102b117cd09SVijay Mahadevan 
103b117cd09SVijay Mahadevan   PetscFunctionBegin;
104b117cd09SVijay Mahadevan 
1059566063dSJacob Faibussowitsch   PetscCall(DMRefine(dm, PetscObjectComm((PetscObject)dm), &dmf[0]));
10648a46eb9SPierre Jolivet   for (i = 1; i < nlevels; i++) PetscCall(DMRefine(dmf[i - 1], PetscObjectComm((PetscObject)dm), &dmf[i]));
1073ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
108b117cd09SVijay Mahadevan }
109b117cd09SVijay Mahadevan 
110*a4e35b19SJacob Faibussowitsch // PetscClangLinter pragma ignore: -fdoc-*
111*a4e35b19SJacob Faibussowitsch /*
11220f4b53cSBarry Smith   DMCoarsenHierarchy_Moab - Generate a multi-level `DM` hierarchy
113e882eb38SVijay Mahadevan   by succesively coarsening a refined mesh.
114b117cd09SVijay Mahadevan 
115d083f849SBarry Smith   Collective
116b117cd09SVijay Mahadevan 
117b117cd09SVijay Mahadevan   Input Parameter:
11820f4b53cSBarry Smith . dm - The `DMMOAB` object
119b117cd09SVijay Mahadevan 
120d8d19677SJose E. Roman   Output Parameters:
121e882eb38SVijay Mahadevan + nlevels - The number of levels of refinement needed to generate the hierarchy
12220f4b53cSBarry Smith - dmc     - The `DM` objects after successive coarsening of the hierarchy
123b117cd09SVijay Mahadevan 
124b117cd09SVijay Mahadevan   Level: beginner
125*a4e35b19SJacob Faibussowitsch */
126d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmc[])
127d71ae5a4SJacob Faibussowitsch {
128b117cd09SVijay Mahadevan   PetscInt i;
129b117cd09SVijay Mahadevan 
130b117cd09SVijay Mahadevan   PetscFunctionBegin;
131b117cd09SVijay Mahadevan 
1329566063dSJacob Faibussowitsch   PetscCall(DMCoarsen(dm, PetscObjectComm((PetscObject)dm), &dmc[0]));
13348a46eb9SPierre Jolivet   for (i = 1; i < nlevels; i++) PetscCall(DMCoarsen(dmc[i - 1], PetscObjectComm((PetscObject)dm), &dmc[i]));
1343ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
135b117cd09SVijay Mahadevan }
136b117cd09SVijay Mahadevan 
1372417220eSVijay Mahadevan PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM, PetscInt *, PetscInt *, PetscInt *, PetscInt *, PetscBool);
138b117cd09SVijay Mahadevan 
139*a4e35b19SJacob Faibussowitsch // PetscClangLinter pragma ignore: -fdoc-*
140*a4e35b19SJacob Faibussowitsch /*
141e882eb38SVijay Mahadevan   DMCreateInterpolation_Moab - Generate the interpolation operators to transform
142e882eb38SVijay Mahadevan   operators (matrices, vectors) from parent level to child level as defined by
14320f4b53cSBarry Smith   the `DM` inputs provided by the user.
144b117cd09SVijay Mahadevan 
145d083f849SBarry Smith   Collective
146b117cd09SVijay Mahadevan 
147d8d19677SJose E. Roman   Input Parameters:
14860225df5SJacob Faibussowitsch + dmp - The `DMMOAB` object
14960225df5SJacob Faibussowitsch - dmc - the second, finer `DMMOAB` object
150b117cd09SVijay Mahadevan 
151d8d19677SJose E. Roman   Output Parameters:
152e882eb38SVijay Mahadevan + interpl - The interpolation operator for transferring data between the levels
153e882eb38SVijay Mahadevan - vec     - The scaling vector (optional)
154b117cd09SVijay Mahadevan 
155e882eb38SVijay Mahadevan   Level: developer
156*a4e35b19SJacob Faibussowitsch */
157d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmp, DM dmc, Mat *interpl, Vec *vec)
158d71ae5a4SJacob Faibussowitsch {
159755f3dfbSVijay Mahadevan   DM_Moab        *dmbp, *dmbc;
160b117cd09SVijay Mahadevan   moab::ErrorCode merr;
161e882eb38SVijay Mahadevan   PetscInt        dim;
1623f1c6e43SVijay Mahadevan   PetscReal       factor;
163ce27a4eeSVijay Mahadevan   PetscInt        innz, *nnz, ionz, *onz;
164755f3dfbSVijay Mahadevan   PetscInt        nlsizp, nlsizc, nlghsizp, ngsizp, ngsizc;
165a86ed394SVijay Mahadevan   const PetscBool use_consistent_bases = PETSC_TRUE;
166b117cd09SVijay Mahadevan 
167b117cd09SVijay Mahadevan   PetscFunctionBegin;
168755f3dfbSVijay Mahadevan   PetscValidHeaderSpecific(dmp, DM_CLASSID, 1);
169755f3dfbSVijay Mahadevan   PetscValidHeaderSpecific(dmc, DM_CLASSID, 2);
170755f3dfbSVijay Mahadevan   dmbp     = (DM_Moab *)(dmp)->data;
171755f3dfbSVijay Mahadevan   dmbc     = (DM_Moab *)(dmc)->data;
172755f3dfbSVijay Mahadevan   nlsizp   = dmbp->nloc;                  // *dmb1->numFields;
173755f3dfbSVijay Mahadevan   nlsizc   = dmbc->nloc;                  // *dmb2->numFields;
174755f3dfbSVijay Mahadevan   ngsizp   = dmbp->n;                     // *dmb1->numFields;
175755f3dfbSVijay Mahadevan   ngsizc   = dmbc->n;                     // *dmb2->numFields;
176755f3dfbSVijay Mahadevan   nlghsizp = (dmbp->nloc + dmbp->nghost); // *dmb1->numFields;
177b117cd09SVijay Mahadevan 
1782417220eSVijay Mahadevan   // Columns = Parent DoFs ;  Rows = Child DoFs
179755f3dfbSVijay Mahadevan   // Interpolation matrix: \sum_{i=1}^P Owned(Child) * (Owned(Parent) + Ghosted(Parent))
180755f3dfbSVijay Mahadevan   // Size: nlsizc * nlghsizp
1813ba16761SJacob Faibussowitsch   PetscCall(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));
182b117cd09SVijay Mahadevan 
1839566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(dmp, &dim));
184a86ed394SVijay Mahadevan 
185941e0cffSVijay Mahadevan   /* allocate the nnz, onz arrays based on block size and local nodes */
1869566063dSJacob Faibussowitsch   PetscCall(PetscCalloc2(nlsizc, &nnz, nlsizc, &onz));
187941e0cffSVijay Mahadevan 
188941e0cffSVijay Mahadevan   /* Loop through the local elements and compute the relation between the current parent and the refined_level. */
1892417220eSVijay Mahadevan   for (moab::Range::iterator iter = dmbc->vowned->begin(); iter != dmbc->vowned->end(); iter++) {
1902417220eSVijay Mahadevan     const moab::EntityHandle vhandle = *iter;
1912417220eSVijay Mahadevan     /* define local variables */
1922417220eSVijay Mahadevan     moab::EntityHandle              parent;
1932417220eSVijay Mahadevan     std::vector<moab::EntityHandle> adjs;
1942417220eSVijay Mahadevan     moab::Range                     found;
1952417220eSVijay Mahadevan 
1962417220eSVijay Mahadevan     /* store the vertex DoF number */
1972417220eSVijay Mahadevan     const int ldof = dmbc->lidmap[vhandle - dmbc->seqstart];
1982417220eSVijay Mahadevan 
1992417220eSVijay Mahadevan     /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects
2002417220eSVijay Mahadevan        to the current vertex. We can then decipher if a vertex is ghosted or not and compute the
2012417220eSVijay Mahadevan        non-zero pattern accordingly. */
2029371c9d4SSatish Balay     merr = dmbc->hierarchy->get_adjacencies(vhandle, dmbc->dim, adjs);
2039371c9d4SSatish Balay     MBERRNM(merr);
2042417220eSVijay Mahadevan 
2052417220eSVijay Mahadevan     /* loop over vertices and update the number of connectivity */
2062417220eSVijay Mahadevan     for (unsigned jter = 0; jter < adjs.size(); jter++) {
2072417220eSVijay Mahadevan       const moab::EntityHandle jhandle = adjs[jter];
208941e0cffSVijay Mahadevan 
209941e0cffSVijay Mahadevan       /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */
2109371c9d4SSatish Balay       merr = dmbc->hierarchy->child_to_parent(jhandle, dmbc->hlevel, dmbp->hlevel, &parent);
2119371c9d4SSatish Balay       MBERRNM(merr);
212941e0cffSVijay Mahadevan 
2132417220eSVijay Mahadevan       /* Get connectivity information in canonical ordering for the local element */
2142417220eSVijay Mahadevan       std::vector<moab::EntityHandle> connp;
2159371c9d4SSatish Balay       merr = dmbp->hierarchy->get_connectivity(parent, dmbp->hlevel, connp);
2169371c9d4SSatish Balay       MBERRNM(merr);
217a86ed394SVijay Mahadevan 
2182417220eSVijay Mahadevan       for (unsigned ic = 0; ic < connp.size(); ++ic) {
2192417220eSVijay Mahadevan         /* loop over each element connected to the adjacent vertex and update as needed */
2202417220eSVijay Mahadevan         /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */
2212417220eSVijay Mahadevan         if (found.find(connp[ic]) != found.end()) continue;                    /* make sure we don't double count shared vertices */
2222417220eSVijay Mahadevan         if (dmbp->vghost->find(connp[ic]) != dmbp->vghost->end()) onz[ldof]++; /* update out-of-proc onz */
2232417220eSVijay Mahadevan         else nnz[ldof]++;                                                      /* else local vertex */
2242417220eSVijay Mahadevan         found.insert(connp[ic]);
2252417220eSVijay Mahadevan       }
2262417220eSVijay Mahadevan     }
227941e0cffSVijay Mahadevan   }
228941e0cffSVijay Mahadevan 
2299371c9d4SSatish Balay   for (int i = 0; i < nlsizc; i++) nnz[i] += 1; /* self count the node */
230941e0cffSVijay Mahadevan 
231941e0cffSVijay Mahadevan   ionz = onz[0];
232941e0cffSVijay Mahadevan   innz = nnz[0];
233755f3dfbSVijay Mahadevan   for (int tc = 0; tc < nlsizc; tc++) {
234941e0cffSVijay Mahadevan     // check for maximum allowed sparsity = fully dense
235755f3dfbSVijay Mahadevan     nnz[tc] = std::min(nlsizp, nnz[tc]);
2362417220eSVijay Mahadevan     onz[tc] = std::min(ngsizp - nlsizp, onz[tc]);
2372417220eSVijay Mahadevan 
2383ba16761SJacob Faibussowitsch     PetscCall(PetscInfo(NULL, "  %d: NNZ = %d, ONZ = %d\n", tc, nnz[tc], onz[tc]));
239941e0cffSVijay Mahadevan 
240941e0cffSVijay Mahadevan     innz = (innz < nnz[tc] ? nnz[tc] : innz);
241941e0cffSVijay Mahadevan     ionz = (ionz < onz[tc] ? onz[tc] : ionz);
242941e0cffSVijay Mahadevan   }
24364e1c140SVijay Mahadevan 
24464e1c140SVijay Mahadevan   /* create interpolation matrix */
2459566063dSJacob Faibussowitsch   PetscCall(MatCreate(PetscObjectComm((PetscObject)dmc), interpl));
2469566063dSJacob Faibussowitsch   PetscCall(MatSetSizes(*interpl, nlsizc, nlsizp, ngsizc, ngsizp));
2479566063dSJacob Faibussowitsch   PetscCall(MatSetType(*interpl, MATAIJ));
2489566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(*interpl));
24964e1c140SVijay Mahadevan 
2509566063dSJacob Faibussowitsch   PetscCall(MatSeqAIJSetPreallocation(*interpl, innz, nnz));
2519566063dSJacob Faibussowitsch   PetscCall(MatMPIAIJSetPreallocation(*interpl, innz, nnz, ionz, onz));
252941e0cffSVijay Mahadevan 
253941e0cffSVijay Mahadevan   /* clean up temporary memory */
2549566063dSJacob Faibussowitsch   PetscCall(PetscFree2(nnz, onz));
255b117cd09SVijay Mahadevan 
256b117cd09SVijay Mahadevan   /* set up internal matrix data-structures */
2579566063dSJacob Faibussowitsch   PetscCall(MatSetUp(*interpl));
258b117cd09SVijay Mahadevan 
259a86ed394SVijay Mahadevan   /* Define variables for assembly */
260a86ed394SVijay Mahadevan   std::vector<moab::EntityHandle> children;
261a86ed394SVijay Mahadevan   std::vector<moab::EntityHandle> connp, connc;
262a86ed394SVijay Mahadevan   std::vector<PetscReal>          pcoords, ccoords, values_phi;
263a86ed394SVijay Mahadevan 
264a86ed394SVijay Mahadevan   if (use_consistent_bases) {
2652417220eSVijay Mahadevan     const moab::EntityHandle ehandle = dmbp->elocal->front();
266a86ed394SVijay Mahadevan 
2679371c9d4SSatish Balay     merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children);
2689371c9d4SSatish Balay     MBERRNM(merr);
269a86ed394SVijay Mahadevan 
270a86ed394SVijay Mahadevan     /* Get connectivity and coordinates of the parent vertices */
2719371c9d4SSatish Balay     merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp);
2729371c9d4SSatish Balay     MBERRNM(merr);
2739371c9d4SSatish Balay     merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc);
2749371c9d4SSatish Balay     MBERRNM(merr);
275a86ed394SVijay Mahadevan 
276a86ed394SVijay Mahadevan     std::vector<PetscReal> natparam(3 * connc.size(), 0.0);
277a86ed394SVijay Mahadevan     pcoords.resize(connp.size() * 3);
278a86ed394SVijay Mahadevan     ccoords.resize(connc.size() * 3);
279a86ed394SVijay Mahadevan     values_phi.resize(connp.size() * connc.size());
280a86ed394SVijay Mahadevan     /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */
2819371c9d4SSatish Balay     merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]);
2829371c9d4SSatish Balay     MBERRNM(merr);
2839371c9d4SSatish Balay     merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]);
2849371c9d4SSatish Balay     MBERRNM(merr);
285a86ed394SVijay Mahadevan 
286a86ed394SVijay Mahadevan     /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */
287a86ed394SVijay Mahadevan     for (unsigned tc = 0; tc < connc.size(); tc++) {
288a86ed394SVijay Mahadevan       const PetscInt offset = tc * 3;
289a86ed394SVijay Mahadevan 
290a86ed394SVijay Mahadevan       /* Scale ccoords relative to pcoords */
2919566063dSJacob Faibussowitsch       PetscCall(DMMoabPToRMapping(dim, connp.size(), &pcoords[0], &ccoords[offset], &natparam[offset], &values_phi[connp.size() * tc]));
292a86ed394SVijay Mahadevan     }
2931baa6e33SBarry Smith   } else {
294a86ed394SVijay Mahadevan     factor = std::pow(2.0 /*degree_P_for_refinement*/, (dmbc->hlevel - dmbp->hlevel) * dmbp->dim * 1.0);
295a86ed394SVijay Mahadevan   }
296a86ed394SVijay Mahadevan 
297755f3dfbSVijay Mahadevan   /* TODO: Decipher the correct non-zero pattern. There is still some issue with onz allocation */
2989566063dSJacob Faibussowitsch   PetscCall(MatSetOption(*interpl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
299b117cd09SVijay Mahadevan 
300e882eb38SVijay Mahadevan   /* Loop through the remaining vertices. These vertices appear only on the current refined_level. */
3012417220eSVijay Mahadevan   for (moab::Range::iterator iter = dmbp->elocal->begin(); iter != dmbp->elocal->end(); iter++) {
302b117cd09SVijay Mahadevan     const moab::EntityHandle ehandle = *iter;
303b117cd09SVijay Mahadevan 
304b117cd09SVijay Mahadevan     /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */
305a86ed394SVijay Mahadevan     children.clear();
306a86ed394SVijay Mahadevan     connc.clear();
3079371c9d4SSatish Balay     merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children);
3089371c9d4SSatish Balay     MBERRNM(merr);
309b117cd09SVijay Mahadevan 
310b117cd09SVijay Mahadevan     /* Get connectivity and coordinates of the parent vertices */
3119371c9d4SSatish Balay     merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp);
3129371c9d4SSatish Balay     MBERRNM(merr);
3139371c9d4SSatish Balay     merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc);
3149371c9d4SSatish Balay     MBERRNM(merr);
315b117cd09SVijay Mahadevan 
316a86ed394SVijay Mahadevan     pcoords.resize(connp.size() * 3);
317a86ed394SVijay Mahadevan     ccoords.resize(connc.size() * 3);
318b117cd09SVijay Mahadevan     /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */
3199371c9d4SSatish Balay     merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]);
3209371c9d4SSatish Balay     MBERRNM(merr);
3219371c9d4SSatish Balay     merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]);
3229371c9d4SSatish Balay     MBERRNM(merr);
323b117cd09SVijay Mahadevan 
324b117cd09SVijay Mahadevan     std::vector<int> dofsp(connp.size()), dofsc(connc.size());
325b117cd09SVijay Mahadevan     /* TODO: specific to scalar system - use GetDofs */
3269566063dSJacob Faibussowitsch     PetscCall(DMMoabGetDofsBlocked(dmp, connp.size(), &connp[0], &dofsp[0]));
3279566063dSJacob Faibussowitsch     PetscCall(DMMoabGetDofsBlocked(dmc, connc.size(), &connc[0], &dofsc[0]));
328b117cd09SVijay Mahadevan 
329a86ed394SVijay Mahadevan     /* Compute the actual interpolation weights when projecting solution/residual between levels */
330a86ed394SVijay Mahadevan     if (use_consistent_bases) {
331145b44c9SPierre Jolivet       /* Use the cached values of natural parametric coordinates and basis pre-evaluated.
332a86ed394SVijay Mahadevan          We are making an assumption here that UMR used in GMG to generate the hierarchy uses
333a86ed394SVijay Mahadevan          the same template for all elements; This will fail for mixed element meshes (TRI/QUAD).
334a86ed394SVijay Mahadevan 
3352417220eSVijay Mahadevan          TODO: Fix the above assumption by caching data for families (especially for Tets and mixed meshes)
336a86ed394SVijay Mahadevan       */
337a86ed394SVijay Mahadevan 
338a86ed394SVijay Mahadevan       /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */
339a86ed394SVijay Mahadevan       for (unsigned tc = 0; tc < connc.size(); tc++) {
340a86ed394SVijay Mahadevan         /* TODO: Check if we should be using INSERT_VALUES instead */
3419566063dSJacob Faibussowitsch         PetscCall(MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi[connp.size() * tc], ADD_VALUES));
342a86ed394SVijay Mahadevan       }
3431baa6e33SBarry Smith     } else {
344e882eb38SVijay Mahadevan       /* Compute the interpolation weights by determining distance of 1-ring
345755f3dfbSVijay Mahadevan          neighbor vertices from current vertex
346755f3dfbSVijay Mahadevan 
347a86ed394SVijay Mahadevan          This should be used only when FEM basis is not used for the discretization.
348a86ed394SVijay Mahadevan          Else, the consistent interface to compute the basis function for interpolation
349a86ed394SVijay Mahadevan          between the levels should be evaluated correctly to preserve convergence of GMG.
3502417220eSVijay Mahadevan          Shephard's basis will be terrible for any unsmooth problems.
351755f3dfbSVijay Mahadevan       */
352755f3dfbSVijay Mahadevan       values_phi.resize(connp.size());
353755f3dfbSVijay Mahadevan       for (unsigned tc = 0; tc < connc.size(); tc++) {
354a86ed394SVijay Mahadevan         PetscReal normsum = 0.0;
355755f3dfbSVijay Mahadevan         for (unsigned tp = 0; tp < connp.size(); tp++) {
356755f3dfbSVijay Mahadevan           values_phi[tp] = 0.0;
3579371c9d4SSatish Balay           for (unsigned k = 0; k < 3; k++) values_phi[tp] += std::pow(pcoords[tp * 3 + k] - ccoords[k + tc * 3], dim);
358755f3dfbSVijay Mahadevan           if (values_phi[tp] < 1e-12) {
359755f3dfbSVijay Mahadevan             values_phi[tp] = 1e12;
3601baa6e33SBarry Smith           } else {
361755f3dfbSVijay Mahadevan             //values_phi[tp] = std::pow(values_phi[tp], -1.0/dim);
362755f3dfbSVijay Mahadevan             values_phi[tp] = std::pow(values_phi[tp], -1.0);
363755f3dfbSVijay Mahadevan             normsum += values_phi[tp];
364b117cd09SVijay Mahadevan           }
365b117cd09SVijay Mahadevan         }
366755f3dfbSVijay Mahadevan         for (unsigned tp = 0; tp < connp.size(); tp++) {
3679371c9d4SSatish Balay           if (values_phi[tp] > 1e11) values_phi[tp] = factor * 0.5 / connp.size();
3689371c9d4SSatish Balay           else values_phi[tp] = factor * values_phi[tp] * 0.5 / (connp.size() * normsum);
369b117cd09SVijay Mahadevan         }
3709566063dSJacob Faibussowitsch         PetscCall(MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi[0], ADD_VALUES));
371b117cd09SVijay Mahadevan       }
372a86ed394SVijay Mahadevan     }
373b117cd09SVijay Mahadevan   }
374304006b3SVijay Mahadevan   if (vec) *vec = NULL;
3759566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(*interpl, MAT_FINAL_ASSEMBLY));
3769566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(*interpl, MAT_FINAL_ASSEMBLY));
3773ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
378b117cd09SVijay Mahadevan }
379b117cd09SVijay Mahadevan 
380*a4e35b19SJacob Faibussowitsch // PetscClangLinter pragma ignore: -fdoc-*
381*a4e35b19SJacob Faibussowitsch /*
382b117cd09SVijay Mahadevan   DMCreateInjection_Moab - Generate a multi-level uniform refinement hierarchy
38320f4b53cSBarry Smith   by succesively refining a coarse mesh, already defined in the `DM` object
384b117cd09SVijay Mahadevan   provided by the user.
385b117cd09SVijay Mahadevan 
386d083f849SBarry Smith   Collective
387b117cd09SVijay Mahadevan 
388b117cd09SVijay Mahadevan   Input Parameter:
38920f4b53cSBarry Smith . dmb  - The `DMMOAB` object
390b117cd09SVijay Mahadevan 
391d8d19677SJose E. Roman   Output Parameters:
392a2b725a8SWilliam Gropp + nlevels   - The number of levels of refinement needed to generate the hierarchy
393a2b725a8SWilliam Gropp - ldegrees  - The degree of refinement at each level in the hierarchy
394b117cd09SVijay Mahadevan 
395b117cd09SVijay Mahadevan   Level: beginner
396*a4e35b19SJacob Faibussowitsch */
397d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dm1, DM dm2, VecScatter *ctx)
398d71ae5a4SJacob Faibussowitsch {
399e882eb38SVijay Mahadevan   //DM_Moab        *dmmoab;
400b117cd09SVijay Mahadevan 
401b117cd09SVijay Mahadevan   PetscFunctionBegin;
402b117cd09SVijay Mahadevan   PetscValidHeaderSpecific(dm1, DM_CLASSID, 1);
403b117cd09SVijay Mahadevan   PetscValidHeaderSpecific(dm2, DM_CLASSID, 2);
404e882eb38SVijay Mahadevan   //dmmoab = (DM_Moab*)(dm1)->data;
405b117cd09SVijay Mahadevan 
4063ba16761SJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[DMCreateInjection_Moab] :: Placeholder\n"));
4073ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
408b117cd09SVijay Mahadevan }
409b117cd09SVijay Mahadevan 
410d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMMoab_UMR_Private(DM dm, MPI_Comm comm, PetscBool refine, DM *dmref)
411d71ae5a4SJacob Faibussowitsch {
412e882eb38SVijay Mahadevan   PetscInt        i, dim;
413b117cd09SVijay Mahadevan   DM              dm2;
414b117cd09SVijay Mahadevan   moab::ErrorCode merr;
415b117cd09SVijay Mahadevan   DM_Moab        *dmb = (DM_Moab *)dm->data, *dd2;
416b117cd09SVijay Mahadevan 
417b117cd09SVijay Mahadevan   PetscFunctionBegin;
418b117cd09SVijay Mahadevan   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
4194f572ea9SToby Isaac   PetscAssertPointer(dmref, 4);
420b117cd09SVijay Mahadevan 
421b117cd09SVijay Mahadevan   if ((dmb->hlevel == dmb->nhlevels && refine) || (dmb->hlevel == 0 && !refine)) {
4223ba16761SJacob Faibussowitsch     if (dmb->hlevel + 1 > dmb->nhlevels && refine) {
4233ba16761SJacob Faibussowitsch       PetscCall(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));
4243ba16761SJacob Faibussowitsch     }
4253ba16761SJacob Faibussowitsch     if (dmb->hlevel - 1 < 0 && !refine) PetscCall(PetscInfo(NULL, "Invalid multigrid coarsen hierarchy level specified (%" PetscInt_FMT "). Creating a NULL object.\n", dmb->hlevel - 1));
426f3fa974cSJacob Faibussowitsch     *dmref = NULL;
4273ba16761SJacob Faibussowitsch     PetscFunctionReturn(PETSC_SUCCESS);
428b117cd09SVijay Mahadevan   }
429b117cd09SVijay Mahadevan 
4309566063dSJacob Faibussowitsch   PetscCall(DMMoabCreate(PetscObjectComm((PetscObject)dm), &dm2));
431b117cd09SVijay Mahadevan   dd2 = (DM_Moab *)dm2->data;
432b117cd09SVijay Mahadevan 
433b117cd09SVijay Mahadevan   dd2->mbiface = dmb->mbiface;
4349daf19fdSVijay Mahadevan #ifdef MOAB_HAVE_MPI
435b117cd09SVijay Mahadevan   dd2->pcomm = dmb->pcomm;
4369daf19fdSVijay Mahadevan #endif
437b117cd09SVijay Mahadevan   dd2->icreatedinstance = PETSC_FALSE;
43864e1c140SVijay Mahadevan   dd2->nghostrings      = dmb->nghostrings;
439b117cd09SVijay Mahadevan 
440b117cd09SVijay Mahadevan   /* set the new level based on refinement/coarsening */
441b117cd09SVijay Mahadevan   if (refine) {
442b117cd09SVijay Mahadevan     dd2->hlevel = dmb->hlevel + 1;
4431baa6e33SBarry Smith   } else {
444b117cd09SVijay Mahadevan     dd2->hlevel = dmb->hlevel - 1;
445b117cd09SVijay Mahadevan   }
446b117cd09SVijay Mahadevan 
447b117cd09SVijay Mahadevan   /* Copy the multilevel hierarchy pointers in MOAB */
448b117cd09SVijay Mahadevan   dd2->hierarchy = dmb->hierarchy;
449b117cd09SVijay Mahadevan   dd2->nhlevels  = dmb->nhlevels;
4509566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(dd2->nhlevels + 1, &dd2->hsets));
451ad540459SPierre Jolivet   for (i = 0; i <= dd2->nhlevels; i++) dd2->hsets[i] = dmb->hsets[i];
452b117cd09SVijay Mahadevan   dd2->fileset = dd2->hsets[dd2->hlevel];
453b117cd09SVijay Mahadevan 
454b117cd09SVijay Mahadevan   /* do the remaining initializations for DMMoab */
455b117cd09SVijay Mahadevan   dd2->bs                = dmb->bs;
456b117cd09SVijay Mahadevan   dd2->numFields         = dmb->numFields;
457b117cd09SVijay Mahadevan   dd2->rw_dbglevel       = dmb->rw_dbglevel;
458b117cd09SVijay Mahadevan   dd2->partition_by_rank = dmb->partition_by_rank;
459c6a7a370SJeremy L Thompson   PetscCall(PetscStrncpy(dd2->extra_read_options, dmb->extra_read_options, sizeof(dd2->extra_read_options)));
460c6a7a370SJeremy L Thompson   PetscCall(PetscStrncpy(dd2->extra_write_options, dmb->extra_write_options, sizeof(dd2->extra_write_options)));
461b117cd09SVijay Mahadevan   dd2->read_mode  = dmb->read_mode;
462b117cd09SVijay Mahadevan   dd2->write_mode = dmb->write_mode;
463b117cd09SVijay Mahadevan 
464b117cd09SVijay Mahadevan   /* set global ID tag handle */
4659566063dSJacob Faibussowitsch   PetscCall(DMMoabSetLocalToGlobalTag(dm2, dmb->ltog_tag));
466b117cd09SVijay Mahadevan 
4679371c9d4SSatish Balay   merr = dd2->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dd2->material_tag);
4689371c9d4SSatish Balay   MBERRNM(merr);
469b117cd09SVijay Mahadevan 
4709566063dSJacob Faibussowitsch   PetscCall(DMSetOptionsPrefix(dm2, ((PetscObject)dm)->prefix));
4719566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(dm, &dim));
4729566063dSJacob Faibussowitsch   PetscCall(DMSetDimension(dm2, dim));
473b117cd09SVijay Mahadevan 
474b117cd09SVijay Mahadevan   /* allow overloaded (user replaced) operations to be inherited by refinement clones */
475b117cd09SVijay Mahadevan   dm2->ops->creatematrix = dm->ops->creatematrix;
476b117cd09SVijay Mahadevan 
477b117cd09SVijay Mahadevan   /* copy fill information if given */
4789566063dSJacob Faibussowitsch   PetscCall(DMMoabSetBlockFills(dm2, dmb->dfill, dmb->ofill));
479b117cd09SVijay Mahadevan 
480b117cd09SVijay Mahadevan   /* copy vector type information */
4819566063dSJacob Faibussowitsch   PetscCall(DMSetMatType(dm2, dm->mattype));
4829566063dSJacob Faibussowitsch   PetscCall(DMSetVecType(dm2, dm->vectype));
4833f1c6e43SVijay Mahadevan   dd2->numFields = dmb->numFields;
48448a46eb9SPierre Jolivet   if (dmb->numFields) PetscCall(DMMoabSetFieldNames(dm2, dmb->numFields, dmb->fieldNames));
485b117cd09SVijay Mahadevan 
4869566063dSJacob Faibussowitsch   PetscCall(DMSetFromOptions(dm2));
487b117cd09SVijay Mahadevan 
488b117cd09SVijay Mahadevan   /* recreate Dof numbering for the refined DM and make sure the distribution is correctly populated */
4899566063dSJacob Faibussowitsch   PetscCall(DMSetUp(dm2));
490b117cd09SVijay Mahadevan 
491b117cd09SVijay Mahadevan   *dmref = dm2;
4923ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
493b117cd09SVijay Mahadevan }
494b117cd09SVijay Mahadevan 
495*a4e35b19SJacob Faibussowitsch // PetscClangLinter pragma ignore: -fdoc-*
496*a4e35b19SJacob Faibussowitsch /*
497b117cd09SVijay Mahadevan   DMRefine_Moab - Generate a multi-level uniform refinement hierarchy
49820f4b53cSBarry Smith   by succesively refining a coarse mesh, already defined in the `DM` object
499b117cd09SVijay Mahadevan   provided by the user.
500b117cd09SVijay Mahadevan 
50120f4b53cSBarry Smith   Collective
502b117cd09SVijay Mahadevan 
503d8d19677SJose E. Roman   Input Parameters:
50420f4b53cSBarry Smith + dm   - The `DMMOAB` object
50520f4b53cSBarry Smith - comm - the communicator to contain the new DM object (or `MPI_COMM_NULL`)
506b117cd09SVijay Mahadevan 
507b117cd09SVijay Mahadevan   Output Parameter:
50820f4b53cSBarry Smith . dmf - the refined `DM`, or `NULL`
509e882eb38SVijay Mahadevan 
510e882eb38SVijay Mahadevan   Level: developer
511b117cd09SVijay Mahadevan 
51220f4b53cSBarry Smith   Note:
51320f4b53cSBarry Smith   If no refinement was done, the return value is `NULL`
514*a4e35b19SJacob Faibussowitsch */
515d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM *dmf)
516d71ae5a4SJacob Faibussowitsch {
517b117cd09SVijay Mahadevan   PetscFunctionBegin;
518b117cd09SVijay Mahadevan   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
519b117cd09SVijay Mahadevan 
5209566063dSJacob Faibussowitsch   PetscCall(DMMoab_UMR_Private(dm, comm, PETSC_TRUE, dmf));
5213ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
522b117cd09SVijay Mahadevan }
523b117cd09SVijay Mahadevan 
524*a4e35b19SJacob Faibussowitsch // PetscClangLinter pragma ignore: -fdoc-*
525*a4e35b19SJacob Faibussowitsch /*
526b117cd09SVijay Mahadevan   DMCoarsen_Moab - Generate a multi-level uniform refinement hierarchy
52720f4b53cSBarry Smith   by succesively refining a coarse mesh, already defined in the `DM` object
528b117cd09SVijay Mahadevan   provided by the user.
529b117cd09SVijay Mahadevan 
53020f4b53cSBarry Smith   Collective
531b117cd09SVijay Mahadevan 
532d8d19677SJose E. Roman   Input Parameters:
53320f4b53cSBarry Smith + dm   - The `DMMOAB` object
53420f4b53cSBarry Smith - comm - the communicator to contain the new `DM` object (or `MPI_COMM_NULL`)
535b117cd09SVijay Mahadevan 
536b117cd09SVijay Mahadevan   Output Parameter:
53760225df5SJacob Faibussowitsch . dmc - the coarsened `DM`, or `NULL`
538b117cd09SVijay Mahadevan 
539e882eb38SVijay Mahadevan   Level: developer
540e882eb38SVijay Mahadevan 
54120f4b53cSBarry Smith   Note:
54220f4b53cSBarry Smith   If no coarsening was done, the return value is `NULL`
543*a4e35b19SJacob Faibussowitsch */
544d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM *dmc)
545d71ae5a4SJacob Faibussowitsch {
546b117cd09SVijay Mahadevan   PetscFunctionBegin;
547b117cd09SVijay Mahadevan   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
5489566063dSJacob Faibussowitsch   PetscCall(DMMoab_UMR_Private(dm, comm, PETSC_FALSE, dmc));
5493ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
550b117cd09SVijay Mahadevan }
551