1: #include <../src/ts/characteristic/impls/da/slda.h>
  2: #include <petscdmda.h>
  3: #include <petscviewer.h>

  5: PetscErrorCode CharacteristicView_DA(Characteristic c, PetscViewer viewer)
  6: {
  7:   Characteristic_DA *da = (Characteristic_DA*) c->data;
  8:   PetscBool         iascii, isstring;
  9:   PetscErrorCode    ierr;

 12:   /* Pull out field names from DM */
 13:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
 14:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERSTRING, &isstring);
 15:   if (iascii) {
 16:     PetscViewerASCIIPrintf(viewer,"  DMDA: dummy=%D\n", da->dummy);
 17:   } else if (isstring) {
 18:     PetscViewerStringSPrintf(viewer,"dummy %D", da->dummy);
 19:   }
 20:   return(0);
 21: }

 23: PetscErrorCode CharacteristicDestroy_DA(Characteristic c)
 24: {
 25:   Characteristic_DA *da = (Characteristic_DA*) c->data;
 26:   PetscErrorCode    ierr;

 29:   PetscFree(da);
 30:   return(0);
 31: }

 33: PetscErrorCode CharacteristicSetUp_DA(Characteristic c)
 34: {
 35:   PetscMPIInt    blockLen[2];
 36:   MPI_Aint       indices[2];
 37:   MPI_Datatype   oldtypes[2];
 38:   PetscInt       dim, numValues;

 41:   DMDAGetInfo(c->velocityDA, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
 42:   if (c->structured) c->numIds = dim;
 43:   else c->numIds = 3;
 44:   if (c->numFieldComp > MAX_COMPONENTS) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "The maximum number of fields allowed is %d, you have %d. You can recompile after increasing MAX_COMPONENTS.", MAX_COMPONENTS, c->numFieldComp);
 45:   numValues = 4 + MAX_COMPONENTS;

 47:   /* Create new MPI datatype for communication of characteristic point structs */
 48:   blockLen[0] = 1+c->numIds; indices[0] = 0;                              oldtypes[0] = MPIU_INT;
 49:   blockLen[1] = numValues;   indices[1] = (1+c->numIds)*sizeof(PetscInt); oldtypes[1] = MPIU_SCALAR;
 50:   MPI_Type_create_struct(2, blockLen, indices, oldtypes, &c->itemType);
 51:   MPI_Type_commit(&c->itemType);

 53:   /* Initialize the local queue for char foot values */
 54:   VecGetLocalSize(c->velocity, &c->queueMax);
 55:   PetscMalloc1(c->queueMax, &c->queue);
 56:   c->queueSize = 0;

 58:   /* Allocate communication structures */
 59:   if (c->numNeighbors <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Invalid number of neighbors %d. Call CharactersiticSetNeighbors() before setup.", c->numNeighbors);
 60:   PetscMalloc1(c->numNeighbors, &c->needCount);
 61:   PetscMalloc1(c->numNeighbors, &c->localOffsets);
 62:   PetscMalloc1(c->numNeighbors, &c->fillCount);
 63:   PetscMalloc1(c->numNeighbors, &c->remoteOffsets);
 64:   PetscMalloc1(c->numNeighbors-1, &c->request);
 65:   PetscMalloc1(c->numNeighbors-1,  &c->status);
 66:   return(0);
 67: }

 69: PETSC_EXTERN PetscErrorCode CharacteristicCreate_DA(Characteristic c)
 70: {
 71:   Characteristic_DA *da;
 72:   PetscErrorCode    ierr;

 75:   PetscNew(&da);
 76:   PetscMemzero(da, sizeof(Characteristic_DA));
 77:   PetscLogObjectMemory((PetscObject)c, sizeof(Characteristic_DA));
 78:   c->data = (void*) da;

 80:   c->ops->setup   = CharacteristicSetUp_DA;
 81:   c->ops->destroy = CharacteristicDestroy_DA;
 82:   c->ops->view    = CharacteristicView_DA;

 84:   da->dummy = 0;
 85:   return(0);
 86: }

 88: /* -----------------------------------------------------------------------------
 89:    Checks for periodicity of a DM and Maps points outside of a domain back onto the domain
 90:    using appropriate periodicity. At the moment assumes only a 2-D DMDA.
 91:    ----------------------------------------------------------------------------------------*/
 92: PetscErrorCode DMDAMapCoordsToPeriodicDomain(DM da, PetscScalar *x, PetscScalar *y)
 93: {
 94:   DMBoundaryType bx, by;
 95:   PetscInt       dim, gx, gy;

 99:   DMDAGetInfo(da, &dim, &gx, &gy, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL);

101:   if (bx == DM_BOUNDARY_PERIODIC) {
102:       while (*x >= (PetscScalar)gx) *x -= (PetscScalar)gx;
103:       while (*x < 0.0)              *x += (PetscScalar)gx;
104:     }
105:     if (by == DM_BOUNDARY_PERIODIC) {
106:       while (*y >= (PetscScalar)gy) *y -= (PetscScalar)gy;
107:       while (*y < 0.0)              *y += (PetscScalar)gy;
108:     }
109:   return(0);
110: }