1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
 void
 ReadBackElementHistory(
    DbInfo_t *dbInfo,    /* database info object */
    int myProcNum,       /* processor rank in MPI_COMM_WORLD */
    int histElem,        /* the element for which history */
    int *numReadBack,    /* number of dump for which history was read back */
    ElementHistory_t **hist /* the resulting history buffer */
 )
 {
    int i, numHist;
    ElementHistory_t *histBuf;

    SAF_Db *db = dbInfo->db;
    SAF_Cat procs = dbInfo->blocks;
    SAF_Cat elems = dbInfo->elems;
    SAF_Set invariantSet = dbInfo->mainMesh;
    SAF_Set *dumpSets = NULL;

    /* describe the self collection on the aggregate set */
    saf_describe_collection(SAF_ALL, &invariantSet, SAF_SELF(db), NULL, &numHist, NULL, NULL, &dumpSets);

    /* allocate some space for the history buffer */
    histBuf = (ElementHistory_t *) malloc(numHist * sizeof(ElementHistory_t));

    /* for each member of this collection, describe the set and the proc's collection on that set */
    for (i = 0; i < numHist; i++)
    {
       int numPieces, numRels, numFields, numElems, localElemID;
       float dofVal, *dofVals=&dofVal;
       SAF_Set *procSets = NULL;
       SAF_Set myPiece;
       SAF_Rel *ssRels=NULL;
       SAF_Field *theFields = NULL;
       RelIndex_t *elemList = NULL;

       /* get the processor collection at this dump */
       saf_describe_collection(SAF_ALL, &(dumpSets[i]), &procs, NULL, &numPieces, NULL, NULL, &procSets);

       /* this processor will take responsibility for the member of this collection whose index in the collection is the
          same as its rank */
       myPiece = procSets[myProcNum];

       /* how many elements are on this processor's piece? */
       saf_describe_collection(SAF_EACH, &myPiece, &elems, NULL, &numElems, NULL, NULL, NULL);

       /* everybody reads their respective subset relations */
       saf_find_subset_relations(SAF_EACH, db, &(dumpSets[i]), &myPiece, &elems, &elems, SAF_BOUNDARY_FALSE, SAF_BOUNDARY_FALSE,
          &numRels, &ssRels);
       if (numRels != 1)
          AbortThisMess("found more than 1 dumpSet to procSet subset relation on elems");

       /* read the subset relation data */
       saf_read_subset_relation(SAF_EACH, &(ssRels[0]), NULL, (void**) &elemList, NULL);

       /* everybody finds the pressure field on thier respective pieces */
       saf_find_fields(SAF_EACH, db, &myPiece, "pressure", SAF_ANY_QUANTITY, SAF_ALGTYPE_ANY, SAF_ANY_BASIS,
          SAF_ANY_UNIT, SAF_ANY_CAT, SAF_ANY_RATIO, SAF_ANY_CAT, SAF_ANY_EFUNC, &numFields, &theFields);
       if (numFields != 1)
          AbortThisMess("found more than 1 field named \"pressure\" on myPiece set");

       /* everybody checks to see if they own the hist elem. This assumes elemList is sorted in ascending order. */
       localElemID = IndexOfValue(elemList, numElems, (RelIndex_t) histElem);

       /* everybody reads (partially) the field (most procs attempt to read 0 dofs) */
       saf_read_field(SAF_EACH, &(theFields[0]), NULL, localElemID==-1?0:1, SAF_TUPLES, &localElemID, (void**) &dofVals);

       AccumulateElementHistory(myProcNum, localElemID, dofVals[0], &histBuf[i]);

       free(theFields); theFields = NULL;
       free(elemList); elemList = NULL;
       free(ssRels); ssRels = NULL;
       free(procSets); procSets = NULL;
    }

    free(dumpSets);

    /* setup return values */
    *numReadBack = numHist;
    *hist = histBuf;

    return;

 }