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; 31b117cd09SVijay Mahadevan PetscInt *pdegrees, i; 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); 40b117cd09SVijay Mahadevan for (i = 0; i < nlevels; i++) pdegrees[i] = 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 */ 5764e1c140SVijay Mahadevan merr = dmmoab->hierarchy->generate_mesh_hierarchy(nlevels, pdegrees, hsets); 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 */ 6664e1c140SVijay Mahadevan for (i = 0; i <= nlevels; i++) 6749d66b22SVijay Mahadevan dmmoab->hsets[i] = hsets[i]; 6864e1c140SVijay Mahadevan 6964e1c140SVijay Mahadevan hsets.clear(); 70b117cd09SVijay Mahadevan if (!ldegrees) { 71b117cd09SVijay Mahadevan ierr = PetscFree(pdegrees);CHKERRQ(ierr); 72b117cd09SVijay Mahadevan } 73b117cd09SVijay Mahadevan PetscFunctionReturn(0); 74b117cd09SVijay Mahadevan } 75b117cd09SVijay Mahadevan 76b117cd09SVijay Mahadevan #undef __FUNCT__ 77b117cd09SVijay Mahadevan #define __FUNCT__ "DMRefineHierarchy_Moab" 78b117cd09SVijay Mahadevan /*@ 79e882eb38SVijay Mahadevan DMRefineHierarchy_Moab - Generate a multi-level DM hierarchy 80e882eb38SVijay Mahadevan by succesively refining a coarse mesh. 81b117cd09SVijay Mahadevan 82b117cd09SVijay Mahadevan Collective on MPI_Comm 83b117cd09SVijay Mahadevan 84b117cd09SVijay Mahadevan Input Parameter: 85e882eb38SVijay Mahadevan + dm - The DMMoab object 86b117cd09SVijay Mahadevan 87b117cd09SVijay Mahadevan Output Parameter: 88e882eb38SVijay Mahadevan + nlevels - The number of levels of refinement needed to generate the hierarchy 89e882eb38SVijay Mahadevan . dmf - The DM objects after successive refinement of the hierarchy 90b117cd09SVijay Mahadevan 91b117cd09SVijay Mahadevan Level: beginner 92b117cd09SVijay Mahadevan 93e882eb38SVijay Mahadevan .keywords: DMMoab, generate, refinement 94b117cd09SVijay Mahadevan @*/ 95304006b3SVijay Mahadevan PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmf[]) 96b117cd09SVijay Mahadevan { 97b117cd09SVijay Mahadevan PetscErrorCode ierr; 98b117cd09SVijay Mahadevan PetscInt i; 99b117cd09SVijay Mahadevan 100b117cd09SVijay Mahadevan PetscFunctionBegin; 101b117cd09SVijay Mahadevan 102e882eb38SVijay Mahadevan ierr = DMRefine(dm, PetscObjectComm((PetscObject)dm), &dmf[0]);CHKERRQ(ierr); 103e882eb38SVijay Mahadevan for (i = 1; i < nlevels; i++) { 104e882eb38SVijay Mahadevan ierr = DMRefine(dmf[i - 1], PetscObjectComm((PetscObject)dm), &dmf[i]);CHKERRQ(ierr); 105b117cd09SVijay Mahadevan } 106b117cd09SVijay Mahadevan PetscFunctionReturn(0); 107b117cd09SVijay Mahadevan } 108b117cd09SVijay Mahadevan 109b117cd09SVijay Mahadevan #undef __FUNCT__ 110b117cd09SVijay Mahadevan #define __FUNCT__ "DMCoarsenHierarchy_Moab" 111b117cd09SVijay Mahadevan /*@ 112e882eb38SVijay Mahadevan DMCoarsenHierarchy_Moab - Generate a multi-level DM hierarchy 113e882eb38SVijay Mahadevan by succesively coarsening a refined mesh. 114b117cd09SVijay Mahadevan 115b117cd09SVijay Mahadevan Collective on MPI_Comm 116b117cd09SVijay Mahadevan 117b117cd09SVijay Mahadevan Input Parameter: 118e882eb38SVijay Mahadevan + dm - The DMMoab object 119b117cd09SVijay Mahadevan 120b117cd09SVijay Mahadevan Output Parameter: 121e882eb38SVijay Mahadevan + nlevels - The number of levels of refinement needed to generate the hierarchy 122e882eb38SVijay Mahadevan . dmc - The DM objects after successive coarsening of the hierarchy 123b117cd09SVijay Mahadevan 124b117cd09SVijay Mahadevan Level: beginner 125b117cd09SVijay Mahadevan 126e882eb38SVijay Mahadevan .keywords: DMMoab, generate, coarsening 127b117cd09SVijay Mahadevan @*/ 128304006b3SVijay Mahadevan PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmc[]) 129b117cd09SVijay Mahadevan { 130b117cd09SVijay Mahadevan PetscErrorCode ierr; 131b117cd09SVijay Mahadevan PetscInt i; 132b117cd09SVijay Mahadevan 133b117cd09SVijay Mahadevan PetscFunctionBegin; 134b117cd09SVijay Mahadevan 135e882eb38SVijay Mahadevan ierr = DMCoarsen(dm, PetscObjectComm((PetscObject)dm), &dmc[0]);CHKERRQ(ierr); 136e882eb38SVijay Mahadevan for (i = 1; i < nlevels; i++) { 137e882eb38SVijay Mahadevan ierr = DMCoarsen(dmc[i - 1], PetscObjectComm((PetscObject)dm), &dmc[i]);CHKERRQ(ierr); 138b117cd09SVijay Mahadevan } 139b117cd09SVijay Mahadevan PetscFunctionReturn(0); 140b117cd09SVijay Mahadevan } 141b117cd09SVijay Mahadevan 142b117cd09SVijay Mahadevan 143b117cd09SVijay Mahadevan #undef __FUNCT__ 144b117cd09SVijay Mahadevan #define __FUNCT__ "DMCreateInterpolation_Moab" 145b117cd09SVijay Mahadevan /*@ 146e882eb38SVijay Mahadevan DMCreateInterpolation_Moab - Generate the interpolation operators to transform 147e882eb38SVijay Mahadevan operators (matrices, vectors) from parent level to child level as defined by 148e882eb38SVijay Mahadevan the DM inputs provided by the user. 149b117cd09SVijay Mahadevan 150b117cd09SVijay Mahadevan Collective on MPI_Comm 151b117cd09SVijay Mahadevan 152b117cd09SVijay Mahadevan Input Parameter: 153e882eb38SVijay Mahadevan + dm1 - The DMMoab object 154e882eb38SVijay Mahadevan - dm2 - the second, finer DMMoab object 155b117cd09SVijay Mahadevan 156b117cd09SVijay Mahadevan Output Parameter: 157e882eb38SVijay Mahadevan + interpl - The interpolation operator for transferring data between the levels 158e882eb38SVijay Mahadevan - vec - The scaling vector (optional) 159b117cd09SVijay Mahadevan 160e882eb38SVijay Mahadevan Level: developer 161b117cd09SVijay Mahadevan 162b117cd09SVijay Mahadevan .keywords: DMMoab, create, refinement 163b117cd09SVijay Mahadevan @*/ 164304006b3SVijay Mahadevan PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmp, DM dmc, Mat* interpl, Vec* vec) 165b117cd09SVijay Mahadevan { 166755f3dfbSVijay Mahadevan DM_Moab *dmbp, *dmbc; 167b117cd09SVijay Mahadevan PetscErrorCode ierr; 168b117cd09SVijay Mahadevan moab::ErrorCode merr; 169e882eb38SVijay Mahadevan PetscInt dim; 1703f1c6e43SVijay Mahadevan PetscReal factor; 171ce27a4eeSVijay Mahadevan PetscInt innz, *nnz, ionz, *onz; 172755f3dfbSVijay Mahadevan PetscInt nlsizp, nlsizc, nlghsizp, ngsizp, ngsizc; 173755f3dfbSVijay Mahadevan moab::Range eowned; 174*a86ed394SVijay Mahadevan const PetscBool use_consistent_bases=PETSC_TRUE; 175b117cd09SVijay Mahadevan 176b117cd09SVijay Mahadevan PetscFunctionBegin; 177755f3dfbSVijay Mahadevan PetscValidHeaderSpecific(dmp, DM_CLASSID, 1); 178755f3dfbSVijay Mahadevan PetscValidHeaderSpecific(dmc, DM_CLASSID, 2); 179755f3dfbSVijay Mahadevan dmbp = (DM_Moab*)(dmp)->data; 180755f3dfbSVijay Mahadevan dmbc = (DM_Moab*)(dmc)->data; 181755f3dfbSVijay Mahadevan nlsizp = dmbp->nloc;// *dmb1->numFields; 182755f3dfbSVijay Mahadevan nlsizc = dmbc->nloc;// *dmb2->numFields; 183755f3dfbSVijay Mahadevan ngsizp = dmbp->n; // *dmb1->numFields; 184755f3dfbSVijay Mahadevan ngsizc = dmbc->n; // *dmb2->numFields; 185755f3dfbSVijay Mahadevan nlghsizp = (dmbp->nloc + dmbp->nghost); // *dmb1->numFields; 186b117cd09SVijay Mahadevan 187755f3dfbSVijay Mahadevan // Interpolation matrix: \sum_{i=1}^P Owned(Child) * (Owned(Parent) + Ghosted(Parent)) 188755f3dfbSVijay Mahadevan // Size: nlsizc * nlghsizp 189755f3dfbSVijay Mahadevan PetscInfo4(NULL, "Creating interpolation matrix %D X %D to apply transformation between levels %D -> %D.\n", ngsizc, nlghsizp, dmbp->hlevel, dmbc->hlevel); 190b117cd09SVijay Mahadevan 191*a86ed394SVijay Mahadevan ierr = DMGetDimension(dmp, &dim);CHKERRQ(ierr); 192*a86ed394SVijay Mahadevan 193941e0cffSVijay Mahadevan /* allocate the nnz, onz arrays based on block size and local nodes */ 194755f3dfbSVijay Mahadevan ierr = PetscCalloc2(nlsizc, &nnz, nlsizc, &onz);CHKERRQ(ierr); 195941e0cffSVijay Mahadevan 196755f3dfbSVijay Mahadevan eowned = *dmbp->elocal; 1979daf19fdSVijay Mahadevan #ifdef MOAB_HAVE_MPI 198755f3dfbSVijay Mahadevan merr = dmbp->pcomm->filter_pstatus(eowned, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr); 1999daf19fdSVijay Mahadevan #endif 200941e0cffSVijay Mahadevan /* Loop through the local elements and compute the relation between the current parent and the refined_level. */ 201755f3dfbSVijay Mahadevan for (moab::Range::iterator iter = eowned.begin(); iter != eowned.end(); iter++) { 202941e0cffSVijay Mahadevan 203941e0cffSVijay Mahadevan const moab::EntityHandle ehandle = *iter; 204941e0cffSVijay Mahadevan std::vector<moab::EntityHandle> children; 205941e0cffSVijay Mahadevan std::vector<moab::EntityHandle> connp, connc; 206755f3dfbSVijay Mahadevan moab::Range connc_owned; 207941e0cffSVijay Mahadevan 208941e0cffSVijay Mahadevan /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */ 209755f3dfbSVijay Mahadevan merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children); MBERRNM(merr); 210941e0cffSVijay Mahadevan 211755f3dfbSVijay Mahadevan /* Get connectivity of the parent elements */ 212755f3dfbSVijay Mahadevan //merr = dmb1->mbiface->get_connectivity(ehandle,connect,vpere,false);MBERRNM(merr); 213755f3dfbSVijay Mahadevan merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp); MBERRNM(merr); 214755f3dfbSVijay Mahadevan merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc_owned); MBERRNM(merr); 215*a86ed394SVijay Mahadevan 216755f3dfbSVijay Mahadevan for (unsigned tc = 0; tc < connc_owned.size(); tc++) { 217755f3dfbSVijay Mahadevan connc.push_back(connc_owned[tc]); 218941e0cffSVijay Mahadevan } 219941e0cffSVijay Mahadevan 220941e0cffSVijay Mahadevan std::vector<int> dofsp(connp.size()), dofsc(connc.size()); 221755f3dfbSVijay Mahadevan /* get DoFs for both coarse and fine scale grid */ 222755f3dfbSVijay Mahadevan ierr = DMMoabGetDofsBlockedLocal(dmp, connp.size(), &connp[0], &dofsp[0]);CHKERRQ(ierr); 223755f3dfbSVijay Mahadevan ierr = DMMoabGetDofsBlockedLocal(dmc, connc.size(), &connc[0], &dofsc[0]);CHKERRQ(ierr); 224941e0cffSVijay Mahadevan 225941e0cffSVijay Mahadevan for (unsigned tp = 0; tp < connp.size(); tp++) { 226755f3dfbSVijay Mahadevan if (dmbp->vowned->find(connp[tp]) != dmbp->vowned->end()) { 22764e1c140SVijay Mahadevan for (unsigned tc = 0; tc < connc.size(); tc++) { 228755f3dfbSVijay Mahadevan if (dmbc->vowned->find(connc[tc]) != dmbc->vowned->end()) { 229755f3dfbSVijay Mahadevan /* FIXME: This will over-allocate since we multi-count shared nodes. 230755f3dfbSVijay Mahadevan Disadvantage of using element traversal for nnz computation. */ 23164e1c140SVijay Mahadevan nnz[dofsc[tc]]++; 23264e1c140SVijay Mahadevan } 233941e0cffSVijay Mahadevan } 234755f3dfbSVijay Mahadevan } 235755f3dfbSVijay Mahadevan else if (dmbp->vghost->find(connp[tp]) != dmbp->vghost->end()) { 23664e1c140SVijay Mahadevan for (unsigned tc = 0; tc < connc.size(); tc++) { 237755f3dfbSVijay Mahadevan if (dmbc->vowned->find(connc[tc]) != dmbc->vowned->end()) { 238755f3dfbSVijay Mahadevan /* FIXME: This will over-allocate since we multi-count shared nodes. 239755f3dfbSVijay Mahadevan Disadvantage of using element traversal for nnz computation. */ 24064e1c140SVijay Mahadevan onz[dofsc[tc]]++; 241941e0cffSVijay Mahadevan } 242941e0cffSVijay Mahadevan } 24364e1c140SVijay Mahadevan } 244755f3dfbSVijay Mahadevan else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid entity in child level %D\n", connp[tp]); 24549d66b22SVijay Mahadevan } 246941e0cffSVijay Mahadevan } 247941e0cffSVijay Mahadevan 248941e0cffSVijay Mahadevan ionz = onz[0]; 249941e0cffSVijay Mahadevan innz = nnz[0]; 250755f3dfbSVijay Mahadevan for (int tc = 0; tc < nlsizc; tc++) { 251941e0cffSVijay Mahadevan // check for maximum allowed sparsity = fully dense 252755f3dfbSVijay Mahadevan nnz[tc] = std::min(nlsizp, nnz[tc]); 253755f3dfbSVijay Mahadevan onz[tc] = std::min(nlghsizp - nlsizp, onz[tc]); 254941e0cffSVijay Mahadevan 255941e0cffSVijay Mahadevan innz = (innz < nnz[tc] ? nnz[tc] : innz); 256941e0cffSVijay Mahadevan ionz = (ionz < onz[tc] ? onz[tc] : ionz); 257941e0cffSVijay Mahadevan } 25864e1c140SVijay Mahadevan 25964e1c140SVijay Mahadevan /* create interpolation matrix */ 260755f3dfbSVijay Mahadevan ierr = MatCreate(PetscObjectComm((PetscObject)dmc), interpl);CHKERRQ(ierr); 261755f3dfbSVijay Mahadevan ierr = MatSetSizes(*interpl, nlsizc, nlsizp, ngsizc, ngsizp);CHKERRQ(ierr); 26264e1c140SVijay Mahadevan //ierr = MatSetType(*interpl, dm1->mattype);CHKERRQ(ierr); 26364e1c140SVijay Mahadevan ierr = MatSetType(*interpl, MATAIJ);CHKERRQ(ierr); 26464e1c140SVijay Mahadevan ierr = MatSetFromOptions(*interpl);CHKERRQ(ierr); 26564e1c140SVijay Mahadevan 266941e0cffSVijay Mahadevan ierr = MatSeqAIJSetPreallocation(*interpl, innz, nnz);CHKERRQ(ierr); 267941e0cffSVijay Mahadevan ierr = MatMPIAIJSetPreallocation(*interpl, innz, nnz, ionz, onz);CHKERRQ(ierr); 268941e0cffSVijay Mahadevan 269941e0cffSVijay Mahadevan /* clean up temporary memory */ 270941e0cffSVijay Mahadevan ierr = PetscFree2(nnz, onz);CHKERRQ(ierr); 271b117cd09SVijay Mahadevan 272b117cd09SVijay Mahadevan /* set up internal matrix data-structures */ 273b117cd09SVijay Mahadevan ierr = MatSetUp(*interpl);CHKERRQ(ierr); 274b117cd09SVijay Mahadevan 275*a86ed394SVijay Mahadevan /* Define variables for assembly */ 276*a86ed394SVijay Mahadevan std::vector<moab::EntityHandle> children; 277*a86ed394SVijay Mahadevan std::vector<moab::EntityHandle> connp, connc; 278*a86ed394SVijay Mahadevan std::vector<PetscReal> pcoords, ccoords, values_phi; 279*a86ed394SVijay Mahadevan moab::Range connc_owned; 280*a86ed394SVijay Mahadevan 281*a86ed394SVijay Mahadevan if (use_consistent_bases) { 282*a86ed394SVijay Mahadevan const moab::EntityHandle ehandle = eowned.front(); 283*a86ed394SVijay Mahadevan 284*a86ed394SVijay Mahadevan merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children); MBERRNM(merr); 285*a86ed394SVijay Mahadevan 286*a86ed394SVijay Mahadevan /* Get connectivity and coordinates of the parent vertices */ 287*a86ed394SVijay Mahadevan connc_owned.clear(); 288*a86ed394SVijay Mahadevan merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp); MBERRNM(merr); 289*a86ed394SVijay Mahadevan merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc_owned); MBERRNM(merr); 290*a86ed394SVijay Mahadevan 291*a86ed394SVijay Mahadevan for (unsigned tc = 0; tc < connc_owned.size(); tc++) { 292*a86ed394SVijay Mahadevan connc.push_back(connc_owned[tc]); 293*a86ed394SVijay Mahadevan } 294*a86ed394SVijay Mahadevan 295*a86ed394SVijay Mahadevan std::vector<PetscReal> natparam(3*connc.size(), 0.0); 296*a86ed394SVijay Mahadevan pcoords.resize(connp.size() * 3); 297*a86ed394SVijay Mahadevan ccoords.resize(connc.size() * 3); 298*a86ed394SVijay Mahadevan values_phi.resize(connp.size()*connc.size()); 299*a86ed394SVijay Mahadevan /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */ 300*a86ed394SVijay Mahadevan merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]); MBERRNM(merr); 301*a86ed394SVijay Mahadevan merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]); MBERRNM(merr); 302*a86ed394SVijay Mahadevan 303*a86ed394SVijay Mahadevan /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */ 304*a86ed394SVijay Mahadevan for (unsigned tc = 0; tc < connc.size(); tc++) { 305*a86ed394SVijay Mahadevan const PetscInt offset = tc * 3; 306*a86ed394SVijay Mahadevan 307*a86ed394SVijay Mahadevan /* Scale ccoords relative to pcoords */ 308*a86ed394SVijay Mahadevan ierr = DMMoabPToRMapping(dim, connp.size(), &pcoords[0], &ccoords[offset], &natparam[offset], &values_phi[connp.size()*tc]);CHKERRQ(ierr); 309*a86ed394SVijay Mahadevan } 310*a86ed394SVijay Mahadevan } 311*a86ed394SVijay Mahadevan else { 312*a86ed394SVijay Mahadevan factor = std::pow(2.0 /*degree_P_for_refinement*/, (dmbc->hlevel - dmbp->hlevel) * dmbp->dim * 1.0); 313*a86ed394SVijay Mahadevan } 314*a86ed394SVijay Mahadevan 315755f3dfbSVijay Mahadevan /* TODO: Decipher the correct non-zero pattern. There is still some issue with onz allocation */ 316755f3dfbSVijay Mahadevan // ierr = MatSetOption(*interpl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);CHKERRQ(ierr); 317b117cd09SVijay Mahadevan 318e882eb38SVijay Mahadevan /* Loop through the remaining vertices. These vertices appear only on the current refined_level. */ 319755f3dfbSVijay Mahadevan for (moab::Range::iterator iter = eowned.begin(); iter != eowned.end(); iter++) { 320b117cd09SVijay Mahadevan 321b117cd09SVijay Mahadevan const moab::EntityHandle ehandle = *iter; 322b117cd09SVijay Mahadevan 323b117cd09SVijay Mahadevan /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */ 324*a86ed394SVijay Mahadevan children.clear(); 325*a86ed394SVijay Mahadevan connc.clear(); 326755f3dfbSVijay Mahadevan merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children); MBERRNM(merr); 327b117cd09SVijay Mahadevan 328b117cd09SVijay Mahadevan /* Get connectivity and coordinates of the parent vertices */ 329*a86ed394SVijay Mahadevan connc_owned.clear(); 330755f3dfbSVijay Mahadevan merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp); MBERRNM(merr); 331755f3dfbSVijay Mahadevan merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc_owned); MBERRNM(merr); 332755f3dfbSVijay Mahadevan for (unsigned tc = 0; tc < connc_owned.size(); tc++) { 333755f3dfbSVijay Mahadevan connc.push_back(connc_owned[tc]); 334b117cd09SVijay Mahadevan } 335b117cd09SVijay Mahadevan 336*a86ed394SVijay Mahadevan pcoords.resize(connp.size() * 3); 337*a86ed394SVijay Mahadevan ccoords.resize(connc.size() * 3); 338b117cd09SVijay Mahadevan /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */ 339755f3dfbSVijay Mahadevan merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]); MBERRNM(merr); 340755f3dfbSVijay Mahadevan merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]); MBERRNM(merr); 341b117cd09SVijay Mahadevan 342b117cd09SVijay Mahadevan std::vector<int> dofsp(connp.size()), dofsc(connc.size()); 343b117cd09SVijay Mahadevan /* TODO: specific to scalar system - use GetDofs */ 344755f3dfbSVijay Mahadevan ierr = DMMoabGetDofsBlocked(dmp, connp.size(), &connp[0], &dofsp[0]);CHKERRQ(ierr); 345755f3dfbSVijay Mahadevan ierr = DMMoabGetDofsBlocked(dmc, connc.size(), &connc[0], &dofsc[0]);CHKERRQ(ierr); 346b117cd09SVijay Mahadevan 347*a86ed394SVijay Mahadevan /* Compute the actual interpolation weights when projecting solution/residual between levels */ 348*a86ed394SVijay Mahadevan if (use_consistent_bases) { 349*a86ed394SVijay Mahadevan 350*a86ed394SVijay Mahadevan /* Use the cached values of natural parameteric coordinates and basis pre-evaluated. 351*a86ed394SVijay Mahadevan We are making an assumption here that UMR used in GMG to generate the hierarchy uses 352*a86ed394SVijay Mahadevan the same template for all elements; This will fail for mixed element meshes (TRI/QUAD). 353*a86ed394SVijay Mahadevan 354*a86ed394SVijay Mahadevan TODO: Fix the above assumption by caching data for families 355*a86ed394SVijay Mahadevan */ 356*a86ed394SVijay Mahadevan 357*a86ed394SVijay Mahadevan /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */ 358*a86ed394SVijay Mahadevan for (unsigned tc = 0; tc < connc.size(); tc++) { 359*a86ed394SVijay Mahadevan /* TODO: Check if we should be using INSERT_VALUES instead */ 360*a86ed394SVijay Mahadevan ierr = MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi[connp.size()*tc], ADD_VALUES);CHKERRQ(ierr); 361*a86ed394SVijay Mahadevan } 362*a86ed394SVijay Mahadevan } 363*a86ed394SVijay Mahadevan else { 364e882eb38SVijay Mahadevan /* Compute the interpolation weights by determining distance of 1-ring 365755f3dfbSVijay Mahadevan neighbor vertices from current vertex 366755f3dfbSVijay Mahadevan 367*a86ed394SVijay Mahadevan This should be used only when FEM basis is not used for the discretization. 368*a86ed394SVijay Mahadevan Else, the consistent interface to compute the basis function for interpolation 369*a86ed394SVijay Mahadevan between the levels should be evaluated correctly to preserve convergence of GMG. 370755f3dfbSVijay Mahadevan RBF basis will be terrible for any unsmooth problems.. 371755f3dfbSVijay Mahadevan */ 372755f3dfbSVijay Mahadevan values_phi.resize(connp.size()); 373755f3dfbSVijay Mahadevan for (unsigned tc = 0; tc < connc.size(); tc++) { 374755f3dfbSVijay Mahadevan 375*a86ed394SVijay Mahadevan PetscReal normsum = 0.0; 376755f3dfbSVijay Mahadevan for (unsigned tp = 0; tp < connp.size(); tp++) { 377755f3dfbSVijay Mahadevan values_phi[tp] = 0.0; 378e882eb38SVijay Mahadevan for (unsigned k = 0; k < 3; k++) 379755f3dfbSVijay Mahadevan values_phi[tp] += std::pow(pcoords[tp * 3 + k] - ccoords[k + tc * 3], dim); 380755f3dfbSVijay Mahadevan if (values_phi[tp] < 1e-12) { 381755f3dfbSVijay Mahadevan values_phi[tp] = 1e12; 382b117cd09SVijay Mahadevan } 383b117cd09SVijay Mahadevan else { 384755f3dfbSVijay Mahadevan //values_phi[tp] = std::pow(values_phi[tp], -1.0/dim); 385755f3dfbSVijay Mahadevan values_phi[tp] = std::pow(values_phi[tp], -1.0); 386755f3dfbSVijay Mahadevan normsum += values_phi[tp]; 387b117cd09SVijay Mahadevan } 388b117cd09SVijay Mahadevan } 389755f3dfbSVijay Mahadevan for (unsigned tp = 0; tp < connp.size(); tp++) { 390755f3dfbSVijay Mahadevan if (values_phi[tp] > 1e11) 391755f3dfbSVijay Mahadevan values_phi[tp] = factor * 0.5 / connp.size(); 392b117cd09SVijay Mahadevan else 393755f3dfbSVijay Mahadevan values_phi[tp] = factor * values_phi[tp] * 0.5 / (connp.size() * normsum); 394b117cd09SVijay Mahadevan } 395755f3dfbSVijay Mahadevan ierr = MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi[0], ADD_VALUES);CHKERRQ(ierr); 396b117cd09SVijay Mahadevan } 397*a86ed394SVijay Mahadevan } 398b117cd09SVijay Mahadevan } 399304006b3SVijay Mahadevan if (vec) *vec = NULL; 400b117cd09SVijay Mahadevan ierr = MatAssemblyBegin(*interpl, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 401b117cd09SVijay Mahadevan ierr = MatAssemblyEnd(*interpl, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 402b117cd09SVijay Mahadevan PetscFunctionReturn(0); 403b117cd09SVijay Mahadevan } 404b117cd09SVijay Mahadevan 405b117cd09SVijay Mahadevan #undef __FUNCT__ 406b117cd09SVijay Mahadevan #define __FUNCT__ "DMCreateInjection_Moab" 407b117cd09SVijay Mahadevan /*@ 408b117cd09SVijay Mahadevan DMCreateInjection_Moab - Generate a multi-level uniform refinement hierarchy 409b117cd09SVijay Mahadevan by succesively refining a coarse mesh, already defined in the DM object 410b117cd09SVijay Mahadevan provided by the user. 411b117cd09SVijay Mahadevan 412b117cd09SVijay Mahadevan Collective on MPI_Comm 413b117cd09SVijay Mahadevan 414b117cd09SVijay Mahadevan Input Parameter: 415b117cd09SVijay Mahadevan . dmb - The DMMoab object 416b117cd09SVijay Mahadevan 417b117cd09SVijay Mahadevan Output Parameter: 418b117cd09SVijay Mahadevan . nlevels - The number of levels of refinement needed to generate the hierarchy 419b117cd09SVijay Mahadevan + ldegrees - The degree of refinement at each level in the hierarchy 420b117cd09SVijay Mahadevan 421b117cd09SVijay Mahadevan Level: beginner 422b117cd09SVijay Mahadevan 423b117cd09SVijay Mahadevan .keywords: DMMoab, create, refinement 424b117cd09SVijay Mahadevan @*/ 425304006b3SVijay Mahadevan PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dm1, DM dm2, VecScatter* ctx) 426b117cd09SVijay Mahadevan { 427e882eb38SVijay Mahadevan //DM_Moab *dmmoab; 428b117cd09SVijay Mahadevan 429b117cd09SVijay Mahadevan PetscFunctionBegin; 430b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm1, DM_CLASSID, 1); 431b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm2, DM_CLASSID, 2); 432e882eb38SVijay Mahadevan //dmmoab = (DM_Moab*)(dm1)->data; 433b117cd09SVijay Mahadevan 434b117cd09SVijay Mahadevan PetscPrintf(PETSC_COMM_WORLD, "[DMCreateInjection_Moab] :: Placeholder\n"); 435b117cd09SVijay Mahadevan PetscFunctionReturn(0); 436b117cd09SVijay Mahadevan } 437b117cd09SVijay Mahadevan 438b117cd09SVijay Mahadevan 439b117cd09SVijay Mahadevan #undef __FUNCT__ 440b117cd09SVijay Mahadevan #define __FUNCT__ "DM_UMR_Moab_Private" 441b117cd09SVijay Mahadevan PetscErrorCode DM_UMR_Moab_Private(DM dm, MPI_Comm comm, PetscBool refine, DM *dmref) 442b117cd09SVijay Mahadevan { 443b117cd09SVijay Mahadevan PetscErrorCode ierr; 444e882eb38SVijay Mahadevan PetscInt i, dim; 445b117cd09SVijay Mahadevan DM dm2; 446b117cd09SVijay Mahadevan moab::ErrorCode merr; 447b117cd09SVijay Mahadevan DM_Moab *dmb = (DM_Moab*)dm->data, *dd2; 448b117cd09SVijay Mahadevan 449b117cd09SVijay Mahadevan PetscFunctionBegin; 450b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 451b117cd09SVijay Mahadevan PetscValidPointer(dmref, 4); 452b117cd09SVijay Mahadevan 453b117cd09SVijay Mahadevan if ( (dmb->hlevel == dmb->nhlevels && refine) || (dmb->hlevel == 0 && !refine) ) { 454e882eb38SVijay 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); 455e882eb38SVijay Mahadevan if (dmb->hlevel - 1 < 0 && !refine) PetscInfo1(NULL, "Invalid multigrid coarsen hierarchy level specified (%D). Creating a NULL object.\n", dmb->hlevel - 1); 456e882eb38SVijay Mahadevan *dmref = PETSC_NULL; 457b117cd09SVijay Mahadevan PetscFunctionReturn(0); 458b117cd09SVijay Mahadevan } 459b117cd09SVijay Mahadevan 460b117cd09SVijay Mahadevan ierr = DMMoabCreate(PetscObjectComm((PetscObject)dm), &dm2);CHKERRQ(ierr); 461b117cd09SVijay Mahadevan dd2 = (DM_Moab*)dm2->data; 462b117cd09SVijay Mahadevan 463b117cd09SVijay Mahadevan dd2->mbiface = dmb->mbiface; 4649daf19fdSVijay Mahadevan #ifdef MOAB_HAVE_MPI 465b117cd09SVijay Mahadevan dd2->pcomm = dmb->pcomm; 4669daf19fdSVijay Mahadevan #endif 467b117cd09SVijay Mahadevan dd2->icreatedinstance = PETSC_FALSE; 46864e1c140SVijay Mahadevan dd2->nghostrings = dmb->nghostrings; 469b117cd09SVijay Mahadevan 470b117cd09SVijay Mahadevan /* set the new level based on refinement/coarsening */ 471b117cd09SVijay Mahadevan if (refine) { 472b117cd09SVijay Mahadevan dd2->hlevel = dmb->hlevel + 1; 473b117cd09SVijay Mahadevan } 474b117cd09SVijay Mahadevan else { 475b117cd09SVijay Mahadevan dd2->hlevel = dmb->hlevel - 1; 476b117cd09SVijay Mahadevan } 477b117cd09SVijay Mahadevan 478b117cd09SVijay Mahadevan /* Copy the multilevel hierarchy pointers in MOAB */ 479b117cd09SVijay Mahadevan dd2->hierarchy = dmb->hierarchy; 480b117cd09SVijay Mahadevan dd2->nhlevels = dmb->nhlevels; 481b117cd09SVijay Mahadevan ierr = PetscMalloc1(dd2->nhlevels + 1, &dd2->hsets);CHKERRQ(ierr); 482b117cd09SVijay Mahadevan for (i = 0; i <= dd2->nhlevels; i++) { 483b117cd09SVijay Mahadevan dd2->hsets[i] = dmb->hsets[i]; 484b117cd09SVijay Mahadevan } 485b117cd09SVijay Mahadevan dd2->fileset = dd2->hsets[dd2->hlevel]; 486b117cd09SVijay Mahadevan 487b117cd09SVijay Mahadevan /* do the remaining initializations for DMMoab */ 488b117cd09SVijay Mahadevan dd2->bs = dmb->bs; 489b117cd09SVijay Mahadevan dd2->numFields = dmb->numFields; 490b117cd09SVijay Mahadevan dd2->rw_dbglevel = dmb->rw_dbglevel; 491b117cd09SVijay Mahadevan dd2->partition_by_rank = dmb->partition_by_rank; 492b117cd09SVijay Mahadevan ierr = PetscStrcpy(dd2->extra_read_options, dmb->extra_read_options);CHKERRQ(ierr); 493b117cd09SVijay Mahadevan ierr = PetscStrcpy(dd2->extra_write_options, dmb->extra_write_options);CHKERRQ(ierr); 494b117cd09SVijay Mahadevan dd2->read_mode = dmb->read_mode; 495b117cd09SVijay Mahadevan dd2->write_mode = dmb->write_mode; 496b117cd09SVijay Mahadevan 497b117cd09SVijay Mahadevan /* set global ID tag handle */ 498b117cd09SVijay Mahadevan ierr = DMMoabSetLocalToGlobalTag(dm2, dmb->ltog_tag);CHKERRQ(ierr); 499b117cd09SVijay Mahadevan 500b117cd09SVijay Mahadevan merr = dd2->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dd2->material_tag); MBERRNM(merr); 501b117cd09SVijay Mahadevan 502b117cd09SVijay Mahadevan ierr = DMSetOptionsPrefix(dm2, ((PetscObject)dm)->prefix);CHKERRQ(ierr); 503b117cd09SVijay Mahadevan ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 504b117cd09SVijay Mahadevan ierr = DMSetDimension(dm2, dim);CHKERRQ(ierr); 505b117cd09SVijay Mahadevan 506b117cd09SVijay Mahadevan /* allow overloaded (user replaced) operations to be inherited by refinement clones */ 507b117cd09SVijay Mahadevan dm2->ops->creatematrix = dm->ops->creatematrix; 508b117cd09SVijay Mahadevan 509b117cd09SVijay Mahadevan /* copy fill information if given */ 510b117cd09SVijay Mahadevan ierr = DMMoabSetBlockFills(dm2, dmb->dfill, dmb->ofill);CHKERRQ(ierr); 511b117cd09SVijay Mahadevan 512b117cd09SVijay Mahadevan /* copy vector type information */ 513b117cd09SVijay Mahadevan ierr = DMSetMatType(dm2, dm->mattype);CHKERRQ(ierr); 514b117cd09SVijay Mahadevan ierr = DMSetVecType(dm2, dm->vectype);CHKERRQ(ierr); 5153f1c6e43SVijay Mahadevan dd2->numFields = dmb->numFields; 5163f1c6e43SVijay Mahadevan if (dmb->numFields) { 517b117cd09SVijay Mahadevan ierr = DMMoabSetFieldNames(dm2, dmb->numFields, dmb->fieldNames);CHKERRQ(ierr); 5183f1c6e43SVijay Mahadevan } 519b117cd09SVijay Mahadevan 520b117cd09SVijay Mahadevan ierr = DMSetFromOptions(dm2);CHKERRQ(ierr); 521b117cd09SVijay Mahadevan 522b117cd09SVijay Mahadevan /* recreate Dof numbering for the refined DM and make sure the distribution is correctly populated */ 523b117cd09SVijay Mahadevan ierr = DMSetUp(dm2);CHKERRQ(ierr); 524b117cd09SVijay Mahadevan 525b117cd09SVijay Mahadevan *dmref = dm2; 526b117cd09SVijay Mahadevan PetscFunctionReturn(0); 527b117cd09SVijay Mahadevan } 528b117cd09SVijay Mahadevan 529b117cd09SVijay Mahadevan 530b117cd09SVijay Mahadevan #undef __FUNCT__ 531b117cd09SVijay Mahadevan #define __FUNCT__ "DMRefine_Moab" 532b117cd09SVijay Mahadevan /*@ 533b117cd09SVijay Mahadevan DMRefine_Moab - Generate a multi-level uniform refinement hierarchy 534b117cd09SVijay Mahadevan by succesively refining a coarse mesh, already defined in the DM object 535b117cd09SVijay Mahadevan provided by the user. 536b117cd09SVijay Mahadevan 537e882eb38SVijay Mahadevan Collective on DM 538b117cd09SVijay Mahadevan 539b117cd09SVijay Mahadevan Input Parameter: 540e882eb38SVijay Mahadevan + dm - The DMMoab object 541e882eb38SVijay Mahadevan - comm - the communicator to contain the new DM object (or MPI_COMM_NULL) 542b117cd09SVijay Mahadevan 543b117cd09SVijay Mahadevan Output Parameter: 544e882eb38SVijay Mahadevan . dmf - the refined DM, or NULL 545b117cd09SVijay Mahadevan 546e882eb38SVijay Mahadevan Note: If no refinement was done, the return value is NULL 547e882eb38SVijay Mahadevan 548e882eb38SVijay Mahadevan Level: developer 549b117cd09SVijay Mahadevan 550b117cd09SVijay Mahadevan .keywords: DMMoab, create, refinement 551b117cd09SVijay Mahadevan @*/ 552304006b3SVijay Mahadevan PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM* dmf) 553b117cd09SVijay Mahadevan { 554b117cd09SVijay Mahadevan PetscErrorCode ierr; 555b117cd09SVijay Mahadevan 556b117cd09SVijay Mahadevan PetscFunctionBegin; 557b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 558b117cd09SVijay Mahadevan 559b117cd09SVijay Mahadevan ierr = DM_UMR_Moab_Private(dm, comm, PETSC_TRUE, dmf);CHKERRQ(ierr); 560b117cd09SVijay Mahadevan PetscFunctionReturn(0); 561b117cd09SVijay Mahadevan } 562b117cd09SVijay Mahadevan 563b117cd09SVijay Mahadevan #undef __FUNCT__ 564b117cd09SVijay Mahadevan #define __FUNCT__ "DMCoarsen_Moab" 565b117cd09SVijay Mahadevan /*@ 566b117cd09SVijay Mahadevan DMCoarsen_Moab - Generate a multi-level uniform refinement hierarchy 567b117cd09SVijay Mahadevan by succesively refining a coarse mesh, already defined in the DM object 568b117cd09SVijay Mahadevan provided by the user. 569b117cd09SVijay Mahadevan 570e882eb38SVijay Mahadevan Collective on DM 571b117cd09SVijay Mahadevan 572b117cd09SVijay Mahadevan Input Parameter: 573e882eb38SVijay Mahadevan . dm - The DMMoab object 574e882eb38SVijay Mahadevan - comm - the communicator to contain the new DM object (or MPI_COMM_NULL) 575b117cd09SVijay Mahadevan 576b117cd09SVijay Mahadevan Output Parameter: 577e882eb38SVijay Mahadevan . dmf - the coarsened DM, or NULL 578b117cd09SVijay Mahadevan 579e882eb38SVijay Mahadevan Note: If no coarsening was done, the return value is NULL 580b117cd09SVijay Mahadevan 581e882eb38SVijay Mahadevan Level: developer 582e882eb38SVijay Mahadevan 583e882eb38SVijay Mahadevan .keywords: DMMoab, create, coarsening 584b117cd09SVijay Mahadevan @*/ 585304006b3SVijay Mahadevan PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM* dmc) 586b117cd09SVijay Mahadevan { 587b117cd09SVijay Mahadevan PetscErrorCode ierr; 588b117cd09SVijay Mahadevan 589b117cd09SVijay Mahadevan PetscFunctionBegin; 590b117cd09SVijay Mahadevan PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 591b117cd09SVijay Mahadevan 592b117cd09SVijay Mahadevan ierr = DM_UMR_Moab_Private(dm, comm, PETSC_FALSE, dmc);CHKERRQ(ierr); 593b117cd09SVijay Mahadevan PetscFunctionReturn(0); 594b117cd09SVijay Mahadevan } 595