Actual source code: vsection.c

      1: /*
      2:    This file contains routines for section object operations on Vecs
      3: */
      4: #include <petsc/private/sectionimpl.h>
      5: #include <petsc/private/vecimpl.h>

      7: static PetscErrorCode PetscSectionVecView_ASCII(PetscSection s, Vec v, PetscViewer viewer)
      8: {
      9:   PetscScalar *array;
     10:   PetscInt     p, i;
     11:   PetscMPIInt  rank;

     13:   PetscFunctionBegin;
     14:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
     15:   PetscCall(VecGetArray(v, &array));
     16:   PetscCall(PetscViewerASCIIPushSynchronized(viewer));
     17:   PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
     18:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
     19:     if (s->bc && (s->bc->atlasDof[p] > 0)) {
     20:       PetscInt b;

     22:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT, p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
     23:       for (i = s->atlasOff[p]; i < s->atlasOff[p] + s->atlasDof[p]; ++i) {
     24:         PetscScalar v = array[i];
     25: #if defined(PETSC_USE_COMPLEX)
     26:         if (PetscImaginaryPart(v) > 0.0) {
     27:           PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v)));
     28:         } else if (PetscImaginaryPart(v) < 0.0) {
     29:           PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g - %g i", (double)PetscRealPart(v), (double)(-PetscImaginaryPart(v))));
     30:         } else {
     31:           PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v)));
     32:         }
     33: #else
     34:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v));
     35: #endif
     36:       }
     37:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " constrained"));
     38:       for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
     39:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
     40:     } else {
     41:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT, p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
     42:       for (i = s->atlasOff[p]; i < s->atlasOff[p] + s->atlasDof[p]; ++i) {
     43:         PetscScalar v = array[i];
     44: #if defined(PETSC_USE_COMPLEX)
     45:         if (PetscImaginaryPart(v) > 0.0) {
     46:           PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v)));
     47:         } else if (PetscImaginaryPart(v) < 0.0) {
     48:           PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g - %g i", (double)PetscRealPart(v), (double)(-PetscImaginaryPart(v))));
     49:         } else {
     50:           PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v)));
     51:         }
     52: #else
     53:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v));
     54: #endif
     55:       }
     56:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
     57:     }
     58:   }
     59:   PetscCall(PetscViewerFlush(viewer));
     60:   PetscCall(PetscViewerASCIIPopSynchronized(viewer));
     61:   PetscCall(VecRestoreArray(v, &array));
     62:   PetscFunctionReturn(PETSC_SUCCESS);
     63: }

     65: /*@
     66:   PetscSectionVecView - View a vector, using the section to structure the values

     68:   Not Collective

     70:   Input Parameters:
     71: + s      - the organizing `PetscSection`
     72: . v      - the `Vec`
     73: - viewer - the `PetscViewer`

     75:   Level: developer

     77: .seealso: `PetscSection`, `PetscViewer`, `PetscSectionCreate()`, `VecSetValuesSection()`
     78: @*/
     79: PetscErrorCode PetscSectionVecView(PetscSection s, Vec v, PetscViewer viewer)
     80: {
     81:   PetscBool isascii;
     82:   PetscInt  f;

     84:   PetscFunctionBegin;
     87:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)v), &viewer));
     89:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
     90:   if (isascii) {
     91:     const char *name;

     93:     PetscCall(PetscObjectGetName((PetscObject)v, &name));
     94:     if (s->numFields) {
     95:       PetscCall(PetscViewerASCIIPrintf(viewer, "%s with %" PetscInt_FMT " fields\n", name, s->numFields));
     96:       for (f = 0; f < s->numFields; ++f) {
     97:         PetscCall(PetscViewerASCIIPrintf(viewer, "  field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]));
     98:         PetscCall(PetscSectionVecView_ASCII(s->field[f], v, viewer));
     99:       }
    100:     } else {
    101:       PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", name));
    102:       PetscCall(PetscSectionVecView_ASCII(s, v, viewer));
    103:     }
    104:   }
    105:   PetscFunctionReturn(PETSC_SUCCESS);
    106: }

    108: /*@C
    109:   VecGetValuesSection - Gets all the values associated with a given point, according to the section, in the given `Vec`

    111:   Not Collective

    113:   Input Parameters:
    114: + v     - the `Vec`
    115: . s     - the organizing `PetscSection`
    116: - point - the point

    118:   Output Parameter:
    119: . values - the array of output values

    121:   Level: developer

    123: .seealso: `PetscSection`, `PetscSectionCreate()`, `VecSetValuesSection()`
    124: @*/
    125: PetscErrorCode VecGetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar *values[])
    126: {
    127:   PetscScalar   *baseArray;
    128:   const PetscInt p = point - s->pStart;

    130:   PetscFunctionBegin;
    133:   PetscCall(VecGetArray(v, &baseArray));
    134:   *values = &baseArray[s->atlasOff[p]];
    135:   PetscCall(VecRestoreArray(v, &baseArray));
    136:   PetscFunctionReturn(PETSC_SUCCESS);
    137: }

    139: /*@C
    140:   VecSetValuesSection - Sets all the values associated with a given point, according to the section, in the given `Vec`

    142:   Not Collective

    144:   Input Parameters:
    145: + v      - the `Vec`
    146: . s      - the organizing `PetscSection`
    147: . point  - the point
    148: . values - the array of input values
    149: - mode   - the insertion mode, either `ADD_VALUES` or `INSERT_VALUES`

    151:   Level: developer

    153:   Fortran Notes:
    154:   This is similar to `MatSetValuesStencil()`. The binding is
    155: $   VecSetValuesSectionF90(vec, section, point, values, mode, ierr)

    157: .seealso: `PetscSection`, `PetscSectionCreate()`, `VecGetValuesSection()`
    158: @*/
    159: PetscErrorCode VecSetValuesSection(Vec v, PetscSection s, PetscInt point, const PetscScalar values[], InsertMode mode)
    160: {
    161:   PetscScalar    *baseArray, *array;
    162:   const PetscBool doInsert    = mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
    163:   const PetscBool doInterior  = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_VALUES || mode == ADD_VALUES ? PETSC_TRUE : PETSC_FALSE;
    164:   const PetscBool doBC        = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_BC_VALUES || mode == ADD_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
    165:   const PetscInt  p           = point - s->pStart;
    166:   const PetscInt  orientation = 0; /* Needs to be included for use in closure operations */
    167:   PetscInt        cDim        = 0;

    169:   PetscFunctionBegin;
    172:   PetscCall(PetscSectionGetConstraintDof(s, point, &cDim));
    173:   PetscCall(VecGetArray(v, &baseArray));
    174:   array = &baseArray[s->atlasOff[p]];
    175:   if (!cDim && doInterior) {
    176:     if (orientation >= 0) {
    177:       const PetscInt dim = s->atlasDof[p];
    178:       PetscInt       i;

    180:       if (doInsert) {
    181:         for (i = 0; i < dim; ++i) array[i] = values[i];
    182:       } else {
    183:         for (i = 0; i < dim; ++i) array[i] += values[i];
    184:       }
    185:     } else {
    186:       PetscInt offset = 0;
    187:       PetscInt j      = -1, field, i;

    189:       for (field = 0; field < s->numFields; ++field) {
    190:         const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */

    192:         for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
    193:         offset += dim;
    194:       }
    195:     }
    196:   } else if (cDim) {
    197:     if (orientation >= 0) {
    198:       const PetscInt  dim  = s->atlasDof[p];
    199:       PetscInt        cInd = 0, i;
    200:       const PetscInt *cDof;

    202:       PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
    203:       if (doInsert) {
    204:         for (i = 0; i < dim; ++i) {
    205:           if ((cInd < cDim) && (i == cDof[cInd])) {
    206:             if (doBC) array[i] = values[i]; /* Constrained update */
    207:             ++cInd;
    208:             continue;
    209:           }
    210:           if (doInterior) array[i] = values[i]; /* Unconstrained update */
    211:         }
    212:       } else {
    213:         for (i = 0; i < dim; ++i) {
    214:           if ((cInd < cDim) && (i == cDof[cInd])) {
    215:             if (doBC) array[i] += values[i]; /* Constrained update */
    216:             ++cInd;
    217:             continue;
    218:           }
    219:           if (doInterior) array[i] += values[i]; /* Unconstrained update */
    220:         }
    221:       }
    222:     } else {
    223:       /* TODO This is broken for add and constrained update */
    224:       const PetscInt *cDof;
    225:       PetscInt        offset  = 0;
    226:       PetscInt        cOffset = 0;
    227:       PetscInt        j       = 0, field;

    229:       PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
    230:       for (field = 0; field < s->numFields; ++field) {
    231:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
    232:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
    233:         const PetscInt sDim = dim - tDim;
    234:         PetscInt       cInd = 0, i, k;

    236:         for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
    237:           if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
    238:             ++cInd;
    239:             continue;
    240:           }
    241:           if (doInterior) array[j] = values[k]; /* Unconstrained update */
    242:         }
    243:         offset += dim;
    244:         cOffset += dim - tDim;
    245:       }
    246:     }
    247:   }
    248:   PetscCall(VecRestoreArray(v, &baseArray));
    249:   PetscFunctionReturn(PETSC_SUCCESS);
    250: }

    252: PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
    253: {
    254:   PetscInt *subIndices;
    255:   PetscInt  Nc, subSize = 0, subOff = 0, p;

    257:   PetscFunctionBegin;
    258:   PetscCall(PetscSectionGetFieldComponents(section, field, &Nc));
    259:   for (p = pStart; p < pEnd; ++p) {
    260:     PetscInt gdof, fdof = 0;

    262:     PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
    263:     if (gdof > 0) PetscCall(PetscSectionGetFieldDof(section, p, field, &fdof));
    264:     subSize += fdof;
    265:   }
    266:   PetscCall(PetscMalloc1(subSize, &subIndices));
    267:   for (p = pStart; p < pEnd; ++p) {
    268:     PetscInt gdof, goff;

    270:     PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
    271:     if (gdof > 0) {
    272:       PetscInt fdof, fc, f2, poff = 0;

    274:       PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
    275:       /* Can get rid of this loop by storing field information in the global section */
    276:       for (f2 = 0; f2 < field; ++f2) {
    277:         PetscCall(PetscSectionGetFieldDof(section, p, f2, &fdof));
    278:         poff += fdof;
    279:       }
    280:       PetscCall(PetscSectionGetFieldDof(section, p, field, &fdof));
    281:       for (fc = 0; fc < fdof; ++fc, ++subOff) subIndices[subOff] = goff + poff + fc;
    282:     }
    283:   }
    284:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)v), subSize, subIndices, PETSC_OWN_POINTER, is));
    285:   PetscCall(VecGetSubVector(v, *is, subv));
    286:   PetscCall(VecSetBlockSize(*subv, Nc));
    287:   PetscFunctionReturn(PETSC_SUCCESS);
    288: }

    290: PetscErrorCode PetscSectionRestoreField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
    291: {
    292:   PetscFunctionBegin;
    293:   PetscCall(VecRestoreSubVector(v, *is, subv));
    294:   PetscCall(ISDestroy(is));
    295:   PetscFunctionReturn(PETSC_SUCCESS);
    296: }

    298: /*@
    299:   PetscSectionVecNorm - Computes the vector norm of each field

    301:   Input Parameters:
    302: + s    - the local Section
    303: . gs   - the global section
    304: . x    - the vector
    305: - type - one of `NORM_1`, `NORM_2`, `NORM_INFINITY`.

    307:   Output Parameter:
    308: . val - the array of norms

    310:   Level: intermediate

    312: .seealso: `VecNorm()`, `PetscSectionCreate()`
    313: @*/
    314: PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
    315: {
    316:   PetscInt Nf, f, pStart, pEnd;

    318:   PetscFunctionBegin;
    322:   PetscAssertPointer(val, 5);
    323:   PetscCall(PetscSectionGetNumFields(s, &Nf));
    324:   if (Nf < 2) PetscCall(VecNorm(x, type, val));
    325:   else {
    326:     PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
    327:     for (f = 0; f < Nf; ++f) {
    328:       Vec subv;
    329:       IS  is;

    331:       PetscCall(PetscSectionGetField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv));
    332:       PetscCall(VecNorm(subv, type, &val[f]));
    333:       PetscCall(PetscSectionRestoreField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv));
    334:     }
    335:   }
    336:   PetscFunctionReturn(PETSC_SUCCESS);
    337: }