Actual source code: swarmpic_da.c

  1: #include <petscdm.h>
  2: #include <petscdmda.h>
  3: #include <petsc/private/dmswarmimpl.h>
  4: #include "../src/dm/impls/swarm/data_bucket.h"

  6: static PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(PetscInt dim, PetscInt np[], PetscInt *_npoints, PetscReal **_xi)
  7: {
  8:   PetscReal *xi;
  9:   PetscInt   d, npoints = 0, cnt;
 10:   PetscReal  ds[] = {0.0, 0.0, 0.0};
 11:   PetscInt   ii, jj, kk;

 13:   PetscFunctionBegin;
 14:   switch (dim) {
 15:   case 1:
 16:     npoints = np[0];
 17:     break;
 18:   case 2:
 19:     npoints = np[0] * np[1];
 20:     break;
 21:   case 3:
 22:     npoints = np[0] * np[1] * np[2];
 23:     break;
 24:   }
 25:   for (d = 0; d < dim; d++) ds[d] = 2.0 / ((PetscReal)np[d]);

 27:   PetscCall(PetscMalloc1(dim * npoints, &xi));
 28:   switch (dim) {
 29:   case 1:
 30:     cnt = 0;
 31:     for (ii = 0; ii < np[0]; ii++) {
 32:       xi[dim * cnt + 0] = -1.0 + 0.5 * ds[d] + ii * ds[0];
 33:       cnt++;
 34:     }
 35:     break;

 37:   case 2:
 38:     cnt = 0;
 39:     for (jj = 0; jj < np[1]; jj++) {
 40:       for (ii = 0; ii < np[0]; ii++) {
 41:         xi[dim * cnt + 0] = -1.0 + 0.5 * ds[0] + ii * ds[0];
 42:         xi[dim * cnt + 1] = -1.0 + 0.5 * ds[1] + jj * ds[1];
 43:         cnt++;
 44:       }
 45:     }
 46:     break;

 48:   case 3:
 49:     cnt = 0;
 50:     for (kk = 0; kk < np[2]; kk++) {
 51:       for (jj = 0; jj < np[1]; jj++) {
 52:         for (ii = 0; ii < np[0]; ii++) {
 53:           xi[dim * cnt + 0] = -1.0 + 0.5 * ds[0] + ii * ds[0];
 54:           xi[dim * cnt + 1] = -1.0 + 0.5 * ds[1] + jj * ds[1];
 55:           xi[dim * cnt + 2] = -1.0 + 0.5 * ds[2] + kk * ds[2];
 56:           cnt++;
 57:         }
 58:       }
 59:     }
 60:     break;
 61:   }
 62:   *_npoints = npoints;
 63:   *_xi      = xi;
 64:   PetscFunctionReturn(PETSC_SUCCESS);
 65: }

 67: static PetscErrorCode private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(PetscInt dim, PetscInt np_1d, PetscInt *_npoints, PetscReal **_xi)
 68: {
 69:   PetscQuadrature  quadrature;
 70:   const PetscReal *quadrature_xi;
 71:   PetscReal       *xi;
 72:   PetscInt         d, q, npoints_q;

 74:   PetscFunctionBegin;
 75:   PetscCall(PetscDTGaussTensorQuadrature(dim, 1, np_1d, -1.0, 1.0, &quadrature));
 76:   PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &quadrature_xi, NULL));
 77:   PetscCall(PetscMalloc1(dim * npoints_q, &xi));
 78:   for (q = 0; q < npoints_q; q++) {
 79:     for (d = 0; d < dim; d++) xi[dim * q + d] = quadrature_xi[dim * q + d];
 80:   }
 81:   PetscCall(PetscQuadratureDestroy(&quadrature));
 82:   *_npoints = npoints_q;
 83:   *_xi      = xi;
 84:   PetscFunctionReturn(PETSC_SUCCESS);
 85: }

 87: static PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA_Q1(DM dm, DM dmc, PetscInt npoints, DMSwarmPICLayoutType layout)
 88: {
 89:   DMSwarmCellDM      celldm;
 90:   PetscInt           dim, npoints_q;
 91:   PetscInt           nel, npe, e, q, k, d;
 92:   const PetscInt    *element_list;
 93:   PetscReal        **basis;
 94:   PetscReal         *xi;
 95:   Vec                coor;
 96:   const PetscScalar *_coor;
 97:   PetscReal         *elcoor;
 98:   PetscReal         *swarm_coor;
 99:   PetscInt          *swarm_cellid;
100:   PetscInt           pcnt, Nfc;
101:   const char       **coordFields, *cellid;

103:   PetscFunctionBegin;
104:   PetscCall(DMGetDimension(dm, &dim));
105:   switch (layout) {
106:   case DMSWARMPIC_LAYOUT_REGULAR: {
107:     PetscInt np_dir[3];
108:     np_dir[0] = np_dir[1] = np_dir[2] = npoints;
109:     PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim, np_dir, &npoints_q, &xi));
110:   } break;
111:   case DMSWARMPIC_LAYOUT_GAUSS:
112:     PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Gauss(dim, npoints, &npoints_q, &xi));
113:     break;

115:   case DMSWARMPIC_LAYOUT_SUBDIVISION: {
116:     PetscInt s, nsub;
117:     PetscInt np_dir[3];
118:     nsub      = npoints;
119:     np_dir[0] = 1;
120:     for (s = 0; s < nsub; s++) np_dir[0] *= 2;
121:     np_dir[1] = np_dir[0];
122:     np_dir[2] = np_dir[0];
123:     PetscCall(private_DMSwarmCreateCellLocalCoords_DA_Q1_Regular(dim, np_dir, &npoints_q, &xi));
124:   } break;
125:   default:
126:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "A valid DMSwarmPIC layout must be provided");
127:   }

129:   PetscCall(DMDAGetElements(dmc, &nel, &npe, &element_list));
130:   PetscCall(PetscMalloc1(dim * npe, &elcoor));
131:   PetscCall(PetscMalloc1(npoints_q, &basis));
132:   for (q = 0; q < npoints_q; q++) {
133:     PetscCall(PetscMalloc1(npe, &basis[q]));

135:     switch (dim) {
136:     case 1:
137:       basis[q][0] = 0.5 * (1.0 - xi[dim * q + 0]);
138:       basis[q][1] = 0.5 * (1.0 + xi[dim * q + 0]);
139:       break;
140:     case 2:
141:       basis[q][0] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
142:       basis[q][1] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
143:       basis[q][2] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
144:       basis[q][3] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
145:       break;

147:     case 3:
148:       basis[q][0] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
149:       basis[q][1] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
150:       basis[q][2] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
151:       basis[q][3] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
152:       basis[q][4] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
153:       basis[q][5] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
154:       basis[q][6] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
155:       basis[q][7] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
156:       break;
157:     }
158:   }

160:   PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
161:   PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
162:   PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
163:   PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));

165:   PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1));
166:   PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
167:   PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));

169:   PetscCall(DMGetCoordinatesLocal(dmc, &coor));
170:   PetscCall(VecGetArrayRead(coor, &_coor));
171:   pcnt = 0;
172:   for (e = 0; e < nel; e++) {
173:     const PetscInt *element = &element_list[npe * e];

175:     for (k = 0; k < npe; k++) {
176:       for (d = 0; d < dim; d++) elcoor[dim * k + d] = PetscRealPart(_coor[dim * element[k] + d]);
177:     }

179:     for (q = 0; q < npoints_q; q++) {
180:       for (d = 0; d < dim; d++) swarm_coor[dim * pcnt + d] = 0.0;
181:       for (k = 0; k < npe; k++) {
182:         for (d = 0; d < dim; d++) swarm_coor[dim * pcnt + d] += basis[q][k] * elcoor[dim * k + d];
183:       }
184:       swarm_cellid[pcnt] = e;
185:       pcnt++;
186:     }
187:   }
188:   PetscCall(VecRestoreArrayRead(coor, &_coor));
189:   PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
190:   PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
191:   PetscCall(DMDARestoreElements(dmc, &nel, &npe, &element_list));

193:   PetscCall(PetscFree(xi));
194:   PetscCall(PetscFree(elcoor));
195:   for (q = 0; q < npoints_q; q++) PetscCall(PetscFree(basis[q]));
196:   PetscCall(PetscFree(basis));
197:   PetscFunctionReturn(PETSC_SUCCESS);
198: }

200: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param)
201: {
202:   DMDAElementType etype;
203:   PetscInt        dim;

205:   PetscFunctionBegin;
206:   PetscCall(DMDAGetElementType(celldm, &etype));
207:   PetscCall(DMGetDimension(celldm, &dim));
208:   switch (etype) {
209:   case DMDA_ELEMENT_P1:
210:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DA support is not currently available for DMDA_ELEMENT_P1");
211:   case DMDA_ELEMENT_Q1:
212:     PetscCheck(dim != 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Support only available for dim = 2, 3");
213:     PetscCall(private_DMSwarmInsertPointsUsingCellDM_DA_Q1(dm, celldm, layout_param, layout));
214:     break;
215:   }
216:   PetscFunctionReturn(PETSC_SUCCESS);
217: }