Actual source code: sectionhdf5.c

  1: #include <petsc/private/sectionimpl.h>
  2: #include <petscsf.h>
  3: #include <petscis.h>
  4: #include <petscviewerhdf5.h>
  5: #include <petsclayouthdf5.h>

  7: #if defined(PETSC_HAVE_HDF5)
  8: static PetscErrorCode PetscSectionView_HDF5_SingleField(PetscSection s, PetscViewer viewer)
  9: {
 10:   MPI_Comm  comm;
 11:   PetscInt  pStart, pEnd, p, n;
 12:   PetscBool hasConstraints, includesConstraints;
 13:   IS        dofIS, offIS, cdofIS, coffIS, cindIS;
 14:   PetscInt *dofs, *offs, *cdofs, *coffs, *cinds, dof, cdof, m, moff, i;

 16:   PetscObjectGetComm((PetscObject)s, &comm);
 17:   PetscSectionGetChart(s, &pStart, &pEnd);
 18:   hasConstraints = (s->bc) ? PETSC_TRUE : PETSC_FALSE;
 19:   MPIU_Allreduce(MPI_IN_PLACE, &hasConstraints, 1, MPIU_BOOL, MPI_LOR, comm);
 20:   for (p = pStart, n = 0, m = 0; p < pEnd; ++p) {
 21:     PetscSectionGetDof(s, p, &dof);
 22:     if (dof >= 0) {
 23:       if (hasConstraints) {
 24:         PetscSectionGetConstraintDof(s, p, &cdof);
 25:         m += cdof;
 26:       }
 27:       n++;
 28:     }
 29:   }
 30:   PetscMalloc1(n, &dofs);
 31:   PetscMalloc1(n, &offs);
 32:   if (hasConstraints) {
 33:     PetscMalloc1(n, &cdofs);
 34:     PetscMalloc1(n, &coffs);
 35:     PetscMalloc1(m, &cinds);
 36:   }
 37:   for (p = pStart, n = 0, m = 0; p < pEnd; ++p) {
 38:     PetscSectionGetDof(s, p, &dof);
 39:     if (dof >= 0) {
 40:       dofs[n] = dof;
 41:       PetscSectionGetOffset(s, p, &offs[n]);
 42:       if (hasConstraints) {
 43:         const PetscInt *cpinds;

 45:         PetscSectionGetConstraintDof(s, p, &cdof);
 46:         PetscSectionGetConstraintIndices(s, p, &cpinds);
 47:         cdofs[n] = cdof;
 48:         coffs[n] = m;
 49:         for (i = 0; i < cdof; ++i) cinds[m++] = cpinds[i];
 50:       }
 51:       n++;
 52:     }
 53:   }
 54:   if (hasConstraints) {
 55:     MPI_Scan(&m, &moff, 1, MPIU_INT, MPI_SUM, comm);
 56:     moff -= m;
 57:     for (p = 0; p < n; ++p) coffs[p] += moff;
 58:   }
 59:   PetscViewerHDF5WriteAttribute(viewer, NULL, "hasConstraints", PETSC_BOOL, (void *)&hasConstraints);
 60:   PetscSectionGetIncludesConstraints(s, &includesConstraints);
 61:   PetscViewerHDF5WriteAttribute(viewer, NULL, "includesConstraints", PETSC_BOOL, (void *)&includesConstraints);
 62:   ISCreateGeneral(comm, n, dofs, PETSC_OWN_POINTER, &dofIS);
 63:   PetscObjectSetName((PetscObject)dofIS, "atlasDof");
 64:   ISView(dofIS, viewer);
 65:   ISDestroy(&dofIS);
 66:   ISCreateGeneral(comm, n, offs, PETSC_OWN_POINTER, &offIS);
 67:   PetscObjectSetName((PetscObject)offIS, "atlasOff");
 68:   ISView(offIS, viewer);
 69:   ISDestroy(&offIS);
 70:   if (hasConstraints) {
 71:     PetscViewerHDF5PushGroup(viewer, "bc");
 72:     ISCreateGeneral(comm, n, cdofs, PETSC_OWN_POINTER, &cdofIS);
 73:     PetscObjectSetName((PetscObject)cdofIS, "atlasDof");
 74:     ISView(cdofIS, viewer);
 75:     ISDestroy(&cdofIS);
 76:     ISCreateGeneral(comm, n, coffs, PETSC_OWN_POINTER, &coffIS);
 77:     PetscObjectSetName((PetscObject)coffIS, "atlasOff");
 78:     ISView(coffIS, viewer);
 79:     ISDestroy(&coffIS);
 80:     PetscViewerHDF5PopGroup(viewer);
 81:     ISCreateGeneral(comm, m, cinds, PETSC_OWN_POINTER, &cindIS);
 82:     PetscObjectSetName((PetscObject)cindIS, "bcIndices");
 83:     ISView(cindIS, viewer);
 84:     ISDestroy(&cindIS);
 85:   }
 86:   return 0;
 87: }

 89: PetscErrorCode PetscSectionView_HDF5_Internal(PetscSection s, PetscViewer viewer)
 90: {
 91:   PetscInt numFields, f;

 93:   PetscViewerHDF5PushGroup(viewer, "section");
 94:   PetscSectionGetNumFields(s, &numFields);
 95:   PetscViewerHDF5WriteAttribute(viewer, NULL, "numFields", PETSC_INT, (void *)&numFields);
 96:   PetscSectionView_HDF5_SingleField(s, viewer);
 97:   for (f = 0; f < numFields; ++f) {
 98:     char        fname[PETSC_MAX_PATH_LEN];
 99:     const char *fieldName;
100:     PetscInt    fieldComponents, c;

102:     PetscSNPrintf(fname, sizeof(fname), "field%" PetscInt_FMT, f);
103:     PetscViewerHDF5PushGroup(viewer, fname);
104:     PetscSectionGetFieldName(s, f, &fieldName);
105:     PetscViewerHDF5WriteAttribute(viewer, NULL, "fieldName", PETSC_STRING, fieldName);
106:     PetscSectionGetFieldComponents(s, f, &fieldComponents);
107:     PetscViewerHDF5WriteAttribute(viewer, NULL, "fieldComponents", PETSC_INT, (void *)&fieldComponents);
108:     for (c = 0; c < fieldComponents; ++c) {
109:       char        cname[PETSC_MAX_PATH_LEN];
110:       const char *componentName;

112:       PetscSNPrintf(cname, sizeof(cname), "component%" PetscInt_FMT, c);
113:       PetscViewerHDF5PushGroup(viewer, cname);
114:       PetscSectionGetComponentName(s, f, c, &componentName);
115:       PetscViewerHDF5WriteAttribute(viewer, NULL, "componentName", PETSC_STRING, componentName);
116:       PetscViewerHDF5PopGroup(viewer);
117:     }
118:     PetscSectionView_HDF5_SingleField(s->field[f], viewer);
119:     PetscViewerHDF5PopGroup(viewer);
120:   }
121:   PetscViewerHDF5PopGroup(viewer);
122:   return 0;
123: }

125: static PetscErrorCode PetscSectionLoad_HDF5_SingleField_SetConstraintIndices(PetscSection s, IS cindIS, IS coffIS)
126: {
127:   MPI_Comm        comm;
128:   PetscInt        pStart, pEnd, p, M, m, i, cdof;
129:   const PetscInt *data;
130:   PetscInt       *cinds;
131:   const PetscInt *coffs;
132:   PetscInt       *coffsets;
133:   PetscSF         sf;
134:   PetscLayout     layout;

136:   PetscObjectGetComm((PetscObject)s, &comm);
137:   PetscSectionGetChart(s, &pStart, &pEnd);
138:   ISGetSize(cindIS, &M);
139:   ISGetLocalSize(cindIS, &m);
140:   PetscMalloc1(m, &coffsets);
141:   ISGetIndices(coffIS, &coffs);
142:   for (p = pStart, m = 0; p < pEnd; ++p) {
143:     PetscSectionGetConstraintDof(s, p, &cdof);
144:     for (i = 0; i < cdof; ++i) coffsets[m++] = coffs[p - pStart] + i;
145:   }
146:   ISRestoreIndices(coffIS, &coffs);
147:   PetscSFCreate(comm, &sf);
148:   PetscLayoutCreate(comm, &layout);
149:   PetscLayoutSetSize(layout, M);
150:   PetscLayoutSetLocalSize(layout, m);
151:   PetscLayoutSetBlockSize(layout, 1);
152:   PetscLayoutSetUp(layout);
153:   PetscSFSetGraphLayout(sf, layout, m, NULL, PETSC_OWN_POINTER, coffsets);
154:   PetscLayoutDestroy(&layout);
155:   PetscFree(coffsets);
156:   PetscMalloc1(m, &cinds);
157:   ISGetIndices(cindIS, &data);
158:   PetscSFBcastBegin(sf, MPIU_INT, data, cinds, MPI_REPLACE);
159:   PetscSFBcastEnd(sf, MPIU_INT, data, cinds, MPI_REPLACE);
160:   ISRestoreIndices(cindIS, &data);
161:   PetscSFDestroy(&sf);
162:   PetscSectionSetUpBC(s);
163:   for (p = pStart, m = 0; p < pEnd; ++p) {
164:     PetscSectionGetConstraintDof(s, p, &cdof);
165:     PetscSectionSetConstraintIndices(s, p, &cinds[m]);
166:     m += cdof;
167:   }
168:   PetscFree(cinds);
169:   return 0;
170: }

172: static PetscErrorCode PetscSectionLoad_HDF5_SingleField(PetscSection s, PetscViewer viewer)
173: {
174:   MPI_Comm comm;
175:   PetscInt pStart, pEnd, p, N, n, M, m;
176:   #if defined(PETSC_USE_DEBUG)
177:   PetscInt N1, M1;
178:   #endif
179:   PetscBool       hasConstraints, includesConstraints;
180:   IS              dofIS, offIS, cdofIS, coffIS, cindIS;
181:   const PetscInt *dofs, *offs, *cdofs;
182:   PetscLayout     map;

184:   PetscObjectGetComm((PetscObject)s, &comm);
185:   PetscViewerHDF5ReadAttribute(viewer, NULL, "includesConstraints", PETSC_BOOL, NULL, (void *)&includesConstraints);
186:   PetscSectionSetIncludesConstraints(s, includesConstraints);
187:   PetscSectionGetChart(s, &pStart, &pEnd);
188:   n = pEnd - pStart;
189:   #if defined(PETSC_USE_DEBUG)
190:   MPIU_Allreduce(&n, &N1, 1, MPIU_INT, MPI_SUM, comm);
191:   #endif
192:   ISCreate(comm, &dofIS);
193:   PetscObjectSetName((PetscObject)dofIS, "atlasDof");
194:   PetscViewerHDF5ReadSizes(viewer, "atlasDof", NULL, &N);
195:   #if defined(PETSC_USE_DEBUG)
197:   #endif
198:   ISGetLayout(dofIS, &map);
199:   PetscLayoutSetSize(map, N);
200:   PetscLayoutSetLocalSize(map, n);
201:   ISLoad(dofIS, viewer);
202:   ISCreate(comm, &offIS);
203:   PetscObjectSetName((PetscObject)offIS, "atlasOff");
204:   PetscViewerHDF5ReadSizes(viewer, "atlasOff", NULL, &N);
205:   #if defined(PETSC_USE_DEBUG)
207:   #endif
208:   ISGetLayout(offIS, &map);
209:   PetscLayoutSetSize(map, N);
210:   PetscLayoutSetLocalSize(map, n);
211:   ISLoad(offIS, viewer);
212:   ISGetIndices(dofIS, &dofs);
213:   ISGetIndices(offIS, &offs);
214:   for (p = pStart, n = 0; p < pEnd; ++p, ++n) {
215:     PetscSectionSetDof(s, p, dofs[n]);
216:     PetscSectionSetOffset(s, p, offs[n]);
217:   }
218:   ISRestoreIndices(dofIS, &dofs);
219:   ISRestoreIndices(offIS, &offs);
220:   ISDestroy(&dofIS);
221:   ISDestroy(&offIS);
222:   PetscViewerHDF5ReadAttribute(viewer, NULL, "hasConstraints", PETSC_BOOL, NULL, (void *)&hasConstraints);
223:   if (hasConstraints) {
224:     PetscViewerHDF5PushGroup(viewer, "bc");
225:     ISCreate(comm, &cdofIS);
226:     PetscObjectSetName((PetscObject)cdofIS, "atlasDof");
227:     PetscViewerHDF5ReadSizes(viewer, "atlasDof", NULL, &N);
228:   #if defined(PETSC_USE_DEBUG)
230:   #endif
231:     ISGetLayout(cdofIS, &map);
232:     PetscLayoutSetSize(map, N);
233:     PetscLayoutSetLocalSize(map, n);
234:     ISLoad(cdofIS, viewer);
235:     ISGetIndices(cdofIS, &cdofs);
236:     for (p = pStart, n = 0; p < pEnd; ++p, ++n) PetscSectionSetConstraintDof(s, p, cdofs[n]);
237:     ISRestoreIndices(cdofIS, &cdofs);
238:     ISDestroy(&cdofIS);
239:     ISCreate(comm, &coffIS);
240:     PetscObjectSetName((PetscObject)coffIS, "atlasOff");
241:     PetscViewerHDF5ReadSizes(viewer, "atlasOff", NULL, &N);
242:   #if defined(PETSC_USE_DEBUG)
244:   #endif
245:     ISGetLayout(coffIS, &map);
246:     PetscLayoutSetSize(map, N);
247:     PetscLayoutSetLocalSize(map, n);
248:     ISLoad(coffIS, viewer);
249:     PetscViewerHDF5PopGroup(viewer);
250:     ISCreate(comm, &cindIS);
251:     PetscObjectSetName((PetscObject)cindIS, "bcIndices");
252:     PetscViewerHDF5ReadSizes(viewer, "bcIndices", NULL, &M);
253:     if (!s->bc) m = 0;
254:     else PetscSectionGetStorageSize(s->bc, &m);
255:   #if defined(PETSC_USE_DEBUG)
256:     MPIU_Allreduce(&m, &M1, 1, MPIU_INT, MPI_SUM, comm);
258:   #endif
259:     ISGetLayout(cindIS, &map);
260:     PetscLayoutSetSize(map, M);
261:     PetscLayoutSetLocalSize(map, m);
262:     ISLoad(cindIS, viewer);
263:     PetscSectionLoad_HDF5_SingleField_SetConstraintIndices(s, cindIS, coffIS);
264:     ISDestroy(&coffIS);
265:     ISDestroy(&cindIS);
266:   }
267:   return 0;
268: }

270: PetscErrorCode PetscSectionLoad_HDF5_Internal(PetscSection s, PetscViewer viewer)
271: {
272:   MPI_Comm comm;
273:   PetscInt N, n, numFields, f;

275:   PetscObjectGetComm((PetscObject)s, &comm);
276:   PetscViewerHDF5PushGroup(viewer, "section");
277:   PetscViewerHDF5ReadAttribute(viewer, NULL, "numFields", PETSC_INT, NULL, (void *)&numFields);
278:   if (s->pStart < 0 && s->pEnd < 0) n = PETSC_DECIDE;
279:   else {
282:     n = s->pEnd;
283:   }
284:   if (numFields > 0) PetscSectionSetNumFields(s, numFields);
285:   PetscViewerHDF5ReadSizes(viewer, "atlasDof", NULL, &N);
286:   if (n == PETSC_DECIDE) PetscSplitOwnership(comm, &n, &N);
287:   PetscSectionSetChart(s, 0, n);
288:   PetscSectionLoad_HDF5_SingleField(s, viewer);
289:   for (f = 0; f < numFields; ++f) {
290:     char     fname[PETSC_MAX_PATH_LEN];
291:     char    *fieldName;
292:     PetscInt fieldComponents, c;

294:     PetscSNPrintf(fname, sizeof(fname), "field%" PetscInt_FMT, f);
295:     PetscViewerHDF5PushGroup(viewer, fname);
296:     PetscViewerHDF5ReadAttribute(viewer, NULL, "fieldName", PETSC_STRING, NULL, &fieldName);
297:     PetscSectionSetFieldName(s, f, fieldName);
298:     PetscFree(fieldName);
299:     PetscViewerHDF5ReadAttribute(viewer, NULL, "fieldComponents", PETSC_INT, NULL, (void *)&fieldComponents);
300:     PetscSectionSetFieldComponents(s, f, fieldComponents);
301:     for (c = 0; c < fieldComponents; ++c) {
302:       char  cname[PETSC_MAX_PATH_LEN];
303:       char *componentName;

305:       PetscSNPrintf(cname, sizeof(cname), "component%" PetscInt_FMT, c);
306:       PetscViewerHDF5PushGroup(viewer, cname);
307:       PetscViewerHDF5ReadAttribute(viewer, NULL, "componentName", PETSC_STRING, NULL, &componentName);
308:       PetscSectionSetComponentName(s, f, c, componentName);
309:       PetscFree(componentName);
310:       PetscViewerHDF5PopGroup(viewer);
311:     }
312:     PetscSectionLoad_HDF5_SingleField(s->field[f], viewer);
313:     PetscViewerHDF5PopGroup(viewer);
314:   }
315:   PetscViewerHDF5PopGroup(viewer);
316:   return 0;
317: }

319: #endif