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 6b117cd09SVijay Mahadevan #undef __FUNCT__ 7b117cd09SVijay Mahadevan #define __FUNCT__ "DMMoabGenerateHierarchy" 8b117cd09SVijay Mahadevan /*@ 9b117cd09SVijay Mahadevan DMMoabGenerateHierarchy - Generate a multi-level uniform refinement hierarchy 10b117cd09SVijay Mahadevan by succesively refining a coarse mesh, already defined in the DM object 11b117cd09SVijay Mahadevan provided by the user. 12b117cd09SVijay Mahadevan 13b117cd09SVijay Mahadevan Collective on MPI_Comm 14b117cd09SVijay Mahadevan 15b117cd09SVijay Mahadevan Input Parameter: 16b117cd09SVijay Mahadevan + dmb - The DMMoab object 17b117cd09SVijay Mahadevan 18b117cd09SVijay Mahadevan Output Parameter: 19b117cd09SVijay Mahadevan + nlevels - The number of levels of refinement needed to generate the hierarchy 20b117cd09SVijay Mahadevan . ldegrees - The degree of refinement at each level in the hierarchy 21b117cd09SVijay Mahadevan 22b117cd09SVijay Mahadevan Level: beginner 23b117cd09SVijay Mahadevan 24b117cd09SVijay Mahadevan .keywords: DMMoab, create, refinement 25b117cd09SVijay Mahadevan @*/ 26b117cd09SVijay Mahadevan PetscErrorCode DMMoabGenerateHierarchy(DM dm, PetscInt nlevels, PetscInt *ldegrees) 27b117cd09SVijay Mahadevan { 28b117cd09SVijay Mahadevan DM_Moab *dmmoab; 29b117cd09SVijay Mahadevan PetscErrorCode ierr; 30b117cd09SVijay Mahadevan moab::ErrorCode merr; 312417220eSVijay Mahadevan PetscInt *pdegrees, ilevel; 32e882eb38SVijay Mahadevan std::vector<moab::EntityHandle> hsets; 33b117cd09SVijay Mahadevan 34b117cd09SVijay Mahadevan PetscFunctionBegin; 35b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 36b117cd09SVijay Mahadevan dmmoab = (DM_Moab*)(dm)->data; 37b117cd09SVijay Mahadevan 38b117cd09SVijay Mahadevan if (!ldegrees) { 39b117cd09SVijay Mahadevan ierr = PetscMalloc1(nlevels, &pdegrees);CHKERRQ(ierr); 402417220eSVijay Mahadevan for (ilevel = 0; ilevel < nlevels; ilevel++) pdegrees[ilevel] = 2; /* default = Degree 2 refinement */ 41b117cd09SVijay Mahadevan } 42b117cd09SVijay Mahadevan else pdegrees = ldegrees; 43b117cd09SVijay Mahadevan 44b117cd09SVijay Mahadevan /* initialize set level refinement data for hierarchy */ 45b117cd09SVijay Mahadevan dmmoab->nhlevels = nlevels; 46b117cd09SVijay Mahadevan 47b117cd09SVijay Mahadevan /* Instantiate the nested refinement class */ 489daf19fdSVijay Mahadevan #ifdef MOAB_HAVE_MPI 493f1c6e43SVijay Mahadevan dmmoab->hierarchy = new moab::NestedRefine(dynamic_cast<moab::Core*>(dmmoab->mbiface), dmmoab->pcomm, dmmoab->fileset); 509daf19fdSVijay Mahadevan #else 519daf19fdSVijay Mahadevan dmmoab->hierarchy = new moab::NestedRefine(dynamic_cast<moab::Core*>(dmmoab->mbiface), NULL, dmmoab->fileset); 529daf19fdSVijay Mahadevan #endif 53b117cd09SVijay Mahadevan 54e882eb38SVijay Mahadevan ierr = PetscMalloc1(nlevels + 1, &dmmoab->hsets);CHKERRQ(ierr); 55e882eb38SVijay Mahadevan 56b117cd09SVijay Mahadevan /* generate the mesh hierarchy */ 57*9c368985SVijay Mahadevan merr = dmmoab->hierarchy->generate_mesh_hierarchy(nlevels, pdegrees, hsets, false);MBERRNM(merr); 58e882eb38SVijay Mahadevan 599daf19fdSVijay Mahadevan #ifdef MOAB_HAVE_MPI 60755f3dfbSVijay Mahadevan if (dmmoab->pcomm->size() > 1) { 6164e1c140SVijay Mahadevan merr = dmmoab->hierarchy->exchange_ghosts(hsets, dmmoab->nghostrings);MBERRNM(merr); 62755f3dfbSVijay Mahadevan } 639daf19fdSVijay Mahadevan #endif 6449d66b22SVijay Mahadevan 6564e1c140SVijay Mahadevan /* copy the mesh sets for nested refinement hierarchy */ 662417220eSVijay Mahadevan for (ilevel = 0; ilevel <= nlevels; ilevel++) 672417220eSVijay Mahadevan { 682417220eSVijay Mahadevan if (ilevel) { /* Update material and other geometric tags from parent to child sets */ 692417220eSVijay Mahadevan merr = dmmoab->hierarchy->update_special_tags(ilevel, hsets[ilevel]);MBERRNM(merr); 702417220eSVijay Mahadevan } 712417220eSVijay Mahadevan 722417220eSVijay Mahadevan dmmoab->hsets[ilevel] = hsets[ilevel]; 73*9c368985SVijay Mahadevan #ifdef MOAB_HAVE_MPI 74*9c368985SVijay Mahadevan if (ilevel) { 75*9c368985SVijay Mahadevan merr = dmmoab->pcomm->assign_global_ids(hsets[ilevel], dmmoab->dim, 1, false, true, false);MBERRNM(merr); 76*9c368985SVijay Mahadevan } 77*9c368985SVijay Mahadevan #endif 782417220eSVijay Mahadevan } 7964e1c140SVijay Mahadevan 8064e1c140SVijay Mahadevan hsets.clear(); 81b117cd09SVijay Mahadevan if (!ldegrees) { 82b117cd09SVijay Mahadevan ierr = PetscFree(pdegrees);CHKERRQ(ierr); 83b117cd09SVijay Mahadevan } 84b117cd09SVijay Mahadevan PetscFunctionReturn(0); 85b117cd09SVijay Mahadevan } 86b117cd09SVijay Mahadevan 87b117cd09SVijay Mahadevan #undef __FUNCT__ 88b117cd09SVijay Mahadevan #define __FUNCT__ "DMRefineHierarchy_Moab" 89b117cd09SVijay Mahadevan /*@ 90e882eb38SVijay Mahadevan DMRefineHierarchy_Moab - Generate a multi-level DM hierarchy 91e882eb38SVijay Mahadevan by succesively refining a coarse mesh. 92b117cd09SVijay Mahadevan 93b117cd09SVijay Mahadevan Collective on MPI_Comm 94b117cd09SVijay Mahadevan 95b117cd09SVijay Mahadevan Input Parameter: 96e882eb38SVijay Mahadevan + dm - The DMMoab object 97b117cd09SVijay Mahadevan 98b117cd09SVijay Mahadevan Output Parameter: 99e882eb38SVijay Mahadevan + nlevels - The number of levels of refinement needed to generate the hierarchy 100e882eb38SVijay Mahadevan . dmf - The DM objects after successive refinement of the hierarchy 101b117cd09SVijay Mahadevan 102b117cd09SVijay Mahadevan Level: beginner 103b117cd09SVijay Mahadevan 104e882eb38SVijay Mahadevan .keywords: DMMoab, generate, refinement 105b117cd09SVijay Mahadevan @*/ 106304006b3SVijay Mahadevan PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmf[]) 107b117cd09SVijay Mahadevan { 108b117cd09SVijay Mahadevan PetscErrorCode ierr; 109b117cd09SVijay Mahadevan PetscInt i; 110b117cd09SVijay Mahadevan 111b117cd09SVijay Mahadevan PetscFunctionBegin; 112b117cd09SVijay Mahadevan 113e882eb38SVijay Mahadevan ierr = DMRefine(dm, PetscObjectComm((PetscObject)dm), &dmf[0]);CHKERRQ(ierr); 114e882eb38SVijay Mahadevan for (i = 1; i < nlevels; i++) { 115e882eb38SVijay Mahadevan ierr = DMRefine(dmf[i - 1], PetscObjectComm((PetscObject)dm), &dmf[i]);CHKERRQ(ierr); 116b117cd09SVijay Mahadevan } 117b117cd09SVijay Mahadevan PetscFunctionReturn(0); 118b117cd09SVijay Mahadevan } 119b117cd09SVijay Mahadevan 120b117cd09SVijay Mahadevan #undef __FUNCT__ 121b117cd09SVijay Mahadevan #define __FUNCT__ "DMCoarsenHierarchy_Moab" 122b117cd09SVijay Mahadevan /*@ 123e882eb38SVijay Mahadevan DMCoarsenHierarchy_Moab - Generate a multi-level DM hierarchy 124e882eb38SVijay Mahadevan by succesively coarsening a refined mesh. 125b117cd09SVijay Mahadevan 126b117cd09SVijay Mahadevan Collective on MPI_Comm 127b117cd09SVijay Mahadevan 128b117cd09SVijay Mahadevan Input Parameter: 129e882eb38SVijay Mahadevan + dm - The DMMoab object 130b117cd09SVijay Mahadevan 131b117cd09SVijay Mahadevan Output Parameter: 132e882eb38SVijay Mahadevan + nlevels - The number of levels of refinement needed to generate the hierarchy 133e882eb38SVijay Mahadevan . dmc - The DM objects after successive coarsening of the hierarchy 134b117cd09SVijay Mahadevan 135b117cd09SVijay Mahadevan Level: beginner 136b117cd09SVijay Mahadevan 137e882eb38SVijay Mahadevan .keywords: DMMoab, generate, coarsening 138b117cd09SVijay Mahadevan @*/ 139304006b3SVijay Mahadevan PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmc[]) 140b117cd09SVijay Mahadevan { 141b117cd09SVijay Mahadevan PetscErrorCode ierr; 142b117cd09SVijay Mahadevan PetscInt i; 143b117cd09SVijay Mahadevan 144b117cd09SVijay Mahadevan PetscFunctionBegin; 145b117cd09SVijay Mahadevan 146e882eb38SVijay Mahadevan ierr = DMCoarsen(dm, PetscObjectComm((PetscObject)dm), &dmc[0]);CHKERRQ(ierr); 147e882eb38SVijay Mahadevan for (i = 1; i < nlevels; i++) { 148e882eb38SVijay Mahadevan ierr = DMCoarsen(dmc[i - 1], PetscObjectComm((PetscObject)dm), &dmc[i]);CHKERRQ(ierr); 149b117cd09SVijay Mahadevan } 150b117cd09SVijay Mahadevan PetscFunctionReturn(0); 151b117cd09SVijay Mahadevan } 152b117cd09SVijay Mahadevan 1532417220eSVijay Mahadevan PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM, PetscInt*, PetscInt*, PetscInt*, PetscInt*, PetscBool); 154b117cd09SVijay Mahadevan 155b117cd09SVijay Mahadevan #undef __FUNCT__ 156b117cd09SVijay Mahadevan #define __FUNCT__ "DMCreateInterpolation_Moab" 157b117cd09SVijay Mahadevan /*@ 158e882eb38SVijay Mahadevan DMCreateInterpolation_Moab - Generate the interpolation operators to transform 159e882eb38SVijay Mahadevan operators (matrices, vectors) from parent level to child level as defined by 160e882eb38SVijay Mahadevan the DM inputs provided by the user. 161b117cd09SVijay Mahadevan 162b117cd09SVijay Mahadevan Collective on MPI_Comm 163b117cd09SVijay Mahadevan 164b117cd09SVijay Mahadevan Input Parameter: 165e882eb38SVijay Mahadevan + dm1 - The DMMoab object 166e882eb38SVijay Mahadevan - dm2 - the second, finer DMMoab object 167b117cd09SVijay Mahadevan 168b117cd09SVijay Mahadevan Output Parameter: 169e882eb38SVijay Mahadevan + interpl - The interpolation operator for transferring data between the levels 170e882eb38SVijay Mahadevan - vec - The scaling vector (optional) 171b117cd09SVijay Mahadevan 172e882eb38SVijay Mahadevan Level: developer 173b117cd09SVijay Mahadevan 174b117cd09SVijay Mahadevan .keywords: DMMoab, create, refinement 175b117cd09SVijay Mahadevan @*/ 176304006b3SVijay Mahadevan PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmp, DM dmc, Mat* interpl, Vec* vec) 177b117cd09SVijay Mahadevan { 178755f3dfbSVijay Mahadevan DM_Moab *dmbp, *dmbc; 179b117cd09SVijay Mahadevan PetscErrorCode ierr; 180b117cd09SVijay Mahadevan moab::ErrorCode merr; 181e882eb38SVijay Mahadevan PetscInt dim; 1823f1c6e43SVijay Mahadevan PetscReal factor; 183ce27a4eeSVijay Mahadevan PetscInt innz, *nnz, ionz, *onz; 184755f3dfbSVijay Mahadevan PetscInt nlsizp, nlsizc, nlghsizp, ngsizp, ngsizc; 185a86ed394SVijay Mahadevan const PetscBool use_consistent_bases=PETSC_TRUE; 186b117cd09SVijay Mahadevan 187b117cd09SVijay Mahadevan PetscFunctionBegin; 188755f3dfbSVijay Mahadevan PetscValidHeaderSpecific(dmp, DM_CLASSID, 1); 189755f3dfbSVijay Mahadevan PetscValidHeaderSpecific(dmc, DM_CLASSID, 2); 190755f3dfbSVijay Mahadevan dmbp = (DM_Moab*)(dmp)->data; 191755f3dfbSVijay Mahadevan dmbc = (DM_Moab*)(dmc)->data; 192755f3dfbSVijay Mahadevan nlsizp = dmbp->nloc;// *dmb1->numFields; 193755f3dfbSVijay Mahadevan nlsizc = dmbc->nloc;// *dmb2->numFields; 194755f3dfbSVijay Mahadevan ngsizp = dmbp->n; // *dmb1->numFields; 195755f3dfbSVijay Mahadevan ngsizc = dmbc->n; // *dmb2->numFields; 196755f3dfbSVijay Mahadevan nlghsizp = (dmbp->nloc + dmbp->nghost); // *dmb1->numFields; 197b117cd09SVijay Mahadevan 1982417220eSVijay Mahadevan // Columns = Parent DoFs ; Rows = Child DoFs 199755f3dfbSVijay Mahadevan // Interpolation matrix: \sum_{i=1}^P Owned(Child) * (Owned(Parent) + Ghosted(Parent)) 200755f3dfbSVijay Mahadevan // Size: nlsizc * nlghsizp 201755f3dfbSVijay Mahadevan PetscInfo4(NULL, "Creating interpolation matrix %D X %D to apply transformation between levels %D -> %D.\n", ngsizc, nlghsizp, dmbp->hlevel, dmbc->hlevel); 202b117cd09SVijay Mahadevan 203a86ed394SVijay Mahadevan ierr = DMGetDimension(dmp, &dim);CHKERRQ(ierr); 204a86ed394SVijay Mahadevan 205941e0cffSVijay Mahadevan /* allocate the nnz, onz arrays based on block size and local nodes */ 206755f3dfbSVijay Mahadevan ierr = PetscCalloc2(nlsizc, &nnz, nlsizc, &onz);CHKERRQ(ierr); 207941e0cffSVijay Mahadevan 208941e0cffSVijay Mahadevan /* Loop through the local elements and compute the relation between the current parent and the refined_level. */ 2092417220eSVijay Mahadevan for (moab::Range::iterator iter = dmbc->vowned->begin(); iter != dmbc->vowned->end(); iter++) { 210941e0cffSVijay Mahadevan 2112417220eSVijay Mahadevan const moab::EntityHandle vhandle = *iter; 2122417220eSVijay Mahadevan /* define local variables */ 2132417220eSVijay Mahadevan moab::EntityHandle parent; 2142417220eSVijay Mahadevan std::vector<moab::EntityHandle> adjs; 2152417220eSVijay Mahadevan moab::Range found; 2162417220eSVijay Mahadevan 2172417220eSVijay Mahadevan /* store the vertex DoF number */ 2182417220eSVijay Mahadevan const int ldof = dmbc->lidmap[vhandle - dmbc->seqstart]; 2192417220eSVijay Mahadevan 2202417220eSVijay Mahadevan /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects 2212417220eSVijay Mahadevan to the current vertex. We can then decipher if a vertex is ghosted or not and compute the 2222417220eSVijay Mahadevan non-zero pattern accordingly. */ 2232417220eSVijay Mahadevan merr = dmbc->hierarchy->get_adjacencies(vhandle, dmbc->dim, adjs);MBERRNM(merr); 2242417220eSVijay Mahadevan 2252417220eSVijay Mahadevan /* loop over vertices and update the number of connectivity */ 2262417220eSVijay Mahadevan for (unsigned jter = 0; jter < adjs.size(); jter++) { 2272417220eSVijay Mahadevan 2282417220eSVijay Mahadevan const moab::EntityHandle jhandle = adjs[jter]; 229941e0cffSVijay Mahadevan 230941e0cffSVijay Mahadevan /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */ 2312417220eSVijay Mahadevan merr = dmbc->hierarchy->child_to_parent(jhandle, dmbc->hlevel, dmbp->hlevel, &parent);MBERRNM(merr); 232941e0cffSVijay Mahadevan 2332417220eSVijay Mahadevan /* Get connectivity information in canonical ordering for the local element */ 2342417220eSVijay Mahadevan std::vector<moab::EntityHandle> connp; 2352417220eSVijay Mahadevan merr = dmbp->hierarchy->get_connectivity(parent, dmbp->hlevel, connp);MBERRNM(merr); 236a86ed394SVijay Mahadevan 2372417220eSVijay Mahadevan for (unsigned ic=0; ic < connp.size(); ++ic) { 2382417220eSVijay Mahadevan 2392417220eSVijay Mahadevan /* loop over each element connected to the adjacent vertex and update as needed */ 2402417220eSVijay Mahadevan /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */ 2412417220eSVijay Mahadevan if (found.find(connp[ic]) != found.end()) continue; /* make sure we don't double count shared vertices */ 2422417220eSVijay Mahadevan if (dmbp->vghost->find(connp[ic]) != dmbp->vghost->end()) onz[ldof]++; /* update out-of-proc onz */ 2432417220eSVijay Mahadevan else nnz[ldof]++; /* else local vertex */ 2442417220eSVijay Mahadevan found.insert(connp[ic]); 2452417220eSVijay Mahadevan } 2462417220eSVijay Mahadevan } 247941e0cffSVijay Mahadevan } 248941e0cffSVijay Mahadevan 2492417220eSVijay Mahadevan for (int i = 0; i < nlsizc; i++) 2502417220eSVijay Mahadevan nnz[i] += 1; /* self count the node */ 251941e0cffSVijay Mahadevan 252941e0cffSVijay Mahadevan ionz = onz[0]; 253941e0cffSVijay Mahadevan innz = nnz[0]; 254755f3dfbSVijay Mahadevan for (int tc = 0; tc < nlsizc; tc++) { 255941e0cffSVijay Mahadevan // check for maximum allowed sparsity = fully dense 256755f3dfbSVijay Mahadevan nnz[tc] = std::min(nlsizp, nnz[tc]); 2572417220eSVijay Mahadevan onz[tc] = std::min(ngsizp - nlsizp, onz[tc]); 2582417220eSVijay Mahadevan 2592417220eSVijay Mahadevan PetscInfo3(NULL, " %d: NNZ = %d, ONZ = %d\n", tc, nnz[tc], onz[tc]); 260941e0cffSVijay Mahadevan 261941e0cffSVijay Mahadevan innz = (innz < nnz[tc] ? nnz[tc] : innz); 262941e0cffSVijay Mahadevan ionz = (ionz < onz[tc] ? onz[tc] : ionz); 263941e0cffSVijay Mahadevan } 26464e1c140SVijay Mahadevan 26564e1c140SVijay Mahadevan /* create interpolation matrix */ 266755f3dfbSVijay Mahadevan ierr = MatCreate(PetscObjectComm((PetscObject)dmc), interpl);CHKERRQ(ierr); 267755f3dfbSVijay Mahadevan ierr = MatSetSizes(*interpl, nlsizc, nlsizp, ngsizc, ngsizp);CHKERRQ(ierr); 26864e1c140SVijay Mahadevan ierr = MatSetType(*interpl, MATAIJ);CHKERRQ(ierr); 26964e1c140SVijay Mahadevan ierr = MatSetFromOptions(*interpl);CHKERRQ(ierr); 27064e1c140SVijay Mahadevan 271941e0cffSVijay Mahadevan ierr = MatSeqAIJSetPreallocation(*interpl, innz, nnz);CHKERRQ(ierr); 272941e0cffSVijay Mahadevan ierr = MatMPIAIJSetPreallocation(*interpl, innz, nnz, ionz, onz);CHKERRQ(ierr); 273941e0cffSVijay Mahadevan 274941e0cffSVijay Mahadevan /* clean up temporary memory */ 275941e0cffSVijay Mahadevan ierr = PetscFree2(nnz, onz);CHKERRQ(ierr); 276b117cd09SVijay Mahadevan 277b117cd09SVijay Mahadevan /* set up internal matrix data-structures */ 278b117cd09SVijay Mahadevan ierr = MatSetUp(*interpl);CHKERRQ(ierr); 279b117cd09SVijay Mahadevan 280a86ed394SVijay Mahadevan /* Define variables for assembly */ 281a86ed394SVijay Mahadevan std::vector<moab::EntityHandle> children; 282a86ed394SVijay Mahadevan std::vector<moab::EntityHandle> connp, connc; 283a86ed394SVijay Mahadevan std::vector<PetscReal> pcoords, ccoords, values_phi; 284a86ed394SVijay Mahadevan 285a86ed394SVijay Mahadevan if (use_consistent_bases) { 2862417220eSVijay Mahadevan const moab::EntityHandle ehandle = dmbp->elocal->front(); 287a86ed394SVijay Mahadevan 288a86ed394SVijay Mahadevan merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children);MBERRNM(merr); 289a86ed394SVijay Mahadevan 290a86ed394SVijay Mahadevan /* Get connectivity and coordinates of the parent vertices */ 291a86ed394SVijay Mahadevan merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp);MBERRNM(merr); 2922417220eSVijay Mahadevan merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc);MBERRNM(merr); 293a86ed394SVijay Mahadevan 294a86ed394SVijay Mahadevan std::vector<PetscReal> natparam(3*connc.size(), 0.0); 295a86ed394SVijay Mahadevan pcoords.resize(connp.size() * 3); 296a86ed394SVijay Mahadevan ccoords.resize(connc.size() * 3); 297a86ed394SVijay Mahadevan values_phi.resize(connp.size()*connc.size()); 298a86ed394SVijay Mahadevan /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */ 299a86ed394SVijay Mahadevan merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]);MBERRNM(merr); 300a86ed394SVijay Mahadevan merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]);MBERRNM(merr); 301a86ed394SVijay Mahadevan 302a86ed394SVijay Mahadevan /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */ 303a86ed394SVijay Mahadevan for (unsigned tc = 0; tc < connc.size(); tc++) { 304a86ed394SVijay Mahadevan const PetscInt offset = tc * 3; 305a86ed394SVijay Mahadevan 306a86ed394SVijay Mahadevan /* Scale ccoords relative to pcoords */ 307a86ed394SVijay Mahadevan ierr = DMMoabPToRMapping(dim, connp.size(), &pcoords[0], &ccoords[offset], &natparam[offset], &values_phi[connp.size()*tc]);CHKERRQ(ierr); 308a86ed394SVijay Mahadevan } 309a86ed394SVijay Mahadevan } 310a86ed394SVijay Mahadevan else { 311a86ed394SVijay Mahadevan factor = std::pow(2.0 /*degree_P_for_refinement*/, (dmbc->hlevel - dmbp->hlevel) * dmbp->dim * 1.0); 312a86ed394SVijay Mahadevan } 313a86ed394SVijay Mahadevan 314755f3dfbSVijay Mahadevan /* TODO: Decipher the correct non-zero pattern. There is still some issue with onz allocation */ 3152417220eSVijay Mahadevan ierr = MatSetOption(*interpl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);CHKERRQ(ierr); 316b117cd09SVijay Mahadevan 317e882eb38SVijay Mahadevan /* Loop through the remaining vertices. These vertices appear only on the current refined_level. */ 3182417220eSVijay Mahadevan for (moab::Range::iterator iter = dmbp->elocal->begin(); iter != dmbp->elocal->end(); iter++) { 319b117cd09SVijay Mahadevan 320b117cd09SVijay Mahadevan const moab::EntityHandle ehandle = *iter; 321b117cd09SVijay Mahadevan 322b117cd09SVijay Mahadevan /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */ 323a86ed394SVijay Mahadevan children.clear(); 324a86ed394SVijay Mahadevan connc.clear(); 325755f3dfbSVijay Mahadevan merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children);MBERRNM(merr); 326b117cd09SVijay Mahadevan 327b117cd09SVijay Mahadevan /* Get connectivity and coordinates of the parent vertices */ 328755f3dfbSVijay Mahadevan merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp);MBERRNM(merr); 3292417220eSVijay Mahadevan merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc);MBERRNM(merr); 330b117cd09SVijay Mahadevan 331a86ed394SVijay Mahadevan pcoords.resize(connp.size() * 3); 332a86ed394SVijay Mahadevan ccoords.resize(connc.size() * 3); 333b117cd09SVijay Mahadevan /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */ 334755f3dfbSVijay Mahadevan merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]);MBERRNM(merr); 335755f3dfbSVijay Mahadevan merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]);MBERRNM(merr); 336b117cd09SVijay Mahadevan 337b117cd09SVijay Mahadevan std::vector<int> dofsp(connp.size()), dofsc(connc.size()); 338b117cd09SVijay Mahadevan /* TODO: specific to scalar system - use GetDofs */ 339755f3dfbSVijay Mahadevan ierr = DMMoabGetDofsBlocked(dmp, connp.size(), &connp[0], &dofsp[0]);CHKERRQ(ierr); 340755f3dfbSVijay Mahadevan ierr = DMMoabGetDofsBlocked(dmc, connc.size(), &connc[0], &dofsc[0]);CHKERRQ(ierr); 341b117cd09SVijay Mahadevan 342a86ed394SVijay Mahadevan /* Compute the actual interpolation weights when projecting solution/residual between levels */ 343a86ed394SVijay Mahadevan if (use_consistent_bases) { 344a86ed394SVijay Mahadevan 345a86ed394SVijay Mahadevan /* Use the cached values of natural parameteric coordinates and basis pre-evaluated. 346a86ed394SVijay Mahadevan We are making an assumption here that UMR used in GMG to generate the hierarchy uses 347a86ed394SVijay Mahadevan the same template for all elements; This will fail for mixed element meshes (TRI/QUAD). 348a86ed394SVijay Mahadevan 3492417220eSVijay Mahadevan TODO: Fix the above assumption by caching data for families (especially for Tets and mixed meshes) 350a86ed394SVijay Mahadevan */ 351a86ed394SVijay Mahadevan 352a86ed394SVijay Mahadevan /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */ 353a86ed394SVijay Mahadevan for (unsigned tc = 0; tc < connc.size(); tc++) { 354a86ed394SVijay Mahadevan /* TODO: Check if we should be using INSERT_VALUES instead */ 355a86ed394SVijay Mahadevan ierr = MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi[connp.size()*tc], ADD_VALUES);CHKERRQ(ierr); 356a86ed394SVijay Mahadevan } 357a86ed394SVijay Mahadevan } 358a86ed394SVijay Mahadevan else { 359e882eb38SVijay Mahadevan /* Compute the interpolation weights by determining distance of 1-ring 360755f3dfbSVijay Mahadevan neighbor vertices from current vertex 361755f3dfbSVijay Mahadevan 362a86ed394SVijay Mahadevan This should be used only when FEM basis is not used for the discretization. 363a86ed394SVijay Mahadevan Else, the consistent interface to compute the basis function for interpolation 364a86ed394SVijay Mahadevan between the levels should be evaluated correctly to preserve convergence of GMG. 3652417220eSVijay Mahadevan Shephard's basis will be terrible for any unsmooth problems. 366755f3dfbSVijay Mahadevan */ 367755f3dfbSVijay Mahadevan values_phi.resize(connp.size()); 368755f3dfbSVijay Mahadevan for (unsigned tc = 0; tc < connc.size(); tc++) { 369755f3dfbSVijay Mahadevan 370a86ed394SVijay Mahadevan PetscReal normsum = 0.0; 371755f3dfbSVijay Mahadevan for (unsigned tp = 0; tp < connp.size(); tp++) { 372755f3dfbSVijay Mahadevan values_phi[tp] = 0.0; 373e882eb38SVijay Mahadevan for (unsigned k = 0; k < 3; k++) 374755f3dfbSVijay Mahadevan values_phi[tp] += std::pow(pcoords[tp * 3 + k] - ccoords[k + tc * 3], dim); 375755f3dfbSVijay Mahadevan if (values_phi[tp] < 1e-12) { 376755f3dfbSVijay Mahadevan values_phi[tp] = 1e12; 377b117cd09SVijay Mahadevan } 378b117cd09SVijay Mahadevan else { 379755f3dfbSVijay Mahadevan //values_phi[tp] = std::pow(values_phi[tp], -1.0/dim); 380755f3dfbSVijay Mahadevan values_phi[tp] = std::pow(values_phi[tp], -1.0); 381755f3dfbSVijay Mahadevan normsum += values_phi[tp]; 382b117cd09SVijay Mahadevan } 383b117cd09SVijay Mahadevan } 384755f3dfbSVijay Mahadevan for (unsigned tp = 0; tp < connp.size(); tp++) { 385755f3dfbSVijay Mahadevan if (values_phi[tp] > 1e11) 386755f3dfbSVijay Mahadevan values_phi[tp] = factor * 0.5 / connp.size(); 387b117cd09SVijay Mahadevan else 388755f3dfbSVijay Mahadevan values_phi[tp] = factor * values_phi[tp] * 0.5 / (connp.size() * normsum); 389b117cd09SVijay Mahadevan } 390755f3dfbSVijay Mahadevan ierr = MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi[0], ADD_VALUES);CHKERRQ(ierr); 391b117cd09SVijay Mahadevan } 392a86ed394SVijay Mahadevan } 393b117cd09SVijay Mahadevan } 394304006b3SVijay Mahadevan if (vec) *vec = NULL; 395b117cd09SVijay Mahadevan ierr = MatAssemblyBegin(*interpl, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 396b117cd09SVijay Mahadevan ierr = MatAssemblyEnd(*interpl, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 397b117cd09SVijay Mahadevan PetscFunctionReturn(0); 398b117cd09SVijay Mahadevan } 399b117cd09SVijay Mahadevan 400b117cd09SVijay Mahadevan #undef __FUNCT__ 401b117cd09SVijay Mahadevan #define __FUNCT__ "DMCreateInjection_Moab" 402b117cd09SVijay Mahadevan /*@ 403b117cd09SVijay Mahadevan DMCreateInjection_Moab - Generate a multi-level uniform refinement hierarchy 404b117cd09SVijay Mahadevan by succesively refining a coarse mesh, already defined in the DM object 405b117cd09SVijay Mahadevan provided by the user. 406b117cd09SVijay Mahadevan 407b117cd09SVijay Mahadevan Collective on MPI_Comm 408b117cd09SVijay Mahadevan 409b117cd09SVijay Mahadevan Input Parameter: 410b117cd09SVijay Mahadevan . dmb - The DMMoab object 411b117cd09SVijay Mahadevan 412b117cd09SVijay Mahadevan Output Parameter: 413b117cd09SVijay Mahadevan . nlevels - The number of levels of refinement needed to generate the hierarchy 414b117cd09SVijay Mahadevan + ldegrees - The degree of refinement at each level in the hierarchy 415b117cd09SVijay Mahadevan 416b117cd09SVijay Mahadevan Level: beginner 417b117cd09SVijay Mahadevan 418b117cd09SVijay Mahadevan .keywords: DMMoab, create, refinement 419b117cd09SVijay Mahadevan @*/ 420304006b3SVijay Mahadevan PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dm1, DM dm2, VecScatter* ctx) 421b117cd09SVijay Mahadevan { 422e882eb38SVijay Mahadevan //DM_Moab *dmmoab; 423b117cd09SVijay Mahadevan 424b117cd09SVijay Mahadevan PetscFunctionBegin; 425b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm1, DM_CLASSID, 1); 426b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm2, DM_CLASSID, 2); 427e882eb38SVijay Mahadevan //dmmoab = (DM_Moab*)(dm1)->data; 428b117cd09SVijay Mahadevan 429b117cd09SVijay Mahadevan PetscPrintf(PETSC_COMM_WORLD, "[DMCreateInjection_Moab] :: Placeholder\n"); 430b117cd09SVijay Mahadevan PetscFunctionReturn(0); 431b117cd09SVijay Mahadevan } 432b117cd09SVijay Mahadevan 433b117cd09SVijay Mahadevan 434b117cd09SVijay Mahadevan #undef __FUNCT__ 435b117cd09SVijay Mahadevan #define __FUNCT__ "DM_UMR_Moab_Private" 436b117cd09SVijay Mahadevan PetscErrorCode DM_UMR_Moab_Private(DM dm, MPI_Comm comm, PetscBool refine, DM *dmref) 437b117cd09SVijay Mahadevan { 438b117cd09SVijay Mahadevan PetscErrorCode ierr; 439e882eb38SVijay Mahadevan PetscInt i, dim; 440b117cd09SVijay Mahadevan DM dm2; 441b117cd09SVijay Mahadevan moab::ErrorCode merr; 442b117cd09SVijay Mahadevan DM_Moab *dmb = (DM_Moab*)dm->data, *dd2; 443b117cd09SVijay Mahadevan 444b117cd09SVijay Mahadevan PetscFunctionBegin; 445b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 446b117cd09SVijay Mahadevan PetscValidPointer(dmref, 4); 447b117cd09SVijay Mahadevan 448b117cd09SVijay Mahadevan if ( (dmb->hlevel == dmb->nhlevels && refine) || (dmb->hlevel == 0 && !refine) ) { 449e882eb38SVijay Mahadevan if (dmb->hlevel + 1 > dmb->nhlevels && refine) PetscInfo2(NULL, "Invalid multigrid refinement hierarchy level specified (%D). MOAB UMR max levels = %D. Creating a NULL object.\n", dmb->hlevel + 1, dmb->nhlevels); 450e882eb38SVijay Mahadevan if (dmb->hlevel - 1 < 0 && !refine) PetscInfo1(NULL, "Invalid multigrid coarsen hierarchy level specified (%D). Creating a NULL object.\n", dmb->hlevel - 1); 451e882eb38SVijay Mahadevan *dmref = PETSC_NULL; 452b117cd09SVijay Mahadevan PetscFunctionReturn(0); 453b117cd09SVijay Mahadevan } 454b117cd09SVijay Mahadevan 455b117cd09SVijay Mahadevan ierr = DMMoabCreate(PetscObjectComm((PetscObject)dm), &dm2);CHKERRQ(ierr); 456b117cd09SVijay Mahadevan dd2 = (DM_Moab*)dm2->data; 457b117cd09SVijay Mahadevan 458b117cd09SVijay Mahadevan dd2->mbiface = dmb->mbiface; 4599daf19fdSVijay Mahadevan #ifdef MOAB_HAVE_MPI 460b117cd09SVijay Mahadevan dd2->pcomm = dmb->pcomm; 4619daf19fdSVijay Mahadevan #endif 462b117cd09SVijay Mahadevan dd2->icreatedinstance = PETSC_FALSE; 46364e1c140SVijay Mahadevan dd2->nghostrings = dmb->nghostrings; 464b117cd09SVijay Mahadevan 465b117cd09SVijay Mahadevan /* set the new level based on refinement/coarsening */ 466b117cd09SVijay Mahadevan if (refine) { 467b117cd09SVijay Mahadevan dd2->hlevel = dmb->hlevel + 1; 468b117cd09SVijay Mahadevan } 469b117cd09SVijay Mahadevan else { 470b117cd09SVijay Mahadevan dd2->hlevel = dmb->hlevel - 1; 471b117cd09SVijay Mahadevan } 472b117cd09SVijay Mahadevan 473b117cd09SVijay Mahadevan /* Copy the multilevel hierarchy pointers in MOAB */ 474b117cd09SVijay Mahadevan dd2->hierarchy = dmb->hierarchy; 475b117cd09SVijay Mahadevan dd2->nhlevels = dmb->nhlevels; 476b117cd09SVijay Mahadevan ierr = PetscMalloc1(dd2->nhlevels + 1, &dd2->hsets);CHKERRQ(ierr); 477b117cd09SVijay Mahadevan for (i = 0; i <= dd2->nhlevels; i++) { 478b117cd09SVijay Mahadevan dd2->hsets[i] = dmb->hsets[i]; 479b117cd09SVijay Mahadevan } 480b117cd09SVijay Mahadevan dd2->fileset = dd2->hsets[dd2->hlevel]; 481b117cd09SVijay Mahadevan 482b117cd09SVijay Mahadevan /* do the remaining initializations for DMMoab */ 483b117cd09SVijay Mahadevan dd2->bs = dmb->bs; 484b117cd09SVijay Mahadevan dd2->numFields = dmb->numFields; 485b117cd09SVijay Mahadevan dd2->rw_dbglevel = dmb->rw_dbglevel; 486b117cd09SVijay Mahadevan dd2->partition_by_rank = dmb->partition_by_rank; 487b117cd09SVijay Mahadevan ierr = PetscStrcpy(dd2->extra_read_options, dmb->extra_read_options);CHKERRQ(ierr); 488b117cd09SVijay Mahadevan ierr = PetscStrcpy(dd2->extra_write_options, dmb->extra_write_options);CHKERRQ(ierr); 489b117cd09SVijay Mahadevan dd2->read_mode = dmb->read_mode; 490b117cd09SVijay Mahadevan dd2->write_mode = dmb->write_mode; 491b117cd09SVijay Mahadevan 492b117cd09SVijay Mahadevan /* set global ID tag handle */ 493b117cd09SVijay Mahadevan ierr = DMMoabSetLocalToGlobalTag(dm2, dmb->ltog_tag);CHKERRQ(ierr); 494b117cd09SVijay Mahadevan 495b117cd09SVijay Mahadevan merr = dd2->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dd2->material_tag);MBERRNM(merr); 496b117cd09SVijay Mahadevan 497b117cd09SVijay Mahadevan ierr = DMSetOptionsPrefix(dm2, ((PetscObject)dm)->prefix);CHKERRQ(ierr); 498b117cd09SVijay Mahadevan ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 499b117cd09SVijay Mahadevan ierr = DMSetDimension(dm2, dim);CHKERRQ(ierr); 500b117cd09SVijay Mahadevan 501b117cd09SVijay Mahadevan /* allow overloaded (user replaced) operations to be inherited by refinement clones */ 502b117cd09SVijay Mahadevan dm2->ops->creatematrix = dm->ops->creatematrix; 503b117cd09SVijay Mahadevan 504b117cd09SVijay Mahadevan /* copy fill information if given */ 505b117cd09SVijay Mahadevan ierr = DMMoabSetBlockFills(dm2, dmb->dfill, dmb->ofill);CHKERRQ(ierr); 506b117cd09SVijay Mahadevan 507b117cd09SVijay Mahadevan /* copy vector type information */ 508b117cd09SVijay Mahadevan ierr = DMSetMatType(dm2, dm->mattype);CHKERRQ(ierr); 509b117cd09SVijay Mahadevan ierr = DMSetVecType(dm2, dm->vectype);CHKERRQ(ierr); 5103f1c6e43SVijay Mahadevan dd2->numFields = dmb->numFields; 5113f1c6e43SVijay Mahadevan if (dmb->numFields) { 512b117cd09SVijay Mahadevan ierr = DMMoabSetFieldNames(dm2, dmb->numFields, dmb->fieldNames);CHKERRQ(ierr); 5133f1c6e43SVijay Mahadevan } 514b117cd09SVijay Mahadevan 515b117cd09SVijay Mahadevan ierr = DMSetFromOptions(dm2);CHKERRQ(ierr); 516b117cd09SVijay Mahadevan 517b117cd09SVijay Mahadevan /* recreate Dof numbering for the refined DM and make sure the distribution is correctly populated */ 518b117cd09SVijay Mahadevan ierr = DMSetUp(dm2);CHKERRQ(ierr); 519b117cd09SVijay Mahadevan 520b117cd09SVijay Mahadevan *dmref = dm2; 521b117cd09SVijay Mahadevan PetscFunctionReturn(0); 522b117cd09SVijay Mahadevan } 523b117cd09SVijay Mahadevan 524b117cd09SVijay Mahadevan 525b117cd09SVijay Mahadevan #undef __FUNCT__ 526b117cd09SVijay Mahadevan #define __FUNCT__ "DMRefine_Moab" 527b117cd09SVijay Mahadevan /*@ 528b117cd09SVijay Mahadevan DMRefine_Moab - Generate a multi-level uniform refinement hierarchy 529b117cd09SVijay Mahadevan by succesively refining a coarse mesh, already defined in the DM object 530b117cd09SVijay Mahadevan provided by the user. 531b117cd09SVijay Mahadevan 532e882eb38SVijay Mahadevan Collective on DM 533b117cd09SVijay Mahadevan 534b117cd09SVijay Mahadevan Input Parameter: 535e882eb38SVijay Mahadevan + dm - The DMMoab object 536e882eb38SVijay Mahadevan - comm - the communicator to contain the new DM object (or MPI_COMM_NULL) 537b117cd09SVijay Mahadevan 538b117cd09SVijay Mahadevan Output Parameter: 539e882eb38SVijay Mahadevan . dmf - the refined DM, or NULL 540b117cd09SVijay Mahadevan 541e882eb38SVijay Mahadevan Note: If no refinement was done, the return value is NULL 542e882eb38SVijay Mahadevan 543e882eb38SVijay Mahadevan Level: developer 544b117cd09SVijay Mahadevan 545b117cd09SVijay Mahadevan .keywords: DMMoab, create, refinement 546b117cd09SVijay Mahadevan @*/ 547304006b3SVijay Mahadevan PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM* dmf) 548b117cd09SVijay Mahadevan { 549b117cd09SVijay Mahadevan PetscErrorCode ierr; 550b117cd09SVijay Mahadevan 551b117cd09SVijay Mahadevan PetscFunctionBegin; 552b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 553b117cd09SVijay Mahadevan 554b117cd09SVijay Mahadevan ierr = DM_UMR_Moab_Private(dm, comm, PETSC_TRUE, dmf);CHKERRQ(ierr); 555b117cd09SVijay Mahadevan PetscFunctionReturn(0); 556b117cd09SVijay Mahadevan } 557b117cd09SVijay Mahadevan 558b117cd09SVijay Mahadevan #undef __FUNCT__ 559b117cd09SVijay Mahadevan #define __FUNCT__ "DMCoarsen_Moab" 560b117cd09SVijay Mahadevan /*@ 561b117cd09SVijay Mahadevan DMCoarsen_Moab - Generate a multi-level uniform refinement hierarchy 562b117cd09SVijay Mahadevan by succesively refining a coarse mesh, already defined in the DM object 563b117cd09SVijay Mahadevan provided by the user. 564b117cd09SVijay Mahadevan 565e882eb38SVijay Mahadevan Collective on DM 566b117cd09SVijay Mahadevan 567b117cd09SVijay Mahadevan Input Parameter: 568e882eb38SVijay Mahadevan . dm - The DMMoab object 569e882eb38SVijay Mahadevan - comm - the communicator to contain the new DM object (or MPI_COMM_NULL) 570b117cd09SVijay Mahadevan 571b117cd09SVijay Mahadevan Output Parameter: 572e882eb38SVijay Mahadevan . dmf - the coarsened DM, or NULL 573b117cd09SVijay Mahadevan 574e882eb38SVijay Mahadevan Note: If no coarsening was done, the return value is NULL 575b117cd09SVijay Mahadevan 576e882eb38SVijay Mahadevan Level: developer 577e882eb38SVijay Mahadevan 578e882eb38SVijay Mahadevan .keywords: DMMoab, create, coarsening 579b117cd09SVijay Mahadevan @*/ 580304006b3SVijay Mahadevan PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM* dmc) 581b117cd09SVijay Mahadevan { 582b117cd09SVijay Mahadevan PetscErrorCode ierr; 583b117cd09SVijay Mahadevan 584b117cd09SVijay Mahadevan PetscFunctionBegin; 585b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 586b117cd09SVijay Mahadevan 587b117cd09SVijay Mahadevan ierr = DM_UMR_Moab_Private(dm, comm, PETSC_FALSE, dmc);CHKERRQ(ierr); 588b117cd09SVijay Mahadevan PetscFunctionReturn(0); 589b117cd09SVijay Mahadevan } 590