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
6ca4445c7SIlya Fursov // A helper function to convert Real vector to Scalar vector (required by MatSetValues)
VecReal_to_VecScalar(const std::vector<PetscReal> & v)7ca4445c7SIlya Fursov static inline std::vector<PetscScalar> VecReal_to_VecScalar(const std::vector<PetscReal> &v)
8ca4445c7SIlya Fursov {
9ca4445c7SIlya Fursov std::vector<PetscScalar> res(v.size());
10ca4445c7SIlya Fursov for (size_t i = 0; i < res.size(); i++) res[i] = v[i];
11ca4445c7SIlya Fursov return res;
12ca4445c7SIlya Fursov }
13ca4445c7SIlya Fursov
14cab5ea25SPierre Jolivet /*@C
15b117cd09SVijay Mahadevan DMMoabGenerateHierarchy - Generate a multi-level uniform refinement hierarchy
1620f4b53cSBarry Smith by succesively refining a coarse mesh, already defined in the `DM` object
17b117cd09SVijay Mahadevan provided by the user.
18b117cd09SVijay Mahadevan
19d083f849SBarry Smith Collective
20b117cd09SVijay Mahadevan
21b117cd09SVijay Mahadevan Input Parameter:
2260225df5SJacob Faibussowitsch . dm - The `DMMOAB` object
23b117cd09SVijay Mahadevan
24d8d19677SJose E. Roman Output Parameters:
25b117cd09SVijay Mahadevan + nlevels - The number of levels of refinement needed to generate the hierarchy
26a2b725a8SWilliam Gropp - ldegrees - The degree of refinement at each level in the hierarchy
27b117cd09SVijay Mahadevan
28b117cd09SVijay Mahadevan Level: beginner
29b117cd09SVijay Mahadevan
30a4e35b19SJacob Faibussowitsch .seealso: `DMMoabCreate()`
31b117cd09SVijay Mahadevan @*/
DMMoabGenerateHierarchy(DM dm,PetscInt nlevels,PetscInt * ldegrees)32d71ae5a4SJacob Faibussowitsch PetscErrorCode DMMoabGenerateHierarchy(DM dm, PetscInt nlevels, PetscInt *ldegrees)
33d71ae5a4SJacob Faibussowitsch {
34b117cd09SVijay Mahadevan DM_Moab *dmmoab;
35b117cd09SVijay Mahadevan moab::ErrorCode merr;
362417220eSVijay Mahadevan PetscInt *pdegrees, ilevel;
37e882eb38SVijay Mahadevan std::vector<moab::EntityHandle> hsets;
38b117cd09SVijay Mahadevan
39b117cd09SVijay Mahadevan PetscFunctionBegin;
40b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
41*57508eceSPierre Jolivet dmmoab = (DM_Moab *)dm->data;
42b117cd09SVijay Mahadevan
43b117cd09SVijay Mahadevan if (!ldegrees) {
449566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(nlevels, &pdegrees));
452417220eSVijay Mahadevan for (ilevel = 0; ilevel < nlevels; ilevel++) pdegrees[ilevel] = 2; /* default = Degree 2 refinement */
469371c9d4SSatish Balay } else pdegrees = ldegrees;
47b117cd09SVijay Mahadevan
48b117cd09SVijay Mahadevan /* initialize set level refinement data for hierarchy */
49b117cd09SVijay Mahadevan dmmoab->nhlevels = nlevels;
50b117cd09SVijay Mahadevan
51b117cd09SVijay Mahadevan /* Instantiate the nested refinement class */
529daf19fdSVijay Mahadevan #ifdef MOAB_HAVE_MPI
533f1c6e43SVijay Mahadevan dmmoab->hierarchy = new moab::NestedRefine(dynamic_cast<moab::Core *>(dmmoab->mbiface), dmmoab->pcomm, dmmoab->fileset);
549daf19fdSVijay Mahadevan #else
559daf19fdSVijay Mahadevan dmmoab->hierarchy = new moab::NestedRefine(dynamic_cast<moab::Core *>(dmmoab->mbiface), NULL, dmmoab->fileset);
569daf19fdSVijay Mahadevan #endif
57b117cd09SVijay Mahadevan
589566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(nlevels + 1, &dmmoab->hsets));
59e882eb38SVijay Mahadevan
60b117cd09SVijay Mahadevan /* generate the mesh hierarchy */
619371c9d4SSatish Balay merr = dmmoab->hierarchy->generate_mesh_hierarchy(nlevels, pdegrees, hsets, false);
629371c9d4SSatish Balay MBERRNM(merr);
63e882eb38SVijay Mahadevan
649daf19fdSVijay Mahadevan #ifdef MOAB_HAVE_MPI
65755f3dfbSVijay Mahadevan if (dmmoab->pcomm->size() > 1) {
669371c9d4SSatish Balay merr = dmmoab->hierarchy->exchange_ghosts(hsets, dmmoab->nghostrings);
679371c9d4SSatish Balay MBERRNM(merr);
68755f3dfbSVijay Mahadevan }
699daf19fdSVijay Mahadevan #endif
7049d66b22SVijay Mahadevan
7164e1c140SVijay Mahadevan /* copy the mesh sets for nested refinement hierarchy */
72e92d1c7cSVijay Mahadevan dmmoab->hsets[0] = hsets[0];
739371c9d4SSatish Balay for (ilevel = 1; ilevel <= nlevels; ilevel++) {
742417220eSVijay Mahadevan dmmoab->hsets[ilevel] = hsets[ilevel];
75e92d1c7cSVijay Mahadevan
769c368985SVijay Mahadevan #ifdef MOAB_HAVE_MPI
779371c9d4SSatish Balay merr = dmmoab->pcomm->assign_global_ids(hsets[ilevel], dmmoab->dim, 0, false, true, false);
789371c9d4SSatish Balay MBERRNM(merr);
799c368985SVijay Mahadevan #endif
80e92d1c7cSVijay Mahadevan
81e92d1c7cSVijay Mahadevan /* Update material and other geometric tags from parent to child sets */
829371c9d4SSatish Balay merr = dmmoab->hierarchy->update_special_tags(ilevel, hsets[ilevel]);
839371c9d4SSatish Balay MBERRNM(merr);
842417220eSVijay Mahadevan }
8564e1c140SVijay Mahadevan
8664e1c140SVijay Mahadevan hsets.clear();
8748a46eb9SPierre Jolivet if (!ldegrees) PetscCall(PetscFree(pdegrees));
883ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
89b117cd09SVijay Mahadevan }
90b117cd09SVijay Mahadevan
91a4e35b19SJacob Faibussowitsch // PetscClangLinter pragma ignore: -fdoc-*
92a4e35b19SJacob Faibussowitsch /*
9320f4b53cSBarry Smith DMRefineHierarchy_Moab - Generate a multi-level `DM` hierarchy
94e882eb38SVijay Mahadevan by succesively refining a coarse mesh.
95b117cd09SVijay Mahadevan
96d083f849SBarry Smith Collective
97b117cd09SVijay Mahadevan
98b117cd09SVijay Mahadevan Input Parameter:
9920f4b53cSBarry Smith . dm - The `DMMOAB` object
100b117cd09SVijay Mahadevan
101d8d19677SJose E. Roman Output Parameters:
102e882eb38SVijay Mahadevan + nlevels - The number of levels of refinement needed to generate the hierarchy
103a2b725a8SWilliam Gropp - dmf - The DM objects after successive refinement of the hierarchy
104b117cd09SVijay Mahadevan
105b117cd09SVijay Mahadevan Level: beginner
106a4e35b19SJacob Faibussowitsch */
DMRefineHierarchy_Moab(DM dm,PetscInt nlevels,DM dmf[])107d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmf[])
108d71ae5a4SJacob Faibussowitsch {
109b117cd09SVijay Mahadevan PetscInt i;
110b117cd09SVijay Mahadevan
111b117cd09SVijay Mahadevan PetscFunctionBegin;
1129566063dSJacob Faibussowitsch PetscCall(DMRefine(dm, PetscObjectComm((PetscObject)dm), &dmf[0]));
11348a46eb9SPierre Jolivet for (i = 1; i < nlevels; i++) PetscCall(DMRefine(dmf[i - 1], PetscObjectComm((PetscObject)dm), &dmf[i]));
1143ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
115b117cd09SVijay Mahadevan }
116b117cd09SVijay Mahadevan
117a4e35b19SJacob Faibussowitsch // PetscClangLinter pragma ignore: -fdoc-*
118a4e35b19SJacob Faibussowitsch /*
11920f4b53cSBarry Smith DMCoarsenHierarchy_Moab - Generate a multi-level `DM` hierarchy
120e882eb38SVijay Mahadevan by succesively coarsening a refined mesh.
121b117cd09SVijay Mahadevan
122d083f849SBarry Smith Collective
123b117cd09SVijay Mahadevan
124b117cd09SVijay Mahadevan Input Parameter:
12520f4b53cSBarry Smith . dm - The `DMMOAB` object
126b117cd09SVijay Mahadevan
127d8d19677SJose E. Roman Output Parameters:
128e882eb38SVijay Mahadevan + nlevels - The number of levels of refinement needed to generate the hierarchy
12920f4b53cSBarry Smith - dmc - The `DM` objects after successive coarsening of the hierarchy
130b117cd09SVijay Mahadevan
131b117cd09SVijay Mahadevan Level: beginner
132a4e35b19SJacob Faibussowitsch */
DMCoarsenHierarchy_Moab(DM dm,PetscInt nlevels,DM dmc[])133d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmc[])
134d71ae5a4SJacob Faibussowitsch {
135b117cd09SVijay Mahadevan PetscInt i;
136b117cd09SVijay Mahadevan
137b117cd09SVijay Mahadevan PetscFunctionBegin;
1389566063dSJacob Faibussowitsch PetscCall(DMCoarsen(dm, PetscObjectComm((PetscObject)dm), &dmc[0]));
13948a46eb9SPierre Jolivet for (i = 1; i < nlevels; i++) PetscCall(DMCoarsen(dmc[i - 1], PetscObjectComm((PetscObject)dm), &dmc[i]));
1403ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
141b117cd09SVijay Mahadevan }
142b117cd09SVijay Mahadevan
1432417220eSVijay Mahadevan PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM, PetscInt *, PetscInt *, PetscInt *, PetscInt *, PetscBool);
144b117cd09SVijay Mahadevan
145a4e35b19SJacob Faibussowitsch // PetscClangLinter pragma ignore: -fdoc-*
146a4e35b19SJacob Faibussowitsch /*
147e882eb38SVijay Mahadevan DMCreateInterpolation_Moab - Generate the interpolation operators to transform
148e882eb38SVijay Mahadevan operators (matrices, vectors) from parent level to child level as defined by
14920f4b53cSBarry Smith the `DM` inputs provided by the user.
150b117cd09SVijay Mahadevan
151d083f849SBarry Smith Collective
152b117cd09SVijay Mahadevan
153d8d19677SJose E. Roman Input Parameters:
15460225df5SJacob Faibussowitsch + dmp - The `DMMOAB` object
15560225df5SJacob Faibussowitsch - dmc - the second, finer `DMMOAB` object
156b117cd09SVijay Mahadevan
157d8d19677SJose E. Roman Output Parameters:
158e882eb38SVijay Mahadevan + interpl - The interpolation operator for transferring data between the levels
159e882eb38SVijay Mahadevan - vec - The scaling vector (optional)
160b117cd09SVijay Mahadevan
161e882eb38SVijay Mahadevan Level: developer
162a4e35b19SJacob Faibussowitsch */
DMCreateInterpolation_Moab(DM dmp,DM dmc,Mat * interpl,Vec * vec)163d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmp, DM dmc, Mat *interpl, Vec *vec)
164d71ae5a4SJacob Faibussowitsch {
165755f3dfbSVijay Mahadevan DM_Moab *dmbp, *dmbc;
166b117cd09SVijay Mahadevan moab::ErrorCode merr;
167e882eb38SVijay Mahadevan PetscInt dim;
1683f1c6e43SVijay Mahadevan PetscReal factor;
169ce27a4eeSVijay Mahadevan PetscInt innz, *nnz, ionz, *onz;
170755f3dfbSVijay Mahadevan PetscInt nlsizp, nlsizc, nlghsizp, ngsizp, ngsizc;
171a86ed394SVijay Mahadevan const PetscBool use_consistent_bases = PETSC_TRUE;
172b117cd09SVijay Mahadevan
173b117cd09SVijay Mahadevan PetscFunctionBegin;
174755f3dfbSVijay Mahadevan PetscValidHeaderSpecific(dmp, DM_CLASSID, 1);
175755f3dfbSVijay Mahadevan PetscValidHeaderSpecific(dmc, DM_CLASSID, 2);
176*57508eceSPierre Jolivet dmbp = (DM_Moab *)dmp->data;
177*57508eceSPierre Jolivet dmbc = (DM_Moab *)dmc->data;
178755f3dfbSVijay Mahadevan nlsizp = dmbp->nloc; // *dmb1->numFields;
179755f3dfbSVijay Mahadevan nlsizc = dmbc->nloc; // *dmb2->numFields;
180755f3dfbSVijay Mahadevan ngsizp = dmbp->n; // *dmb1->numFields;
181755f3dfbSVijay Mahadevan ngsizc = dmbc->n; // *dmb2->numFields;
182755f3dfbSVijay Mahadevan nlghsizp = (dmbp->nloc + dmbp->nghost); // *dmb1->numFields;
183b117cd09SVijay Mahadevan
1842417220eSVijay Mahadevan // Columns = Parent DoFs ; Rows = Child DoFs
185755f3dfbSVijay Mahadevan // Interpolation matrix: \sum_{i=1}^P Owned(Child) * (Owned(Parent) + Ghosted(Parent))
186755f3dfbSVijay Mahadevan // Size: nlsizc * nlghsizp
1873ba16761SJacob 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));
188b117cd09SVijay Mahadevan
1899566063dSJacob Faibussowitsch PetscCall(DMGetDimension(dmp, &dim));
190a86ed394SVijay Mahadevan
191941e0cffSVijay Mahadevan /* allocate the nnz, onz arrays based on block size and local nodes */
1929566063dSJacob Faibussowitsch PetscCall(PetscCalloc2(nlsizc, &nnz, nlsizc, &onz));
193941e0cffSVijay Mahadevan
194941e0cffSVijay Mahadevan /* Loop through the local elements and compute the relation between the current parent and the refined_level. */
1952417220eSVijay Mahadevan for (moab::Range::iterator iter = dmbc->vowned->begin(); iter != dmbc->vowned->end(); iter++) {
1962417220eSVijay Mahadevan const moab::EntityHandle vhandle = *iter;
1972417220eSVijay Mahadevan /* define local variables */
1982417220eSVijay Mahadevan moab::EntityHandle parent;
1992417220eSVijay Mahadevan std::vector<moab::EntityHandle> adjs;
2002417220eSVijay Mahadevan moab::Range found;
2012417220eSVijay Mahadevan
2022417220eSVijay Mahadevan /* store the vertex DoF number */
2032417220eSVijay Mahadevan const int ldof = dmbc->lidmap[vhandle - dmbc->seqstart];
2042417220eSVijay Mahadevan
2052417220eSVijay Mahadevan /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects
2062417220eSVijay Mahadevan to the current vertex. We can then decipher if a vertex is ghosted or not and compute the
2072417220eSVijay Mahadevan non-zero pattern accordingly. */
2089371c9d4SSatish Balay merr = dmbc->hierarchy->get_adjacencies(vhandle, dmbc->dim, adjs);
2099371c9d4SSatish Balay MBERRNM(merr);
2102417220eSVijay Mahadevan
2112417220eSVijay Mahadevan /* loop over vertices and update the number of connectivity */
2122417220eSVijay Mahadevan for (unsigned jter = 0; jter < adjs.size(); jter++) {
2132417220eSVijay Mahadevan const moab::EntityHandle jhandle = adjs[jter];
214941e0cffSVijay Mahadevan
215941e0cffSVijay Mahadevan /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */
2169371c9d4SSatish Balay merr = dmbc->hierarchy->child_to_parent(jhandle, dmbc->hlevel, dmbp->hlevel, &parent);
2179371c9d4SSatish Balay MBERRNM(merr);
218941e0cffSVijay Mahadevan
2192417220eSVijay Mahadevan /* Get connectivity information in canonical ordering for the local element */
2202417220eSVijay Mahadevan std::vector<moab::EntityHandle> connp;
2219371c9d4SSatish Balay merr = dmbp->hierarchy->get_connectivity(parent, dmbp->hlevel, connp);
2229371c9d4SSatish Balay MBERRNM(merr);
223a86ed394SVijay Mahadevan
2242417220eSVijay Mahadevan for (unsigned ic = 0; ic < connp.size(); ++ic) {
2252417220eSVijay Mahadevan /* loop over each element connected to the adjacent vertex and update as needed */
2262417220eSVijay Mahadevan /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */
2272417220eSVijay Mahadevan if (found.find(connp[ic]) != found.end()) continue; /* make sure we don't double count shared vertices */
2282417220eSVijay Mahadevan if (dmbp->vghost->find(connp[ic]) != dmbp->vghost->end()) onz[ldof]++; /* update out-of-proc onz */
2292417220eSVijay Mahadevan else nnz[ldof]++; /* else local vertex */
2302417220eSVijay Mahadevan found.insert(connp[ic]);
2312417220eSVijay Mahadevan }
2322417220eSVijay Mahadevan }
233941e0cffSVijay Mahadevan }
234941e0cffSVijay Mahadevan
2359371c9d4SSatish Balay for (int i = 0; i < nlsizc; i++) nnz[i] += 1; /* self count the node */
236941e0cffSVijay Mahadevan
237941e0cffSVijay Mahadevan ionz = onz[0];
238941e0cffSVijay Mahadevan innz = nnz[0];
239755f3dfbSVijay Mahadevan for (int tc = 0; tc < nlsizc; tc++) {
240941e0cffSVijay Mahadevan // check for maximum allowed sparsity = fully dense
241755f3dfbSVijay Mahadevan nnz[tc] = std::min(nlsizp, nnz[tc]);
2422417220eSVijay Mahadevan onz[tc] = std::min(ngsizp - nlsizp, onz[tc]);
2432417220eSVijay Mahadevan
2443ba16761SJacob Faibussowitsch PetscCall(PetscInfo(NULL, " %d: NNZ = %d, ONZ = %d\n", tc, nnz[tc], onz[tc]));
245941e0cffSVijay Mahadevan
246941e0cffSVijay Mahadevan innz = (innz < nnz[tc] ? nnz[tc] : innz);
247941e0cffSVijay Mahadevan ionz = (ionz < onz[tc] ? onz[tc] : ionz);
248941e0cffSVijay Mahadevan }
24964e1c140SVijay Mahadevan
25064e1c140SVijay Mahadevan /* create interpolation matrix */
2519566063dSJacob Faibussowitsch PetscCall(MatCreate(PetscObjectComm((PetscObject)dmc), interpl));
2529566063dSJacob Faibussowitsch PetscCall(MatSetSizes(*interpl, nlsizc, nlsizp, ngsizc, ngsizp));
2539566063dSJacob Faibussowitsch PetscCall(MatSetType(*interpl, MATAIJ));
2549566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(*interpl));
25564e1c140SVijay Mahadevan
2569566063dSJacob Faibussowitsch PetscCall(MatSeqAIJSetPreallocation(*interpl, innz, nnz));
2579566063dSJacob Faibussowitsch PetscCall(MatMPIAIJSetPreallocation(*interpl, innz, nnz, ionz, onz));
258941e0cffSVijay Mahadevan
259941e0cffSVijay Mahadevan /* clean up temporary memory */
2609566063dSJacob Faibussowitsch PetscCall(PetscFree2(nnz, onz));
261b117cd09SVijay Mahadevan
262b117cd09SVijay Mahadevan /* set up internal matrix data-structures */
2639566063dSJacob Faibussowitsch PetscCall(MatSetUp(*interpl));
264b117cd09SVijay Mahadevan
265a86ed394SVijay Mahadevan /* Define variables for assembly */
266a86ed394SVijay Mahadevan std::vector<moab::EntityHandle> children;
267a86ed394SVijay Mahadevan std::vector<moab::EntityHandle> connp, connc;
268a86ed394SVijay Mahadevan std::vector<PetscReal> pcoords, ccoords, values_phi;
269ca4445c7SIlya Fursov std::vector<PetscScalar> values_phi_scalar;
270a86ed394SVijay Mahadevan
271a86ed394SVijay Mahadevan if (use_consistent_bases) {
2722417220eSVijay Mahadevan const moab::EntityHandle ehandle = dmbp->elocal->front();
273a86ed394SVijay Mahadevan
2749371c9d4SSatish Balay merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children);
2759371c9d4SSatish Balay MBERRNM(merr);
276a86ed394SVijay Mahadevan
277a86ed394SVijay Mahadevan /* Get connectivity and coordinates of the parent vertices */
2789371c9d4SSatish Balay merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp);
2799371c9d4SSatish Balay MBERRNM(merr);
2809371c9d4SSatish Balay merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc);
2819371c9d4SSatish Balay MBERRNM(merr);
282a86ed394SVijay Mahadevan
283a86ed394SVijay Mahadevan std::vector<PetscReal> natparam(3 * connc.size(), 0.0);
284a86ed394SVijay Mahadevan pcoords.resize(connp.size() * 3);
285a86ed394SVijay Mahadevan ccoords.resize(connc.size() * 3);
286a86ed394SVijay Mahadevan values_phi.resize(connp.size() * connc.size());
287a86ed394SVijay Mahadevan /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */
2889371c9d4SSatish Balay merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]);
2899371c9d4SSatish Balay MBERRNM(merr);
2909371c9d4SSatish Balay merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]);
2919371c9d4SSatish Balay MBERRNM(merr);
292a86ed394SVijay Mahadevan
293a86ed394SVijay Mahadevan /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */
294a86ed394SVijay Mahadevan for (unsigned tc = 0; tc < connc.size(); tc++) {
295a86ed394SVijay Mahadevan const PetscInt offset = tc * 3;
296a86ed394SVijay Mahadevan
297a86ed394SVijay Mahadevan /* Scale ccoords relative to pcoords */
2989566063dSJacob Faibussowitsch PetscCall(DMMoabPToRMapping(dim, connp.size(), &pcoords[0], &ccoords[offset], &natparam[offset], &values_phi[connp.size() * tc]));
299a86ed394SVijay Mahadevan }
3001baa6e33SBarry Smith } else {
301a86ed394SVijay Mahadevan factor = std::pow(2.0 /*degree_P_for_refinement*/, (dmbc->hlevel - dmbp->hlevel) * dmbp->dim * 1.0);
302a86ed394SVijay Mahadevan }
303a86ed394SVijay Mahadevan
304755f3dfbSVijay Mahadevan /* TODO: Decipher the correct non-zero pattern. There is still some issue with onz allocation */
3059566063dSJacob Faibussowitsch PetscCall(MatSetOption(*interpl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
306b117cd09SVijay Mahadevan
307e882eb38SVijay Mahadevan /* Loop through the remaining vertices. These vertices appear only on the current refined_level. */
308ca4445c7SIlya Fursov values_phi_scalar = VecReal_to_VecScalar(values_phi);
3092417220eSVijay Mahadevan for (moab::Range::iterator iter = dmbp->elocal->begin(); iter != dmbp->elocal->end(); iter++) {
310b117cd09SVijay Mahadevan const moab::EntityHandle ehandle = *iter;
311b117cd09SVijay Mahadevan
312b117cd09SVijay Mahadevan /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */
313a86ed394SVijay Mahadevan children.clear();
314a86ed394SVijay Mahadevan connc.clear();
3159371c9d4SSatish Balay merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children);
3169371c9d4SSatish Balay MBERRNM(merr);
317b117cd09SVijay Mahadevan
318b117cd09SVijay Mahadevan /* Get connectivity and coordinates of the parent vertices */
3199371c9d4SSatish Balay merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp);
3209371c9d4SSatish Balay MBERRNM(merr);
3219371c9d4SSatish Balay merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc);
3229371c9d4SSatish Balay MBERRNM(merr);
323b117cd09SVijay Mahadevan
324a86ed394SVijay Mahadevan pcoords.resize(connp.size() * 3);
325a86ed394SVijay Mahadevan ccoords.resize(connc.size() * 3);
326b117cd09SVijay Mahadevan /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */
3279371c9d4SSatish Balay merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]);
3289371c9d4SSatish Balay MBERRNM(merr);
3299371c9d4SSatish Balay merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]);
3309371c9d4SSatish Balay MBERRNM(merr);
331b117cd09SVijay Mahadevan
332b117cd09SVijay Mahadevan std::vector<int> dofsp(connp.size()), dofsc(connc.size());
333b117cd09SVijay Mahadevan /* TODO: specific to scalar system - use GetDofs */
3349566063dSJacob Faibussowitsch PetscCall(DMMoabGetDofsBlocked(dmp, connp.size(), &connp[0], &dofsp[0]));
3359566063dSJacob Faibussowitsch PetscCall(DMMoabGetDofsBlocked(dmc, connc.size(), &connc[0], &dofsc[0]));
336b117cd09SVijay Mahadevan
337a86ed394SVijay Mahadevan /* Compute the actual interpolation weights when projecting solution/residual between levels */
338a86ed394SVijay Mahadevan if (use_consistent_bases) {
339145b44c9SPierre Jolivet /* Use the cached values of natural parametric coordinates and basis pre-evaluated.
340a86ed394SVijay Mahadevan We are making an assumption here that UMR used in GMG to generate the hierarchy uses
341a86ed394SVijay Mahadevan the same template for all elements; This will fail for mixed element meshes (TRI/QUAD).
342a86ed394SVijay Mahadevan
3432417220eSVijay Mahadevan TODO: Fix the above assumption by caching data for families (especially for Tets and mixed meshes)
344a86ed394SVijay Mahadevan */
345a86ed394SVijay Mahadevan
346a86ed394SVijay Mahadevan /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */
347a86ed394SVijay Mahadevan for (unsigned tc = 0; tc < connc.size(); tc++) {
348a86ed394SVijay Mahadevan /* TODO: Check if we should be using INSERT_VALUES instead */
349ca4445c7SIlya Fursov PetscCall(MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi_scalar[connp.size() * tc], ADD_VALUES));
350a86ed394SVijay Mahadevan }
3511baa6e33SBarry Smith } else {
352e882eb38SVijay Mahadevan /* Compute the interpolation weights by determining distance of 1-ring
353755f3dfbSVijay Mahadevan neighbor vertices from current vertex
354755f3dfbSVijay Mahadevan
355a86ed394SVijay Mahadevan This should be used only when FEM basis is not used for the discretization.
356a86ed394SVijay Mahadevan Else, the consistent interface to compute the basis function for interpolation
357a86ed394SVijay Mahadevan between the levels should be evaluated correctly to preserve convergence of GMG.
3582417220eSVijay Mahadevan Shephard's basis will be terrible for any unsmooth problems.
359755f3dfbSVijay Mahadevan */
360755f3dfbSVijay Mahadevan values_phi.resize(connp.size());
361755f3dfbSVijay Mahadevan for (unsigned tc = 0; tc < connc.size(); tc++) {
362a86ed394SVijay Mahadevan PetscReal normsum = 0.0;
363ca4445c7SIlya Fursov std::vector<PetscScalar> values_phi_scalar2;
364755f3dfbSVijay Mahadevan for (unsigned tp = 0; tp < connp.size(); tp++) {
365755f3dfbSVijay Mahadevan values_phi[tp] = 0.0;
3669371c9d4SSatish Balay for (unsigned k = 0; k < 3; k++) values_phi[tp] += std::pow(pcoords[tp * 3 + k] - ccoords[k + tc * 3], dim);
367755f3dfbSVijay Mahadevan if (values_phi[tp] < 1e-12) {
368755f3dfbSVijay Mahadevan values_phi[tp] = 1e12;
3691baa6e33SBarry Smith } else {
370755f3dfbSVijay Mahadevan //values_phi[tp] = std::pow(values_phi[tp], -1.0/dim);
371755f3dfbSVijay Mahadevan values_phi[tp] = std::pow(values_phi[tp], -1.0);
372755f3dfbSVijay Mahadevan normsum += values_phi[tp];
373b117cd09SVijay Mahadevan }
374b117cd09SVijay Mahadevan }
375755f3dfbSVijay Mahadevan for (unsigned tp = 0; tp < connp.size(); tp++) {
3769371c9d4SSatish Balay if (values_phi[tp] > 1e11) values_phi[tp] = factor * 0.5 / connp.size();
3779371c9d4SSatish Balay else values_phi[tp] = factor * values_phi[tp] * 0.5 / (connp.size() * normsum);
378b117cd09SVijay Mahadevan }
379ca4445c7SIlya Fursov values_phi_scalar2 = VecReal_to_VecScalar(values_phi);
380ca4445c7SIlya Fursov PetscCall(MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi_scalar2[0], ADD_VALUES));
381b117cd09SVijay Mahadevan }
382a86ed394SVijay Mahadevan }
383b117cd09SVijay Mahadevan }
384ca4445c7SIlya Fursov if (vec) *vec = NULL; // TODO: <-- is it safe/appropriate?
3859566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(*interpl, MAT_FINAL_ASSEMBLY));
3869566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(*interpl, MAT_FINAL_ASSEMBLY));
3873ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
388b117cd09SVijay Mahadevan }
389b117cd09SVijay Mahadevan
390a4e35b19SJacob Faibussowitsch // PetscClangLinter pragma ignore: -fdoc-*
391a4e35b19SJacob Faibussowitsch /*
392b117cd09SVijay Mahadevan DMCreateInjection_Moab - Generate a multi-level uniform refinement hierarchy
39320f4b53cSBarry Smith by succesively refining a coarse mesh, already defined in the `DM` object
394b117cd09SVijay Mahadevan provided by the user.
395b117cd09SVijay Mahadevan
396d083f849SBarry Smith Collective
397b117cd09SVijay Mahadevan
398b117cd09SVijay Mahadevan Input Parameter:
39920f4b53cSBarry Smith . dmb - The `DMMOAB` object
400b117cd09SVijay Mahadevan
401d8d19677SJose E. Roman Output Parameters:
402a2b725a8SWilliam Gropp + nlevels - The number of levels of refinement needed to generate the hierarchy
403a2b725a8SWilliam Gropp - ldegrees - The degree of refinement at each level in the hierarchy
404b117cd09SVijay Mahadevan
405b117cd09SVijay Mahadevan Level: beginner
406a4e35b19SJacob Faibussowitsch */
DMCreateInjection_Moab(DM dm1,DM dm2,VecScatter * ctx)407d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dm1, DM dm2, VecScatter *ctx)
408d71ae5a4SJacob Faibussowitsch {
409b117cd09SVijay Mahadevan PetscFunctionBegin;
410b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm1, DM_CLASSID, 1);
411b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm2, DM_CLASSID, 2);
4123ba16761SJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[DMCreateInjection_Moab] :: Placeholder\n"));
4133ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
414b117cd09SVijay Mahadevan }
415b117cd09SVijay Mahadevan
DMMoab_UMR_Private(DM dm,MPI_Comm comm,PetscBool refine,DM * dmref)416d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMMoab_UMR_Private(DM dm, MPI_Comm comm, PetscBool refine, DM *dmref)
417d71ae5a4SJacob Faibussowitsch {
418e882eb38SVijay Mahadevan PetscInt i, dim;
419b117cd09SVijay Mahadevan DM dm2;
420b117cd09SVijay Mahadevan moab::ErrorCode merr;
421b117cd09SVijay Mahadevan DM_Moab *dmb = (DM_Moab *)dm->data, *dd2;
422b117cd09SVijay Mahadevan
423b117cd09SVijay Mahadevan PetscFunctionBegin;
424b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
4254f572ea9SToby Isaac PetscAssertPointer(dmref, 4);
426b117cd09SVijay Mahadevan
427b117cd09SVijay Mahadevan if ((dmb->hlevel == dmb->nhlevels && refine) || (dmb->hlevel == 0 && !refine)) {
4283ba16761SJacob Faibussowitsch if (dmb->hlevel + 1 > dmb->nhlevels && refine) {
4293ba16761SJacob 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));
4303ba16761SJacob Faibussowitsch }
4313ba16761SJacob 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));
432f3fa974cSJacob Faibussowitsch *dmref = NULL;
4333ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
434b117cd09SVijay Mahadevan }
435b117cd09SVijay Mahadevan
4369566063dSJacob Faibussowitsch PetscCall(DMMoabCreate(PetscObjectComm((PetscObject)dm), &dm2));
437b117cd09SVijay Mahadevan dd2 = (DM_Moab *)dm2->data;
438b117cd09SVijay Mahadevan
439b117cd09SVijay Mahadevan dd2->mbiface = dmb->mbiface;
4409daf19fdSVijay Mahadevan #ifdef MOAB_HAVE_MPI
441b117cd09SVijay Mahadevan dd2->pcomm = dmb->pcomm;
4429daf19fdSVijay Mahadevan #endif
443b117cd09SVijay Mahadevan dd2->icreatedinstance = PETSC_FALSE;
44464e1c140SVijay Mahadevan dd2->nghostrings = dmb->nghostrings;
445b117cd09SVijay Mahadevan
446b117cd09SVijay Mahadevan /* set the new level based on refinement/coarsening */
447b117cd09SVijay Mahadevan if (refine) {
448b117cd09SVijay Mahadevan dd2->hlevel = dmb->hlevel + 1;
4491baa6e33SBarry Smith } else {
450b117cd09SVijay Mahadevan dd2->hlevel = dmb->hlevel - 1;
451b117cd09SVijay Mahadevan }
452b117cd09SVijay Mahadevan
453b117cd09SVijay Mahadevan /* Copy the multilevel hierarchy pointers in MOAB */
454b117cd09SVijay Mahadevan dd2->hierarchy = dmb->hierarchy;
455b117cd09SVijay Mahadevan dd2->nhlevels = dmb->nhlevels;
4569566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(dd2->nhlevels + 1, &dd2->hsets));
457ad540459SPierre Jolivet for (i = 0; i <= dd2->nhlevels; i++) dd2->hsets[i] = dmb->hsets[i];
458b117cd09SVijay Mahadevan dd2->fileset = dd2->hsets[dd2->hlevel];
459b117cd09SVijay Mahadevan
460b117cd09SVijay Mahadevan /* do the remaining initializations for DMMoab */
461b117cd09SVijay Mahadevan dd2->bs = dmb->bs;
462b117cd09SVijay Mahadevan dd2->numFields = dmb->numFields;
463b117cd09SVijay Mahadevan dd2->rw_dbglevel = dmb->rw_dbglevel;
464b117cd09SVijay Mahadevan dd2->partition_by_rank = dmb->partition_by_rank;
465c6a7a370SJeremy L Thompson PetscCall(PetscStrncpy(dd2->extra_read_options, dmb->extra_read_options, sizeof(dd2->extra_read_options)));
466c6a7a370SJeremy L Thompson PetscCall(PetscStrncpy(dd2->extra_write_options, dmb->extra_write_options, sizeof(dd2->extra_write_options)));
467b117cd09SVijay Mahadevan dd2->read_mode = dmb->read_mode;
468b117cd09SVijay Mahadevan dd2->write_mode = dmb->write_mode;
469b117cd09SVijay Mahadevan
470b117cd09SVijay Mahadevan /* set global ID tag handle */
4719566063dSJacob Faibussowitsch PetscCall(DMMoabSetLocalToGlobalTag(dm2, dmb->ltog_tag));
472b117cd09SVijay Mahadevan
4739371c9d4SSatish Balay merr = dd2->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dd2->material_tag);
4749371c9d4SSatish Balay MBERRNM(merr);
475b117cd09SVijay Mahadevan
4769566063dSJacob Faibussowitsch PetscCall(DMSetOptionsPrefix(dm2, ((PetscObject)dm)->prefix));
4779566063dSJacob Faibussowitsch PetscCall(DMGetDimension(dm, &dim));
4789566063dSJacob Faibussowitsch PetscCall(DMSetDimension(dm2, dim));
479b117cd09SVijay Mahadevan
480b117cd09SVijay Mahadevan /* allow overloaded (user replaced) operations to be inherited by refinement clones */
481b117cd09SVijay Mahadevan dm2->ops->creatematrix = dm->ops->creatematrix;
482b117cd09SVijay Mahadevan
483b117cd09SVijay Mahadevan /* copy fill information if given */
4849566063dSJacob Faibussowitsch PetscCall(DMMoabSetBlockFills(dm2, dmb->dfill, dmb->ofill));
485b117cd09SVijay Mahadevan
486b117cd09SVijay Mahadevan /* copy vector type information */
4879566063dSJacob Faibussowitsch PetscCall(DMSetMatType(dm2, dm->mattype));
4889566063dSJacob Faibussowitsch PetscCall(DMSetVecType(dm2, dm->vectype));
4893f1c6e43SVijay Mahadevan dd2->numFields = dmb->numFields;
49048a46eb9SPierre Jolivet if (dmb->numFields) PetscCall(DMMoabSetFieldNames(dm2, dmb->numFields, dmb->fieldNames));
491b117cd09SVijay Mahadevan
4929566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(dm2));
493b117cd09SVijay Mahadevan
494b117cd09SVijay Mahadevan /* recreate Dof numbering for the refined DM and make sure the distribution is correctly populated */
4959566063dSJacob Faibussowitsch PetscCall(DMSetUp(dm2));
496b117cd09SVijay Mahadevan
497b117cd09SVijay Mahadevan *dmref = dm2;
4983ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
499b117cd09SVijay Mahadevan }
500b117cd09SVijay Mahadevan
501a4e35b19SJacob Faibussowitsch // PetscClangLinter pragma ignore: -fdoc-*
502a4e35b19SJacob Faibussowitsch /*
503b117cd09SVijay Mahadevan DMRefine_Moab - Generate a multi-level uniform refinement hierarchy
50420f4b53cSBarry Smith by succesively refining a coarse mesh, already defined in the `DM` object
505b117cd09SVijay Mahadevan provided by the user.
506b117cd09SVijay Mahadevan
50720f4b53cSBarry Smith Collective
508b117cd09SVijay Mahadevan
509d8d19677SJose E. Roman Input Parameters:
51020f4b53cSBarry Smith + dm - The `DMMOAB` object
51120f4b53cSBarry Smith - comm - the communicator to contain the new DM object (or `MPI_COMM_NULL`)
512b117cd09SVijay Mahadevan
513b117cd09SVijay Mahadevan Output Parameter:
51420f4b53cSBarry Smith . dmf - the refined `DM`, or `NULL`
515e882eb38SVijay Mahadevan
516e882eb38SVijay Mahadevan Level: developer
517b117cd09SVijay Mahadevan
51820f4b53cSBarry Smith Note:
51920f4b53cSBarry Smith If no refinement was done, the return value is `NULL`
520a4e35b19SJacob Faibussowitsch */
DMRefine_Moab(DM dm,MPI_Comm comm,DM * dmf)521d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM *dmf)
522d71ae5a4SJacob Faibussowitsch {
523b117cd09SVijay Mahadevan PetscFunctionBegin;
524b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
525b117cd09SVijay Mahadevan
5269566063dSJacob Faibussowitsch PetscCall(DMMoab_UMR_Private(dm, comm, PETSC_TRUE, dmf));
5273ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
528b117cd09SVijay Mahadevan }
529b117cd09SVijay Mahadevan
530a4e35b19SJacob Faibussowitsch // PetscClangLinter pragma ignore: -fdoc-*
531a4e35b19SJacob Faibussowitsch /*
532b117cd09SVijay Mahadevan DMCoarsen_Moab - Generate a multi-level uniform refinement hierarchy
53320f4b53cSBarry Smith by succesively refining a coarse mesh, already defined in the `DM` object
534b117cd09SVijay Mahadevan provided by the user.
535b117cd09SVijay Mahadevan
53620f4b53cSBarry Smith Collective
537b117cd09SVijay Mahadevan
538d8d19677SJose E. Roman Input Parameters:
53920f4b53cSBarry Smith + dm - The `DMMOAB` object
54020f4b53cSBarry Smith - comm - the communicator to contain the new `DM` object (or `MPI_COMM_NULL`)
541b117cd09SVijay Mahadevan
542b117cd09SVijay Mahadevan Output Parameter:
54360225df5SJacob Faibussowitsch . dmc - the coarsened `DM`, or `NULL`
544b117cd09SVijay Mahadevan
545e882eb38SVijay Mahadevan Level: developer
546e882eb38SVijay Mahadevan
54720f4b53cSBarry Smith Note:
54820f4b53cSBarry Smith If no coarsening was done, the return value is `NULL`
549a4e35b19SJacob Faibussowitsch */
DMCoarsen_Moab(DM dm,MPI_Comm comm,DM * dmc)550d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM *dmc)
551d71ae5a4SJacob Faibussowitsch {
552b117cd09SVijay Mahadevan PetscFunctionBegin;
553b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
5549566063dSJacob Faibussowitsch PetscCall(DMMoab_UMR_Private(dm, comm, PETSC_FALSE, dmc));
5553ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
556b117cd09SVijay Mahadevan }
557