Actual source code: plexnatural.c

      1: #include <petsc/private/dmpleximpl.h>

      3: /*@
      4:   DMPlexSetMigrationSF - Sets the SF for migrating from a parent DM into this DM

      6:   Input Parameters:
      7: + dm        - The DM
      8: - naturalSF - The PetscSF

     10:   Note: It is necessary to call this in order to have DMCreateSubDM() or DMCreateSuperDM() build the Global-To-Natural map

     12:   Level: intermediate

     14: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateMigrationSF(), DMPlexGetMigrationSF()
     15: @*/
     16: PetscErrorCode DMPlexSetMigrationSF(DM dm, PetscSF migrationSF)
     17: {
     20:   dm->sfMigration = migrationSF;
     21:   PetscObjectReference((PetscObject) migrationSF);
     22:   return(0);
     23: }

     25: /*@
     26:   DMPlexGetMigrationSF - Gets the SF for migrating from a parent DM into this DM

     28:   Input Parameter:
     29: . dm          - The DM

     31:   Output Parameter:
     32: . migrationSF - The PetscSF

     34:   Level: intermediate

     36: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateMigrationSF(), DMPlexSetMigrationSF
     37: @*/
     38: PetscErrorCode DMPlexGetMigrationSF(DM dm, PetscSF *migrationSF)
     39: {
     41:   *migrationSF = dm->sfMigration;
     42:   return(0);
     43: }

     45: /*@
     46:   DMPlexSetGlobalToNaturalSF - Sets the SF for mapping Global Vec to the Natural Vec

     48:   Input Parameters:
     49: + dm          - The DM
     50: - naturalSF   - The PetscSF

     52:   Level: intermediate

     54: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexGetGlobaltoNaturalSF()
     55: @*/
     56: PetscErrorCode DMPlexSetGlobalToNaturalSF(DM dm, PetscSF naturalSF)
     57: {
     60:   dm->sfNatural = naturalSF;
     61:   PetscObjectReference((PetscObject) naturalSF);
     62:   dm->useNatural = PETSC_TRUE;
     63:   return(0);
     64: }

     66: /*@
     67:   DMPlexGetGlobalToNaturalSF - Gets the SF for mapping Global Vec to the Natural Vec

     69:   Input Parameter:
     70: . dm          - The DM

     72:   Output Parameter:
     73: . naturalSF   - The PetscSF

     75:   Level: intermediate

     77: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexSetGlobaltoNaturalSF
     78: @*/
     79: PetscErrorCode DMPlexGetGlobalToNaturalSF(DM dm, PetscSF *naturalSF)
     80: {
     82:   *naturalSF = dm->sfNatural;
     83:   return(0);
     84: }

     86: /*@
     87:   DMPlexCreateGlobalToNaturalSF - Creates the SF for mapping Global Vec to the Natural Vec

     89:   Input Parameters:
     90: + dm          - The DM
     91: . section     - The PetscSection describing the Vec before the mesh was distributed,
     92:                 or NULL if not available
     93: - sfMigration - The PetscSF used to distribute the mesh, or NULL if it cannot be computed

     95:   Output Parameter:
     96: . sfNatural   - PetscSF for mapping the Vec in PETSc ordering to the canonical ordering

     98:   Note: This is not typically called by the user.

    100:   Level: intermediate

    102: .seealso: DMPlexDistribute(), DMPlexDistributeField()
    103:  @*/
    104: PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural)
    105: {
    106:   MPI_Comm       comm;
    107:   Vec            gv, tmpVec;
    108:   PetscSF        sf, sfEmbed, sfSeqToNatural, sfField, sfFieldInv;
    109:   PetscSection   gSection, sectionDist, gLocSection;
    110:   PetscInt      *spoints, *remoteOffsets;
    111:   PetscInt       ssize, pStart, pEnd, p, globalSize;
    112:   PetscLayout    map;
    113:   PetscBool      destroyFlag = PETSC_FALSE;

    117:   PetscObjectGetComm((PetscObject) dm, &comm);
    118:   if (!sfMigration) {
    119:     /* If sfMigration is missing,
    120:     sfNatural cannot be computed and is set to NULL */
    121:     *sfNatural = NULL;
    122:     return(0);
    123:   } else if (!section) {
    124:     /* If the sequential section is not provided (NULL),
    125:     it is reconstructed from the parallel section */
    126:     PetscSF sfMigrationInv;
    127:     PetscSection localSection;

    129:     DMGetLocalSection(dm, &localSection);
    130:     PetscSFCreateInverseSF(sfMigration, &sfMigrationInv);
    131:     PetscSectionCreate(PetscObjectComm((PetscObject) dm), &section);
    132:     PetscSFDistributeSection(sfMigrationInv, localSection, NULL, section);
    133:     PetscSFDestroy(&sfMigrationInv);
    134:     destroyFlag = PETSC_TRUE;
    135:   }
    136:   /* PetscPrintf(comm, "Point migration SF\n");
    137:    PetscSFView(sfMigration, 0); */
    138:   /* Create a new section from distributing the original section */
    139:   PetscSectionCreate(comm, &sectionDist);
    140:   PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);
    141:   /* PetscPrintf(comm, "Distributed Section\n");
    142:    PetscSectionView(sectionDist, PETSC_VIEWER_STDOUT_WORLD); */
    143:   DMSetLocalSection(dm, sectionDist);
    144:   /* If a sequential section is provided but no dof is affected,
    145:   sfNatural cannot be computed and is set to NULL */
    146:   DMCreateGlobalVector(dm, &tmpVec);
    147:   VecGetSize(tmpVec, &globalSize);
    148:   DMRestoreGlobalVector(dm, &tmpVec);
    149:   if (globalSize) {
    150:   /* Get a pruned version of migration SF */
    151:     DMGetGlobalSection(dm, &gSection);
    152:     PetscSectionGetChart(gSection, &pStart, &pEnd);
    153:     for (p = pStart, ssize = 0; p < pEnd; ++p) {
    154:       PetscInt dof, off;

    156:       PetscSectionGetDof(gSection, p, &dof);
    157:       PetscSectionGetOffset(gSection, p, &off);
    158:       if ((dof > 0) && (off >= 0)) ++ssize;
    159:     }
    160:     PetscMalloc1(ssize, &spoints);
    161:     for (p = pStart, ssize = 0; p < pEnd; ++p) {
    162:       PetscInt dof, off;

    164:       PetscSectionGetDof(gSection, p, &dof);
    165:       PetscSectionGetOffset(gSection, p, &off);
    166:       if ((dof > 0) && (off >= 0)) spoints[ssize++] = p;
    167:     }
    168:     PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);
    169:     PetscFree(spoints);
    170:     /* PetscPrintf(comm, "Embedded SF\n");
    171:     PetscSFView(sfEmbed, 0); */
    172:     /* Create the SF for seq to natural */
    173:     DMGetGlobalVector(dm, &gv);
    174:     VecGetLayout(gv,&map);
    175:     /* Note that entries of gv are leaves in sfSeqToNatural, entries of the seq vec are roots */
    176:     PetscSFCreate(comm, &sfSeqToNatural);
    177:     PetscSFSetGraphWithPattern(sfSeqToNatural, map, PETSCSF_PATTERN_GATHER);
    178:     DMRestoreGlobalVector(dm, &gv);
    179:     /* PetscPrintf(comm, "Seq-to-Natural SF\n");
    180:     PetscSFView(sfSeqToNatural, 0); */
    181:     /* Create the SF associated with this section */
    182:     DMGetPointSF(dm, &sf);
    183:     PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection);
    184:     PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);
    185:     PetscSFDestroy(&sfEmbed);
    186:     PetscSectionDestroy(&gLocSection);
    187:     /* PetscPrintf(comm, "Field SF\n");
    188:     PetscSFView(sfField, 0); */
    189:     /* Invert the field SF so it's now from distributed to sequential */
    190:     PetscSFCreateInverseSF(sfField, &sfFieldInv);
    191:     PetscSFDestroy(&sfField);
    192:     /* PetscPrintf(comm, "Inverse Field SF\n");
    193:     PetscSFView(sfFieldInv, 0); */
    194:     /* Multiply the sfFieldInv with the */
    195:     PetscSFComposeInverse(sfFieldInv, sfSeqToNatural, sfNatural);
    196:     PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view");
    197:     /* Clean up */
    198:     PetscSFDestroy(&sfFieldInv);
    199:     PetscSFDestroy(&sfSeqToNatural);
    200:   } else {
    201:     *sfNatural = NULL;
    202:   }
    203:   PetscSectionDestroy(&sectionDist);
    204:   PetscFree(remoteOffsets);
    205:   if (destroyFlag) {PetscSectionDestroy(&section);}
    206:   return(0);
    207: }

    209: /*@
    210:   DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.

    212:   Collective on dm

    214:   Input Parameters:
    215: + dm - The distributed DMPlex
    216: - gv - The global Vec

    218:   Output Parameters:
    219: . nv - Vec in the canonical ordering distributed over all processors associated with gv

    221:   Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().

    223:   Level: intermediate

    225: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalEnd()
    226: @*/
    227: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
    228: {
    229:   const PetscScalar *inarray;
    230:   PetscScalar       *outarray;
    231:   PetscMPIInt        size;
    232:   PetscErrorCode     ierr;

    235:   PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
    236:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
    237:   if (dm->sfNatural) {
    238:     VecGetArray(nv, &outarray);
    239:     VecGetArrayRead(gv, &inarray);
    240:     PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray,MPI_REPLACE);
    241:     VecRestoreArrayRead(gv, &inarray);
    242:     VecRestoreArray(nv, &outarray);
    243:   } else if (size == 1) {
    244:     VecCopy(gv, nv);
    245:   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.\n");
    246:   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
    247:   PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
    248:   return(0);
    249: }

    251: /*@
    252:   DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.

    254:   Collective on dm

    256:   Input Parameters:
    257: + dm - The distributed DMPlex
    258: - gv - The global Vec

    260:   Output Parameters:
    261: . nv - The natural Vec

    263:   Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().

    265:   Level: intermediate

    267:  .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
    268:  @*/
    269: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
    270: {
    271:   const PetscScalar *inarray;
    272:   PetscScalar       *outarray;
    273:   PetscMPIInt        size;
    274:   PetscErrorCode     ierr;

    277:   PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
    278:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
    279:   if (dm->sfNatural) {
    280:     VecGetArrayRead(gv, &inarray);
    281:     VecGetArray(nv, &outarray);
    282:     PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray,MPI_REPLACE);
    283:     VecRestoreArrayRead(gv, &inarray);
    284:     VecRestoreArray(nv, &outarray);
    285:   } else if (size == 1) {
    286:   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.\n");
    287:   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
    288:   PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
    289:   return(0);
    290: }

    292: /*@
    293:   DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order.

    295:   Collective on dm

    297:   Input Parameters:
    298: + dm - The distributed DMPlex
    299: - nv - The natural Vec

    301:   Output Parameters:
    302: . gv - The global Vec

    304:   Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().

    306:   Level: intermediate

    308: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(),DMPlexGlobalToNaturalEnd()
    309: @*/
    310: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
    311: {
    312:   const PetscScalar *inarray;
    313:   PetscScalar       *outarray;
    314:   PetscMPIInt        size;
    315:   PetscErrorCode     ierr;

    318:   PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
    319:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
    320:   if (dm->sfNatural) {
    321:     /* We only have acces to the SF that goes from Global to Natural.
    322:        Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
    323:        Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
    324:     VecZeroEntries(gv);
    325:     VecGetArray(gv, &outarray);
    326:     VecGetArrayRead(nv, &inarray);
    327:     PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
    328:     VecRestoreArrayRead(nv, &inarray);
    329:     VecRestoreArray(gv, &outarray);
    330:   } else if (size == 1) {
    331:     VecCopy(nv, gv);
    332:   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.\n");
    333:   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
    334:   PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
    335:   return(0);
    336: }

    338: /*@
    339:   DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order.

    341:   Collective on dm

    343:   Input Parameters:
    344: + dm - The distributed DMPlex
    345: - nv - The natural Vec

    347:   Output Parameters:
    348: . gv - The global Vec

    350:   Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().

    352:   Level: intermediate

    354: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
    355:  @*/
    356: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
    357: {
    358:   const PetscScalar *inarray;
    359:   PetscScalar       *outarray;
    360:   PetscMPIInt        size;
    361:   PetscErrorCode     ierr;

    364:   PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
    365:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
    366:   if (dm->sfNatural) {
    367:     VecGetArrayRead(nv, &inarray);
    368:     VecGetArray(gv, &outarray);
    369:     PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
    370:     VecRestoreArrayRead(nv, &inarray);
    371:     VecRestoreArray(gv, &outarray);
    372:   } else if (size == 1) {
    373:   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov.\n");
    374:   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
    375:   PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
    376:   return(0);
    377: }