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