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;

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

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

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

     71:   Not collective

     73:   Input Parameters:
     74: + s      - the organizing PetscSection
     75: . v      - the Vec
     76: - viewer - the PetscViewer

     78:   Level: developer

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

     91:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)v), &viewer);}
     93:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
     94:   if (isascii) {
     95:     const char *name;

     97:     PetscObjectGetName((PetscObject) v, &name);
     98:     if (s->numFields) {
     99:       PetscViewerASCIIPrintf(viewer, "%s with %D fields\n", name, s->numFields);
    100:       for (f = 0; f < s->numFields; ++f) {
    101:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
    102:         PetscSectionVecView_ASCII(s->field[f], v, viewer);
    103:       }
    104:     } else {
    105:       PetscViewerASCIIPrintf(viewer, "%s\n", name);
    106:       PetscSectionVecView_ASCII(s, v, viewer);
    107:     }
    108:   }
    109:   return(0);
    110: }

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

    115:   Not collective

    117:   Input Parameters:
    118: + v - the Vec
    119: . s - the organizing PetscSection
    120: - point - the point

    122:   Output Parameter:
    123: . values - the array of output values

    125:   Level: developer

    127: .seealso: PetscSection, PetscSectionCreate(), VecSetValuesSection()
    128: @*/
    129: PetscErrorCode VecGetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar **values)
    130: {
    131:   PetscScalar    *baseArray;
    132:   const PetscInt p = point - s->pStart;

    138:   VecGetArray(v, &baseArray);
    139:   *values = &baseArray[s->atlasOff[p]];
    140:   VecRestoreArray(v, &baseArray);
    141:   return(0);
    142: }

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

    147:   Not collective

    149:   Input Parameters:
    150: + v - the Vec
    151: . s - the organizing PetscSection
    152: . point - the point
    153: . values - the array of input values
    154: - mode - the insertion mode, either ADD_VALUES or INSERT_VALUES

    156:   Level: developer

    158:   Note: This is similar to MatSetValuesStencil(). The Fortran binding is
    159: $
    160: $   VecSetValuesSectionF90(vec, section, point, values, mode, ierr)
    161: $

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

    179:   PetscSectionGetConstraintDof(s, point, &cDim);
    180:   VecGetArray(v, &baseArray);
    181:   array = &baseArray[s->atlasOff[p]];
    182:   if (!cDim && doInterior) {
    183:     if (orientation >= 0) {
    184:       const PetscInt dim = s->atlasDof[p];
    185:       PetscInt       i;

    187:       if (doInsert) {
    188:         for (i = 0; i < dim; ++i) array[i] = values[i];
    189:       } else {
    190:         for (i = 0; i < dim; ++i) array[i] += values[i];
    191:       }
    192:     } else {
    193:       PetscInt offset = 0;
    194:       PetscInt j      = -1, field, i;

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

    199:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
    200:         offset += dim;
    201:       }
    202:     }
    203:   } else if (cDim) {
    204:     if (orientation >= 0) {
    205:       const PetscInt dim  = s->atlasDof[p];
    206:       PetscInt       cInd = 0, i;
    207:       const PetscInt *cDof;

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

    236:       PetscSectionGetConstraintIndices(s, point, &cDof);
    237:       for (field = 0; field < s->numFields; ++field) {
    238:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
    239:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
    240:         const PetscInt sDim = dim - tDim;
    241:         PetscInt       cInd = 0, i ,k;

    243:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
    244:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
    245:           if (doInterior) array[j] = values[k];   /* Unconstrained update */
    246:         }
    247:         offset  += dim;
    248:         cOffset += dim - tDim;
    249:       }
    250:     }
    251:   }
    252:   VecRestoreArray(v, &baseArray);
    253:   return(0);
    254: }

    256: PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
    257: {
    258:   PetscInt      *subIndices;
    259:   PetscInt       Nc, subSize = 0, subOff = 0, p;

    263:   PetscSectionGetFieldComponents(section, field, &Nc);
    264:   for (p = pStart; p < pEnd; ++p) {
    265:     PetscInt gdof, fdof = 0;

    267:     PetscSectionGetDof(sectionGlobal, p, &gdof);
    268:     if (gdof > 0) {PetscSectionGetFieldDof(section, p, field, &fdof);}
    269:     subSize += fdof;
    270:   }
    271:   PetscMalloc1(subSize, &subIndices);
    272:   for (p = pStart; p < pEnd; ++p) {
    273:     PetscInt gdof, goff;

    275:     PetscSectionGetDof(sectionGlobal, p, &gdof);
    276:     if (gdof > 0) {
    277:       PetscInt fdof, fc, f2, poff = 0;

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

    295: PetscErrorCode PetscSectionRestoreField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
    296: {

    300:   VecRestoreSubVector(v, *is, subv);
    301:   ISDestroy(is);
    302:   return(0);
    303: }

    305: /*@C
    306:   PetscSectionVecNorm - Computes the vector norm, separated into field components.

    308:   Input Parameters:
    309: + s    - the local Section
    310: . gs   - the global section
    311: . x    - the vector
    312: - type - one of NORM_1, NORM_2, NORM_INFINITY.

    314:   Output Parameter:
    315: . val  - the array of norms

    317:   Level: intermediate

    319: .seealso: VecNorm(), PetscSectionCreate()
    320: @*/
    321: PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
    322: {
    323:   PetscInt       Nf, f, pStart, pEnd;

    331:   PetscSectionGetNumFields(s, &Nf);
    332:   if (Nf < 2) {VecNorm(x, type, val);}
    333:   else {
    334:     PetscSectionGetChart(s, &pStart, &pEnd);
    335:     for (f = 0; f < Nf; ++f) {
    336:       Vec subv;
    337:       IS  is;

    339:       PetscSectionGetField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
    340:       VecNorm(subv, type, &val[f]);
    341:       PetscSectionRestoreField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
    342:     }
    343:   }
    344:   return(0);
    345: }